Field Quantization for Radiative Decay of Plasmons in Finite and Infinite Geometries

TL;DR

This paper develops a quantum field theory framework for plasmons in nanostructures, enabling calculation of their radiative decay rates across various geometries relevant to nanotechnology and quantum sensing.

Contribution

It introduces a quasistatic quantization method for plasmons in complex geometries and derives analytical expressions for their radiative decay rates.

Findings

Quantized plasmon fields are obtained for multiple geometries.

Analytical decay rate formulas are derived using perturbation theory.

The approach applies to both finite and infinite nanostructures.

Abstract

We investigate field quantization in high-curvature geometries. The models and calculations can help with understanding the elastic and inelastic scattering of photons and electrons in nanostructures and probe-like metallic domains. The results find important applications in high-resolution photonic and electronic modalities of scanning probe microscopy, nano-optics, plasmonics, and quantum sensing. Quasistatic formulation, leading to nonretarded quantities, is employed and justified on the basis of the nanoscale, here subwavelength, dimensions of the considered domains of interest. Within the quasistatic framework, we represent the nanostructure material domains with frequency-dependent dielectric functions. Quantities associated with the normal modes of the electronic systems, the nonretarded plasmon dispersion relations, eigenmodes, and fields are then calculated for several geometric entities of use in nanoscience and nanotechnology. From the classical energy of the charge density oscillations in the modeled nanoparticle, we then derive the Hamiltonian of the system, which is used for quantization. The quantized plasmon field is obtained and, employing an interaction Hamiltonian derived from the first-order perturbation theory within the hydrodynamic model of the electron gas, we obtain an analytical expression for the radiative decay rate of the plasmons. The established treatment is applied to multiple geometries to investigate the quantized charge density oscillations on their bounding surfaces. Specifically, using one sheet of a two-sheeted hyperboloid of revolution, paraboloid of revolution, and cylindrical domains, all with one infinite dimension, and the finite spheroidal and toroidal domains are treated. ...

Equations

$$F({\mathbf{r}}_i)=m\frac{d^2{\mathbf{r}}_i}{d{t}^2},$$

$$F({\mathbf{r}}_i)=-k\mathbf{r},$$

$$m\frac{d^2{\mathbf{r}}_i}{d{t}^2}=-k{\mathbf{r}}_i,$$

$${\mathbf{r}}_i(t)=A\cos ({\omega }_{0}t)+B\sin ({\omega }_{0}t),$$

$$\frac{d^2{\mathbf{r}}_i}{d{t}^2}+\gamma \frac{d{\mathbf{r}}_i}{dt}+{\omega }_0^2{\mathbf{r}}_i=-\frac{e}{m_e}\vec{E},$$

$$\vec{E}(t)=\cos (\omega t+\theta ),$$

$$\vec{E}(t)={\vec{E}}_0\text{Re}(e^{-(i\omega t+\theta )}),$$

$${\mathbf{r}}_i(t)={\vec{R}}_0\text{Re}(e^{-(i\omega t+\tilde{\theta })}),$$

$$\left(-{\omega }^2-i\gamma \omega +{\omega }_0^2\right){\vec{R}}_0=-\frac{e}{m_0}{\vec{E}}_0,$$

$${\vec{R}}_0=\frac{-e{\vec{E}}_0}{m_e\left(-{\omega }^2-i\gamma \omega +{\omega }_0^2\right)},$$

$$\vec{P}=n_{0}p(t),$$

$$n_{0}p(t)=-n_{0}q{\mathbf{r}}_i=\frac{n_0{e}^2{\vec{E}}_0}{m_e\left[{\omega }_0^2-\left({\omega }^2+i\gamma \omega \right)\right]}.$$

$$\vec{D}=\vec{E}+4\pi \vec{P}.$$

$$\vec{P}=-\frac{{\chi }_e(\omega )}{4\pi }\vec{E},$$

$$\vec{D}=\frac{1}{4\pi }\epsilon (\omega )\vec{E},$$

$$\epsilon (\omega )=1+{\chi }_e(\omega ),$$

$$\epsilon (\omega )=1+\frac{{\omega }_p^2}{{\omega }_0^2-{\omega }^2-i\gamma \omega },$$

$${\omega }_p^2=\frac{4\pi n_0{e}^2}{m_e},$$

$$\epsilon (\omega )=1-\frac{{\omega }_p^2}{\omega (\omega +i\gamma )}.$$

$$\epsilon (\omega )=1-\frac{{\omega }_p^2}{{\omega }^2}.$$

$$\vec{\nabla }\cdot \vec{E}$$

$$\vec{\nabla }\times \vec{E}$$

$$\vec{\nabla }\cdot \vec{B}$$

$$c^2\vec{\nabla }\times \vec{B}$$

$$\vec{E}=-\left[\frac{1}{c}\frac{\partial \mathbf{A}}{\partial t}+\vec{\nabla }\Phi (\mathbf{r},t)\right],$$

$$\vec{E}=-\vec{\nabla }\Phi (\mathbf{r},t).$$

$${\vec{\nabla }}^2\Phi (\mathbf{r},t)=-4\pi \rho .$$

$${\vec{\nabla }}^2\Phi (\mathbf{r},t)=0.$$

$${\oint }_{\partial \Pi }\vec{D}\cdot \mathbf{n}d\mathcal{A}={\int }_{\Pi }\rho d\mathcal{V}.$$

$${\oint }_{\partial \Pi }\vec{D}\cdot \mathbf{n}d\mathcal{A}=({\vec{D}}_2-{\vec{D}}_1)\cdot \mathbf{n}\Delta \mathcal{A},$$

$${\int }_{\Pi }\rho d\mathcal{V}={\int }_{\partial \Pi }\sigma d\mathcal{A}=\sigma \Delta \mathcal{A}.$$

$$({\vec{D}}_2-{\vec{D}}_1)\cdot \mathbf{n}=\sigma ,$$

Full text

1. 1 INTRODUCTION

2. 1.1 Introduction

3. 1.2 Preliminaries

4. 1.3 General Approach

5. 1.4 Harmonic Oscillator Model

6. 1.5 Quantum Field & Second Quantization

7. 1.6 Quantization of the Electromagnetic field

8. 1.7 Radiative Decay Rate

9. 2 RADIATIVE DECAY ON PARABOLOIDAL, HYPERBOLOIDAL,                   & PROLATE SPHEROIDAL SURFACES

10. 2.1 Introduction

11. 2.2 Paraboloidal Surfaces

12. 2.2.1 Nonretarded Potential & Dispersion Relations

13. 2.2.2 Classical Energy

14. 2.2.3 Field Quantization & Interaction of Plasmons with Photons

15. 2.2.4 Emission of Light & Radiative Decay

16. 2.3 Hyperboloidal Surfaces

17. 2.3.1 Nonretarded Potential & Dispersion Relations

18. 2.3.2 Classical Energy & Decay Rate

19. 2.4 Prolate Spheroidal Surfaces

20. 2.4.1 Nonretarded Potential & Dispersion Relations

21. 2.4.2 Classical Energy & Radiative Decay Rate

22. 2.5 Comparison

23. 3 AN APPROXIMATION FOR RADIATIVE DECAY RATE ON T-                   HE NON-SIMPLY CONNECTED SURFACE OF A TORUS

24. 3.1 Introduction

25. 3.2 Nonretarded Potential & Dispersion Relations

26. 3.3 Classical Energy

27. 3.4 Interaction Hamiltonian & Radiative Decay Rate

28. 4 CONCLUSION

29. A.1 Radiative Decay on Spherical Surfaces

30. A.2 Radiative Decay on Cylindrical Surfaces

31. B.3 Collected Facts about Coordinate Systems

32. B.3.1 An Arbitrary Coordinate System

33. B.3.2 Spherical Coordinate

34. B.3.3 Cylindrical Coordinate

35. B.3.4 Parabaloidal Coordinate

36. B.3.5 Prolate Spheroidal Coordinate

37. B.3.6 Toroidal Coordinate

38. B.4 Collected Facts about Special Functions

39. B.4.1 Bessel Functions

40. B.4.2 Associated Legendre Functions

41. B.4.3 Conical Functions

42. B.5 Sturm-Liouville Problem

1. 1.1 A tiny pillbox of volume $\mathcal{V}$ at the interface between two
2. 1.2 The solid surface of revolution region defined by fixing $\zeta=\zeta_{0}$,
3. 2.1 Modeling systems: One sheet of a two-sheeted hyperboloid
4. 2.2 Modeling systems for the spatial distribution of the lowest
5. 2.3 Paraboloidal nonretarded surface plasmon dispersion relations
6. 2.4 Modeling systems for the spatial distribution
7. 2.5 Hyperboloidal nonretarded surface plasmon dispersion relations
8. 2.6 Figure shows the spatial distribution of the lowest
9. 2.7 Prolate spheroidal nonretarded surface plasmon dispersion relations
10. 2.8 Here shows the calculations of the wavelength against
11. 2.9 The contour plots are obtained by plotting decay rate given in
12. 2.10 Here we consider a sphere as a limiting case for prolate spheroid.
13. 2.11 Modeling systems and their potential distributions
14. 2.12 The spatial distribution of the one low and one high lying
15. 2.13 Paraboloidal, hyperboloidal, and prolate spheroidal nonretarded
16. 2.14 Paraboloidal and hyperboloidal nonretarded surface
17. 2.15 The radiative decay rate of plasmons engendered
18. 2.16 Comparison of the higher order modes’ radiative
19. 2.17 Curvature induced shift in the radiation pattern

ABSTRACT

We investigate field quantization in high-curvature geometries. The models and calculations can help with understanding the elastic and inelastic scattering of photons and electrons in nanostructures and probe-like metallic domains. The results find important applications in high-resolution photonic and electronic modalities of scanning probe microscopy, nano-optics, plasmonics, and quantum sensing. Quasistatic formulation, leading to nonretarded quantities, is employed and justified on the basis of the nanoscale, here subwavelength, dimensions of the considered domains of interest. Within the quasistatic framework, we represent the nanostructure material domains with frequency-dependent dielectric functions. Quantities associated with the normal modes of the electronic systems, the nonretarded plasmon dispersion relations, eigenmodes, and fields are then calculated for several geometric entities of use in nanoscience and nanotechnology. From the classical energy of the charge density oscillations in the modeled nanoparticle, we then derive the Hamiltonian of the system, which is used for quantization. The quantized plasmon field is obtained and, employing an interaction Hamiltonian derived from the first-order perturbation theory within the hydrodynamic model of an electron gas, we obtain an analytical expression for the radiative decay rate of the plasmons. The established treatment is applied to multiple geometries to investigate the quantized charge density oscillations on their bounding surfaces. Specifically, using one sheet of a two-sheeted hyperboloid of revolution, paraboloid of revolution, and cylindrical domains, all with one infinite dimension, and the finite spheroidal and toroidal domains are treated. In addition to a comparison of the paraboloidal and hyperboloidal results, interesting similarities are observed for the paraboloidal domains with respect to the surface modes and radiation patterns of a prolate spheroid, a finite geometric domain highly suitable for modeling of nanoparticles such as quantum dots. The prolate and oblate spheroidal calculations are validated by comparison to the spherical case, which is obtained as a special case of a spheroid. In addition to calculating the potential and field distributions, and dispersion relations, we study the angular intensity and the relation between the emission angle with the rate of radiative decay. The various morphologies are compared for their plasmon dispersion properties, field distributions, and radiative decay rates, which are shown to be consistent. For the specific case of a nanoring, modeled in the toroidal geometry, significant complexity arises due to an inherent coupling among the various modes. Within reasonable approximations to decouple the modes, we study the radiative decay channel for a vacuum bounded single solid nanoring by quantizing the fields associated with charge density oscillations on the nanoring surface. Further suggestions are made for future studies. The obtained results are relevant to other material domains that model a nanostructure such as a probe tip, quantum dot, or nanoantenna.

CHAPTER 1INTRODUCTION

Introduction

Interaction of photos with surface plasmons has been studied extensively during the past five decades, particularly in relation to their application in scanning probe microscopy (SPM), they present a high potential for emerging applications in fields such as quantum sensing [1, 2, 3, 4, 5]. It is known as a branch of microscopy that forms images of surfaces using a physical probe. Throughout these studies, various geometries have been investigated. In particular, Ritchie together with Crowell, Little, Ashley and Ferell [6, 7, 8, 19, 10] have studied the interaction of general and special cases of surface of a metallic sphere and oblate spheroid with photons and derives the very first relations describing the field quntization of surface plasmons. Inspired by these results, quantization relations for a special case of a long string-like cylinder has been studied by Burmistrova [11]. In this work, we wish to investigate the interaction in more detail as well as to extend the obtained results to both new finite and infinite geometries such as cylinder, paraboloid, hyperboloid, prolate spheroid and the non-simply connected case of a torus. It should be pointed out that the surface plasmon-photon interaction (SPP) can be analyzed for various cases such as absorption, emission, Thomson and Rayleigh scattering, and scattering cross-sections. Our main focus in this study, however, is the emission case and its application to the decay rate. As we shall see later in this chapter, there is a close relation between both cases of emission and absorption in terms of the interaction matrix element. As a result, our study also addresses the absorption phenomena [8, 12].

Throughout this dissertation, we consider a metallic particle confined in vacuum whose shape is described by one of the above mentioned geometries. The conduction electrons and ion centers in a metal together to form a plasma. The authors in [13], while studying the collective oscillations in the metallic plasma, showed that through the plasma system exhibits resonances along the directions that charge density waves propagate, with no effect of damping. Ritchie [10] in 1957 showed that the electron density can be supported at the surface of a bounded plane free-electron gas. In [14, 15], it has been discussed both the property of plasmons and their interaction with charged particles. The properties of plasmons and details on how they interact with light or charges particles, such as photons, are also discussed and reviewed in [16, 17, 18, 9, 20]. To describe the properties of a metal, mostly from a dynamical point of view, the complex-valued dielectric function $\varepsilon(\omega)$ was used. Generally speaking, the dielectric function specifies the relation between the frequency and the wave-vector in a certain metal. It is also known as permitivity. The main assumption here is that in the incident of large wave-vectors, the dielectric function would be no longer wave-vector-dependent [13].

Throughout this chapter, we adopt the conventions of Ritchie mostly presented in [7], unless otherwise specified. The dielectric function, $\varepsilon(\omega)$, is mostly used to express certain properties of a metal. It is a complex and frequency-dependent function which in general also depends on wave vector of the plasmon field, $k_{p}\equiv 2\pi/\lambda_{p}$ where $\lambda_{p}$ denotes the plasmon wavelength. However, under certain suitable assumptions, the wave vector dependency could be ignored [13]. Under the assumption that the wavelength is large enough, the dielectric function can be treated as independent from the wave vector. This requirement is met for all cases which have been investigated in this work.

Free-electron-gas dispersion relation

We start by considering the case of a spring oscillator, with no damping forces. The equation of motion, using Newton’s second law, can be written as [21]:

[TABLE]

where $m$ denotes the mass and $\mathbf{r}_{i}$ is the displacement vector of $i$-th electron with respect to the equilibrium position. Hooke’s law, on the other hand, gives:

[TABLE]

where $k$ denotes the spring constant. From Eqs. (1.1) and (1.2), one could write the second order differential equation:

[TABLE]

whose solutions are given by:

[TABLE]

for some constants $A$ and $B$, where $\omega_{0}=\sqrt{\frac{k}{m}}$ is the resonant frequency. The equation of motion for the $i$-th electron, including both the driving forces and the damping $\gamma$ is written as [22]:

[TABLE]

where $e$, $m_{e}$ denote the electric charge and the mass of an electron, respectively, and $\vec{E}$ represents the electric field. The terms $\omega^{2}_{0}\mathbf{r}_{i}$ and $-\frac{e}{m_{0}}\,\vec{E}$ represent the restoring and electric forces, respectively. The time dependent electric field of a light wave inducing oscillations is given by:

[TABLE]

where $\omega$ is the frequency of the system and $\theta$ is denotes the phase of the wave. Eq. (1.6) could be written as:

[TABLE]

where $\vec{E}_{0}$ denotes the amplitude. We are mostly interested in those solutions of Eq. (1.5) which have the form similar to Eq. (1.7), namely,

[TABLE]

where $\vec{R}_{0}$ and $\tilde{\theta}$ represent the amplitude and phase of the oscillations, respectively. Substituting Eqs.  (1.7) and (1.8) in Eq. (1.5), one may write:

[TABLE]

and hence:

[TABLE]

where $\omega$ is the frequency of the electric field. The electric polarization vector $\vec{P}$ is given by [23]:

[TABLE]

where $p(t)$ denotes the dipole moment per unit volume and $n_{0}$ is the number of effective electrons per unit volume. One can write:

[TABLE]

The displacement vector $\vec{D}$ throughout the volume of the gas is written as [24]:

[TABLE]

Assuming the substance is isotropic, the relation between polarization vector and electric field is given by:

[TABLE]

where $\chi_{e}(\omega)$ denotes the electrical susceptibility of a dielectric material. Using Eq. (1.14), one may write Eq. (1.13) as:

[TABLE]

where we put:

[TABLE]

denoting dielectric function for a materiel. Substituting Eqs. (1.12) and (1.16) in Eq. (1.15), one finds the dielectric function for a material as:

[TABLE]

where $\omega_{p}$ denotes the bulk plasma frequency also known as plasmon frequency, and is given by:

[TABLE]

in Gaussian’s unit. The Drude model for the dielectric function of metals is obtained in the assumption that the free electrons in metals are not bound to any atoms and are not subject to restoring forces. This implies the spring constant $k$ to be zero and therefore $\omega_{0}=0$ [25, 22]. Using this in Eq. (1.17) leads to:

[TABLE]

In the absence of damping forces, the free-electron-gas dielectric function reduces to the well-known expression:

[TABLE]

In the generalized Drude model one may assume $\gamma$ is frequency-independent and that it is complex valued for $\varepsilon$ to undergo Kramers Kroing dispersion relation, in which the real part stays constant with respect to plasma energy, while the imaginary part could be neglected. Thus the simple Drude model stays valid throughout all the steps up to plasma energy. If one ignores the damping force, then the real dielectric function of a free electron gas could be expressed by Eq. (1.20). In order to determine the value of dielectric function for surface plasmon oscillations, we impose the requirement of continuity of electric potential and normal component of the displacement vector across the surface of the metal. This is where the geometric effects of the surface plays a role in finding the values of dielectric function $\varepsilon(\omega)$ and hence the value of frequency $\omega$ itself. We will discuss these geometric effects in the coming chapters thoroughly.

Throughout the whole discussion, we took on the assumption that all the interactions occur without a delay, i.e. they are instantaneous. In other words, the system is considered in a way that electrostatic solutions are applicable. This gives us a relatively fair approximation as long as our assumption of having large wavevectors for plasmons (comparing to that of light) is valid. This large-wavevector region is known in the literature the as nonretarted region. This allows us to consider Poisson and Laplace’s equations instead of wave equation to determine the allowed frequencies. What one might loose under this assumption is thoroughly discussed by Ritchie in [9]. We devote this chapter to describing our approach regarding the radiative decay of the surface plasmons of an excited metallic surface which mostly follows the one given in [8]. Furthermore, the geometric effect of the curvatures for such metallic surface has been investigated thoroughly in [27].

Preliminaries

The electrical size is determined by comparing physical dimension of an electronic structure and the signal wavelength. Let us consider a field which varies at low enough frequency when compared to its wavelength. Therefore, considering a small lapse of time, one could assume that the system remains in an internal equilibrium. In other words, the time delay due to wave propagation from one point to any other point can be ignored if the circuit geometry is relatively small compared to the signal wavelength (which prevents the field distribution to vary significantly). Under these circumstances, quantities like potential and surface current (which will be defined later in this chapter) will not be dependent on position and time [28]. We start this section with elaborating the mathematical formulation in greater detail.

Let us recall Maxwell’s equations in their most standard formulations as:

[TABLE]

where $\vec{E}$, $\vec{B}$, $\mathbf{J}$, and $\rho$ denote electric field, magnetic field, current and charge density, respectively. As customary, $c$ indicates the speed of light. In the absence of any magnetic field $\vec{B}$, one could set:

[TABLE]

where $\Phi(\mathbf{r},t)$ and $\mathbf{A}$ are scalar potential and vector potential field, respectively, and $\mathbf{r}$ denotes the position vector (see any classical reference on this topic, also [24, 29, 23, 30, 31, 32]). Furthermore, if the wave speed is relatively small compared to $c$, one may ignore the first term in the right-hand side of Eq. (1.25). Under this assumption, we can write Eq. (1.25) as:

[TABLE]

By means of Eqs. (1.21) and (1.26), we may assume that the scalar potential $\Phi$ satisfies Poisson’s equation:

[TABLE]

In a charge-free system, Poisson equation is known as Laplace’s equation:

[TABLE]

It is of fundamental importance that under proper boundary conditions the solution to either Poisson or Laplace’s equations is unique [33]. In terms of boundary conditions, they are discussed in the next section, for the convenience of the reader and the sake of self consistency.

General Approach

It is well-known that the retarded potentials are due to time-varying charge densities which are associated with surface plasmons [34]. In the non-retarded limit, the instantaneous potentials are calculated assuming the source is a static charge density. Retardation for the fields can be ignored near the surface of a particle if the dimensions of the particle are much less than the wavelength of light at frequency $\omega$.

We start by considering a pillbox with volume $\mathcal{V}$ at the interface of two media with dielectric function $\varepsilon_{1}$ and $\varepsilon_{2}$, such that the areas filled with each material are the same (see Fig. 1.1). Maxwell’s equations imply the continuity of the tangential component of the electric field $\vec{E}$ and the normal component of displacement vector $\vec{D}$. We determine these boundary conditions using Gauss and Stoke’s laws [35].

In view of Eq. (1.15), for any external sources in vacuum with $\varepsilon(\omega)=1$, using Gauss’s law for a continuous charge density $\rho$, we can write:

[TABLE]

Utilizing the identity:

[TABLE]

where $\vec{D}_{1}$ and $\vec{D}_{2}$ correspond to displacement vectors of the inside and outside regions and $\mathbf{n}$ is the surface unit normal vector. Since the polarization charge is only confined to the surface, we may write:

[TABLE]

It follows from Eqs. (1.30) and (1.31):

[TABLE]

where $\sigma$ notes the surface charge density on the surface boundary $\partial\Pi$. Eq. (1.32) presents no polarization charge and therefore indicates that the normal component of displacement vector $\vec{D}$ has a shift equals to surface charge density $\sigma$. This is the first boundary condition. Using Eq. (1.15), Eq. (1.32) may be also written as:

[TABLE]

Making use of Eq. (1.26), we may also write:

[TABLE]

the second boundary condition is widely used to determine the dielectric function and, as a consequence, normal mode frequencies for the given geometry. Using the definition $\oint_{C}\vec{E}\cdot dl=0$ on a closed path $C$, it follows that:

[TABLE]

In virtue of Stoke’s theorem, we have:

[TABLE]

Under electrostatic condition, we derive the identity $\vec{\nabla}\times E=0$. Considering a closed loop as shown in Fig. 1.1 and applying the line-integral along the loop, we obtain:

[TABLE]

having used the fact that $dl_{1}=-dl_{2}$. In vector notation, since Eq. (1.37) holds for any component $dl$, it makes the tangential component of $\vec{E}$ to be continuous along the boundary. Therefore,

[TABLE]

By means of Eq. (1.15), it follows from Eq. (1.38) that:

[TABLE]

Moreover, using Eq. (1.26), we can obtain:

[TABLE]

where $db$ is shown in Fig. 1.1. Once $db\to 0$, the right hand side of Eq. (1.40) vanishes. Hence at the boundary we have:

[TABLE]

Given a coordinate system $(\zeta,\xi,\varphi)$ with the scaling factors :

[TABLE]

we consider a body $\Pi$ with dielectric function $\varepsilon_{1}$ immersed in a medium with dielectric function $\varepsilon_{2}$. Boundary surface $\partial\Pi$, finite or infinite, is described via revolution through an angle $\varphi$ about the $z$–axis resulting in the usual azimuthal symmetry and is defined by fixing one of the coordinates, say $\zeta=\zeta_{0}$ [36]. Furthermore, we denote the inside and outside regions of $\Pi$ by $\zeta<\zeta_{0}$ and $\zeta>\zeta_{0}$ as well as their respective potentials as $\Phi_{\text{i}}$ and $\Phi_{\text{o}}$, respectively, noting that:

[TABLE]

while at the boundary $\partial\Pi$,

[TABLE]

[TABLE]

Considering the quotient $\varepsilon_{1}/\varepsilon_{2}$, we may assume with no loss of generality, $\varepsilon_{2}=1$ which corresponds the dielectric constant of vacuum, and $\varepsilon_{1}=\varepsilon$. Thus Neumann boundary condition given in Eq. (1.45) can be rewritten as:

[TABLE]

An arbitrary surface of revolution along with inside and outside potentials and dielectric constants are illustrated in Fig. 1.2.

Laplace equation

Since the non-retarded potential $\Phi$’s inside and outside parts satisfy the Laplace equation (see Eq. (1.43)), we consider general form of the Laplacian in a coordinate system $(\zeta,\xi,\varphi)$ :

[TABLE]

A coordinate system is called separable if it obeys Robertson condition (see [37], p. 655). Except the case of a nanoring in which a quasi-separation of variables holds, all considered coordinate systems in this manuscript satisfy Robertson condition and are separable. The case of a nanoring is treated separately in Chapter 3. Inspired by the azimuthal symmetry assumption and hence the $2\pi$-periodicity of the potential $\Phi$, an ansatz of the form:

[TABLE]

where $m=0,\pm 1,\pm 2,\dots$ is used. As a consequence, we obtain two Sturm-Liouville problems in $\zeta$ and $\xi$ (see [37] p. 656 and Appendix A, Eqs. (B-62)–(B-64)). Denoting by $F_{mn}(\zeta)$, $\tilde{F}_{mn}(\zeta)$, and $G_{mn}(\xi)$, $\tilde{G}_{mn}(\xi)$ the linearly independent pair of solutions to the corresponding obtained Sturm-Liouville problems in $Z(\zeta)$ and $E(\xi)$, where $n$ stands for the second separation constant, the general form of eigenfunctions for the Laplacian may be written as:

[TABLE]

For simplicity, we may assume that $n$ takes discrete values, i.e. $n=0,\pm 1,\pm 2,\dots$. However, as we shall see in Chapter 2, we will also deal with real continuous spectrum with respect to $n$. Utilizing the superposition principle, we may write the general solution as:

[TABLE]

where $\mathcal{D}_{mn}(t)$ is a general Fourier component with respect to time $t$ and they are considered as the non-retarded, time-dependence amplitudes at time $t$. Consequently, the Inside and outside potentials, $\Phi_{\text{i}}$ and $\Phi_{\text{o}}$, could be defined by considering the asymptotic behavior of the eigenfunctions $F_{mn}(\zeta)$, $\tilde{F}_{mn}(\zeta)$ and $G_{mn}(\xi)$, $\tilde{G}_{mn}(\xi)$ as the potential has to satisfy the following conditions:

[TABLE]

Without loss of generality, we may assume the sets $F_{mn}(\zeta)\,G_{mn}(\xi)$ and $\tilde{F}_{mn}(\zeta)\,G_{mn}(\xi)$ determine finite potentials for the inside and the outside of the geometry, respectively. Using Eq. (1.51) together with boundary condition given in Eq. (1.44), we may write:

[TABLE]

Harmonic Oscillator Model

In this section, we will show that amplitudes $\mathcal{D}_{mn}(t)$ undergo harmonic oscillator motion. We start by deriving two relations between surface charge density and quasi-static potential using the total charge for the whole space. Equating the two relations gives us an equation in amplitudes $\mathcal{D}_{mn}(t)$ and their second time-derivative $\ddot{\mathcal{D}}_{mn}(t)$ known as the harmonic oscillator equation of motion, see Eq. (1.69). The total charge $Q$ is defined by:

[TABLE]

where $\rho$ denotes the volume-charge density. Since the polarization charge is only confined to the surface, Eq. (1.54) may be written as:

[TABLE]

where $\sigma$ denotes the surface charge density and $d\mathcal{A}$ is the surface element. Equating Eq. (1.54) with Eq. (1.55), we have:

[TABLE]

It follows now from Eq. (1.56) that

[TABLE]

where $\delta$ denotes the Dirac delta function. In view of Poisson’s equation given in Eq. (1.27), one finds:

[TABLE]

On the other hand, using the electric polarization vector given by [23]:

[TABLE]

where $\vec{u}$ denotes the charge displacement of the system. Assuming there is no free-charge density on surface, then using Eq. (1.13), we can write:

[TABLE]

It follows from Eq. (1.21) that $\vec{\nabla}\cdot\vec{P}=-\rho$, which in view of Eq. (1.54) and Divergence Theorem gives:

[TABLE]

where $\mathbf{n}=-\mathbf{\vec{e}}_{\zeta}$. Comparing the above equation with Eq. (1.55) and using Eq. (1.59), one can find:

[TABLE]

The surface charge due to polarization is also derived as the motion of free electrons in the metal surface due to a time-harmonic external electric field, following standard descriptions of anomalous dispersion [30]. We can thus also describe this surface charge as $\sigma=-n_{0}\,e\,\vec{u}.$ The electrostatic force on the electron current per unit volume is $F=-n_{0}e\nabla\Phi_{\text{i}}.$ Putting these together through Newton’s force law, since the acceleration is the force per unit mass, we have:

[TABLE]

Differentiating both sides of Eq. (1.62) twice with respect to time gives:

[TABLE]

which together with Eq. (1.63) implies:

[TABLE]

Utilizing the identity:

[TABLE]

one may write:

[TABLE]

This equation is also provided in [8] without a detailed proof. Differentiating Eq. (1.58) twice with respect to time $t$ gives:

[TABLE]

By equating Eq. (1.67), in terms of ${\mathcal{D}}_{mn}(t)$, and Eq. (1.68), in terms of second time-derivatives of amplitudes $\ddot{\mathcal{D}}_{mn}(t)$, we obtain the harmonic oscillator equation with the frequencies $\omega_{mn}$ of the form:

[TABLE]

The importance of Eq. (1.69) is the fact that the (allowed) resonant values of the dielectric function $\varepsilon$ can also be independently calculated from this transcendental equation obtained by imposing the quasistatic boundary conditions at the bounding surfaces of the domain within which the scalar electric field satisfies the Laplace equation [38].

Classical energy

Once the scalar potential and charge density are known, we shall be able to calculate the total energy of the system. In classical mechanics, the total energy $E$, also known as total mechanical energy, of a position dependent force, in our case surface plasmons, is given by:

[TABLE]

where $V$ and $T$ denote potential and kinetic energy, respectively. This is nothing more than an expression of the conservation of mechanical energy principle [39, 40]. We shall start by calculating potential energy.

The potential energy is given by the volume integral [8]:

[TABLE]

where $\mathcal{V}$ is the volume of the geometry of interest. Using Eq. (1.57), one could write:

[TABLE]

Since inside and outside scalar potential $\Phi_{\text{i}}$ and $\Phi_{\text{o}}$ coincide on the surface, namely the first boundary condition given in Eq. (1.44), it is a direct computation using any of the relations of inside or outside scalar potentials. Utilizing the $\delta$ function, we can write:

[TABLE]

Using the identity $d\mathcal{A}=h_{\xi}h_{\varphi}\,d\xi d\varphi$, we find:

[TABLE]

The kinetic energy $T$ of the oscillations is given by the volume integral [8]:

[TABLE]

Depending on whether time-dependent amplitudes obey pure harmonic oscillator or not, we may take different steps in order to find kinetic energy. In all the cases investigated here, amplitudes satisfy the relation given in Eq. (1.69) but the case for nano-ring. We proceed by assuming the general case and we defer the discussion for nano-ring to Chapter 3. We may use the relation given in Eq. (1.63) to write:

[TABLE]

where we put:

[TABLE]

Now, using Eq. (1.69), we may write:

[TABLE]

Putting this in Eq. (1.77), we could write:

[TABLE]

Replacing this is Eq. (1.76), we obtain:

[TABLE]

Integrating both sides of the above equation with respect to time $t$ gives:

[TABLE]

Using the latter expression for $\dot{\vec{u}}$ in Eq. (1.75), we may write:

[TABLE]

where $\overline{\dot{\Phi}_{mn}}$ denotes the conjugate of ${\dot{\Phi}}_{mn}$. The fact that scalar potential is real-valued allows us to replace it with its complex conjugate. Here, Fubini’s Theorem [41] allows us to change the order of the summation and integration, assuming the integral is bounded. Using above argument as well as plasmon-bulk frequency relation given in Eq. (1.18), Eq. (1.82) could be rewritten as:

[TABLE]

Since Laplace equation is independent from time $t$, it is easy to see that if $\Phi_{\text{i}}$ satisfies the Laplace equation, so does its time-derivative $\dot{\Phi}_{\text{i}}$. We can use the following general identity:

[TABLE]

to write the kinetic energy as:

[TABLE]

The volume integral in Eq. (1.85), using Divergence Theorem could be transformed to a surface integral through the following relation:

[TABLE]

Using the identity given in Eq. (1.66), we find the final expression for kinetic energy as:

[TABLE]

From Eqs. (1.74) and (1.87), the classical energy is obtained by using Eq. (1.70). Depending on amplitudes being real or complex valued, they could be written in the following form:

[TABLE]

with time-dependent operators $d_{mn}$ directly proportional to $e^{i\omega_{mn}t}$, where $d^{*}_{mn}$ is its respective complex conjugate. The purely imaginary case could be assumed by multiplying the real case of Eq. (1.88) by $i$. For the simplicity of the discussion, we may consider the case for which amplitudes are real-valued. The goal is now to determine the values of $\mathbf{\alpha_{mn}}$. First find the time-derivative of amplitudes $\mathcal{D}_{mn}(t)$ using Eq. (1.88), i.e.:

[TABLE]

Using Eqs. (1.88) and (1.89) to write the total energy $E$ in terms of $d_{mn}$ and $d^{*}_{mn}$ and comparing the new expression with Hamiltonian operator for the quantum mechanical harmonic oscillator given as [13]:

[TABLE]

where $\hbar$ is the Planck’s constant . The following substitution allows us to find the values for $\mathbf{\delta_{mn}}$:

[TABLE]

The operators $\left({\hat{d}}_{mn}\,,\,\hat{d}^{\dagger}_{mn}\right)$ are interpreted as the non-commuting boson creation and annihilation operators, respectively. Calculations of the interactions between surface plasmons and other excitations are often time complicated, yet have been simplified by quantizing the fields (second quantization) [10]. The total classic energy was calculated to help us finding coefficients in Hamiltonian operator. In the next section, we will use the Hamiltonian operator to calculate matrix element and probability amplitudes for interactions.

Quantum Field & Second Quantization

Among the many quantum numbers that distinguish particles, we can count their spin [42, 43, 44]. There are two main families of particles: bosons, with integer spin and fermions, with half-integer spin. Photons, gluons and vector bosons are examples of bosonic particles and electrons, neutrinos, and quarks are classified as fermionic particles.

There have been two main approaches to quantum field theory and on how to study the interaction between particles and surfaces. One approach is due to Feynman: His celebrated path-integral approach based on Feynman’s diagram. Two detailed accounts can be found in [42, 45]. This approach does not rely on operators and the theory behind the operator in Hilbert spaces. Another approach, by means of operator theory, can be traced back to Heisenberg, Born, Jordan, Dirac, Pauli, and von Neumann. This method, more formal than the previous one, is more difficult to apply to practical computations, due to the nonlinearities of the interactions. In [42], the crucial aspect of a quantum field is that it can be treated as a system with infinitely many quantum particles that can be created and annihilated. The main purpose of second quantization is to describe quantum particles, either bosons or fermions, in terms of operators. The number of particles is assumed to be infinite. We first set our full system as the combined surface plasmon and photon field. Then we consider the ground state or zero-population for both plasmon and photon as:

[TABLE]

A state consisting of precisely $N$ plasmon particles is shown by the notation:

[TABLE]

where $\nu_{mn}$’s are the number of plamons in the state $(mn)$. Similarly, the state consisting of exactly $\tilde{N}$ photons can be represented by:

[TABLE]

where $\nu_{sq}$’s show the number of photons with wavevector $s$ and polarization $q=1,2$. Polarization vectors are chosen based on the fact that they are perpendicular to each other and to the direction of travel. It is customary to represent them by s-polarization and p-polarization vectors. We shall discuss it in more details in Chapter 2. Our main purpose concerns the investigation of interacting quantum fields in rigorous mathematical terms. We shall use the fact that the operators describing photons only act on and hence change the photon state. The same situation holds for the plasmon operators while acting over the full system. The states being orthonormal results in the Kronecker delta function as:

[TABLE]

for some arbitrary states $1,\cdots,k$. Considering these facts, we are now able to write a general state for the combined plasmon and photon field in the form:

[TABLE]

Next, we shall discuss how creation and annihilation operators, $\hat{d}_{mn}$ and $\hat{d}^{\dagger}_{mn}$ acting on an arbitrary state, including ground state of the field. A prototype of commutation relations for these operators when act on bosonic particles is given in [42, 43]:

[TABLE]

for all plasmon states $mn$ and $m^{\prime}n^{\prime}$ and

[TABLE]

with $\hat{a}$ and $\hat{a}^{*}$, creation and annihilation operators for photon and for all photon states $sq$ and $s^{\prime}q^{\prime}$.

Quantization of the Electromagnetic field

In view of Maxwell’s equations (recall Eqs. (1.21)–(1.24)), both electric and magnetic fields can be expressed in terms of vector potential $\mathbf{A}$ as:

[TABLE]

In case of absence of any sources of the field, it is possible to choose the Coulomb gauge condition as:

[TABLE]

This condition is also known as transversality condition since the vector potential under this condition is purely transverse and could be expressed as [29, 46]:

[TABLE]

where $\hbar\omega_{s}$ is the photon energy, the polarization vector $\mathbf{\hat{e}}_{q}$ is perpendicular to $s$ for both values of polarization on index $q=1,2$ and $a_{sq}(t)$ and its conjugate $a^{*}_{sq}(t)$ are the photon operators, such that:

[TABLE]

therefore:

[TABLE]

which is considered as the equations for motion of the field for all $s$. To write the photon field as a sum with discrete momentum eigenstates as opposed to the continuous representation, we have to make some specific choices. We consider our electromagnetic field to be confined to a volume $\mathcal{V},$ which is normally taken to be represented by a cube over which we impose periodic boundary conditions. Since the electromagnetic energy confined to this volume is independent of the shape of the volume [47, 48, 44], we take as our quantization volume as a box $\mathcal{B}$, with side $L$. The transition to the discrete sum then follows the lattice strategy in quantum electrodynamics [42] which we summarize as follows:

Periodic free field: Classical free field is represented in terms of finite Fourier coefficients. The first step is to confine the photons and electrons within the above-mentioned box $\mathcal{B}$. 2. 2.

We impose periodic boundary conditions by assuming that fields have the period $L$ with respect to the position of particles. Observe that this new condition on boundary by no means contradicts the boundary conditions we first defined in Eqs. (1.44) and (1.46). 3. 3.

We introduce a momentum space which yields finite Fourier series. 4. 4.

The Fourier coefficients of the finite Fourier series are replaced by creation and annihilation operators we introduced in this section. This is the bridge from classical mechanics to quantum mechanics. 5. 5.

Using finite number of creation and annihilation operators, the matrix elements could be formed. 6. 6.

Finally, we let the finite box $\mathcal{B}$ tend to infinity, in other words, we undertake the limit of $L\to\infty$. We let the period $L$ become infinite.

In practice, we perform the transition:

[TABLE]

to write the discretized vector potential as:

[TABLE]

where the photon operators have been replaced by $\hat{a}_{sq}\,,\hat{a}^{\dagger}_{sq}$, annihilation and creation operators for photon. Following Eq. (1.98), we replace Dirac delta function $\delta(s-s^{\prime})$ by the discrete Kronecker delta function $\delta_{ss^{\prime}}$ and obtain new the commutation relations as:

[TABLE]

Given the explicit form of vector potential $\mathbf{A}$, the goal is to calculate the probability amplitude and interaction Hamiltonian.

Probability amplitude and matrix elements

The probability amplitude that a surface plasmon in a given state, defined by $m^{\prime},n^{\prime}$ will interact with a photon, described by $s^{\prime},q^{\prime}$ is expressed as:

[TABLE]

where $H_{\text{int}}$ denotes the interaction Hamiltonian. We shall call:

[TABLE]

as the interaction matrix element. We could define various matrix elements of the interaction Hamiltonian based on the type of interaction. A collection of matrix elements are given in [8] and we list them here:

Direct scattering (elastic, including the Thomson limit):

[TABLE]

where index $f$ shows final state while index $i$ indicates initial state and $H^{(0)}_{RP}$ denotes the direct scattering interaction Hamiltonian, given by:

[TABLE]

where $\zeta_{0}$ is the shape parameter. This Hamiltonian is also known as the zeroth-term of Hamiltonian. 2. 2.

Emission (radiative decay of surface plasmons):

[TABLE]

where

[TABLE]

The term $H^{(1)}_{RP}$ represents the interaction of one photon and one plasmon, which can be used to predicct the creation of a plasmon by a photon or the decay of a plasmon into a photon. 3. 3.

Absorption:

[TABLE]

which is the Hermitian dual to Eq. (1.112), since it represents the inverse process. 4. 4.

Total elastic scattering: Here $\kappa$ is the plasmon damping factor and the matrix element is obtained as shown below:

[TABLE] 5. 5.

Inelastic scattering:

[TABLE]

where

[TABLE]

This term involves two photons and one plasmon and can describe, for example, the inelastic scattering of a photon in creating a plasmon.

Interaction Hamiltonian

Here, to describe the plasmon-photon interaction, we resort to the hydrodynamical formulation of the electron gas by Crowell and Ritchie [6] as:

[TABLE]

where $\Phi(\bf r,t)$, $\Psi(\bf r,t)$ and $n(\bf r,t)$ are the electric potential, velocity potential, and electronic density, respectively, in the electron gas, while $n_{0}$ is the electronic density in the undisturbed state of the electron gas and $\xi$ is the propagation speed of the disturbance through electron gas. By linearizing these equations, employing perturbation theory, Ritchie obtained the first order interaction Hamiltonian given as:

[TABLE]

where $\mathbf{J}$ denotes the current density. This interaction has been used previously to describe the emission of photons via plasmon decay on finite surfaces of an oblate spheroid modeling silver nanoparticles evaporated on a dielectric substrate [7]. Here, we will apply this Hamiltonian to modeling the creation of a plasmon on the surface by a photon or the decay of a plasmon and emission of a photon [8]. This application requires the explicit determination of the current density operator. Let us first consider the current given by:

[TABLE]

We define:

[TABLE]

in which we make the replacement of $\dot{D}_{mn}(t)$ with its equivalent given in Eq. (1.89), and use Eq.  (1.81) to write:

[TABLE]

Thus one can write:

[TABLE]

Using the vector identity:

[TABLE]

Imposing the Coulomb gauge condition given in Eq. (1.101), and the divergence theorem, recalling that the current only exists on the surface of the body, we can find:

[TABLE]

which gives the interaction Hamiltonian as:

[TABLE]

Once $H_{\text{em}}$ is known, can calculate the emission matrix element given in Eq. (1.112). Using Eqs. (1.106) and (1.123), we write:

[TABLE]

Using Fubini’s Theorem in Eq. (1.129), we obtain:

[TABLE]

We use this equation in the relevant transition matrix element for photon emission, given by Eq. (1.112) where $\hat{a}_{q_{f}s_{f}}$ and $\hat{d}^{\dagger}_{m_{i}n_{i}}$ describe the final photon state and the initial plasmon state, respectively. It is noteworthy to emphasize again that photon operators act only on the photon state and plasmon operators act only on the plasmon state [44, 49, 42, 43, 50]. The question of which particular operators, either creation or annihilation, give rise to a non vanishing matrix element depends entirely on initial and final states. The following commutation relation is our main reference to calculate the probability amplitude and therefore, the matrix element:

[TABLE]

Using all these commutation relations, we could discard the summations in Eq. (1.130) and can proceed to calculate the matrix element. Having calculated the emission matrix element, we have all the ingredients to discuss and calculate the radiative decay in the following section.

Radiative Decay Rate

In the study of the phenomenon of radiation in the quantum domain, the starting point is to talk about vector potential $\mathbf{A}$ [31]. Imposing guage condition given in Eq. (1.101), the fields, either electric or magnetic, which are derived from vector potentials are called radiation fields. This term is often used to describe, improperly, the vector potential itself. Fermi’s formalism based on emission interaction Hamiltonian, derived in the previous section, is called the radiation or Coulomb gauge method. A theory once quantized, provides transparent details about the process of bosonic particles, here photons, being emitted, absorbed and scattered. Let us now specify our discussion with respect to emission process. We consider a solid angle element $d\Omega$ through which a photon is emitted. The number of allowed states in an energy interval can be written as:

[TABLE]

with $s$ pointing into the solid angle $d\Omega$, This equation also is known as the the density of a single photon [31]. Fermi’s Golden rule [23] gives the transition probability in the form:

[TABLE]

where $\mathcal{M}_{em}$ represents the matrix element given in Eq. (1.112) and $\rho_{f}$ is the density of final states. The total radiative decay rate, $\gamma_{mn}$, is obtained by summing over final photon states as:

[TABLE]

If we let the quantization volume, $\mathcal{V}$, become large enough so that the sum over $s$ changes into an integral, then we could write:

[TABLE]

We could now write, using Eq. (1.132):

[TABLE]

Lastly, using Eq. (1.133), we find:

[TABLE]

which gives the radiative decay rate per solid angle for two polarization vectors $q=1,2$. Eq.  (1.137), as one can easily see, becomes independent from the quantization volume $\mathcal{V}$, also known as normalized volume. In calculating Eq. (1.137), we also require to expand the plane wave $e^{\pm i\,s\mathbf{\cdot r}},$ using the transition expressions from the given coordinate system, here $(\zeta,\xi,\varphi)$, to Cartesian coordinates. Often time, the resulting wave functions are difficult to calculate. Depending on the case, we either avoid direct calculations of wave functions using rotating and scaling the wave vector $s$ with algebraic substitutions, or we used numerical calculations to compute integrals due to the absence of analytical solutions.

In the next chapter, we apply the whole procedure given in the current chapter to two infinite geometries as solid paraboloid and hyperboloid and one finite geometry, prolate spheroid. The results for these two cases have been provided in [51, 52].

CHAPTER 2RADIATIVE DECAY ON PARABOLOIDAL, HYPERBOLOIDAL, & PROLATE SPHEROIDAL SURFACES

11footnotetext: Portions of this chapter were reprinted from: Bagherian, M. and Kouchekian, S. and Rothstein, I. and Passian, A., Quantization of surface charge density on hyperboloidal and paraboloidal domains with application to plasmon decay rate on nanoprobes, 125413, vol. 98, Sep 2018, with permission from Phys. Rev. B

Permission is included in Appendix C.

Introduction

Materials confined to microscopic elongated "probe-like" domains, in addition to having been tremendously enabling in various forms of scanning probe microscopy (SPM), hold great potential for emerging applications in fields such as quantum sensing [1, 2, 3, 4]. Recent interest in the properties of such tip-shaped material domains is reflected in the demonstration of laser pulse induced electron emission from a gold solid tip under grating coupled plasmon excitation [53, 54, 55]. These applications make use of the excitation and resonant properties of surface modes on bounding surfaces and interfaces of metallic, dielectric, and metallo-dielectric domains that take the form of a tip [27]. Examples of systems that use a probe tip include scanning tunneling microscope, photon scanning tunneling microscope, apertured and apertureless nearfield scanning optical microscope, nanoantennas, and processes such as tip-enhanced spectroscopy and lithography. The intriguing excitations, typically studied near the tip apex, are expected to receive contributions not only from stationary modes, such as those occurring at the surfaces of finite nanoparticles, but also from non-stationary modes propagating at the infinite interfaces. The theoretical and modeling tools for investigating the response of these systems and their dependence specifically upon the geometric characteristics have been indispensable in the development of these applications. In particular, analytical techniques that lend themselves to provide complete or partial information on the system are often regarded as necessary not only for obtaining the system response (e.g., energy distribution in the nearfield of the nanoparticles), but also for elucidating the inner working of the systems (e.g., the contributing eigenstates and eigenvalues). Calculation of geometric and material dependencies of surface mode excitation, decay and scattering on the bounding surfaces of nanoscale domains are both instructive and necessary for better design and fabrication.

Here, we investigate the radiative decay rate of plasmons by quantizing the surface modes engendered on the surface of a metallic probe modeled as one sheet of a two-sheeted hyperboloid of revolution, shown in Fig. 2.1. This geometry offers an elegant adaptability not only for the description of the local curvature of a fabricated probe but also for the modeling of nearly planar interfaces [27]. In addition, it has the property that the hyperboloidal domain translates along its symmetry axis when changing the opening angle $\theta_{0}$, that is, smaller $\mu_{0}=\cos\theta_{0},$ yields smaller gap $z_{\text{min}}$, the apex distance to the origin $o$ in Fig. 2.1. To provide a basis for comparison, we quantize the surface charge density oscillations on the useful system of a paraboloid of revolution, which offers a similar apex morphology but a different asymptotic behavior away from the apex. Importantly, the apex and off-apex curvature of a paraboloid of revolution presents a more natural topology for comparison of its spectral and scattering properties with that of a finite body of similar curvature, e.g., a spheroidal domain. Therefore, for the sake of validation, we extend our investigations to study the radiative decay of plasmons excited on a prolate spheroid, which owing to its finite volume presents a more tangible system.

Following the steps given in Chapter 1, we organize this chapter as follows: In section 2.2, we treat the paraboloidal plasmons. Here, within the quasi-static framework, representing the material domain with a frequency-dependent dielectric function, we derive the nonretarded plasmon dispersion relations, eigenmodes, and fields. From the classical energy of paraboloidal charge density oscillations, we then derive the Hamiltonian of the system. We then proceed to quantize the plasmon field and, employing an interaction Hamiltonian derived from the first order perturbation theory within the hydrodynamic model of an electron gas [10, 6], obtain an analytical expression for the radiative decay rate of the plasmons. Having established the full treatment of the paraboloidal system, in section 2.3, we proceed to investigate the quantized charge density oscillations on the surface of one sheet of a two-sheeted hyperboloidal of revolution. In both these sections the use of non-retarded potentials and dispersion relations is justified due to the sub-wavelength dimension of the tip. We would also discuss our findings and compare both the paraboloidal and hyperboloidal results. An interesting comparison of the paraboloidal domain can be made with respect to the surface modes and radiation patterns of a prolate spheroid, a finite geometric domain highly suitable for modeling of nanoparticles such as a quantum dot. In specific cases, we further validate the results using computational techniques to obtain the lower energy eigenmodes and farfield radiation patterns. Quantum calculations that take into account the geometric effects of the bounding surfaces of the material domains are important to corroborate experimental observations in nanophysics. Here, our goal is to calculate the probability amplitude that a surface plasmon in a given initial state, engendered near the apex region of a probe-shaped material domain, will emit a photon into a given final state. To model a tip-like domain, we will employ three specific cases of surfaces of revolution: a single sheet of a two-sheeted hyperboloid, a paraboloid, and a prolate spheroid. Fig. 2.1 shows an example of a hyperboloidal domain and its relation with solid angle $d\Omega$ (see Eq. (1.137)).

Paraboloidal Surfaces

In the following, we will first treat a paraboloid, for which we begin by seeking pertinent classical quantities. The paraboloidal domain allows the parametric study of the various scattering processes as functions of the local curvature without displacement of the domain. The quasi-static solution of the electric scalar potential near a paraboloidal domain has been reported in [56].

The parabaloidal coordinates $(\xi,\eta,\varphi)$, are related to the rectangular coordinates by

[TABLE]

with the corresponding scale factors

[TABLE]

where $0\leq\varphi<2\pi$ denotes the usual azimuthal angle, $a>0$ is a dimensionless constant to be determined later, and the two coordinates $\eta,\xi\geq 0$ are such that the surfaces of constant $\eta>0$ and $\xi>0$ describe upward and downward paraboloids of revolution about the $z$-axis, respectively.

Nonretarded Potential & Dispersion Relations

We consider a vacuum-bounded solid paraboloid of revolution defined by $\eta=\eta_{0}$, via the coordinates $(\xi,\eta,\varphi)$ given in Eqs. (2.1)–(2.4). Denoting the frequency $\omega$ dependent dielectric function of the paraboloidal material domain with $\varepsilon(\omega)$, and we consider the outside medium to be vacuum with dielectric constant $1$ (see explanation after Eq. (1.46)). Following the procedure outlined through Eqs. (1.48)–(1.53), assuming the ansatz $\Phi(\xi,\eta,\varphi)=F(\lambda\xi)G(\lambda\eta)e^{im\varphi}$ results in two Sturm-Liouville problems, where $F$ satisfies the Bessel equation:

[TABLE]

and $G$ satisfies the modified Bessel equation:

[TABLE]

The solutions to Eq. (2.5) are the Bessel functions of the first and second kind $J_{m}(\lambda\xi)$ and $Y_{m}(\lambda\xi)$, whereas solutions to Eq. (2.6) are the modified Bessel Functions $I_{m}(\lambda\eta)$ and $K_{m}(\lambda\eta)$. These eigenfunctions form a basis for the solution of corresponding equation (see Appendix B) [57, 58, 59]. Considering the two resulting Sturm-Liouville problems Eqs. (2.5) and (2.6), with unbounded domain $\eta,\xi\in[0,\infty)$, lead to a continuous spectrum of real eigenvalues and eigenfunctions [60] in terms of Bessel and modified Bessel functions given by Eqs. (10.3.64) and (10.3.65) [37], the reader is also encouraged to read [61, 62]. Using the fact that the potential is bounded on the $z$-axis and vanishes as $r\to\infty$ (see Eq. (1.51)), together with the asymptotic behavior of the Bessel functions given in Appendix B, Eq. (B-40), we denote the potentials with $\Phi_{\text{i}}$ and $\Phi_{\text{o}}$, for the interior and the exterior domains, respectively, and utilize the Heaviside function $\Theta$, given in Eq. (A-16), with the half-maximum convention $\Theta(0)=\frac{1}{2}$ to write the total potential as:

[TABLE]

or explicitly:

[TABLE]

where $m=0,1,2,\cdots,$ and $p=0,1,$ while $A_{m\lambda p}(t)$ are the time $t$ dependent amplitudes to be determined by the boundary conditions, and $\{S^{p}_{m}(\varphi)\}$ indicate the azimuthal symmetry of the eigenmodes, explicitly:

[TABLE]

satisfying the orthogonality relation given:

[TABLE]

where

[TABLE]

Fig. 2.2 shows the spatial distribution of the lowest lying eigenmodes of the quasi-static electric potential for the paraboloid.

The Laplacian of $\Phi$ in paraboloidal coordinates is given by (see Eq. (1.47)):

[TABLE]

Since the Heaviside function depends on the coordinate $\eta,$ we only need to consider $\frac{\partial\Phi}{\partial\eta}$ and $\frac{\partial^{2}\Phi}{\partial\eta^{2}}$; The partial derivatives of $\Phi$ with respect to the other coordinates, $\xi$ and $\varphi$, can be trivially expressed in terms the corresponding partials of $\Phi_{\text{i}}$ and $\Phi_{\text{o}},$ respectively. We can write:

[TABLE]

where the last two terms of the above expression cancel out due to the fact that $\Phi_{\text{i}}=\Phi_{\text{o}}$ at the boundary $\eta=\eta_{0}.$ Thus,

[TABLE]

and as a result:

[TABLE]

Substituting Eqs. (2.14) and (2.15), together with partials derivatives of $\Phi$ with respect to $\xi$ and $\varphi,$ into the Laplacian Eq. (2.12) gives

[TABLE]

Treating the above expression in the sense of a distribution in $\eta$ only and noting that $\vec{\nabla}^{2}\Phi_{\text{i}}$ and $\vec{\nabla}^{2}\Phi_{\text{o}}$ vanish for $\eta<\eta_{0}$ and $\eta>\eta_{0},$ respectively, we obtain the following identity for the Laplacian of $\Phi$:

[TABLE]

where $h_{\eta}$, given by Eq. (2.4), is a scale factor of the paraboloidal system.

Surface charge density

From Eq. (2.17) and using the relation given in Eq. (1.68), one can solve for the charge density to obtain:

[TABLE]

The expression on the right-hand side of Eq. (2.18) can be easily calculated from Eq. (2.8) as:

[TABLE]

where $\mathcal{W}_{m}(\cdot)$ denotes the Wronskian given by:

[TABLE]

In view of the identity $\mathcal{W}_{m}(z)=-\frac{1}{z}$ ($z\neq 0$) for the modified Bessel functions [57], and from Eqs. (2.18) and (2.20) we obtain:

[TABLE]

From the definition of the polarization, the surface charge $\sigma$ can be written as relation provided in Eq. (1.65) and it follows that:

[TABLE]

Differentiating the charge density in Eq. (2.21) twice with respect to time $t$ and equating it with Eq. (2.22), it follows from the orthogonality of system $\{S_{m}^{p}\}_{m,p}$ given in Eq. (2.10), that for each fixed $m$ and $p$:

[TABLE]

where

[TABLE]

Utilizing the orthogonality relation for Bessel functions given by (see Eq. (11.59) [63]):

[TABLE]

for $m\geq-1,$ and $\lambda,\lambda^{\prime}>0,$ it follows from Eq. (2.23) that the amplitudes $A_{m\lambda p}(t)$ undergo harmonic oscillations at continuous frequencies $\omega_{m\lambda}$ given by Eq. (2.24), that is:

[TABLE]

The harmonic behavior of the field amplitudes will play an essential role in the possibility of analytically calculating the energy of the paraboloidal charge system prior to quantization.

Moreover, as it has already been discussed in Section 1.4, the obtained $\omega_{m\lambda}$ in Eq. (2.24) satisfies the Drude model given in Eq. (1.20), namely:

[TABLE]

The resonance values of the dielectric function $\varepsilon$ given in Eq. (2.27) for a few modes are shown in Fig. 2.3.

Classical Energy

Having obtained closed form expressions for the potential and induced surface charge density, we are now in the position to calculate the potential energy $V$ due to the polarization charge distribution $\rho$ given by Eq. (1.71). Since $\rho$ is only confined to the paraboloidal surface $\eta=\eta_{0}$ and vanishes elsewhere in the space, one may integrate over an infinitesimal thin cover across the boundary: $\eta_{0}-\epsilon\leq\eta\leq\eta_{0}+\epsilon.$ Letting $\epsilon\to 0^{+}$ and making the observation $\rho h_{\eta}=\delta(\eta-\eta_{0})\sigma$, the potential energy can be expressed as a surface integral given in Eq. (1.74). From Eq. (2.8), the potential \Phi_{\text{i}}\bigg{|}_{\eta=\eta_{0}} is given by:

[TABLE]

Substituting Eq.  (2.28) and the expressions for $\sigma,$ $h_{\xi},$ and $h_{\varphi}$ from Eqs. (2.23) and (2.4) into Eq. (1.74) gives

[TABLE]

From the orthogonality relations for the trigonometric system $\{S^{p}_{m}(\varphi)\}$ given in Eq. (2.10), we can write Eq. (2.29) as:

[TABLE]

Finally, in view of the orthogonality relation for Bessel functions given by Eq. (2.25), we obtain we obtain the potential:

[TABLE]

where $\hat{\delta}_{mp}$ is given in Eq. (2.11).

Note that in the above calculations, we have changed the order of integration to simplify the obtained expressions. This is done in view of Fubini’s theorem based on the absolute integrability of the above expressions. Our argument is based on the asymptotic behavior of Bessel functions (see Appendix B, Eqs. (B-36)–(B-38)) and the fact that the potential, and thus Eq. (2.8), is finite in the entire space. The same argument is also used in the other cases using the asymptotic behavior of the corresponding eigenfunctions.

We will now seek the kinetic energy $T$ of the paraboloidal charge system. Employing the charge displacement vector, kinetic energy is given in Eq. (1.75). To obtain an expansion for $\dot{\vec{u}},$ we note that in the expression for the potential $\Phi_{\text{i}}$ given by Eq. (2.8), one may use the harmonic oscillator equation Eq. (2.26) for the amplitudes to get:

[TABLE]

As a result, integrating the charge displacement vector with respect to time, we obtain: $m_{e}\dot{\vec{u}}=-e\vec{\nabla}\dot{\Psi},$ where

[TABLE]

With this expression in the integral for kinetic energy, making use of the Gauss theorem [35] and orthogonality relations for $\{S^{p}_{m}\}_{m,p}$ and Bessel functions, given in Eqs. (2.10) and (2.25), we thus calculate Kinetic energy by rewriting Eq. (1.87) for the case of paraboloid as:

[TABLE]

Using the definition of $\Psi$ given by Eq. (2.33) and the expressions for the scale factors Eq. (2.4), the right-hand side of Eq. (2.34) can be written as:

[TABLE]

Invoking the orthogonality relations, given in Eqs. (2.10) and (2.25), to the Eq. (2.35) by using a similar argument as the one given for the case of the potential, it follows that

[TABLE]

Now in view of Eq. (2.24), we have :

[TABLE]

This substitution in Eq. (2.36) gives the kinetic energy as:

[TABLE]

In view of Eqs. (2.31) and (2.37), we find that the total classical energy $E$ of the paraboloidal system takes the following form:

[TABLE]

which is the total energy of the surface plasmon field, suitable for quantization.

Field Quantization & Interaction of Plasmons with Photons

To obtain the quantized plasmon field, the expression for the classical field $E$ needs to be rewritten in a suitable form. We begin by noting that since the potential together with the paraboloidal harmonics $I_{m}(\lambda\eta),$ $J_{m}(\lambda\xi),$ $K_{m}(\lambda\eta),$ and $S^{p}_{m}(\varphi)$ are all real-valued, we have that the amplitudes $A_{m\lambda p}(t)$ are real-valued and satisfy the harmonic oscillator equation Eq. (2.26). In analogous to Eq. (1.88), to pursue the conversion of $E$ into the Hamiltonian operator for a scalar boson field (plasmons are spinless quasi-particles) and we write:

[TABLE]

where $a_{m\lambda p}$ are complex-valued time dependent functions proportional to $e^{-i\omega_{m\lambda}t},$ and $\alpha_{mp\lambda}$ are some real-valued constants to be determined later. The time derivative of Eq. (2.39) can now be expressed as:

[TABLE]

Substituting Eqs. (2.40) and (2.39) in the expression for classical energy given in Eq. (2.38), we may write:

[TABLE]

Performing field quantization (see Section 1.6), we replace the amplitudes $a_{m\lambda p}(t)$ and $a^{*}_{m\lambda p}(t)$ with operator valued distributions $\hat{a}_{m\lambda p}$ and its conjugate $\hat{a}^{\dagger}_{m\lambda p},$ and following Eq. (1.131), we note the commutation relations:

[TABLE]

owing to the continuous spectrum of the eigenvalues $\lambda$ that originates from the infinite axial dimension of the paraboloid. Taking the normal ordered expansion of the non-commuting boson creation $\hat{a}^{\dagger}_{m\lambda p}$ and annihilation $\hat{a}_{m\lambda p}$ operators, a comparison with the normal ordered expression of the Hamiltonian operator for a scalar boson field [8], yields the Hamiltonian:

[TABLE]

upon choosing:

[TABLE]

For the interaction of the plasmon system with a photon field with Hamiltonian $\mathcal{H}_{p}$, we require the plasmon-photon interaction Hamiltonian $\mathcal{H}_{i}$, which is required to determine the radiative decay rate of surface plasmons excited on the paraboloidal surface.

Our full system, plasmon field + photon field is described by the tensor product of the Fock spaces for the two constituent fields.

This interaction has been used previously to describe the emission of photons via plasmon decay on finite surfaces of an oblate spheroid modeling silver nanoparticles vacuum evaporated on a dielectric substrate [7]. Here, we will apply this Hamiltonian to modeling the creation of a plasmon on the surface of the paraboloid by a photon or the decay of a paraboloidal plasmon and emission of a photon [8]. This application requires the explicit determination of the current density operator, which in light of the displacement operator: $\dot{\vec{u}}=-(e/m_{e})\vec{\nabla}\dot{\Psi},$ can be written as $\mathbf{J}=-(n_{0}e^{2}/m_{e})\vec{\nabla}\dot{\Psi}$ (see Eqs. (1.122)–(1.124)).

To write the photon field as a sum with discrete momentum eigenstates as opposed to the continuous representation, we consider our electromagnetic field to be confined in a volume $\mathcal{V},$ which is normally taken to be represented by a cube. Taking the electromagnetic energy confined to a volume to be independent of the shape of the volume [47, 48], we take as our quantization volume a paraboloidal structure given by $\xi=\eta=L$ with $L\gg\eta_{0}$. The volume is then found to be $\frac{\pi}{2}L^{6}.$ In the paraboloidal coordinates, $\xi$ and $\eta$ have dimensions of length1/2 and so our volume has dimension $L^{3}$. Hence, the discretized vector potential given in Eq. (1.106) for paraboloidal case may be written as:

[TABLE]

with $\mathbf{s}$ being the three-dimensional photon wavevector $\omega_{s}=\mathbf{s}c$, $\mathbf{\hat{e}}_{j}$ the polarization vector being perpendicular to $\mathbf{s}$ for both values of $j$ and $\hat{c}^{\dagger}_{\mathbf{s}j}$ and $\hat{c}_{\mathbf{s}j}$, the photon creation and annihilation time-dependent operators, following Eq. (1.107), satisfying the commutation relations:

[TABLE]

For the physical quantities of interest in our work, we find that the quantization volume $\mathcal{V}$ cancels out in the calculations. The photon field Hamiltonian corresponding to the vector potential above is

[TABLE]

Taking $\mathbf{A}$ to be in the radiation gauge with both $\Phi=0$ and $\vec{\nabla}\cdot\mathbf{A}=0$ (see Section 1.6), we note

[TABLE]

since the current is confined to the surface of the paraboloid, we have

[TABLE]

Differentiating $\Psi$ given by Eq. (2.33) with respect to time and replacing $\dot{A}_{m\lambda p}(t)$ using Eq. (2.40), we can now write the interaction Hamiltonian as:

[TABLE]

where we have taken the integral in $\xi$ to be from 0 to $L$, since the space is bounded by $\xi=L$.

Emission of Light & Radiative Decay

Having obtained the explicit form of the interaction Hamiltonian, we now aim to calculate the probability amplitude that a surface plasmon in a given initial state, defined by $m,\lambda,p$ will emit a $j$-polarized photon in $\mathbf{s}$ state with momentum $\hbar\mathbf{s}$. Following Section 1.7, the total radiative decay rate for a given initial state, $\gamma_{m\lambda p}$, is obtained by summing over final photon states, $\gamma_{m\lambda p}=\sum_{s}\sum_{j=1,2}\gamma^{(j\mathbf{s})}_{m\lambda p},$ where by the Fermi Golden rule, the transition probabilities are:

[TABLE]

with $\mathcal{M}^{(j\mathbf{s})}_{m\lambda p}=|\bra{0}\hat{c}_{j_{f}\textbf{s}_{f}}\mathcal{H}_{i}~{}\hat{a}^{\dagger}_{m_{i}\lambda_{i}p_{i}}\ket{0}|$ denoting the emission matrix elements, $f$ the final state, and $i$ the initial state, by means of the commutation properties:

[TABLE]

the non-vanishing terms yield:

[TABLE]

where we have omitted $i$ and $f$ for convenience and we have set:

[TABLE]

For an arbitrary photon wavevector of the form:

[TABLE]

we may consider the polarization vectors $\hat{\mathbf{e}}_{j}$ for s-polarization and p-polarization, being perpendicular and parallel to $\hat{z}\mathbf{s}$-plane, respectively, noting that $\mathbf{s}\cdot\hat{\mathbf{e}}_{j}=0$. Thus, we consider the polarization vectors $\hat{\mathbf{e}}_{1}=(0,1,0)$ and $\hat{\mathbf{e}}_{2}=(\sin\psi,0,-\cos\psi).$ Under these polarization conditions, Eq. (2.55) leads to two different integrals:

[TABLE]

and

[TABLE]

where $E_{m\lambda}(\eta,\xi)=e^{-i\mathcal{A}}$ with $\mathcal{A}$ given by:

[TABLE]

Thus, depending on the choice of the polarization vector $\hat{\mathbf{e}}_{j}$, each integral contributes to a polarization state. In other words, $\mathcal{I}^{(1)}_{m\lambda}$ represents the s-polarization and $\mathcal{I}^{(2)}_{m\lambda}$ corresponds to the p-polarization. For a photon emitted pointing to the solid angle $d\Omega$ as depicted in Fig. 2.1, final photon states result in a continuous energy spectrum. Hence, we let the quantization volume $\mathcal{V}$ go to infinity so that we recover a continuum set of $\mathbf{s}$-states, that is:

[TABLE]

where $\mathcal{V}/(2\pi)^{3}$ is the density of photon states per polarization. In view of transition probability given in Eq. (2.51), the radiative decay rate per unit solid angle may be written as:

[TABLE]

in which using $\omega_{s}=sc$, we obtain the final expression:

[TABLE]

which is observed to be independent of the volume $\mathcal{V}$ as $|\mathcal{M}^{(j\mathbf{s})}_{m\lambda p}|^{2}$ is proportional to $\mathcal{V}^{-1}$ in . Thus, we can take our quantization volume to be infinite and convert the sum over the photon states to an integral. Hence, the paraboloidal decay rate per solid angle is given by:

[TABLE]

where $\mathcal{I}^{(1)}_{m\lambda}$ and $\mathcal{I}^{(2)}_{m\lambda}$ are given by Eqs. (2.57) and (2.58), respectively, and must be computed numerically.

Hyperboloidal Surfaces

Similarly to the paraboloidal case, the hyperboloidal domain, shown in Fig. 2.1, is highly relevant for modeling of electronic and photonic response of probe-like nanostructures [PPRB]. Under the same assumptions as the paraboloidal case, we proceed as follow: An arbitrary point in the Cartesian space can be expressed in terms of the prolate spheroidal coordinates $(\eta,\mu,\varphi)$ [57, 64] as:

[TABLE]

in the domain defined by :

[TABLE]

or alternatively,

[TABLE]

with the corresponding scale factors:

[TABLE]

where $z_{0}>0$ is an overall scale factor, $1\leq\eta<\infty,$ $-1\leq\mu\leq 1,$ and $0\leq\varphi<2\pi$ denote the usual azimuthal angle. The surfaces of constant $\mu$ define hyperboloids of revolution about the $z$-axis, while the surfaces of constant $\eta$ correspond to prolate spheroids. It is customary to use the substitutions $\sinh\zeta=\eta$ and $\sin\theta=\mu$ with $1\leq\eta<\infty$, and $-1\leq\mu\leq 1;$ however, in this presentation the surfaces of constant coordinates are no longer geometrically meaningful. The radius vector is given by:

[TABLE]

where $(\hat{i},\hat{j},\hat{k})$ are the Cartesian unit vectors. The unit vectors for the prolate spheroidal coordinates are:

[TABLE]

Nonretarded Potential & Dispersion Relations

Considering a solid hyperboloid of revolution $\mu\geq\mu_{0}$ ($\mu_{0}>0$) with dielectric constant $\varepsilon_{1}$ immersed in a medium whose dielectric constant $\varepsilon_{2}$ is assumed to be $1,$ it follows that the potential of the electric $\Phi(\mu,\eta,\varphi)$ has to satisfy the Laplace equation $\vec{\nabla}^{2}\Phi=0$ everywhere except on the boundary surface $\mu=\mu_{0},$ where the Laplacian is given by [57]:

[TABLE]

Assuming the ansatz [57] $\Phi(\mu,\eta,\varphi)=F(\eta)G(\mu)e^{im\varphi},$ one finds that $F$ and $G$ satisfy the differential equations:

[TABLE]

Letting $c=\nu(\nu+1)$ (see [60, 65, 57]), where $\nu=-\frac{1}{2}+iq$ and $q\in[0,\infty),$ one can obtain a continuous spectrum of real eigenvalues and eigenfunctions in terms of the associated Legendre functions or conical functions $P^{m}_{-\frac{1}{2}+iq}(z)$ (denoted by $P_{mq}(z)$) of complex lower index with $-\infty<z<\infty$ satisfying the orthogonality relation [60, 66]:

[TABLE]

with

[TABLE]

A bounded satisfactory real-valued solution to the first equation in Eq. (2.69) is given by $P^{m}_{-\frac{1}{2}+iq}(\eta)$ (see [60, 65, 57]), whereas a pair of satisfactory real-valued solutions to the second equation in Eq. (2.69) are given by $P^{m}_{-\frac{1}{2}+iq}(\mu)$ and $P^{m}_{-\frac{1}{2}+iq}(-\mu)$ [57, 67, 68]. Avoiding the singularity at $\mu=-1$ [57], we shall choose the solution $P^{m}_{-\frac{1}{2}+iq}(\mu)$ for the inside region $\mu_{0}\leq\mu\leq 1$ and $P^{m}_{-\frac{1}{2}+iq}(-\mu)$ for the outside region $-1\leq\mu\leq\mu_{0}.$ Thus, the inside and outside potentials $\Phi_{\text{i}}$ and $\Phi_{\text{o}}$ are expressed as sums over the Fourier harmonics

[TABLE]

and

[TABLE]

respectively, where $A_{mq}(t)$ and $B_{mq}(t)$ are complex valued amplitudes and $m\in\mathbb{Z}$. Using the fact [69] that

[TABLE]

for all $z\in\mathbb{R}$ and $m\geq 0,$ it follows that one only needs to expand $\Phi_{\text{i}}$ and $\Phi_{\text{o}}$ in terms of $\cos m\varphi,$ for $m=0,1,2,\dots.$ More precisely, we use the system $H_{m}(\varphi)=(2-\delta_{0m})\cos\varphi$ satisfying the orthogonality relation:

[TABLE]

where

[TABLE]

Finally, imposing the continuity of the potential across the boundary $\mu=\mu_{0},$ we find the inside and outside potentials, which are expressed as:

[TABLE]

or explicitly [PPRB]:

[TABLE]

where $A_{mq}(t)$ are real time dependent amplitudes and $H_{m}(\varphi)=(2-\delta_{0m})\cos m\varphi,$ and, for convenience, we use $P_{mq}(\cdot)$ to denote $P^{m}_{-\frac{1}{2}+iq}(\cdot)$. Fig. 2.4 shows the spatial distribution of the lowest lying eigenmodes of the quasi-static electric potential for the hyperboloid. The values of conical functions $P^{m}_{-\frac{1}{2}+iq}(\cdot)$ at any given point were calculated through the algorithms provided in [70, 71].

Applying the Laplacian in Eq. (2.68) to the potential Eq. (2.79) and using a similar argument given in Eqs. (2.13)–(2.16), we can write:

[TABLE]

hence

[TABLE]

we have

[TABLE]

by imposing the first boundary condition $\mu=\mu_{0}$. Second derivatives could be calculated as:

[TABLE]

hence the Laplace equation becomes

[TABLE]

since Laplacian vanishes everywhere but on the surface of the hyperboloid, one finds

[TABLE]

Using the expressions for $\Phi_{\text{i}}$ and $\Phi_{\text{o}}$ in Eq. (2.79), we get

[TABLE]

where the Wronskian for $\mathcal{W}_{mq}(\mu)$ for the conical functions is given by:

[TABLE]

The exact value of Eq. (2.93) is given by [60]:

[TABLE]

where $Z_{q}^{m}$ is given by Eq. (2.72).

It follows from the orthogonality relation for the system $\{H_{m}(\varphi)\}$, given in Eq. (2.76), that for each fixed $m=0,1,2,\dots,$ we have:

[TABLE]

Applying the Van-Nostrand orthogonality relation for the conical functions given in Eq. (2.71) and the exact expression for the Wronskian given in Eq. (2.94), it follows from Eq. (2.95) that for each fixed $m$ and $q\geq 1$ the amplitudes $A_{mq}(t)$ satisfy the harmonic oscillator model :

[TABLE]

where the frequencies $\omega^{2}_{mq}$ are given by:

[TABLE]

Again, these resonant values can also be independently calculated from the transcendental equation generated by satisfying the boundary conditions. Dispersion relations $\varepsilon_{mq}$, using Drude model, are shown in Fig. 2.5.

Following Eqs. (2.13)–(2.16), the surface charge density $\sigma$ on $\mu=\mu_{0}$ is found to be:

[TABLE]

Classical Energy & Decay Rate

To calculate the classical energy $E,$ we follow the exact same argument as in the case of a paraboloid outlined in previous section. We hence find potential and kinetic energy in hyperboloid as:

[TABLE]

and

[TABLE]

Consequently, using the orthogonality relations in the hyperboloidal case given by Eqs. (2.76) and (2.71), one finds the hyperboloidal energy:

[TABLE]

Using similar ansatz for $A_{mq}(t)$ and $\dot{A}_{mq}(t)$ as in Eqs. (2.39) and (2.40), the coefficients $\alpha^{2}_{mq}$ can be obtained by a comparison of Eq. (2.101) with the Hamiltonian in Eq. (2.43) as:

[TABLE]

where $\hat{\delta}_{0m}$ is given by Eq. (2.77) and thus:

[TABLE]

Utilizing Eq. (2.45) and the hyperboloidal analogue of Eq. (2.49), setting up the plasmon current density, the interaction Hamiltonian for the photon emission can be expressed as:

[TABLE]

Thus, we may calculate the hyperboloidal emission matrix element as:

[TABLE]

where for different polarization directions, we obtain:

[TABLE]

and

[TABLE]

where we have set $E_{mq}(\eta,\varphi)=e^{-i\mathcal{B}}$, with:

[TABLE]

The wavevector $\mathbf{s}$ and the polarization vectors $\hat{\mathbf{e}}_{1},~{}\hat{\mathbf{e}}_{2}$ remain the same as before. Moreover, similar to the case for paraboloid, the choice of the polarization vector $\hat{\mathbf{e}}_{j}$ in each integral determines the polarization state. Hence, $\mathcal{I}^{(1)}_{mq}$ represents the s-polarization and $\mathcal{I}^{(2)}_{mq}$ corresponds to the p-polarization. Using Eq. (2.102) in Eq. (2.60), and now the position vector in hyperboloidal coordinates given by Eq. (2.3), we arrive at:

[TABLE]

Prolate Spheroidal Surfaces

In order to provide a geometric basis for comparison, we now treat the case of plasmon excitation on the surface of a prolate spheroidal domain. This structure presents an almost identical curvature to that of the paraboloid but encompasses a finite domain, making the interpretation of the surface modes and their associated radiation patterns more tangible. The closely related structure of an oblate spheroid has been employed in previous plasmon studies [7]. In both prolate and oblate systems, a limiting case is that of a sphere [8, 6], which can serve to validate the results. To obtain the quantized surface modes of the prolate spheroid, we closely follow the oblate case [7]. A prolate spheroid is defined by fixing the coordinate $\eta=\eta_{0}$ in coordinate system $(\eta,\mu,\varphi)$ given in Eq. (2.3). To avoid confusion, we use $\zeta$ in the case for prolate spheroid and recall the coordinate system $(\zeta,\mu,\varphi)$ as:

[TABLE]

in the domain defined by :

[TABLE]

Scale factors in this coordinate system are:

[TABLE]

The radius vector is

[TABLE]

The unit vectors for the prolate spheroidal coordinates are

[TABLE]

Nonretarded Potential & Dispersion Relations

The quasi-static scalar potential for a spheroid defined by $\zeta=\zeta_{0}$ may then be written as [72, 73]:

[TABLE]

or explicitly:

[TABLE]

for some real coefficients $A_{mlp}(t)$, with $l=1,2,3,\cdots$, $m=0,\cdots,l$, $p=0,1$ and $P^{m}_{l}(\cdot)$ and $Q^{m}_{l}(\cdot)$ being the associated Legendre polynomial of first and second kind, respectively, while $Y^{p}_{lm}(\mu,\varphi)$ are the real spherical harmonics given by [69]:

[TABLE]

with associated Legendre functions $P^{m}_{l}(\cdot)$ and orthogonality relation [74]:

[TABLE]

and $\delta$ is the Kroneker delta symbol. Note that the first argument of the real spherical harmonic, here $\mu$, is also taken to be $\mu=\cos\theta$ for some $-\pi\leq\theta\leq\pi$ is certain references. Since we developed the subsequent in terms of $\mu$, we avoid using conventional angle $\theta$.

In Fig. 2.6, the potential distribution for the fixed shape parameter $\zeta_{0}=65~{}$(nm) and fixed lowest modes $(l,m)=(1,0)$ is shown. The magnitude of the potential could be read by the color bar. Yellow color shows the highest intensity of the potential distribution while the dark blue shows the lowest. As expected, the lowest mode shows the mono-pole approximation for the field distribution.

With Heaviside function $\Theta$, the inside and outside potential could be written as given in Eq. (2.114) as :

[TABLE]

Laplace equation in this coordinate system is given by:

[TABLE]

We define the derivative operators $\mathcal{D}$ in each component $\zeta$, $\mu$ and $\varphi$ as:

[TABLE]

we can put:

[TABLE]

we write:

[TABLE]

Since the Laplace equation vanishes inside and outside but on the surface, then we have

[TABLE]

Using Eq. (2.115) along with Maxwell’s equations give us the charge density as $\rho=-\frac{\nabla^{2}\Phi}{4\pi}$, using the Laplacian we just obtained, we have

[TABLE]

where we have used the Wronskian identity as

[TABLE]

Using Eq. (1.57) one can find

[TABLE]

where $\sigma$ is the surface charge. We have

[TABLE]

On the other hand, Eq. (1.67) gives:

[TABLE]

Putting (2.129) and (2.130) equal, one finds frequency as

[TABLE]

The allowed values of the dielectric function $\varepsilon_{lm}$ are then found to be:

[TABLE]

following similar steps for Eqs. (2.8)–(2.24). In Fig. 2.7, the different plots of dielectric function $\varepsilon_{lm}$ as a function of shape parameter $\zeta$ for different $l$ and $m$, where $m<l$ are obtained.

Classical Energy & Radiative Decay Rate

We proceed to find the classical energy $E=T+V.$ The potential energy, taking the same steps as we did for all other cases, and using the orthogonality relation for spherical harmonics given in Eq. (2.117), can be calculated as:

[TABLE]

Similarly, kinetic energy is given by:

[TABLE]

The classical energy $E$ could be obtained by adding kinetic and potential energies, given in Eqs. (2.134) and (2.133), respectively, as:

[TABLE]

Real coefficients $A_{lmp}(t)$ obey the harmonic oscillators equation for motion, we can choose the form:

[TABLE]

for some complex time dependent functions $a_{lmp}$ which are proportional to $e^{-i\omega_{lm}t}$ and for some constant $\alpha_{lmp}$ which will be chosen later.

[TABLE]

This gives $\alpha_{lmp}$ as:

[TABLE]

Considering now the interaction Hamiltonian, followed by Eq. (1.130), as:

[TABLE]

where

[TABLE]

and $\mathbf{A}$ is given by Eq. (1.106), we could find:

[TABLE]

Emission matrix element in prolate spheroidal case following the same relations as in the paraboloidal and hyperboloidal cases, initiates by considering:

[TABLE]

We set:

[TABLE]

Assuming:

[TABLE]

then

[TABLE]

and

[TABLE]

We use the change of variable $\mu=\cosh\xi,$ to write:

[TABLE]

We consider the following coordinate transformation as:

[TABLE]

using the fact that $\mu_{0}$ is fixed. We define the gradient operator $\vec{\nabla}^{\prime}$ as

[TABLE]

Then comparing with

[TABLE]

we could write

[TABLE]

We use the following integral identity in the calculation of prolate spheroidal matrix element given by Eq. (16.127)[23]:

[TABLE]

where $j_{l}(\mathbf{k}r)$ is the Spherical Bessel function of order $l$ whose relation with regular Bessel function $J_{l}(\cdot)$ is given by:

[TABLE]

We then obtain the matrix elements for photon emission:

[TABLE]

using initial expression as in Eq. (2.142), similar to the oblate spheroidal case in [7], with $\mathbf{\hat{e}}_{j}$ being the unit polarization vector for the emitted light, $j_{l}(\cdot)$ spherical Bessel function of order $l$, (see Eq. (2.157)) and $\vec{\nabla}_{s^{\prime}}$ is the gradient of the wavevector given in spherical coordinate $(s^{\prime},\theta^{\prime},\varphi^{\prime})$. Following the same quantization scheme as before, radiative decay rate of plasmons per unit solid angle $\Omega$, is then calculated to be:

[TABLE]

where

[TABLE]

for the two polarization states s and p [7], with $(s^{\prime\prime},\theta^{\prime\prime},\varphi^{\prime\prime})$ denoting the wavevector in spherical coordinates, for which the transformations are given by taking the following transformations:

[TABLE]

Fig. 2.8 is obtained by using the experimental data for a silver spheroid provided by Hageman in [75, 76], in line with Fig. 2 in [7]. It shows, for the two specific modes $(l,m)=(1,0)$ and $(2,0)$, that the light emission reaches its maximal value at [math] and [math], respectively. Fig. 2.9, however, is a direct contour plot of Eq. (2.159) which illustrates the same result as shown in Fig. 2.8. Lastly, considering a sphere, as a limiting case of a spheroid by letting the shape parameter $\zeta\to 2$, Fig. 2.10 shows that the radiative decay rate becomes independent from emission angle $\theta$ as its shape approaches that of a sphere.

Comparison

The harmonic functions, constituting the eigenmodes of the hyperboloidal and paraboloidal systems, represent the normal modes of the charge density. The infinite axial dimension of these domains results in the corresponding eigenvalues ($\lambda$ and $q$) to form continuous spectra as opposed to the discrete value spectrum ($l$) of finite domains such as spheroidal particles. While, this difference in the nature of the eigenvalues is less noticeable when visualizing the potentials, as seen in Fig. 2.11,(b), (c) and (d), it becomes relevant when performing the quantization since now an integration over the corresponding eigenvalue spectrum enters the Hamiltonian, Eqs. (2.50) and (2.104). However, an inspection of the asymptotic properties of the integrands in the Hamiltonians reveal fast convergence, facilitating the calculation of the needed matrix elements.

To proceed, we note that for a more reasonable comparison of the two systems, their physical dimensions must be made comparable. Therefore, we geometrically adjust the paraboloid and the hyperboloid so as to control and match their apex curvature and further match it with that of a finite prolate spheroidal domain. The latter, due to its finite domain, makes a natural case to validate the findings for the two infinite domain cases investigated.

For settings consistent with nanofabrication and photon scattering experiments involving gold probes, where strong plasmon excitation has been reported [sam_poster, 55], we employ the following settings: $\eta_{0}=60$ nm and $\mu_{0}=86$ nm for the systems in Fig. 2.11(b), (c) and (d). The close relationship of the calculated eigenvalues with the surface plasmon momenta are clearly observable from the ”wavelength” of the charge density oscillations in Fig. 2.11(b), (c), (d) and Fig. 2.13. For proper geometric and modal adjustments, the similarities in the potential distributions are clearly evident from the contour plots. Since for the argument values $\mu\approx 1,-1$, the conical functions become singular, a numerical artifact in the form of a discontinuity in the contour lines appear near the symmetry axis of the hyperboloids, where $\mu$ attains those values. Comparison with the corresponding potential distribution in the spheroidal case can be facilitated by taking $\zeta_{0}=65$ nm matching its curvature with that of the apexes in paraboloidal and hyperboloidal domains. With the dimensions adjusted, using Eqs. (2.8), (2.79) and (2.115) to simulate the spatial distributions of the potentials for the three cases, one clearly observes the analogous role of $\lambda$ and $q$ to the discrete spectrum $l$. The lowest azimuthal mode $m=0$ for the three domains shown in Fig. 2.11(b), (c) and (d) was simulated by taking the second indices as $\lambda=1$, $q=0.2$ and $l=1$, generating relatively same potential distributions. The potential distributions for higher modes are shown in Fig. 2.12, where two higher modes for each domain (shape parameters $\eta_{0}$, $\mu_{0}$ and $\zeta_{0}$) are shown. In doing so, the Legendre functions $P_{mq}(\mu)$ of imaginary order and their derivates have been calculated using the computational algorithm of Gil and Segura [70], and the integral expansion of Kölbig [71]. Analogous to plasmon wavevector in the case of excitations on a planar interface (or a Cartesian thin film), which can be emulated by $\mu_{0}\to\pi/2$ in Fig. 2.11(a), the higher the $\lambda$, $q$, and $l$ for the same $m$, the higher the number of fluctuations for the same spatial domains, as shown in Fig. 2.13. Furthermore, from the spheroidal nearfield distribution, one can readily observe the multipole order so that $(m,l)=(0,1)$ corresponds to a dipolar distributions, whereas $(m,l)=(0,2)$ leads to a quadrupolar behavior, etc. Similarly, for $(m,q)=(0,1.5)$ or $(m,\lambda)=(0,2)$ one obtains the corresponding multipolar nearfield distributions of the apex regions. Thus, guided by these charge density oscillations, controlled by the parameters $(m,\lambda,\eta_{0}),$ $(m,q,\mu_{0})$, and $(m,l,\zeta_{0})$ for the three cases and by a geometric matching of their apex curvatures, one may discuss the eigenmode dependent radiative decay rates.

In addition to the material-specific electronic and optical properties, in as far as the effect of the local curvature is concerned, we expect the considered cases to exhibit similarities in their resonance spectra. The plasmon dispersion relations are good indicators of this resemblance. To study the relation between resonance modes of isolated solid paraboloids and hyperboloids in vacuum, we assume a local dielectric function and calculate their eigenmode dependent allowed values $\varepsilon_{m\lambda}$ and $\varepsilon_{mq}$ using Eqs. (2.24) and (2.97). A comparison of the lowest lying plasmon modes ($m=0,1,2$ and $3$ of fixed probes $\eta_{0}$ and $\mu_{0}$) is given in Fig. 2.14. The modes may alternatively be displayed with reference to bulk plasma frequency $\omega_{p}$. Interestingly, the higher $m$ modes show a higher sensitivity to the morphological differences between the two systems, in particular for lower $\lambda$ and $q$ values, that is, in the long wavelength limit, which in analogy with the planar plasmons would be near the light line, where retardation effects are more pronounced [PPRB]. We also note that in the hyperboloidal case, in the limit $\mu_{t}\to 0$, we have $\varepsilon_{mq}\to-1$; that is, the modes asymptotically approach the surface plasmon resonance ($\omega\to\omega_{p}/\sqrt{2}$), as expected for a Cartesian half-space. This limit is also approached by large $m$ values as seen in Fig. 2.14. In the short wavelength regime $\lambda,q\to\infty,$ the dielectric function reads:

[TABLE]

yielding the same surface plasmon limit [PPRB, 67]. Similar to the paraboloidal case, using the asymptotic behavior of modified Bessel functions $I_{m}(\lambda\eta)$ and $K_{m}(\lambda\eta)$ for large order $m$ and large arguments $\lambda$ and $\eta$ (see Eqs. (9.7.8) and (9.7.9) in [74]), we can write:

[TABLE]

which implies the same limit for large $\lambda$.

Following the field distribution and resonant dielectric values corresponding to the normal modes of the surface charge density oscillations, we may assume that a plasmon has been created in a given eigenmode $(m,o,p)$ where $o=\lambda,q,l$ designate paraboloid, hyperboloid and spheroid, respectively. If the plasmons on the paraboloidal, hyperboloidal or spheroidal surface are in the initial state $\hat{a}^{\dagger}_{(mop)_{i}}\ket{0},$ the probability amplitude for their emission into the final photon state $\hat{c}^{\dagger}_{(\mathbf{s}j)_{f}}\ket{0}$ could be obtained using expressions (2.54), (2.105) and (2.158). The $\ket{0}$ indicates that the fields have been populated with 0 plasmons or photons (noting $a_{mop}\ket{0}=c_{(\mathbf{s}j)}\ket{0}=0$) whereas a general photon-plasmon state is written as $\ket{\Psi}=\ket{n_{\mathbf{s}j}}\otimes\ket{n_{mop}},$ that is, a state with $n_{\mathbf{s}j}$ $j$-polarized photons of momentum $\mathbf{s},$ and $n_{mop}$ plasmons in the $(m,o,p)$ state.

Keeping the mode patterns, here, the $\varphi=0,\,\pi$-plane projection of the relative potential distributions in Fig. 2.11(b), (c) and (d) and dispersion relations in Fig. 2.14 in mind, we now compare the radiative decay rate per solid angle of plasmons engendered on the three domains, using Eq. (2.61) (paraboloid with $\eta_{0}=60$ nm), Eq. (2.108) (hyperboloid with $\mu_{0}=86$ nm) and Eq. (2.159) (prolate spheroid with $\zeta_{0}=65$ nm) for two lowest azimuthal mods $m=0$ and $1$. To calculate the matrix elements given in Eqs. (2.54) and (2.105), the integrations $\mathcal{I}^{(j)}_{m\lambda}$ and $\mathcal{I}^{(j)}_{mq}$ (in case for paraboloid and hyperboloid, respectively) must be carried out numerically as they lack analytical solutions in variable $\xi$ in the case for paraboloid and $\eta$ in the case for hyperboloid. The choice of polarization vectors $\mathbf{\hat{e}}_{j}$ for $j=1,2$ corresponds to $\mathcal{I}^{(1)}_{m\lambda}$ given by Eq. (2.57) to represent the s-polarization, and $\mathcal{I}^{(2)}_{m\lambda}$ given by Eq. (2.58), the p-polarization. Inspecting the integrands for their convergence, we compute the integrals using an iterative numerical integration scheme (Runge-Kutta) due to lack of fast oscillations. The result is shown in Fig. 2.15. Note that to facilitate visual comparison within the same plot window, the radiative decay rate for prolate spheroid for modes $m=0$ and $m=1$ have been multiplied by $0.006$ and $0.01$, respectively. The effect of higher index modes $\lambda$, $q$ and $l$ with the same azimuthal order $m$ on radiative decay rate per solid angle may also be studied, as shown in Fig. 2.16.

As one may expect from the nearfield patterns of higher $\lambda,$ and $q,$ analogous to $l=2$, a quadrupolar pattern appears for the emitted photons. Here, it is understood that an angular segment is occupied by the probe itself, as opposed to the $0-2\pi$ range for the finite volume spheroidal systems. It is further understood that for larger particles or probe apex size retardation effects may modify the higher order modes.

The dependence of the radiative decay rate upon the parameter that sets the boundary, allows for control of the curvature and thus the photon emission patters. As can be seen from Fig. 2.17, for the $m=0$ mode, the higher the curvature of the hyperboloidal apex, the lower the amplitude and the narrower the angular distribution of the emitted photons. Moreover, for $m=1$ the higher curvature while resulting in a lower amplitude does not result in in a higher angular confinement.

CHAPTER 3AN APPROXIMATION FOR RADIATIVE DECAY RATE ON THE NON-SIMPLY CONNECTED SURFACE OF A TORUS

Introduction

Photon scattering and excitation of surface plasmons in single rings, single composite rings and many-ring systems was recently shown to provide useful dispersion and field distributions [1, 2] with specific applications in particle and molecular trapping [77, 78]. Following photon or electron excitation of metallic nanorings, both radiative and nonradiative decay channels are important in considering of the nanoparticle as a photon or phonon source, and in related applications. Here, we investigate the radiative decay channel for a vacuum bounded single solid nanoring by quantizing the fields associated with charge density oscillations on the nanoring surface. In quantizing the fields, the frequency spectrum of the charge density normal modes of the nanoring is obtained and shown to agree with the exact quasi-static plasmon dispersion relations. We calculate the radiative decay rate per unit solid angle and show that both poloidal and toroidal modes contribute to the radiation field.

Particle scattering from material domains that take the form of nanoparticles, lead to useful excitations and surface modes [8, 79, 80]. Both electron and photon bombardment are known to lead to resonant collective electronic effects in thin films and particles, from which emitted light due to the damping of excited plasmons has been observed in many experiments [81]. Following an excitation event of surface plasmon, the radiative decay rate of the associated modes on the bounding surfaces of the underlying material domains, and their radiation pattern can help understand and interpret experimental observations and emission measurements in the nearfield as well as in the far field. In this chapter, we present a calculation of the differential scattering cross section for a nanoring, from quantizing the charge density oscillations. The presented formalism and results can be used to study the decay properties of other photon scattering and plasmon-related processes such as photoabsorption, photoemission, tip-enhanced, and surface-enhanced Raman scattering. By obtaining the angular, polarization and wavelength characteristics of the radiation from the excited surface-plasmon modes on various particles, experimental results may be interpreted.

Our approach to calculate the radiative decay of the toroidal and poloidal plasmons closely follows previous works [7, 11, 82] based on Bloch linearized hydrodynamic model [6], leading to the Hamiltonian

[TABLE]

Specifically, the radiative decay of surface collective excitations in polarizable cylinders by [11], calculation for a system of small metallic spheres embedded at random in a thin dielectric slab by [6], investigation of a plane slab by [10], and oblate spheroids by [7] and its special case of a metallic 2D disc by [82]. In these works, using nonretarded electrodynamics and bulk-metal optical properties, plasmons have been described on various particles with subwavelength dimensions, although the results may be made more accurate if surface-metal optical properties can be used.

This chapter is organized as follows. We begin by modeling a nanoring as a solid ring toroid in toroidal coordinates. We then proceed to determine the surface plasmon field and dispersion relations assuming a local dielectric function and nonretarded electrodynamics. Relevant modes are discussed within the free electron gas model. Subsequently, the energy of the system is obtained and the field is quantized to obtain the surface-plasmon operators. Using the interaction Hamiltonian of Ritchie, and presenting the photon field via quantized vector potential, we determine the radiative decay rate of the plasmons per unit solid angle as a function of the geometric and excitation parameters.

Nonretarded Potential & Dispersion Relations

The toroidal coordinates $(\mu,\eta,\varphi)$ are related to the Cartesian coordinates $(x,y,z)$ via:

[TABLE]

where $a>0$ is an overall scaling factor, $\mu\geq 0,$ and $\eta,\phi\in[0,2\pi)$. The associated toroidal scale factors are given by

[TABLE]

For any fixed value $\mu_{0}>0,$ the equation $\mu=\mu_{0}$ describes the surface of a torus whose inside and outside domains are given by $\mu>\mu_{0}$ and $\mu<\mu_{0},$ respectively. [57].

We consider a solid torus in the toroidal coordinate system $(\mu,\eta,\varphi)$ whose boundary surface is described by $\mu=\mu_{0}$ and its dielectric function is denoted by $\varepsilon.$ The surface modes of a single and a composite nanoring were recently reported in [2]. Specifically, following the quasi-static formulation of the toroidal problem, the graphic presentation, and the plasmon dispersion relations of engendered normal modes in [1], we write the scalar electric potential of a nanoring as:

[TABLE]

where $\Phi_{\text{i}}$ and $\Phi_{\text{o}}$ denote the inside- and the outside-potentials [2, 83]:

[TABLE]

[TABLE]

where $m,n=0,\pm 1,\pm 2,\cdots,$ $C_{mn}(t)$ are the complex-valued time dependent amplitudes,

[TABLE]

and $P_{n}^{m}(\mu),$ $Q_{n}^{m}(\mu),$ $P^{\prime m}_{n}(\mu),$ $Q;^{m}_{n}(\mu)$ stand for the associated Legendre functions of the first and second kind $P^{m}_{n-\frac{1}{2}}(\cosh\mu),$ $Q^{m}_{n-\frac{1}{2}}(\cosh\mu),$ and their derivatives with respect to $\cosh\mu$. Their evaluations at the boundary $\mu=\mu_{0}$, $P^{m}_{n-\frac{1}{2}}(\cosh\mu_{0}),$ $Q^{m}_{n-\frac{1}{2}}(\cosh\mu_{0}),$ are shown by $P_{n}^{m}$ and $Q_{n}^{m},$ respectively. This convention will be used throughout the rest of the paper along with $f$ and $f_{0}$ to denote $f(\mu,\eta)$ and $f(\mu_{0},\eta)$, respectively. Moreover, the format adopted in Eqs. (3.5) and (3.6) implies the continuity of the potential across the toroidal boundary; i.e.,

[TABLE]

Since $\vec{\nabla}^{2}{\Phi}=0$ for the $\mu<\mu_{0}$ and $\mu>\mu_{0},$ the Poisson equation takes the form

[TABLE]

where $\sigma$ denotes the surface charge density on $\mu=\mu_{0}$ and $\delta(x)$ is the Dirac delta function.

In toroidal coordinate system (see [37]), the Laplacian of a function $F$ in the form $F=fg,$ where $f$ is given by Eq. (3.7), may be written as:

[TABLE]

where the operator $D$ is given by:

[TABLE]

In view of Eqs. (3.5) and (3.6), we may write the inside and outside potentials as $\Phi_{\text{i}}=f\Psi_{\text{i}}$ and $\Phi_{\text{o}}=f\Psi_{\text{o}},$ where

[TABLE]

[TABLE]

Using Eq. (3.10), the Laplacian of $\Phi$ given by Eq. (3.4) may be calculated as:

[TABLE]

where

[TABLE]

To find $\mathcal{D}\Psi,$ noting that $\Theta$ only depends on $\mu,$ it suffices to calculate the expressions for $\frac{\partial\Psi}{\partial\mu}$ and $\frac{\partial^{2}\Psi}{\partial\mu^{2}}.$

Using the continuity of the potential $\Phi$ across the boundary surface, it follows that $\Psi_{\text{i}}=\Psi_{\text{o}}$ at $\mu=\mu_{0}$ and thus the term $\delta(\mu-\mu_{0})\Psi_{\text{i}}-\delta(\mu_{0}-\mu)\Psi_{\text{o}}$ vanishes identically yielding

[TABLE]

Differentiating Eq. (3.16) once again with respect to $\mu$ gives

[TABLE]

Using the obtained expressions for $\frac{\partial\Psi}{\partial\mu}$ and $\frac{\partial^{2}\Psi}{\partial\mu^{2}}$ given in Eq. (3.16) and Eq. (3.17), one may write $\mathcal{D}\Psi$ as:

[TABLE]

Now, since $\vec{\nabla}^{2}\Phi_{\text{i}}=0$ for $\mu>\mu_{0}$ and $\vec{\nabla}^{2}\Phi_{\text{o}}=0$ for $\mu<\mu_{0},$ it follows from Eq. (3.10) that $\mathcal{D}\Psi_{\text{i}}=0$ for $\mu>\mu_{0}$ and $\mathcal{D}\Psi_{\text{o}}=0$ for $\mu<\mu_{0},$ and thus the last two terms in Eq. (3.18) both vanish identically. As a result, Eq. (3.18) and Eq. (3.14) imply:

[TABLE]

where from Eq. (3.12) and Eq. (3.13)

[TABLE]

with $\mathcal{W}^{m}_{n}(\mu)$ denoting the Wronskian given as:

[TABLE]

Finally, in view of the Poisson equation, Eq. (3.19) and Eq. (3.20) imply the claimed expression for the surface charge density given as:

[TABLE]

On the other hand, the surface charge density may be expressed using the polarization argument given in Eq. (1.67). To establish $\ddot{\sigma}$, we first observe that [8],

[TABLE]

Letting $\Phi_{\text{i}}=f\Psi_{\text{i}}$ as above, together with the facts that $h_{\mu_{0}}=af_{0}^{-2}$ and

[TABLE]

we may write the right-hand side of Eq. (3.23) as:

[TABLE]

Utilizing the expression for $\Psi_{\text{i}}$ given in Eq. (3.12), we arrive at Eq. (3.26) in view of Eq. (3.25).

A straightforward calculation shows that

[TABLE]

Differentiating Eq. (3.22) twice with respect to time $t$ and comparing the result with Eq. (3.26) , it follows form the orthogonality of the system $\{e^{im\varphi}\}$ that for each fixed $m=0,\pm 1,\dots,$ the amplitudes $C_{mn}$ satisfy

[TABLE]

For $\mu_{0}>0,$ function $\ \frac{1}{\cosh\mu_{0}-\cos\eta}$ can be expanded in a uniformly convergent Fourier series on $[-\pi,\pi]$ as

[TABLE]

where the Fourier coefficients $a_{l}$ are given by:

[TABLE]

Using the substitution: $z=e^{i\eta},$ $\dfrac{1}{iz}\,dz=d\eta,$ and $\cos\eta=\frac{1}{2}(z+z^{-1}),$ the integral in (3.29) may be transformed into a contour integral over the positively oriented unit circle $\partial\mathbb{D}$ in the complex plane as:

[TABLE]

where

[TABLE]

is a rational function having two simple poles at $z=e^{\pm\mu_{0}}$ for all $l\in\mathbb{Z}$ with only $e^{-\mu_{0}}$ lying inside the unit circle and a pole of order $l$ at $z=0$ for $l\geq 1$. It follows from the Residue Theorem [41] that

[TABLE]

where

[TABLE]

and, for $l\geq 1,$

[TABLE]

Noting that:

[TABLE]

together with the Leibniz’s rule:

[TABLE]

it follows form Eq. (3.34) that:

[TABLE]

Finally, in view of Eqs. (3.32)–(3.37), we conclude that:

[TABLE]

Now, replacing $\dfrac{1}{2f_{0}^{2}}$ by its Fourier series expansion:

[TABLE]

and utilizing the substitution $n+l\mapsto n,$ it follows from the orthogonality of the system $\{e^{in\eta}\}$ that for each fixed $n\in\mathbb{Z}$ in Eq. (3.27), we may write:

[TABLE]

where

[TABLE]

and

[TABLE]

Cylindrical limit

In order to use the well-developed theories of a cylindrical model, it is often useful to consider the cylindrical limit of the toroidal coordinate through a transformation from toroidal to cylindrical geometry by keeping the minor radius of a torus, $r$, fixed while the major radius $R$ increases indefinitely with $\mu=\mu_{0}$. Under these assumptions, the ring transforms into the column and one could utilize the theories of a cylindrical model. We refer the reader to [84, 85] for a detailed description of the cylindrical model. To this end, we define the ratio:

[TABLE]

and we let $m$ and $\mu$ increase indefinitely such that $k$ remains fixed. Following Eq. (3.43), we may write:

[TABLE]

Next, we consider the cylindrical limits of $\mathcal{V}^{m}_{n,l}$. Using Eq. (3.41), we can write:

[TABLE]

noting:

[TABLE]

Following Eq. (3.42), we obtain:

[TABLE]

The cylindrical limit of $P_{n}^{m}(\mu_{0})$ and $Q_{n}^{m}(\mu_{0})$ in terms of modified Bessel functions of first and second kind are given by [69]:

[TABLE]

where

[TABLE]

Now, in view of Eqs. (3.48) –(3.50),

[TABLE]

and

[TABLE]

where, in Eq. (3.51), we have utilized the Wronskian identity for modified Bessel functions [57]:

[TABLE]

Finally, noting that

[TABLE]

and

[TABLE]

we obtain:

[TABLE]

where

[TABLE]

which is the same as the cylindrical-limit frequency given in [85].

Classical Energy

Potential energy, following Eq. (1.74), and using Eq. (3.22), may be expressed as:

[TABLE]

where using orthogonality in $n$ and $\eta$ ($m$ and $\varphi$) as:

[TABLE]

and following the same procedure as was taken in Chapters 1 and 2, we may find the potential energy $V$ as:

[TABLE]

Kinetic energy calculation’s steps are different from the general procedure we introduced in Chapter 1 due to the coupling between the modes in case of nanoring. The fact that the geometric structure of a ring is not simply connected and that Laplace’s equation is not separable in this coordinate system results in some changes. Eq. (3.40) shows that amplitudes in this case does not undergo pure harmonic oscillator equation of motion.

The kinetic energy is given by Eq. (1.75). The acceleration is the force per unit mass following Eq. (1.76) gives:

[TABLE]

where

[TABLE]

Hence

[TABLE]

and

[TABLE]

Using the identity:

[TABLE]

in Eq. (3.65), we can write:

[TABLE]

Hence

[TABLE]

The Divergence theorem leads to:

[TABLE]

By means of the expression for inside potential given in Eqs. (3.5) and  (3.6), we write:

[TABLE]

where $\overline{\Phi_{\text{i}}}$ denotes the complex conjugate of inside potential given in Eq. (3.5) as well as using the fact:

[TABLE]

The chain rule gives:

[TABLE]

Hence:

[TABLE]

Using scale factor relations:

[TABLE]

and some simple calculations:

[TABLE]

Using expression for $\frac{1}{f_{0}^{2}}$, as given in Eq. (3.39), we get

[TABLE]

Letting $n+l\rightarrow n$ and isolating the term $l=0$, the orthogonality in $\eta$ and $\varphi$ implies:

[TABLE]

with the bulk frequency identity given in Eq. (1.18). If we multiply and divide each term in $m,n$ by the Wronskian, and use the relations given in Eqs. (3.41) and (3.42), we can write:

[TABLE]

From Eq. (3.40), it follows that:

[TABLE]

Using Eq. (3.79) in Eq. (3.78) and integrating with respect to time $t$, we may write:

[TABLE]

Inspired by cylindrical limit, we may assume amplitudes undergo harmonic oscillator equation of motion (see Eqs. (3.46)–(3.58) and [85]). We therefore set:

[TABLE]

in Eq. (3.80) to obtain:

[TABLE]

One can easily find the total classic energy $E=T+V$, using Eqs. (3.61), (3.82) as:

[TABLE]

Complex coefficients $C_{mn}(t)$ could be written in the form [44]:

[TABLE]

for some complex time dependent functions $c_{mn}$ which are proportional to $e^{-i\omega_{mn}t}$ and for some constants $\gamma_{mn}$ which will be chosen later. The derivative with respect to time gives us:

[TABLE]

Using this, total energy could be written as:

[TABLE]

where $\bar{c}_{mn}$ denotes the complex conjugate of $c_{mn}$.

Interaction Hamiltonian & Radiative Decay Rate

The full classic Hamiltonian could be written as:

[TABLE]

This is understood to act on bosonic Fock space [44]. States of which we will denote using the usual Dirac ket notation. Using Wronskian in terms of Gamma function [74] given in Eq. (3.21), we can choose:

[TABLE]

in Eq. (3.86).

In analogy with other cases, the interaction of these two fields is described via the minimal coupling:

[TABLE]

then the matrix element for emission reduces to

[TABLE]

with

[TABLE]

and

[TABLE]

where

[TABLE]

and

[TABLE]

Lastly, the decay rate per solid angle $\Omega$ is given by:

[TABLE]

CHAPTER 4CONCLUSION

44footnotetext: Portions of this chapter were reprinted from: Bagherian, M. and Kouchekian, S. and Rothstein, I. and Passian, A., Quantization of surface charge density on hyperboloidal and paraboloidal domains with application to plasmon decay rate on nanoprobes, 125413, vol. 98, Sep 2018, with permission from Phys. Rev. B

Permission is included in Appendix C.

On the basis of the simulated results, the presented quantization scheme emerges as a reasonable approach to obtaining the qualitative behavior of extended electronic systems, for which many body calculations may be demanding or unfeasible, analytically or computationally. The obtained analytical results for the charge density oscillations on the surface of geometric entities with local curvature but with an extended dimension constitute first time results with potential for modeling quantum effects in plasmonics. The obtained analytical expressions for the operators associated with surface plasmon quantities, help calculating interactions with other quantized fields, e.g., the interaction of the probe with a nearby quantum emitter, or with the radiation field of a quantum dipole. Owing to their apex symmetry and curvature, the comparison of the calculated quantities show that hyperboloidal and paraboloidal plasmons qualitatively exhibit similar dispersion relations and radiative decay rates. We conclude that, the quantitative differences observed in the allowed resonance values of the dielectric function and the emitted radiation patterns are therefore primarily attributed to the difference in curvature in the asymptotic region, away from the apex. From a comparison with the modes excited on a prolate spheroidal surface, for which experimentally typically only low energy dipolar and quadrupolar eigenmodes or their mixtures have been observed to contribute to far field radiation in the visible spectrum, we expect that also only low lying modes will contribute to the emission spectra of the probe. Unlike the all discrete spectrum of the quantum numbers associated with the spheroidal modes, the eigenvalues characterizing the hyperboloidal or paraboloidal plasmons along their infinite dimensions exhibit a continuous spectrum, making a direct comparison of plasmon wavevectors unclear. However, visualization of the eigenmode field patterns, that is, how the fields fluctuate for given continuous eigenvalues in the case of hyperboloids and paraboloids, can facilitate the comparison with the discrete eigenvalues for the spheroidal case. For a given set of eigenmodes, photon emission from the higher apex curvature tips occurs with a more localized radiation pattern. Following the presented results for the allowed values of the dielectric functions $\varepsilon_{ml}(\zeta_{0})$, $\varepsilon_{m\lambda}(\eta_{0})$ and $\varepsilon_{mq}(\mu_{0})$, the corresponding dispersion relations and radiative decay rates can be obtained for real materials from a comparison with the experimentally determined optical properties of solids (such the compilation by Hageman [75], or by Palik [86]). In summary, the presented results can aid the design and fabrication of tips with specific photonic and plasmonic characteristics. For example for a gold or silver tip, such as those used in electron emission experiments, using the presented results one may determine both the fabrication design parameters and the excitation laser wavelength and polarization. In such applications, a comparison of Fig. 2.14 and 2.13 with the optical properties of, for example silver [75], indicates the availability of several resonance modes in the visible. The results can help the analysis of the radiation emitted from the nanotips as a result of electron or photon scattering, which are of importance to plasmonics in experiments such as EELS (electron energy loss spectroscopy) and SPM. For a specific material typically employed in plasmonics such as gold and silver, the results offer estimates of the polarization, angular, and spectral properties of the emitted radiation. In such instances, the emitted photons may be analyzed for the specific eigenmode of charge density oscillation $(m,\lambda)$ that created them. In light of the obtained multi-parameter dependence, any agreement with the theory would require fabrication control of the geometric features of the nanostructure. Undoubtedly, many features remain to be explored from a mathematical as well as physical point of view. For example, the geometric origin of the mode coupling in the toroidal solutions and its significance within a bigger topological theory or morphological picture warrant further exploration. What are the class of systems for which similar couplings occur? What are the implication of the geometric characteristics for the dynamics of the electronic system within the many-body theory? What are some of the implications of the specifics of the quantization volumes for the unbounded material systems? Devising appropriate truncation criteria or introducing suitable perturbation theoretical approaches to calculating the Hamiltonian of the system are other examples of topics needing further investigation. The connection between the underlying physical processes and the corresponding mathematical formulation, such as presented here, offers a rich spectrum of heretofore unknown mathematics and physics. Examples of the productivity of such connections are many including the discovery of a new orthogonality relation for the MacDonald’s functions [87] and new integral expansions [88] .

Radiative Decay on Spherical Surfaces

Spherical coordinates are given by [57]:

[TABLE]

where $a_{0}$ is the overall scale factor and:

[TABLE]

The scale factors for a parametrization in each component are given by [57]:

[TABLE]

Laplacian in spherical coordinate is given by [64]:

[TABLE]

Using scaling factors given in Eq. (A-3), we can write:

[TABLE]

Assuming ansatz, using separation of variable:

[TABLE]

where $R(r)$, $\Theta(\theta)$ and $\Psi(\varphi)$ are functions of variables $r,\theta$ and $\varphi$, respectively. Therefore,

[TABLE]

and,

[TABLE]

Substituting Eqs. (A-7) and (A-8) into Eq. (A-5), we get:

[TABLE]

Having separated Laplace’s equation into two ordinary differential equations, the general solution given by the suitable ansatz for the potential $\Phi(\mathbf{r},t)$ as:

[TABLE]

for some real time-dependent amplitudes $A_{l}(t)$ and $B_{l}(t)$ and real spherical harmonics $Y^{p}_{lm}(\theta,\varphi)$ given by: [57] :

[TABLE]

where $\delta_{kp}$ is the Kroneker delta symbol and $P^{m}_{l}(\cos\theta)$ denotes real associated Legendre functions of first kind [69, 89]. Inside and outside potential of a sphere can be constructed by imposing proper boundary conditions over the solution of Laplace’s equation given in Eq. (A-10). Inside and outside potential, denoted by $\Phi_{\text{i}}$ and $\Phi_{\text{o}}$ respectively, should satisfy the following conditions along the boundary of a solid sphere defined by $r=r_{0}$ using Eqs. (1.44) and (1.46), we may write:

[TABLE]

Therefore, we could write inside and outside potentials as:

[TABLE]

Utilizing the Heaviside function $\Theta(x)$:

[TABLE]

whose derivative gives:

[TABLE]

where $\delta$ denotes the Dirac Delta function, we define potential function implicitly in terms of inside and outside potentail $\Phi_{\text{i}}$ and $\Phi_{\text{o}}$ given by Eqs. (A-14) and (A-15) as:

[TABLE]

Its derivative with respect to variable $r$ gives:

[TABLE]

Using that fact that on the boundary $r=r_{0}$,

[TABLE]

we can write:

[TABLE]

Similarly,

[TABLE]

and,

[TABLE]

Second derivatives can be calculated as:

[TABLE]

and,

[TABLE]

and,

[TABLE]

Using Eqs. (A-24)–(A-26) in Eq. (A-9), we have:

[TABLE]

Since Laplacian vanishes inside and outside but on the surface of the sphere, then we can write:

[TABLE]

where using (A-14) and (A-15):

[TABLE]

Using this in Eq. (A-28), we can write:

[TABLE]

In order to find surface charge density $\sigma$, we use Eq. (A-30) in Eq. (1.58) to write:

[TABLE]

The only time-dependent component in the latter equation are amplitudes $A_{l}(t)$. Hence, the second time-derivative of Eq. (A-31) gives:

[TABLE]

Amplitudes $A_{l}(t)$ undergo harmonic oscillations at continuous frequencies $\omega_{l}$, that is:

[TABLE]

We utilize the relation given in Eq. (1.67) to write:

[TABLE]

Using the equality of (A-32) and (A-34), we find:

[TABLE]

where,

[TABLE]

which is known as the surface plasmon eigen-frequency and is given by [8], p. 129. It is noteworthy that the frequency in case of a sphere is independent of index $m$. One can use the latter expression for frequency in Drude model and finds dielectric function of a sphere as:

[TABLE]

The potential energy, using Eq. (A-14) and (A-31), could be written:

[TABLE]

as we replaced real-valued potential with its complex conjugate as they are equal. This allows us to utilize the orthogonality relations for the orthogonal system $\{Y^{p}_{lm}\}_{l,m,p}$ known as spherical harmonics system. This relation is given by [57, 69]:

[TABLE]

Applying Eq. (A-39) to Eq. (A-38), we may write:

[TABLE]

using Eq. (A-36), the potential energy could be written as:

[TABLE]

Similar calculations, using Eq. (A-39), give the kinetic energy as:

[TABLE]

By means of Eqs. (A-41) and (A-42), total energy $E$ could be found as:

[TABLE]

As described in Chapter 1, in order to convert $E$ into the Hamiltonian operator for a scalar boson field (plasmons are spinless quasi-particles), we write:

[TABLE]

where $a_{lmp}$ are some complex time-dependent functions, proportional to $e^{-i\omega_{l}t}$, $a^{*}_{lmp}$ its complex conjugate and coefficients $\alpha_{lmp}$ would be determined later. Its time derivative gives:

[TABLE]

Using Eqs. (A-44) and (A-45) in Eq. (A-43), we have:

[TABLE]

Performing field quantization [8, 32], we replace the amplitudes $a_{lmp}(t)$ and $a^{*}_{lmp}(t)$ with operator valued distributions $\hat{a}_{lmp}$ and its conjugate $\hat{a}^{\dagger}_{lmp},$ and note the commutation relations:

[TABLE]

Taking the normal ordered expansion of the noncommuting boson creation $\hat{a}^{\dagger}_{lmp}$ and annihilation $\hat{a}_{lmp}$ operators, a comparison with the normal ordered expression of the Hamiltonian operator for a scalar boson field [8], yields the Hamiltonian:

[TABLE]

Comparing Eq. (A-46) with Eq. (A-48), one finds:

[TABLE]

Next, in order to calculate the current $\mathbf{J}$ as given in Eq. (1.122), following Eq. (1.123), we set:

[TABLE]

Using the inside potential expression given in Eq. (A-14), one can write:

[TABLE]

upon making the replacement $\dot{A}_{l}(t)\to-i\frac{\alpha_{lmp}}{2}~{}(a^{*}_{lmp}-a_{lmp})$ given by Eq. (A-45) in the expression for $\dot{\Psi}$. Given interaction Hamiltonian by Eq. (1.128), using Eqs.(1.102) and (A-51), we can write:

[TABLE]

which, after using scaling factors given in Eq. (A-3), simplifies to:

[TABLE]

The relevant transition matrix element for photon emission, given by Eq. (1.112), using Eq. (A-53), could be obtained as:

[TABLE]

Using the identities:

[TABLE]

we can write:

[TABLE]

hence,

[TABLE]

Let:

[TABLE]

using Eq. (A-61), we write:

[TABLE]

Using the identity (see [23], p. 767):

[TABLE]

where $\bf s$ is a wave vector in spherical coordinate $(r,\theta^{\prime},\varphi^{\prime})$. Replacing $i$ with $-i$, we have:

[TABLE]

Hence the emission matrix element becomes:

[TABLE]

Using Eq. (A-49), we write:

[TABLE]

which for scale factor $a_{0}=1$ is the same as matrix element given in by Eq. (51) in [8], p. 129, as:

[TABLE]

On the other hand, considering $\bf s$ being a vector in spherical coordinate, we have:

[TABLE]

Since the polarization vector must be perpendicular to the direction of travel, thus $(\mathbf{\hat{e}}_{q}\cdot\mathbf{\hat{e}}_{r})=0$, hence:

[TABLE]

where:

[TABLE]

and,

[TABLE]

Taking the unit vectors $\mathbf{\hat{e}}_{\theta}$ and $\mathbf{\hat{e}}_{\varphi}$ as:

[TABLE]

Taking $(\mathbf{\hat{e}}_{q_{1}}=\mathbf{\hat{e}}_{\theta})$ and $(\mathbf{\hat{e}}_{q_{2}}=\mathbf{\hat{e}}_{\varphi})$ for s and p-polarization vectors, respectively, then:

[TABLE]

Letting:

[TABLE]

[TABLE]

Eq. (A-71), by taking $s=\omega_{s}/c$, could be rewritten as:

[TABLE]

Therefore, radiative decay rate, using Eq. (A-82) in Eq. (1.133), could be found as:

[TABLE]

which leads to:

[TABLE]

Using frequency given in Eq. (A-36), we can write:

[TABLE]

where:

[TABLE]

Knowing:

[TABLE]

and,

[TABLE]

hence:

[TABLE]

It simplifies to:

[TABLE]

Let:

[TABLE]

therefore:

[TABLE]

Looking at the fixed modes $m=0$ and $l=1$, one could see:

[TABLE]

which $l=1$ gives the same result as the one:

[TABLE]

given in [8, p. 52]. Next, using the asymptotic behavior for small argument in spherical Bessel functions:

[TABLE]

hence:

[TABLE]

and,

[TABLE]

which is the same as Eq. (14) in [6].

Radiative Decay on Cylindrical Surfaces

The interaction of surface collective excitation, like surface plasmons, in the vicinity of a cylindrical boundary with external probes has been studied in [11] following the general case argued in [8]. Authors in [11] studied one specific case of an extremely long and thin cylinder for which the $z$-axis coincides with the axis of the cylinder. The case investigated in [8] yet is more general and is considered the first infinite geometry being quantized. Throughout this section, we take a general cylinder and we follow the exact steps we have taken in Chapter I, to drive the results in [8].

Cylindrical coordinate is given by:

[TABLE]

where:

[TABLE]

with scale factors:

[TABLE]

Laplace equation in Cylindrical coordinate [57] $(\rho,\varphi,z)$ is:

[TABLE]

which has infinitely many solutions of the form:

[TABLE]

Applying separation of variable techniques, one obtains the equations:

[TABLE]

replacing $k$ with $ik$ gives:

[TABLE]

Scalar inside and outside potentials of an infinite cylinder with base radius $\rho_{0}$ could be written as:

[TABLE]

Utilizing Heaviside function $\Theta$, we define:

[TABLE]

hence Laplacian becomes:

[TABLE]

Since we are only interested on the surface and since Laplacian vanishes inside and outside, then:

[TABLE]

Using (A-109) and (A-110), we have:

[TABLE]

using Wronskian identity [57]:

[TABLE]

On the one hand, using the relation $h_{\rho}\vec{\nabla}^{2}\Phi=-4\pi\delta(\rho-\rho_{0})\sigma,$ we get:

[TABLE]

therefore:

[TABLE]

On the other hand,

[TABLE]

Putting (A-117) and (A-118) equal and using the orthogonality in $m$ and $\varphi$ and also $q$ and $z$, we get:

[TABLE]

where:

[TABLE]

is known as the frequency in cylindrical domain.

Next, we proceed by the intention of calculating classical energy $E$. for which we need to find potential and kinetic energy. We shall follow the steps provided in Chapter 1, through Eqs. (1.70)–(1.87). Rewriting Eq. (1.74) for the specific cylinder, potential energy is given by:

[TABLE]

Using Eq. (A-109) for inside potential $\Phi_{\text{i}}$ and Eq. (A-116) for surface charge density $\sigma$, one finds:

[TABLE]

as we replaced real-valued potential with its complex conjugate. Using the orthogonality relation:

[TABLE]

It follows that (A-122) can be written as:

[TABLE]

We use frequency relation we have obtained in Eq.  (A-120) to write:

[TABLE]

Similar calculations gives kinetic energy as:

[TABLE]

Therefore, total classic energy is found by adding Eqs. (A-125) and (A-126) as:

[TABLE]

Before moving forward,we note the following symmetry relations for modified Bessel functions given by Eq.  (5.7.10) (see [57], p. 110) as:

[TABLE]

where for integers $\nu=n$, $k_{n}(z)=\lim\limits_{\nu\to n}k_{\nu}(z)$, then:

[TABLE]

hence $k_{n}(z)=k_{-n}(z)$ for $n\in\mathbb{Z}$. Considering inside potential given in Eq. (A-109), the conjugate of potential is hence given by:

[TABLE]

letting $m\to-m$ and $k\to-k$, we have:

[TABLE]

using symmetry relations for modified Bessel functions given in (A-128) [57], we may write:

[TABLE]

since potential is real-valued, then $\Phi_{\text{i}}(\mathbf{r},t)=\overline{\Phi_{\text{i}}}(\mathbf{r},t)$, implies:

[TABLE]

for all $m$ and $k$. The complex coefficients $B_{mk}(t)$ could be written as:

[TABLE]

where $b_{mk}$ are some complex function (of time) proportional to $e^{-i\omega_{mk}t}$ in which symmetric (in $m$ and $k$) coefficients $\beta_{mk}$ would be determined later. We could follow the identities:

[TABLE]

noting that according to Eq. (A-133), since $\omega_{mk}=\omega_{-m-k}$, we could replace one with another. The time derivative of coefficients given in Eq. (A-134) gives:

[TABLE]

then we follow similar steps to write:

[TABLE]

Therefore:

[TABLE]

and,

[TABLE]

Comparing Eq. (A-145) with the quantized Hamiltonian in terms of operators as:

[TABLE]

where $\hat{b}_{mk}$ and $\hat{b}^{\dagger}_{mk}$ are annihilation and creation operators respectively. We find $\beta_{mk}$ as:

[TABLE]

Using Eq. (A-147), Eq. (A-134) could be rewritten as:

[TABLE]

Replacing Eq. (A-148) in Eqs.  (A-109) and (A-110), we can write the quantized inside and outside potentials as:

[TABLE]

where:

[TABLE]

also see [11]. Carrying from integral to sum, over a finite box $\mathcal{B}$ with side $L$, we have:

[TABLE]

with one-dimensional wave vector $k$ and:

[TABLE]

Next, we continue by calculating the current given by $\mathbf{J}=-n_{0}\,e\,\dot{\vec{u}}$. For the charge accerelation vector $\dot{\vec{u}}$, we use the relation given in Eq. (1.81) to write:

[TABLE]

using:

[TABLE]

we find the current as:

[TABLE]

We may use the bulk frequency relation, given in Eq. (1.18) to write:

[TABLE]

Here if we follow [8, 49, 42, 43], every formula remains correct if we carry the substitutions:

[TABLE]

using the normalization factor $\frac{1}{\sqrt{\mathcal{V}}}$. Using this in Eq. (A-157), gives:

[TABLE]

For $b_{0}=1$, it gives the Eq.  (70) on [8], p. 134, as:

[TABLE]

Considering Eq (1.128) in case for a cylinder:

[TABLE]

we could find the interaction Hamiltonian as:

[TABLE]

where:

[TABLE]

is the volume of the finite cylindrical box with length $L$. The emission matrix element can be calculate as:

[TABLE]

where the integrals $\mathcal{I}_{mk}$ are given as:

[TABLE]

Taking:

[TABLE]

allows us to take $\hat{\mathbf{e}}_{q}$, for $q=1,2$ as:

[TABLE]

Leaving us with two different integrals:

[TABLE]

and,

[TABLE]

where:

[TABLE]

Radiative decay rate per solid angel $\Omega$, after using Eqs. (A-147) and (A-165) in Eq. (1.137), takes the form:

[TABLE]

The analytic solution for $\mathcal{I}^{(1)}_{mk}$ and $\mathcal{I}^{(2)}_{mk}$ are given as following:

[TABLE]

and,

[TABLE]

where:

[TABLE]

APPENDIX B : COLLECTED FACTS

Collected Facts about Coordinate Systems

An Arbitrary Coordinate System

Let $(\hat{i},\hat{j},\hat{k})$ denote the Cartesian unit vectors and $(\alpha,\beta,\gamma)$ be an arbitrary coordinate system. Let $\vec{r}$ be the radius vector in $(\alpha,\beta,\gamma)$-coordinate. Lastly, let $h_{\alpha},h_{\beta}$ and $h_{\gamma}$ be the scale factors given by:

[TABLE]

with ${\vec{r}}=x\hat{i}+y\hat{j}+z\hat{k}$ with $(\hat{i},\hat{j},\hat{k})$ the Cartesian unit vectors. Let us assume the surface of revolution is defined by fixing $\zeta=\zeta_{0}$, then unit vectors for this coordinate system are defined by [63] :

[TABLE]

Laplace equation in an arbitrary coordinate system is given by:

[TABLE]

Here are some collected facts about the specific coordinate systems we have considered throughout this manuscript.

Spherical Coordinate

[TABLE]

with

[TABLE]

The scale factors for a parametrization in each component are given by [57]:

[TABLE]

Cylindrical Coordinate

[TABLE]

where

[TABLE]

with scale factors

[TABLE]

Parabaloidal Coordinate

are related to the rectangular coordinates by

[TABLE]

with the corresponding scale factors

[TABLE]

where $0\leq\varphi<2\pi$ denotes the usual azimuthal angle and the two coordinates $\eta,\xi\geq 0$ are such that the surfaces of constant $\eta>0$ and $\xi>0$ describe upward and downward paraboloids of revolution about the $z$-axis.

Prolate Spheroidal Coordinate

[TABLE]

in the domain defined by :

[TABLE]

or alternatively,

[TABLE]

with the corresponding scale factors:

[TABLE]

Toroidal Coordinate

[TABLE]

where

[TABLE]

The associated toroidal scale factors are given by

[TABLE]

Collected Facts about Special Functions

Bessel Functions

Asymptotic behaviour for large arguments and fixed order $m$:

For large arguments, the series expansions are Eqs. (5.11.6), (5.11.8) and (5.11.9) [57], p.122, 123.

[TABLE]

[TABLE]

[TABLE]

where $|argz|\leq\pi-\delta$ and

[TABLE]

In general [57, 62]:

[TABLE]

The Bessel functions of the first kind:

[TABLE]

The Bessel functions of the second kind:

[TABLE]

Bessel functions of imaginary argument

Also called the modified Bessel functions of first, $I_{\nu}(z)$ and second, $K_{\nu}(z)$, kind, respectively. Their relation with Bessel functions of first and second kind are given as [57, 62]:

[TABLE]

The asymtotic behaviour for modified Bessel functions for a fixed argument $z$ and large complex order $\nu$ when $\nu\to\infty$ are given by [74]:

[TABLE]

Associated Legendre Functions

Here are some collected facts and indentities about Associated Legendre Functions [57, 74, 69]. The associated Legendre functions are symmetric in $n$:

[TABLE]

Following relations hold upon making $m\to-m$:

[TABLE]

Conical Functions

[90, 60] Integral representation (Mehler integrals) of conical functions:

[TABLE]

with $Z^{m}_{q}$ given in Eq.  (2.72). Moreover,

[TABLE]

Successive relations for Conical functions:

[TABLE]

Asymptotic representations for fixed $0<\xi\leq\alpha<\infty$ and $0<\theta\pi$ can be derived from the asymptotic form of Lagendre functions as [91]:

[TABLE]

[TABLE]

Sturm-Liouville Problem

A Sturm-Liouville system is defined by the self-adjoint operator $\mathcal{L}$ as:

[TABLE]

where $p>0$ and $q$ are two continuous functions on an interval $[a,b]$, (see [92, al2008sturm]). The regular Sturm-Liouville eigenvalue problem has the form:

[TABLE]

subject to the separated homogeneous boundary conditions :

[TABLE]

where such values of $\lambda$, if any, are called the eigenvalues of the system and $\rho(x)$ is a weight function. Since $\mathcal{L}$ is a self-adjoint operator under the boundary conditions given in Eqs. (B-64), eigenvalues of Eq. (B-63) are real and can be ordered in an ascending manner. To each eigenvalue $\lambda_{n}$, there corresponds a unique eigenfunction $u_{n}(x)$ which is considered as the $n$-th fundamental solution to Eq. (B-63). The set of all these fundamental solutions $\{u_{n}(x)\}_{n=0}^{\infty}$, up to a normalization constant, form an orthonormal basis in $L^{2}([a,b],\rho(x)dx)$ as they satisfy the following orthogonality relation:

[TABLE]

The spectrum of an operator, e.g. $\mathcal{L}$, is a set consists of all $\lambda$ (in $\mathbb{C}$) such that the operator (eigenvalue problem) $\lambda I-\mathcal{L}$ fails to be an invertible element of the respective algebra. The spectrum can be partitioned to many disjoint parts one of which is called the point spectrum and consists of all eigenvalues $\lambda_{n}$. Eigenvalue problems over finite intervals lead to discrete spectrum whereas infinite intervals result in a continuous spectrum. Among the geometries investigated here, sphere, prolate spheroid and torus are examples in which the spectrum space is discrete and paraboloid, hyperboloid and infinite cylinder are examples in which a continuous spectrum is assumed [41, 63, 94, 93].

APPENDIX C : COPYRIGHT CLEARANCE FORM

Below is the permission for Chapters 2 and 5.

ABOUT THE AUTHOR

Maryam Bagherian received her Bachelor of Science in pure mathematics from the Shahid Beheshti University in 2006. She started to pursue the Doctor of Philosophy in Applied Mathematics at the University of South Florida in August 2014. In between, she earned a Master of Arts in applied mathematics with mathematical physics concentration in 2016. Her research interests include nanoparticles, modeling and simulation and signal/image processing. Some of her non-academic pursuits are yoga, classical music, photography and traveling.