Impact of Graphene on the Polarizability of a Neighbour Nanoparticle: A Dyadic Green’s Function Study

Abstract

We discuss the renormalization of the polarizability of a nanoparticle in the presence of either: (1) a continuous graphene sheet; or (2) a plasmonic graphene grating, taking into account retardation effects. Our analysis demonstrates that the excitation of surface plasmon polaritons in graphene produces a large enhancement of the real and imaginary parts of the renormalized polarizability. We show that the imaginary part can be changed by a factor of up to 100 relative to its value in the absence of graphene. We also show that the resonance in the case of the grating is narrower than in the continuous sheet. In the case of the grating it is shown that the resonance can be tuned by changing the grating geometric parameters.

1. Introduction

The polarizability of a nanoparticle is a response function which relates the electric dipole moment produced in it to an externally applied electric field. The polarizability is not an intrinsic property of the nanoparticle, but actually depends on the environment which it is embedded in [1,2,3]. As such, a nanoparticle’s polarizability will be modified by the presence of an underlying substrate. The study of this problem is of significant interest, since in most experimental setups the nanoparticle (NP) is placed directly onto a dielectric substrate or at a given distance from it. In previous studies in which the radiation scattered by a dielectric NP was measured using dark-field microscopy, it was shown that the presence of the substrate leads to a redshift of the NP’s resonance with respect to the situation where the NP is in vacuum [4,5,6].

The polarizability of a nanoparticle at a given frequency is a complex quantity, with its real and imaginary parts describing, respectively, the reactive and dissipative responses of a nanoparticle subjected to an electromagnetic field. Therefore, the imaginary part of the polarizability controls the extinction and absorption cross-sections of a nanoparticle subjected to an impinging electromagnetic field [7,8]. These quantities are essential for the understanding of scattering experiments of electromagnetic radiation involving nanoparticles, either isolated or forming clusters. In particular, the former case has been a topic of much interest in the context of single-molecule or single-particle spectroscopies [9,10]. The knowledge of the imaginary part of the polarizability is also essential in order to understand the phenomena of blackbody and thermal friction experienced by a neutral nanoparticle in close proximity to an interface between two media [11]. It is therefore of major importance to understand how the imaginary part of the polarizability is renormalized relative to its value in a vacuum when it is near an interface, the most common setup in experiments.

It would be of particular relevance (from the device engineering viewpoint), if the dielectric properties of the interface, near which the nanoparticle is located, could be tuned. This would provide a route for controlling the value of the nanoparticle polarizability in real time. Such an approach is not viable when we consider the interface between two conventional dielectrics or between a metal and a dielectric, since the interface has fixed properties by construction. Fortunately, there is a possible and technologically feasible route to overcome this limitation. Adding a graphene sheet between an interface involving two different dielectrics provides an additional degree of freedom to the problem. Indeed, the Fermi energy of a graphene sheet can be controlled in real time using a gate. Tuning the Fermi energy of graphene changes the local dielectric environment around the nanoparticle and therefore the value of the imaginary part of the nanoparticle polarizability. This is the opportunity we will explore in this paper.

Incidentally, the problem of nanoparticle polarizability renormalization in the presence of a substrate is also relevant for the characterization of the dielectric properties of a scanning near-field optical microscope (SNOM). Scanning near-field optical microscopy is a technique frequently used to image and characterize surface polaritons in graphene [12] and other two-dimensional materials, such as boron nitride [13]. More recently, exciton-polaritons have also been studied in layered transition metal dichalcogenides using the same method [14]. Indeed, the SNOM tip is frequently modeled as a dipole, as is the nanoparticle [15]. Therefore, understanding how a nanoparticle changes its dielectric properties under illumination allows us to also understand the problem of a SNOM tip illuminated with THz radiation during the excitation of surface polaritons in graphene and other two-dimensional materials.

In this work, we study the renormalization of a nanoparticle polarizability located near the interface between two dielectrics interspaced with a doped graphene sheet, or with an array of graphene ribbons (see Figure 1). One of the dielectrics is the vacuum and the other acts as substrate for the support of the graphene sheet. In order to keep the analysis simple we shall restrict ourselves to the case of a non-dispersive and non-dissipative substrate, characterized by a frequency-independent and real dielectric constant. We explore the imaginary part of the polarizability in the THz range of the electromagnetic spectrum, a spectral region where graphene supports surface plasmon polaritons [16,17,18]. As we will see, the excitation of these polaritons leads to a significant change of the polarizability of both metallic and semiconductor nanoparticles. Indeed, the bare polarizability of a metallic nanoparticle in vacuum is essentially constant in the THz, with a very small imaginary part. However, when located near a graphene sheet, the polarizability undergoes strong renormalization, especially with respect to its imaginary part.

Although the problem of modeling the polarizability of a nanoparticle close to a graphene sheet has been considered before by some of the authors of the present paper [19], that work relied on a electrostatic approximation. The present work goes beyond that, taking into account retardation effects, allowing us to correctly describe the imaginary part of the polarizability. It should be noted that the problem of determining the nanoparticle’s polarizability in the presence of a homogeneous flat dielectric substrate has also been considered previously both in the electrostatic approximation [20] and in the full electrodynamic approach [6,21,22].

The goals of this work are as follows: (1) to extend the study of [19] including retardation effects, thus using a more general formalism; (2) to bring together in a single paper a formalism that is scattered in the literature using many different notations; (3) to introduce a rigorous formulation of the dyadic Green’s function formalism that is absent in many papers; and (4) to extend this approach to the case where a nanoparticle has both dipolar electric and dipolar magnetic moments.

This paper is organized as follows: in Section 2 we introduce the concept of dyadic Green’s function for the electric field as a tool to obtain the electric field in the presence of source currents. In Section 2.1 we study in detail the electric field dyadic Green’s function in free-space (or in a homogeneous medium). The Weyl’s, or angular spectrum, representation of the dyadic Green’s function is introduced in Section 2.2. This representation is well-suited to deal with the problem of radiation scattering at planar interfaces. It is also shown that the dyadic Green’s functions can be expressed in terms of the tensor product of the electric field s- and p-polarization vectors. In Section 2.3, we focus on the problem of scattering at a planar interface and define the reflected and transmitted Green’s functions. In Section 3 we deduce the polarizability of a nanoparticle close to an interface covered by graphene. We start defining and studying the polarizability of a nanoparticle embedded in a vacuum, in Section 3.1. The approach is generalized in Section 3.2 to the case of a nanoparticle close to a planar interface. This general description is then used to describe the renormalization of a nanoparticle’s polarizability close to a continuous graphene sheet and to a graphene grating in Section 3.3 and Section 3.4. In Section 4 we present a generalization of the formalism to the case where the nanoparticle has both electric and magnetic dipole moments. Such a magnetic moment can be generated, even for nanoparticles formed by a non-magnetic material, due to induced currents inside the nanoparticle [23], and can actually be the main contribution for the polarizability in the case of dielectric NPs [4,5,6]. Finally, a set of appendices provides some auxiliary results.

2. Dyadic Green’s Function for the Electric Field

2.1. Free-Space Dyadic Green’s Function

The goal of this section is to introduce the necessary tools to study the polarizability of a nanoparticle in the presence of substrate. In particular, we will introduce the electric field dyadic Green’s function, which allows us to solve the wave equation for the electric field in the presence of source currents. Although the material in this section is relatively well-known, we present it here in some detail both for the sake of completeness and to fix notations used throughout the paper.

In frequency space, the dyadic Green’s function, $\overleftrightarrow{G_0}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )$, allows us to express the total electric field, $\mathbf{E}(\mathbf{r},\omega )$, in the presence of a source current, ${\mathbf{j}}_f(\mathbf{r},\omega )$, as

$$\mathbf{E}(\mathbf{r},\omega )={\mathbf{E}}_0(\mathbf{r},\omega )+i\omega {\mu }_n{\mu }_0\int d^3{\mathbf{r}}^{\prime }\overleftrightarrow{G_0}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )·{\mathbf{j}}_f({\mathbf{r}}^{\prime },\omega ),$$

(1)

where ${\mathbf{E}}_0(\mathbf{r},\omega )$ is the electric field in the absence of the source current, ${\mu }_0$ is vacuum’s permeability, and ${\mu }_n$ is the relative permeability of the medium where the source current is embedded in. If the source current is due to a point dipole located at $\mathbf{r}={\mathbf{r}}_0$, we have that ${\mathbf{j}}_f(\mathbf{r},\omega )=-i\omega {\mathbf{d}}_0\delta \left(\mathbf{r}-{\mathbf{r}}_0\right)$, where ${\mathbf{d}}_0$ is the electric dipole moment. In this case, Equation (1) reduces to (for $\mathbf{r}\ne {\mathbf{r}}_0$)

$$\mathbf{E}(\mathbf{r},\omega )={\mathbf{E}}_0(\mathbf{r},\omega )+{\omega }^2{\mu }_n{\mu }_0{\overleftrightarrow{G}}_0(\mathbf{r},{\mathbf{r}}_0,\omega )·{\mathbf{d}}_0.$$

(2)

In order to determine the dyadic Green’s function we will follow a method originally described in [24]. The first step to determining $\overleftrightarrow{G_0}$ in this approach is noticing that in the presence of a source current, the wave equation for the electric field can be written as an inhomogeneous Helmholtz equation (see Appendix A for a derivation) as follows:

$$-{\nabla }^2\mathbf{E}(\mathbf{r},\omega )-k_n^2\mathbf{E}(\mathbf{r},\omega )=i\omega {\mu }_n{\mu }_0\left[{\mathbf{j}}_f(\mathbf{r},\omega )+\frac{1}{k_n^2}\nabla \left(\nabla ·{\mathbf{j}}_f(\mathbf{r},\omega )\right)\right],$$

(3)

where we have introduce the quantity $k_n=\omega /v_n$, where $v_n=1/\sqrt{{ϵ}_0{ϵ}_n{\mu }_0{\mu }_n}$ is the speed of light in the medium, ${ϵ}_0$ is the vacuum’s permittivity, and ${ϵ}_n$ is the medium relative dielectric constant. The general solution of the Helmholtz equation can be written as (see Appendix B)

$$\mathbf{E}(\mathbf{r},\omega )={\mathbf{E}}_0(\mathbf{r},\omega )+i\omega {\mu }_n{\mu }_0{\int }_{\backslash V_{\delta }\left(\mathbf{r}\right)}{d}^3{\mathbf{r}}^{\prime }{g}_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)\left[\overleftrightarrow{I}+\frac{1}{k_n^2}{\nabla }^{\prime }{\nabla }^{\prime }\right]{\mathbf{j}}_f({\mathbf{r}}^{\prime },\omega ),$$

(4)

where

$$g_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)=\frac{e^{i{k}_n\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|}}{4\pi \left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|},$$

(5)

is the Green’s function for the scalar Helmholtz equation [1,24,25] and ${\int }_{\backslash V_{\delta }\left(\mathbf{r}\right)}$ represents integration in the principal value sense, where an infinitesimal volume, $V_{\delta }\left(\mathbf{r}\right)$, enclosing the point ${\mathbf{r}}^{\prime }=\mathbf{r}$ is excluded. We have written ${\nabla }^{\prime }{\nabla }^{\prime }\equiv {\nabla }^{\prime }\otimes {\nabla }^{\prime }$ with ⊗ denoting the tensor product and the prime indicates that the derivative acts on ${\mathbf{r}}^{\prime }$ variables. The Helmholtz Green’s function is the solution of

$$\left[-{\nabla }^2-k_n^2\right]g_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)=\delta (\mathbf{r}-{\mathbf{r}}^{\prime })$$

(6)

in a way that is clarified in Appendix B. Notice that $g_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)$ is integrable, and therefore, the exclusion of the volume $V_{\delta }\left(\mathbf{r}\right)$ is not usually emphasized. However, it will be important when determining the behaviour of $\overleftrightarrow{G_0}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )$ when $\mathbf{r}={\mathbf{r}}^{\prime }$. Although Equation (4) already allows us to compute the electric field as a function of the current, it is useful to obtain an alternative expression which does not involve derivatives of the current. Such an expression can be obtained by carefully performing integration by parts. It must the noticed that due to the excluded volume surrounding ${\mathbf{r}}^{\prime }=\mathbf{r}$, boundary terms are generated when doing this. In particular we have that

$$\begin{array}{cc}\hfill {\int }_{\backslash V_{\delta }\left(\mathbf{r}\right)}{d}^3{\mathbf{r}}^{\prime }{g}_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right){\nabla }^{\prime }\left({\nabla }^{\prime }·{\mathbf{j}}_f({\mathbf{r}}^{\prime },\omega )\right)=& -{\int }_{\partial V_{\delta }\left(\mathbf{r}\right)}{d}^2{\mathbf{r}}^{\prime }{\mathbf{n}}^{\prime }\left[g_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)\left({\nabla }^{\prime }·{\mathbf{j}}_f({\mathbf{r}}^{\prime },\omega )\right)\right]\hfill \\ & -{\int }_{\backslash V_{\delta }\left(\mathbf{r}\right)}{d}^3{\mathbf{r}}^{\prime }{\nabla }^{\prime }{g}_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)\left({\nabla }^{\prime }·{\mathbf{j}}_f({\mathbf{r}}^{\prime },\omega )\right),\hfill \end{array}$$

(7)

where ${\mathbf{n}}^{\prime }$ is a outward-pointing unit vector, normal to the surface $\partial V_{\delta }\left(\mathbf{r}\right)$ of the enclosing volume $V_{\delta }\left(\mathbf{r}\right)$. In the limit of infinitesimal excluded volume, the first term of the above equation vanishes, since the element of area scales as $d^2{\mathbf{r}}^{\prime }\sim {\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|}^2$, while $g_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)\sim 1/\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|$. For the second term, we perform integration by parts once again. For clarity we explicitly write the tensorial components in a Cartesian basis, with repeated induces being summed over, and obtain

$$\begin{array}{cc}\hfill -{\int }_{\backslash V_{\delta }\left(\mathbf{r}\right)}{d}^3{\mathbf{r}}^{\prime }{\partial }_i^{\prime }{g}_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right){\partial }_k^{\prime }{j}_f^k({\mathbf{r}}^{\prime },\omega )=& {\int }_{\partial V_{\delta }\left(\mathbf{r}\right)}{d}^2{\mathbf{r}}^{\prime }{n}_k^{\prime }\left[{\partial }_i^{\prime }{g}_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)j_f^k({\mathbf{r}}^{\prime },\omega )\right]\hfill \\ & +{\int }_{\backslash V_{\delta }\left(\mathbf{r}\right)}{d}^3{\mathbf{r}}^{\prime }\left[{\partial }_i^{\prime }{\partial }_k^{\prime }{g}_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)j_f^k({\mathbf{r}}^{\prime },\omega )\right].\hfill \end{array}$$

(8)

Now the boundary term is finite. In the limit of an infinitesimal volume, we take ${\mathbf{r}}^{\prime }\to \mathbf{r}$, such that $j_f^k({\mathbf{r}}^{\prime },\omega )\to j_f^k(\mathbf{r},\omega )$ and use the $\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|\to 0$ limit of ${\partial }_i^{\prime }{g}_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)$ Equation (A26). This allows us to write

$$\begin{array}{c}\hfill {\int }_{\backslash V_{\delta }\left(\mathbf{r}\right)}{d}^3{\mathbf{r}}^{\prime }{g}_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right){\nabla }^{\prime }\left({\nabla }^{\prime }·{\mathbf{j}}_f({\mathbf{r}}^{\prime },\omega )\right)={\int }_{\backslash V_{\delta }\left(\mathbf{r}\right)}{d}^3{\mathbf{r}}^{\prime }{\nabla }^{\prime }{\nabla }^{\prime }{g}_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)j_f^k({\mathbf{r}}^{\prime },\omega )-\frac{1}{k_n^2}{\overleftrightarrow{L}}_{V_{\delta }}·{\mathbf{j}}_f(\mathbf{r},\omega )\end{array}$$

where the dyadic ${\overleftrightarrow{L}}_{V_{\delta }}$ is defined as [24]

$${\overleftrightarrow{L}}_{V_{\delta }}={\int }_{\partial V_{\delta }\left(\mathbf{r}\right)}\frac{d^2{\mathbf{r}}^{\prime }}{4\pi }\frac{\left({\mathbf{r}}^{\prime }-\mathbf{r}\right)\otimes {\mathbf{n}}^{\prime }}{{\left|{\mathbf{r}}^{\prime }-\mathbf{r}\right|}^3},$$

(9)

which can be interpreted as a depolarization term. Therefore, Equation (4) can be written in the form of Equation (1) with the dyadic Green’s function, $\overleftrightarrow{G_0}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )$, given by

$$\overleftrightarrow{G_0}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )={\mathrm{P}.\mathrm{V}.}_{V_{\delta }}\left[\overleftrightarrow{I}+\frac{1}{k_n^2}\nabla \nabla \right]g_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)-\frac{1}{k_n^2}{\overleftrightarrow{L}}_{V_{\delta }}\delta \left(\mathbf{r}-{\mathbf{r}}^{\prime }\right),$$

(10)

where ${\mathrm{P}.\mathrm{V}.}_{V_{\delta }}$ indicates that the small volume $V_{\delta }$ centered at ${\mathbf{r}}^{\prime }=\mathbf{r}$ is to be excluded. Notice that ${\overleftrightarrow{L}}_{V_{\delta }}$ depends on the shape of the chosen excluded volume [24]. For a sphere, it is straightforward to show that ${\overleftrightarrow{L}}_{{\mathrm{Sphere}}_{\delta }}=\overleftrightarrow{I}/3$. In this case, the free-space dyadic Green’s function in real space can be written as the sum of four terms [26,27]:

$$\overleftrightarrow{G_0}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )={\overleftrightarrow{G_0}}^{\mathrm{FF}}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )+{\overleftrightarrow{G_0}}^{\mathrm{IF}}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )+{\overleftrightarrow{G_0}}^{\mathrm{NF}}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )+{\overleftrightarrow{G_0}}^{\mathrm{SF}}(\mathbf{r},{\mathbf{r}}^{\prime },\omega ),$$

(11)

respectively, the far-, intermediate-, near- and self-field terms, which are written as

$$\begin{array}{cc}\hfill {\overleftrightarrow{G_0}}^{\mathrm{FF}}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )& =\left(\overleftrightarrow{I}-\widehat{\mathbf{R}}\widehat{\mathbf{R}}\right)\frac{e^{i{k}_n|\mathbf{r}-{\mathbf{r}}^{\prime }|}}{4\pi |\mathbf{r}-{\mathbf{r}}^{\prime }|},\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill {\overleftrightarrow{G_0}}^{\mathrm{IF}}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )& =i\left(\overleftrightarrow{I}-3\widehat{\mathbf{R}}\widehat{\mathbf{R}}\right)\frac{e^{i{k}_n|\mathbf{r}-{\mathbf{r}}^{\prime }|}}{4\pi k_n{|\mathbf{r}-{\mathbf{r}}^{\prime }|}^2},\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill {\overleftrightarrow{G_0}}^{\mathrm{NF}}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )& =-\left(\overleftrightarrow{I}-3\widehat{\mathbf{R}}\widehat{\mathbf{R}}\right)\frac{e^{i{k}_n|\mathbf{r}-{\mathbf{r}}^{\prime }|}}{4\pi k_n^2{|\mathbf{r}-{\mathbf{r}}^{\prime }|}^3}\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill {\overleftrightarrow{G_0}}^{\mathrm{SF}}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )& =-\overleftrightarrow{I}\frac{1}{3k_n^2}\delta (\mathbf{r}-{\mathbf{r}}^{\prime }),\hfill \end{array}$$

(15)

where the terms ${\overleftrightarrow{G_0}}^{\mathrm{FF}}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )$, ${\overleftrightarrow{G_0}}^{\mathrm{IF}}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )$ and ${\overleftrightarrow{G_0}}^{\mathrm{NF}}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )$ are to be understood in the principal value sense, and we have introduced the definitions $\widehat{\mathbf{R}}=\left(\mathbf{r}-{\mathbf{r}}^{\prime }\right)/\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|$ and $\widehat{\mathbf{R}}\widehat{\mathbf{R}}=\widehat{\mathbf{R}}\otimes \widehat{\mathbf{R}}$.

We point out that in more standard derivations of the dyadic Green’s function it is easy to miss the depolarization term Equation (9) [1,28], or fail to recognize its dependence on the shape of the excluded volume [29,30]. It will be shown in Section 3 that this term is crucial to correctly capture self-field contributions to a particle polarizability. The dependence of depolarization term on the shape of the volume is also important when dealing with non-spherical particles. An alternative derivation of the depolarization term Equation (9) based on the vector potential can be found in [31].

2.2. Weyl’s or Angular Spectrum Representation of the Dyadic Green’s Function: A Useful Formulation for Interfaces

Although Equation (10) can be used directly to evaluate $\overleftrightarrow{G_0}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )$, for many applications such formulation might not be the most useful. In the the case of scattering by planar interfaces it is useful to make a (two-dimensional) Fourier transform of the fields in the coordinates parallel to the interface. This representation of the fields and of the Green’s function is generally referred to as Weyl’s or angular spectrum representation. In this representation, the electric field is written as

$$\mathbf{E}(\mathbf{r},\omega )=\int \frac{d^2{\mathbf{p}}_{\parallel }}{{\left(2\pi \right)}^2}{e}^{i{\mathbf{p}}_{\parallel }·\mathit{\rho }}\mathbf{E}({\mathbf{p}}_{\parallel },z,\omega ),$$

(16)

where ${\mathbf{p}}_{\parallel }$ is the in-plane wave-vector and $\mathit{\rho }=\left(x,y\right)$ are in-plane coordinates. In this representation Equation (4) becomes

$$\mathbf{E}({\mathbf{p}}_{\parallel },z,\omega )={\mathbf{E}}_0({\mathbf{p}}_{\parallel },z,\omega )+i\omega {\mu }_n{\mu }_0{\displaystyle ⨏}d{z}^{\prime }{g}_0\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right)\left[\overleftrightarrow{I}-\frac{1}{k_n^2}{𝒟}^{\prime }{𝒟}^{\prime }\right]{\mathbf{j}}_f({\mathbf{p}}_{\parallel },z^{\prime },\omega ),$$

(17)

where ${\mathbf{j}}_f({\mathbf{p}}_{\parallel },z,\omega )$ is the Weyl representation of the current density, defined in analogous way to Equation (16), ${𝒟}^{\prime }={\mathbf{p}}_{\parallel }-i{\widehat{e}}_z{\partial }_z^{\prime }$, ${𝒟}^{\prime }{𝒟}^{\prime }\equiv {𝒟}^{\prime }\otimes {𝒟}^{\prime }$, ⨏ represents the principal value integral in one dimension, excluding the point $z^{\prime }=z$, and $g_0\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right)$ is the Helmholtz Green’s function in the Weyl representation, defined such that

$$g_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)=\int \frac{d^2{\mathbf{p}}_{\parallel }}{{\left(2\pi \right)}^2}{e}^{i{\mathbf{p}}_{\parallel }·\left(\mathit{\rho }-{\mathit{\rho }}^{\prime }\right)}{g}_0\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right).$$

(18)

The function $g_0({\mathbf{p}}_{\parallel },z,z^{\prime },\omega )$ can be easily obtained from the components of the three-dimensional Fourier transform of the Helmholtz Green’s function, $g_0\left(\mathbf{p},\omega \right)={\left(p_{\parallel }^2+p_z^2-k_n^2\right)}^{-1}$, as

$$g_0\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right)=\int \frac{d{p}_z}{2\pi }{e}^{i{p}_z\left(z-z^{\prime }\right)}{g}_0\left(\mathbf{p},\omega \right).$$

(19)

This integral can be easily performed by contour integration yielding

$$g_0\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right)=\frac{i}{2{\beta }_n}{e}^{i{\beta }_n\left|z-z^{\prime }\right|},$$

(20)

where ${\beta }_n$ is defined as

$${\beta }_n=\left\{\begin{array}{cc}\sqrt{k_n^2-p_{\parallel }^2},\hfill & k_n^2>p_{\parallel }^2\hfill \\ i\sqrt{p_{\parallel }^2-k_n^2},\hfill & k_n^2<p_{\parallel }^2\hfill \end{array}\right..$$

(21)

Clearly, Equation (20) is written is terms of both propagating and evanescent waves [32]. Similarly to what we have done in the previous section, we can rewrite Equation (17) by moving the derivatives ${\partial }_z^{\prime }$ from ${\mathbf{j}}_f({\mathbf{p}}_{\parallel },z^{\prime },\omega )$ to $g_0\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right)$. Doing this yields

$$\mathbf{E}({\mathbf{p}}_{\parallel },z,\omega )={\mathbf{E}}_0({\mathbf{p}}_{\parallel },z,\omega )+i\omega {\mu }_n{\mu }_0\int d{z}^{\prime }\overleftrightarrow{G_0}\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right){\mathbf{j}}_f({\mathbf{p}}_{\parallel },z^{\prime },\omega ),$$

(22)

with $\overleftrightarrow{G_0}\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right)$ being the dyadic Green’s function in Weyl’s representation

$$\overleftrightarrow{G_0}\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right)=\mathrm{P}.\mathrm{V}.\left[\overleftrightarrow{I}-\frac{1}{k_n^2}{\mathbf{p}}_n^{\pm }{\mathbf{p}}_n^{\pm }\right]\frac{i}{2{\beta }_n}{e}^{i{\beta }_n\left|z-z^{\prime }\right|}-\frac{1}{k_n^2}{\widehat{e}}_z{\widehat{e}}_z\delta \left(z-z^{\prime }\right),$$

(23)

where we have introduced ${\mathbf{p}}_n^{\pm }={\mathbf{p}}_{\parallel }\pm {\beta }_n{\widehat{e}}_z$, with the ± sign applying for $z\gtrless z^{\prime }$. The last term in the above equation is the depolarization term, that arises from the boundary contributions when performing integration by parts, due to the exclusion of an infinitesimal line element around $z^{\prime }=z$ in ⨏. The principal value in the first term indicates that a small region around $z^{\prime }=z$ is to be excluded. We also notice that this depolarization term could also have been obtained from the general depolarization dyadic in real space, Equation (9), if we choose as excluded volume an infinite slab located at $-\delta <z<\delta$ (with $\delta \to 0$). For this excluded volume, we would obtain ${\overleftrightarrow{L}}_{{\mathrm{Slab}}_{\delta }}={\widehat{e}}_z{\widehat{e}}_z$.

It is possible to write Equation (23) in a more meaningful way by introducing the s- and p-polarization vectors. The s-polarization vector lies in the $xy-$plane and is therefore written as [33]

$${\widehat{e}}_s=\frac{{\mathbf{p}}_{\parallel }\times {\widehat{e}}_z}{p_{\parallel }}.$$

(24)

On the other hand, the p-polarization vector is orthogonal to ${\mathbf{p}}_n^{\pm }$ and ${\widehat{e}}_s$, and therefore we write it as [33]

$${\widehat{e}}_{p,n}^{\pm }=\frac{{\widehat{e}}_s\times {\mathbf{p}}_n^{\pm }}{k_n}=\frac{p_{\parallel }}{k_n}{\widehat{e}}_z\mp \frac{{\beta }_n}{k_n}\frac{{\mathbf{p}}_{\parallel }}{p_{\parallel }},$$

(25)

where ${\widehat{e}}_{p,n}^{\pm }$ is the p-polarization vector for a field propagating in the positive/negative $z-$direction. With these definitions we have that $\overleftrightarrow{I}-{\mathbf{p}}_n^{\pm }{\mathbf{p}}_n^{\pm }/k_n^2={\widehat{e}}_s{\widehat{e}}_s+{\widehat{e}}_{p,n}^{\pm }{\widehat{e}}_{p,n}^{\pm }$ and, therefore, we can write Equation (23) as [34,35]

$$\overleftrightarrow{G_0}\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right)={\widehat{e}}_s{\widehat{e}}_s\frac{i}{2{\beta }_n}{e}^{i{\beta }_n\left|z-z^{\prime }\right|}+{\widehat{e}}_{p,n}^{\pm }{\widehat{e}}_{p,n}^{\pm }\frac{i}{2{\beta }_n}{e}^{i{\beta }_n\left|z-z^{\prime }\right|}-\frac{1}{k_n^2}{\widehat{e}}_z{\widehat{e}}_z\delta \left(z-z^{\prime }\right),$$

(26)

with the first and the second terms corresponding to the s- and p-polarization components of the free-space dyadic Green’s function, respectively. A different derivation of previous two equations has been given in the literature previously [36,37,38]. The same decomposition has been used in the study of an emitter’s life-time near a graphene sheet [39,40] and in the context of the calculation of the electric field of a dipole near graphene [41].

2.3. Source and Scattered Green’s Functions: Scattering at a Planar Interface

We now want to address the problem of determining the Green’s function in a system with a planar interface between media 1 and 2. To that end, we shall evaluate the electric field generated by a point dipole, characterized by an electric dipole moment ${\mathbf{d}}_0$, located at a distance $z_0>0$ from the interface. We assume that medium 1 is located in the half-space $z>0$, whereas medium 2 is located in the complementary space, as represented in Figure 1. Note that in general ${\beta }_1\ne {\beta }_2$ due to the different values of the speed of light in the media. Assuming that $z\ne z_0$, the primary field emitted by the oscillating dipole in the half-space $z>0$ reads

$${\mathbf{E}}_0\left({\mathbf{p}}_{\parallel },z,z_0,\omega \right)={\mu }_1{\mu }_0{\omega }^2\overleftrightarrow{G_0}\left({\mathbf{p}}_{\parallel },z,z_0,\omega \right)·{\mathbf{d}}_0.$$

(27)

We have two different values for the field, depending on whether $z\gtrless z_0$. Recalling Equation (26) we obtain

$${\mathbf{E}}_0\left({\mathbf{p}}_{\parallel },z\gtrless z_0,\omega \right)=E_{0,s}{e}^{i{\beta }_1\left|z-z_0\right|}{\widehat{e}}_s+E_{0,p}{e}^{i{\beta }_1\left|z-z_0\right|}{\widehat{e}}_{p,1}^{\pm },$$

(28)

with s- and p-polarization amplitudes being given by

$$\begin{array}{cc}\hfill E_{0,s}& ={\mu }_1{\mu }_0{\omega }^2\frac{i}{2{\beta }_1}\left({\widehat{e}}_s·{\mathbf{d}}_0\right),\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill E_{0,p}& ={\mu }_1{\mu }_0{\omega }^2\frac{i}{2{\beta }_1}\left({\widehat{e}}_{p,n}^{\pm }·{\mathbf{d}}_0\right).\hfill \end{array}$$

(30)

This field corresponds to waves impinging on the interface at $z=0$, being partially reflected and partially transmitted as depicted in Figure 2. The reflected and transmitted fields can be expressed in terms of the amplitudes of the incident field at $z=0$ and of the reflection, $r_s$ and $r_p$, and transmission, $t_s$ and $t_p$, coefficients at the interface for the s- and p-polarizations as [33,35]

$$\begin{array}{cc}\hfill {\mathbf{E}}_r\left({\mathbf{p}}_{\parallel },z>0,\omega \right)& =r_s{E}_{0,s}{e}^{i{\beta }_1\left(z+z_0\right)}{\widehat{e}}_s+r_p{E}_{0,p}{e}^{i{\beta }_1\left(z+z_0\right)}{\widehat{e}}_{p,1}^{+},\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill {\mathbf{E}}_t\left({\mathbf{p}}_{\parallel },z<0,\omega \right)& =t_s{E}_{0,s}{e}^{i{\beta }_1{z}_0}{e}^{-i{\beta }_{2}z}{\widehat{e}}_s+t_s{E}_{0,s}{e}^{i{\beta }_1{z}_0}{e}^{-i{\beta }_{2}z}{\widehat{e}}_{p,2}^{-}.\hfill \end{array}$$

(32)

The factor $e^{i{\beta }_{1}z}\left(e^{-i{\beta }_{2}z}\right)$ is acquired by the field while propagating along the positive (negative) z direction in medium 1 (2). The p-polarization vector for the reflected field is ${\widehat{e}}_{p,1}^{+}$ since it propagates along the positive z direction. Conversely, we have ${\widehat{e}}_{p,2}^{-}$ for the transmitted field, since it propagates along the negative z direction.

Therefore, the total field for $z>0$ can be written as

$$\mathbf{E}\left({\mathbf{p}}_{\parallel },z>0,z_0,\omega \right)={\mu }_1{\mu }_0{\omega }^2\left[\overleftrightarrow{G_0}({\mathbf{p}}_{\parallel },z-z_0,\omega )+\overleftrightarrow{G_r}({\mathbf{p}}_{\parallel },z,z_0,\omega )\right]·{\mathbf{d}}_0,$$

(33)

where we have introduced the reflected Green’s function

$$\overleftrightarrow{G_r}\left({\mathbf{p}}_{\parallel },z,z_0,\omega \right)=r_s\frac{i}{2{\beta }_1}{\widehat{e}}_s{\widehat{e}}_s{e}^{i{\beta }_1(z+z_0)}+r_p\frac{i}{2{\beta }_1}{\widehat{e}}_{p,1}^{+}{\widehat{e}}_{p,1}^{-}{e}^{i{\beta }_1(z+z_0)}.$$

(34)

Similarly, the transmitted field for $z<0$ can be written as

$$\mathbf{E}({\mathbf{p}}_{\parallel },z<0,z_0,\omega )={\mu }_1{\mu }_0{\omega }^2\overleftrightarrow{G_t}\left({\mathbf{p}}_{\parallel },z,z_0,\omega \right)·{\mathbf{d}}_0.$$

(35)

with the transmitted Green’s function being written as

$$\overleftrightarrow{G_t}\left({\mathbf{p}}_{\parallel },z,z_0,\omega \right)=t_s\frac{i}{2{\beta }_1}{\widehat{e}}_s{\widehat{e}}_s{e}^{-i{\beta }_{2}z}{e}^{i{\beta }_1{z}_0}+t_p\frac{i}{2{\beta }_1}{\widehat{e}}_{p,2}^{-}{\widehat{e}}_{p,1}^{-}{e}^{-i{\beta }_{2}z}{e}^{i{\beta }_1{z}_0}.$$

(36)

At this point, we have now in our possession all the relevant tools to study the renormalization of the polarizability of a nanoparticle in the vicinity of a planar interface.

3. Renormalization of the Polarizability of a Quantum Emitter Near a Graphene Sheet and a Graphene-Based Grating

The dyadic Green’s function method is a powerful tool for describing the modification of the properties of a quantum emitter near interfaces, as it takes into account the change in the density of electromagnetic modes induced by the presence of the interface. Problems such as the calculation of the Purcell factor and Förster energy transfer are two examples [42,43] well-suited for the Green’s function approach. Here we consider another problem that also depends on the density of electromagnetic modes, the calculation of the effective polarizability of a quantum emitter.

3.1. Polarizability of a Quantum Emitter in a Homogeneous Medium

The polarizability of a nanoparticle, $\overleftrightarrow{\alpha }$, treated as a point object, relates the electric dipole moment, $\mathbf{d}$, that is induced in the nanoparticle to the value of the externally applied electric field, ${\mathbf{E}}_{\mathrm{ext}}\left({\mathbf{r}}_0\right)$, at the nanoparticle’s position, ${\mathbf{r}}_0$, via

$$\mathbf{d}=\overleftrightarrow{\alpha }·{\mathbf{E}}_{\mathrm{ext}}\left({\mathbf{r}}_0\right).$$

(37)

Note that ${\mathbf{E}}_{\mathrm{ext}}\left({\mathbf{r}}_0\right)$ does not include self-field effects, that is, the electric field generated by the nanoparticle itself when subjected to ${\mathbf{E}}_{\mathrm{ext}}\left({\mathbf{r}}_0\right)$. Let us consider a homogeneous medium characterized by ${ϵ}_1$ and ${\mu }_1$, in which a nanoparticle with dielectric function ${ϵ}_{\mathrm{np}}\left(\omega \right)$ lives. Then, the electric field obeys Equation (1) with the free current due to the nanoparticle polarization (excluding the current due to the polarization density of the homogeneous medium) being written as

$${\mathbf{j}}_f(\mathbf{r},\omega )=-i\omega \left[{\mathbf{P}}_{\mathrm{np}}(\mathbf{r},\omega )-{\mathbf{P}}_1\left(\omega \right)\right]=-i\omega {ϵ}_0\left({ϵ}_{\mathrm{np}}\left(\omega \right)-{ϵ}_1\right)\mathbf{E}(\mathbf{r},\omega ),$$

(38)

where we have used the usual linear constitutive relation ${\mathbf{P}}_n(\omega ,\mathbf{r})={ϵ}_0\left({ϵ}_n-1\right)\mathbf{E}(\mathbf{r},\omega )$. ${\mathbf{P}}_{\mathrm{np}}(\mathbf{r},\omega )$ is the polarization density due to the nanoparticle, ${\mathbf{P}}_1\left(\omega \right)$ is the polarization density due to the homogeneous medium, and $\mathbf{E}\left(\mathbf{r}\right)$ is the total electric field in the nanoparticle. Therefore, from Equation (1), the electric field obeys a Lippmann–Schwinger equation [44]

$$\mathbf{E}(\mathbf{r},\omega )={\mathbf{E}}_{\mathrm{ext}}(\mathbf{r},\omega )+{\omega }^2{\mu }_1{\mu }_0{\int }_{V}d{\mathbf{r}}^{\prime }\left({ϵ}_{\mathrm{np}}\left(\omega \right)-{ϵ}_1\right)\overleftrightarrow{G_0}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )·\mathbf{E}(\mathbf{r},\omega )$$

(39)

where ${\mathbf{E}}_{\mathrm{ext}}(\mathbf{r},\omega )$ is a solution of the wave equation in the homogeneous medium, and V is the volume of the nanoparticle. We want to solve for the electric field inside the nanoparticle. We will follow the approximate approach of [22]. We shall assume a spherical nanoparticle, with radius R, and assume that $k_{n}R\ll 1$. This allows us to approximate $\mathbf{E}(\mathbf{r},\omega )$ as being constant inside the nanoparticle and to take the limit $\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|\to 0$ for $\overleftrightarrow{G_0}(\mathbf{r}-{\mathbf{r}}^{\prime },\omega )$. Taking into account Equations (12)–(14), in the limit $\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|\to 0$ we can write the regular part (excluding the Dirac $\delta$-function term) of the free-space dyadic Green’s function as

$${\overleftrightarrow{G_0}}^{\mathrm{reg}}\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)\simeq \frac{1}{6\pi }\frac{1}{\left|\mathbf{R}\right|}\overleftrightarrow{I}+i\frac{k_1}{6\pi }\overleftrightarrow{I},$$

(40)

where $\mathbf{R}=\mathbf{r}-{\mathbf{r}}^{\prime }$. We neglect the real part of ${\overleftrightarrow{G_0}}^{\mathrm{reg}}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )$ when compared to ${\overleftrightarrow{G_0}}^{\mathrm{SF}}(\mathbf{r}-{\mathbf{r}}^{\prime },\omega )$, thereby approximating

$$\overleftrightarrow{G_0}\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)\simeq -\frac{1}{3k_1^2}\overleftrightarrow{I}\delta (\mathbf{r}-{\mathbf{r}}^{\prime })+i\frac{k_1}{6\pi }\overleftrightarrow{I}.$$

(41)

Using the above approximation in Equation (39) and assuming that the electric field within the nanoparticle varies slowly, that is, $\mathbf{E}(\mathbf{r},\omega )=\mathbf{E}({\mathbf{r}}_0,\omega )$ throughout V, we can solve the obtained algebraic equation for the total field $\mathbf{E}({\mathbf{r}}_0,\omega )$ as a function of the external field ${\mathbf{E}}_{\mathrm{ext}}({\mathbf{r}}_0,\omega )$, obtaining

$$\mathbf{E}({\mathbf{r}}_0,\omega )=\frac{1}{1-\frac{1}{{ϵ}_1}\left({ϵ}_{\mathrm{np}}\left(\omega \right)-{ϵ}_1\right)\left(-\frac{1}{3}+i\frac{k_1^3}{6\pi }V\right)}{\mathbf{E}}_{\mathrm{ext}}({\mathbf{r}}_0,\omega ).$$

(42)

Therefore, the electric dipole moment follows from

$$\begin{array}{cc}\hfill \mathbf{d}& ={ϵ}_0\left({ϵ}_{\mathrm{np}}\left(\omega \right)-{ϵ}_1\right){\int }_V{d}^3\mathbf{r}\mathbf{E}(\mathbf{r},\omega )\hfill \\ & \simeq \frac{{\alpha }_{\mathrm{CM}}}{1-i\frac{k_1^3}{6\pi {ϵ}_0{ϵ}_1}{\alpha }_{\mathrm{CM}}}{\mathbf{E}}_{\mathrm{ext}}({\mathbf{r}}_0,\omega ),\hfill \end{array}$$

(43)

where we have introduced the Clausius–Mossotti polarizability [45]

$${\alpha }_{\mathrm{CM}}=4\pi {ϵ}_1{ϵ}_0{R}^3\frac{{ϵ}_{\mathrm{np}}\left(\omega \right)-{ϵ}_1}{{ϵ}_{\mathrm{np}}\left(\omega \right)+2{ϵ}_1}.$$

(44)

The polarizability of a nanoparticle embedded in a homogeneous medium with relative permittivity ${ϵ}_1$ can be read from Equation (43)

$$\overleftrightarrow{{\alpha }_0}=\frac{{\alpha }_{\mathrm{CM}}}{1-i\frac{k_1^3}{6\pi {ϵ}_0{ϵ}_1}{\alpha }_{CM}}\overleftrightarrow{I}.$$

(45)

It must be pointed out that the above equation is only approximate. As a matter of fact it is easy to see that if we had kept the term $\propto 1/\left|\mathbf{R}\right|$ in ${\overleftrightarrow{G_0}}^{\mathrm{reg}}(\mathbf{r}-{\mathbf{r}}^{\prime },\omega )$ we would have generated a real term, $\propto k_1^2$, in the denominator of Equation (45), which can be interpreted as a dynamic depolarization effect [46]. The obtained term would still be incorrect, as additional terms of the same order in $k_1$ would appear from taking into account that the electric field inside the nanoparticle is not constant. An exact treatment using Mie’s scattering theory for a spherical particle would lead to [45,47]

$${\overleftrightarrow{\alpha }}_{\mathrm{Mie}}=4\pi {ϵ}_0{ϵ}_1{R}^3{\left[\frac{{ϵ}_{\mathrm{np}}\left(\omega \right)+2{ϵ}_1}{{ϵ}_{\mathrm{np}}\left(\omega \right)-{ϵ}_1}-\frac{3}{5}\frac{{ϵ}_{\mathrm{np}}\left(\omega \right)-2{ϵ}_1}{{ϵ}_{\mathrm{np}}\left(\omega \right)-{ϵ}_1}{R}^2{k}_1^2-i\frac{2}{3}{R}^3{k}_1^3\right]}^{-1}\overleftrightarrow{I}.$$

(46)

There is indeed a term of order $k_1^2$, but the term of order $k_1^3$ is unchanged. The imaginary term of order $k_1^3$ is usually denoted by radiation damping correction [1,48] and is essentially to enforce the optical theorem for electromagnetic scattering in the lowest order [1,49,50,51]. We note in passing that the radiation damping correction is also responsible for the decay rate of the dipole [1]. In the next sections, we will ignore the term of order $k_1^2$ as it plays no significant role. However, it will become clear that it is essential to keep the radiation damping correction.

3.2. Polarizability of a Quantum Emitter in Proximity to a Planar Interface

If the nanoparticle is situated in the vicinity of an interface, it is also possible to write an equation of the Lippmann–Schwinger type for the electric field similar to Equation (39). The only difference is that in order to take into account the interface, the free-space dyadic Green’s function must be replaced by the one which incorporates the reflection from the substrate, $\overleftrightarrow{G_0}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )\to \overleftrightarrow{G}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )=\overleftrightarrow{G_0}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )+\overleftrightarrow{G_r}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )$. Likewise, the external field ${\mathbf{E}}_0(\mathbf{r},\omega )$ must be replaced by a solution of the electric field wave equation in the presence of the substrate, ${\mathbf{E}}_0(\mathbf{r},\omega )\to {\mathbf{E}}_{\mathrm{ext}}(\mathbf{r},\omega )={\mathbf{E}}_0(\mathbf{r},\omega )+{\mathbf{E}}_r(\mathbf{r},\omega )$, where ${\mathbf{E}}_0(\mathbf{r},\omega )$ is the incident field and ${\mathbf{E}}_r(\mathbf{r},\omega )$ is the reflected field. Therefore, the Lippmann–Schwinger equation for the electric field taking into account the substrate is given by

$$\mathbf{E}(\mathbf{r},\omega )={\mathbf{E}}_{\mathrm{ext}}(\mathbf{r},\omega )+{\omega }^2{\mu }_1{\mu }_0{\int }_V{d}^3{\mathbf{r}}^{\prime }\left({ϵ}_{\mathrm{np}}\left(\omega \right)-{ϵ}_1\right)\overleftrightarrow{G}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )·\mathbf{E}(\mathbf{r},\omega ).$$

(47)

We now proceed in the same fashion as before, assuming $k_{1}R\ll 1$, approximating $\mathbf{E}(\mathbf{r},\omega )=\mathbf{E}({\mathbf{r}}_0,\omega )$ as constant inside the nanoparticle, and keeping only the dominant contributions from $\overleftrightarrow{G}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )$ in the limit $|\mathbf{r}-{\mathbf{r}}^{\prime }|\to 0$. Therefore, we write [22]

$$\overleftrightarrow{G}(\mathbf{r},{\mathbf{r}}^{\prime },\omega )\simeq -\frac{1}{3k_1^2}\overleftrightarrow{I}\delta (\mathbf{r}-{\mathbf{r}}^{\prime })+i\frac{k_1}{6\pi }\overleftrightarrow{I}+{\overleftrightarrow{G}}_r({\mathbf{r}}_0,{\mathbf{r}}_0,\omega ),$$

(48)

where we have used the fact that $\overleftrightarrow{G_r}({\mathbf{r}}_0,{\mathbf{r}}_0,\omega )$ is regular. Introducing Equation (48) into Equation (47), and approximating $\mathbf{E}(\mathbf{r},\omega )\simeq \mathbf{E}({\mathbf{r}}_0,\omega )$, we can solve for $\mathbf{E}({\mathbf{r}}_0,\omega )$ as a function of ${\mathbf{E}}_{\mathrm{ext}}({\mathbf{r}}_0,\omega )$, obtaining

$$\begin{array}{cc}\hfill \mathbf{E}({\mathbf{r}}_0,\omega )& =\frac{3{ϵ}_1}{{ϵ}_{\mathrm{np}}\left(\omega \right)+2{ϵ}_1}{\left[\overleftrightarrow{I}-{\alpha }_{\mathrm{CM}}{\omega }^2{\mu }_1{\mu }_0\left(\frac{i{k}_1}{6\pi }\overleftrightarrow{I}+{\overleftrightarrow{G}}_r({\mathbf{r}}_0,{\mathbf{r}}_0,\omega )\right)\right]}^{-1}·{\mathbf{E}}_{\mathrm{ext}}({\mathbf{r}}_0,\omega ).\hfill \end{array}$$

(49)

The electric dipole moment is thus given by

$$\mathbf{d}=V{ϵ}_0\left({ϵ}_{\mathrm{np}}\left(\omega \right)-{ϵ}_1\right)\mathbf{E}({\mathbf{r}}_0,\omega )={\overleftrightarrow{\alpha }}_{\mathrm{eff}}·{\mathbf{E}}_{\mathrm{ext}}({\mathbf{r}}_0,\omega ),$$

(50)

with the effective polarizability defined as

$$\begin{array}{cc}\hfill {\overleftrightarrow{\alpha }}_{\mathrm{eff}}& ={\alpha }_{\mathrm{CM}}{\left[\overleftrightarrow{I}-{\mu }_1{\mu }_0{\omega }^2{\alpha }_{\mathrm{CM}}{\left(i\Im \overleftrightarrow{G_0}({\mathbf{r}}_0,{\mathbf{r}}_0,\omega )+{\overleftrightarrow{G}}_r({\mathbf{r}}_0,{\mathbf{r}}_0,\omega )\right)}^{-1}\right]}^{-1}.\hfill \end{array}$$

(51)

This equation can be expressed in terms of the free-space polarizability Equation (45) as

$${\overleftrightarrow{\alpha }}_{\mathrm{eff}}={\overleftrightarrow{\alpha }}_0{\left[\overleftrightarrow{I}-{\mu }_1{\mu }_0{\omega }^2{\overleftrightarrow{G}}_r({\mathbf{r}}_0,{\mathbf{r}}_0,\omega )·{\overleftrightarrow{\alpha }}_0\right]}^{-1}.$$

(52)

Equation (51) has been derived in the literature before, following a similar argumentation [6,21,22].

The importance of keeping the free-space radiation damping correction, $i\Im \overleftrightarrow{G_0}({\mathbf{r}}_0,{\mathbf{r}}_0,\omega )$, will now become clear. According to Poynting’s theorem, the power dissipated by the nanoparticle is given by

$$P_{\mathrm{dis}}=\frac{\omega }{2}\Im \left[{\mathbf{E}}_{\mathrm{ext}}^{†}({\mathbf{r}}_0,\omega )·{\overleftrightarrow{\alpha }}_{\mathrm{eff}}·{\mathbf{E}}_{\mathrm{ext}}({\mathbf{r}}_0,\omega )\right].$$

(53)

This implies that the imaginary part of the diagonal components of ${\overleftrightarrow{\alpha }}_{\mathrm{eff}}$ must be positive, since the dissipated power must be positive. It is easily checked that if a and g are complex quantities then $\Im \left[a/(1-ag)\right]=(\Im a+{\left|a\right|}^2\Im g)/{\left|1-ag\right|}^2.$ If $\Im a>0$, but a is otherwise arbitrary, the requirement that $\Im \left[a/(1-ag)\right]>0$ implies that $\Im g>0$. Translating this into the problem of the polarizability, since we have that in general $\Im {\alpha }_{\mathrm{CM}}\ge 0$, the requirement that $\Im {\alpha }_{\mathrm{eff}}>0$ demands that $\Im \left[\overleftrightarrow{G}({\mathbf{r}}_0,{\mathbf{r}}_0,\omega )\right]=\Im \left[i\Im \overleftrightarrow{G_0}({\mathbf{r}}_0,{\mathbf{r}}_0,\omega )+{\overleftrightarrow{G}}_r({\mathbf{r}}_0,{\mathbf{r}}_0,\omega )\right]>0$. This is true only for the complete Green’s function, following from the general properties of Green’s functions, and can be understood either classically as the fact that $\Im \left[\overleftrightarrow{G}({\mathbf{r}}_0,{\mathbf{r}}_0,\omega )\right]$ gives the total power emitted by a point dipole, or quantum mechanically, since the diagonal elements or $\Im \left[\overleftrightarrow{G}({\mathbf{r}}_0,{\mathbf{r}}_0,\omega )\right]$ correspond to a spectral function (a density of electromagnetic states), that is always positive. However, it will not be true in general that $\Im \left[\overleftrightarrow{G_r}({\mathbf{r}}_0,{\mathbf{r}}_0,\omega )\right]$ is positive. As a matter of fact, it becomes negative in situations where subradiance of a quantum emitter occurs. Therefore, the requirement that $\Im {\alpha }_{\mathrm{eff}}>0$, forces us to keep the free-space radiation damping correction.

3.3. Renormalized Polarizability of an Isotropic Quantum Emitter Near a Continuous Graphene Sheet

In what follows we shall consider the case of an isotropic quantum emitter in close proximity to a graphene sheet. In the previous sections, we have seen how the effective polarizability of a nanoparticle depends on the reflected Green’s function, ${\overleftrightarrow{G}}_r({\mathbf{r}}_0,{\mathbf{r}}_0,\omega )$, which can be reconstructed from its angular spectrum representation as

$${\overleftrightarrow{G}}_r({\mathbf{r}}_0,{\mathbf{r}}_0,\omega )=\int \frac{d^2{\mathbf{p}}_{\parallel }}{{\left(2\pi \right)}^2}{\overleftrightarrow{G}}_r\left({\mathbf{p}}_{\parallel },z_0,z_0,\omega \right).$$

(54)

As shown in Section 2.3, the reflected Green’s function in the angular spectrum representation can be written in terms of the Fresnel reflection coefficients. For a planar interface covered by graphene, the reflection coefficients are given by [18,52]

$$\begin{array}{cc}\hfill r_s& =\frac{{\beta }_1-{\beta }_2-{\mu }_0\omega {\sigma }_T\left(\omega \right)}{{\beta }_1+{\beta }_2+{\mu }_0\omega {\sigma }_T\left(\omega \right)},\hfill \end{array}$$

(55)

$$\begin{array}{cc}\hfill r_p& =\frac{{\beta }_1{ϵ}_2-{\beta }_2{ϵ}_1+{\beta }_1{\beta }_2{\sigma }_L/\left({ϵ}_0\omega \right)}{{\beta }_1{ϵ}_2+{\beta }_2{ϵ}_1+{\beta }_1{\beta }_2{\sigma }_L/\left({ϵ}_0\omega \right)},\hfill \end{array}$$

(56)

where ${\sigma }_T\left(\omega \right)$ and ${\sigma }_L\left(\omega \right)$ are the transverse and longitudinal optical conductivities of graphene and ${\beta }_{1/2}$ given by Equation (21). Neglecting nonlocal effects in the conductivities we have ${\sigma }_T\left(\omega \right)={\sigma }_L\left(\omega \right)=\sigma \left(\omega \right)$, which we will model with a Drude-like term [18,53,54]

$$\sigma \left(\omega \right)=\frac{e^2}{4\hslash }\frac{4}{\pi }\frac{{ϵ}_F}{\hslash \gamma -i\hslash \omega },$$

(57)

where ${ϵ}_F$ is graphene’s Fermi energy and $\gamma$ is the broadening factor. The transmission coefficients $t_s$ and $t_p$ are related to the reflection coefficients via [18]

$$\begin{array}{cc}\hfill t_s& =1+r_s,\hfill \end{array}$$

(58)

$$\begin{array}{cc}\hfill t_p& =\frac{{\beta }_1}{{\beta }_2}\sqrt{\frac{{ϵ}_2}{{ϵ}_1}}(1-r_p).\hfill \end{array}$$

(59)

After performing the integration over the angular variable in Equation (54), we find that $\overleftrightarrow{G_r}({\mathbf{r}}_0,{\mathbf{r}}_0,\omega )$ is diagonal. Rotational invariance along the z direction imposes that $G_r^{xx}({\mathbf{r}}_0,{\mathbf{r}}_0,\omega )=G_r^{yy}({\mathbf{r}}_0,{\mathbf{r}}_0,\omega )$, which will differ from $G_r^{zz}({\mathbf{r}}_0,{\mathbf{r}}_0,\omega )$. The same will be true for the polarizability of the nanoparticle, which, using Equation (51), we can write as

$$\begin{array}{cc}\hfill {\alpha }_{\mathrm{eff}}^{xx}=& {\alpha }_{\mathrm{eff}}^{yy}=4\pi {ϵ}_0{ϵ}_1{R}^3\frac{{\tilde{\alpha }}_0}{1-{\left(k_{1}R\right)}^3{\mathcal{G}}_r^{\parallel }\left({\mathbf{r}}_0,{\mathbf{r}}_0,\omega \right){\tilde{\alpha }}_0},\hfill \end{array}$$

(60)

$$\begin{array}{c}{\alpha }_{\mathrm{eff}}^{zz}=4\pi {ϵ}_0{ϵ}_1{R}^3\frac{{\tilde{\alpha }}_0}{1-{\left(k_{1}R\right)}^3{\mathcal{G}}_r^{zz}\left({\mathbf{r}}_0,{\mathbf{r}}_0,\omega \right){\tilde{\alpha }}_0},\hfill \end{array}$$

(61)

where we have defined the dimensionless quantities ${\tilde{\alpha }}_0={\alpha }_0/\left(4\pi {ϵ}_0{ϵ}_1{R}^3\right)$ with ${\alpha }_0$ the diagonal element of the nanoparticle polarizability from Equation (45). ${\mathcal{G}}_r^{\parallel }\left({\mathbf{r}}_0,{\mathbf{r}}_0,\omega \right)=\left(4\pi /k_1\right)G_r^{xx}\left({\mathbf{r}}_0,{\mathbf{r}}_0,\omega \right)=\left(4\pi /k_1\right)G_r^{yy}\left({\mathbf{r}}_0,{\mathbf{r}}_0,\omega \right)$ and ${\mathcal{G}}_r^{zz}\left({\mathbf{r}}_0,{\mathbf{r}}_0,\omega \right)=\left(4\pi /k_1\right)G_r^{zz}\left({\mathbf{r}}_0,{\mathbf{r}}_0,\omega \right)$.

More explicitly, these quantities can be evaluated from

$$\begin{array}{cc}\hfill {\mathcal{G}}_r^{\parallel }\left({\mathbf{r}}_0,{\mathbf{r}}_0,\omega \right)& =\frac{i}{2}{\int }_0^{\infty }ds{e}^{i2k_1{z}_0\sqrt{1-s^2}}s\left(\frac{1}{\sqrt{1-s^2}}{r}_s-\sqrt{1-s^2}{r}_p\right),\hfill \end{array}$$

(62)

$$\begin{array}{cc}\hfill {\mathcal{G}}_r^{zz}\left({\mathbf{r}}_0,{\mathbf{r}}_0,\omega \right)& =i{\int }_0^{\infty }ds{e}^{i2k_{1}d\sqrt{1-s^2}}\frac{s^3}{\sqrt{1-s^2}}{r}_p,\hfill \end{array}$$

(63)

where $s=p_{\parallel }/k_1$.

Some insight on the previous expressions can be obtained by estimating them in the electrostatic limit, valid for $k_1{z}_0\ll 1$. In this limit, the main contribution is due to the $r_p$ reflection coefficient. Approximating $\sqrt{1-s^2}\simeq \sqrt{{ϵ}_2/{ϵ}_1-s^2}\simeq is$ we obtain

$${\mathcal{G}}_r^{zz}\left({\mathbf{r}}_0,{\mathbf{r}}_0,\omega \right)\simeq 2{\mathcal{G}}_r^{\parallel }\left({\mathbf{r}}_0,{\mathbf{r}}_0,\omega \right)\simeq {\int }_0^{\infty }ds{e}^{-2k_1{z}_{0}s}{s}^2{r}_p,$$

(64)

with the reflection coefficient being approximated by

$$\begin{array}{cc}\hfill r_p& =1-\frac{2{\beta }_2{ϵ}_1}{{\beta }_1{ϵ}_2+{\beta }_2{ϵ}_1+{\beta }_1{\beta }_2{\sigma }_L/\left({ϵ}_0\omega \right)}\hfill \\ & \simeq 1-\frac{2{ϵ}_1}{{ϵ}_2+{ϵ}_1}\frac{k_{\mathrm{spp}}\left(\omega \right)}{k_{\mathrm{spp}}\left(\omega \right)-k_{1}s},\hfill \end{array}$$

(65)

where we have approximated ${\beta }_1\simeq {\beta }_2\simeq i{k}_{1}s$ and have introduced the graphene’s surface plasmon polariton complex wavenumber

$$k_{\mathrm{spp}}\left(\omega \right)=\frac{\omega }{c}\frac{{ϵ}_1+{ϵ}_2}{4{\alpha }_f}\frac{\hslash \omega +i\hslash \gamma }{{ϵ}_F},$$

(66)

and ${\alpha }_f\simeq 1/137$ is the fine structure constant. From these results we can already estimate when the effect of the graphene substrate on the NP polarizability will be more significant. From Equation (65), $r_p$ is peaked at $s=\Re k_{\mathrm{spp}}\left(\omega \right)/k_1$, while the term $e^{-2k_1{z}_{0}s}{s}^2$ in the integrand of Equation (64) has a maximum at $s={\left(k_1{z}_0\right)}^{-1}$. Therefore, ${\mathcal{G}}_r^{zz}\left({\mathbf{r}}_0,{\mathbf{r}}_0,\omega \right)$ will have a maximum when these two peaks coincide, [43] which occurs for $\Re k_{\mathrm{spp}}\left(\omega \right)z_0\simeq 1$. In the electrostatic limit, Equation (64) can be written in terms of known functions as

$${\mathcal{G}}_r^{zz}\left({\mathbf{r}}_0,{\mathbf{r}}_0,\omega \right)\simeq 2{\mathcal{G}}_r^{\parallel }\left({\mathbf{r}}_0,{\mathbf{r}}_0,\omega \right)\simeq {\left(\frac{k_{\mathrm{spp}}\left(\omega \right)}{k_1}\right)}^{3}f\left(2k_{\mathrm{spp}}\left(\omega \right)z_0\right),$$

(67)

where the function $f\left(z\right)$ is given by

$$f\left(z\right)=\frac{2}{z^3}+\frac{2{ϵ}_1}{{ϵ}_1+{ϵ}_2}\left(\frac{1}{z^2}+\frac{1}{z}+e^{-z}\left[i\pi -\mathrm{Ei}\left(z\right)\right]\right),$$

(68)

with $\mathrm{Ei}\left(z\right)$ the exponential integral function, which for a real positive argument is written as $\mathrm{Ei}\left(x\right)=-{⨏}_{-x}^{\infty }dt{e}^{-t}/t.$ However, we point out that Equation (67) is valid even in the presence of finite broadening $\gamma$ in graphene.

We shall consider both metallic and polar semiconductor nanoparticles, with the relative dielectric function described, respectively, by the Drude and Lorentz models. The Drude model for the dielectric function reads

$${ϵ}_{\mathrm{Drude}}\left(\omega \right)=1-\frac{{\omega }_p^2}{\omega (\omega +i{\mathsf{\Gamma }}_0/\hslash )}$$

(69)

where ${\omega }_p$ is the metal’s plasma frequency and ${\mathsf{\Gamma }}_0$ is the relaxation rate, while the Lorentz model for the dielectric function of a polar material is given by

$${ϵ}_{\mathrm{Lorentz}}\left(\omega \right)={ϵ}_{\infty }\left(1+\frac{{\omega }_{\mathrm{L}0}^2-{\omega }_{\mathrm{T}0}^2}{{\omega }_{\mathrm{T}0}^2-{\omega }^2-i\omega {\mathsf{\Gamma }}_{\mathrm{TO}}}\right),$$

(70)

where ${\omega }_{\mathrm{TO}}$ and ${\omega }_{\mathrm{LO}}$ are the frequencies of the transverse and longitudinal optical phonons, ${\mathsf{\Gamma }}_{\mathrm{TO}}$ is a phonon decay rate, and ${ϵ}_{\infty }$ is the high frequency limit of the dielectric function. As examples of commonly used materials for the production of nanoparticles, we consider gold (metallic) and CdSe (polar semiconducing) nanoparticles. Typical values of the polarizability for different substances are given in [55]. The used values for the Lorentz model of CdSe are taken from [56].

In Figure 3 we depict the real and imaginary parts of the polarizability of a gold nanoparticle near a doped graphene sheet on a substrate with ${ϵ}_2=2$. The figure clearly shows the strong renormalization of the polarizability of the nanoparticle relative to its value in the presence of the interface without graphene (blue dashed line). This is due to the close proximity of the nanoparticle to the graphene sheet, $z_0=151$nm (compared to the typical free-space wavelengths in the THz range). Nowadays, with the ubiquitous use of hexagonal boron nitride (h-BN) for encapsulating graphene, together with the possibility of controlling the number of layers of h-BN, it poses no difficulty to routinely produce structures where nanoparticles are positioned very close to the graphene sheet, that is, at distances much smaller than their radius. Also, the $zz$-component of the polarizability tensor (black dotted line) is renormalized differently from the $xx$-component (red solid line). This is a consequence of breaking the translation symmetry along the z–direction introduced by the graphene sheet and the dielectric change as we cross the $z=0$ plane. We have verified that the broadband resonance seen in the imaginary part of the polarizability tensor is due to the excitation of surface plasmon polaritons in graphene. This was assessed studying the dispersion of the resonance as a function of the Fermi energy (more on this below).

Given the close proximity of the nanoparticle to the graphene sheet, the question of the necessity of a nonlocal description of the graphene conductivity arises. In order to check the correctness of our local description, we have performed simulations (results not shown) using the nonlocal Drude-like conductivity [18] of graphene. We have found that nonlocal effects in the graphene conductivity (its dependence on the wavevector) play no visible role in the position or the intensity of the resonance in the effective polarizability of the nanoparticle (when $z_0=151$ nm). We also expect nonlocal effects akin to the nanoparticle to be negligible provided that $z_0\gtrsim 10$nm, below which nonlocal contributions in metallic nanoparticles usually arise [57,58] (the situation is different for semiconductor nanoparticles [58]).

In Figure 4 we present the spectrum of the polarizability of a CdSe nanoparticle in the presence of graphene on a substrate. As in Figure 3, the observed broad band resonance in the imaginary part of the polarizability tensor components is due to the excitation of surface plasmons in graphene. As discussed previously, the order of magnitude of the plasmonic resonance frequency can be estimated from the relation $k_{\mathrm{spp}}\left(\omega \right)z_0\simeq 1$, which is independent of the nanoparticle’s material. Indeed, we observe the resonance here approximately at the same frequency as for the metallic nanoparticle (see Figure 3). In order to further access the plasmonic nature of the broad band resonance, we have studied its position as a function of the Fermi energy and found a complete agreement with the above relation, that is, in the peak of the resonance scales with the Fermi energy $\sqrt{E_F}$. Interestingly, the intensity of the resonance in Figure 4 is smaller by a factor of 2.5 when compared to the case of the metallic nanoparticle. This can be understood from the following simple consideration. We notice that in the imaginary part of the bare polarizability, $\Im {\tilde{\alpha }}_0$ is negligible for both metallic and CdSe nanoparticles in the considered frequency range. Therefore, we can write for the imaginary part of the renormalized polarizability (e.g., the $xx$ component):

$$\begin{array}{cc}\hfill {\left[\Im \frac{{\alpha }_{\mathrm{eff}}^{xx}}{4\pi {ϵ}_0{ϵ}_1{R}^3}\right]}_{max}=& {\tilde{\alpha }}_0{\left[\Im \frac{1}{1-{\left(k_{1}R\right)}^3{\tilde{\alpha }}_0{\mathcal{G}}_r^{\parallel }\left({\mathbf{r}}_0,{\mathbf{r}}_0,\omega \right)}\right]}_{max}\hfill \\ & \approx {\tilde{\alpha }}_0^2\frac{{\left[{\left(k_{1}R\right)}^3\Im {\mathcal{G}}_r^{\parallel }\left({\mathbf{r}}_0,{\mathbf{r}}_0,\omega \right)\right]}_{max}}{{\left|1-{\left(k_{1}R\right)}^3{\tilde{\alpha }}_0\Re {\mathcal{G}}_r^{\parallel }\left({\mathbf{r}}_0,{\mathbf{r}}_0,\omega \right)\right|}^2}.\hfill \end{array}$$

(71)

Since ${ϵ}_{CdSe}\approx {ϵ}_{\infty }$ for $\omega >$ 15 THz, we have ${\tilde{\alpha }}_0\approx 0.684$ for CdSe, while for the metal particles ${\tilde{\alpha }}_0\approx 1$, because the metal permittivity is negative and large in modulus in the far-infrared region. The denominator is close to unity in both cases, so in the case of a metal nanoparticle, we have a stronger renormalization of the polarizability (by a factor of ≈2.5) in the presence of a graphene sheet. This is the maximum renormalization that one can obtain if the bare polarizability is dispersionless. In contrast to metals, a nanoparticle made of CdSe has a dipolar mode due to optical phonons, which occurs at the so-called Fröhlich frequency (unless the particle’s size is so small that it gives rise to quantum confinement effects, of the order of few nanometers) [59], ${\omega }_F={\left[({ϵ}_{\infty }{\omega }_{\mathrm{LO}}^2+2{\omega }_{\mathrm{TO}}^2)({ϵ}_{\infty }+2)\right]}^{1/2}$ , which is equal to ≈6 THz in this case. When the plasmonic resonance overlaps with ${\omega }_F$, the phonon resonance in the nanoparticle is greatly enhanced because $\Im {\tilde{\alpha }}_0\gg 1$ at $\omega \approx {\omega }_F$. In this case, a non-trivial dependence on the Fermi energy takes place [19]. Note that none of these effects would take place in the presence of a metallic substrate for the same studied spectral range, as plasmons in metals at these frequencies are essentially free radiation.

3.4. Renormalized Polarizability of an Isotropic Quantum Emitter Near a Plasmonic Graphene Grating

In this section we revisit the problem of the renormalization of the polarizability of a quantum emitter now considering it near a plasmonic graphene grating. The used procedure is only approximate, relying on a semi-analytic approach. However, the analysis performed is sufficient to capture the effect of plasmonic resonances of the graphene grating in the nanoparticle polarizability.

3.4.1. Optical Properties of a Plasmonic Graphene Grating

For a metamaterial such as the graphene-based grating depicted in Figure 1b, the description of the interaction of the material with a quantum emitter can be quite complex. One possible method for overcoming such a difficulty is by computing the effective conductivity of the metamaterial, in this case the plasmonic graphene grating. The general method for accomplishing this was given in [60] and was later applied to the problem of tuning total absorption in graphene [61], but no details of its calculation were given. Instrumental to the calculation of the effective conductivity is the knowledge of the reflection and transmission Fresnel coefficients. These coefficients were computed in approximated analytical form in [62] and we give here only the final results:

$$\begin{array}{cc}\hfill r_{p,m}& =-{\delta }_{m,0}+t_{p,m}+{\mu }_0\chi \left(\omega \right)\frac{w}{4}{J}_1(m\pi w/L)\hfill \end{array}$$

(72)

$$\begin{array}{cc}\hfill t_{p,m}& =\frac{{ϵ}_2{\beta }_{1,m}}{{ϵ}_1{\beta }_{2,m}+{ϵ}_2{\beta }_{1,m}}\left(2{\delta }_{m,0}-{\mu }_0\chi \left(\omega \right)\frac{w}{4}{J}_1(m\pi w/L)\right)\hfill \end{array}$$

(73)

where $r_{p,0}$ and $t_{p,0}$ are the reflection and transmission coefficients, respectively, of the zero diffraction-order of the grating (the only propagating order for a sub-wavelength grating), w is the width of the graphene ribbons in the grating, L is the period of the grating, and the function $\chi \left(\omega \right)$ reads

$$\chi \left(\omega \right)=\frac{2{\beta }_{2,0}{\beta }_{1,0}}{{ϵ}_1{\beta }_{2,0}+{ϵ}_2{\beta }_{1,0}}\frac{{\sigma }_L\left(\omega \right)c^2}{\omega \mathsf{\Lambda }\left(\omega \right)}$$

(74)

which encodes the information about the plasmonic resonance in the grating. With $\mathsf{\Lambda }\left(\omega \right)$ given by

$$\mathsf{\Lambda }\left(\omega \right)=\frac{w}{4}\sum _{n=-\infty }^{\infty }\frac{1}{n}{J}_1(n\pi w/L)\left[1+\frac{{\sigma }_L\left(\omega \right)}{\omega {ϵ}_0}\frac{{\beta }_{2,n}{\beta }_{1,n}}{{ϵ}_1{\beta }_{2,n}+{ϵ}_2{\beta }_{1,n}}\right]$$

(75)

where ${\beta }_{1,n}=\sqrt{k_1^2-k_x^2-q_n^2}$ and ${\beta }_{2,n}=\sqrt{k_2^2-k_x^2-q_n^2}$, with $q_n=k_y+n2\pi /L$, $J_1\left(x\right)$ is the Bessel function of order 1. Here, the summation in $\mathsf{\Lambda }\left(\omega \right)$ is delicate due to the oscillatory nature of the Bessel function (see [62]). For simplicity of the calculation, we approximate ${\beta }_{j,n\ne 0}$ by ${\beta }_{j,n\ne 0}\approx \sqrt{k_j^2-p_{\parallel }^2-n^{2}4{\pi }^2/L^2}$. In addition to $r_{p,0}$ and $t_{p,0}$ there is an infinite number of other coefficients associated with higher diffraction order, but they are all evanescent in nature for the parameters chosen in the figures. Therefore, we approximate the optical properties of the grating considering only $r_{p,0}$ and $t_{p,0}$, and $r_{p,1}$ and $t_{p,1}$ (we have checked that introducing more evanescent terms does not change the results). This gives us an analytical description of its optical properties. As noted above, from the knowledge of $r_{p,0}$ and $t_{p,0}$, and $r_{p,1}$ and $t_{p,1}$ we can derive an effective conductivity for the graphene grating along the direction perpendicular to the axis of the graphene ribbon. This effective conductivity shows a maximum in its real part associated with the excitation of surface plasmon polaritons. The same information is encoded in the function $\chi \left(\omega \right)$, as can be seen in Figure 5 and, in fact, for our analysis this latter function is all we need for including plasmonic effects into the calculation.

Notice that the conductivity of the system is no longer isotropic. Therefore, we will introduce this anisotropy in an effective way, choosing different Fermi energies for the $r_s$ and $r_p$ reflection coefficients. Also, while the $r_{p,m}$ coefficients are given by Equation (72), the $r_s$ coefficient is given by Equation (55). This procedure renders our results qualitative and no quantitative agreement is expected with an exact calculation. The exact solution would require us to extend the formalism to the case on a non-isotropic system in the $xy-$plane. Note that this system has broken rotational symmetry around the $z-$axis. Therefore, we expected that the equality ${\alpha }_{xx}={\alpha }_{yy}$ seen in the case of continuous graphene sheet should not hold in the case of grating. Our qualitative results show that this is indeed the case.

3.4.2. Renormalization of the Polarizability of a Quantum Emitter

In this section we study the renormalization of the polarizability of a quantum emitter near a plasmonic graphene-based grating. As explained above, we use the reflection coefficients $r_{p,0}$ and $r_{p,1}$ in the reflected p–Green’s function and an effective Fermi energy, given by $E_F^{\mathrm{eff}}=E_{F}w/L$ in the $r_s$ coefficient in Equation (55), and use this in the reflected s-Green’s function. We consider only the case of a metallic nanoparticle, as the results are qualitatively the same for a semiconductor one. In Figure 6 we depict the real and imaginary parts of the renormalized polarizability of a gold nanoparticle in the proximity of a graphene-based grating. A strong renormalization of the real part of the polarizability can be seen at the same frequency where the grating supports the excitation of surface plasmon polaritons (see Figure 5). The same happens in the imaginary part. However, the relative change of the imaginary part is much larger than for the real part. The results for the imaginary part of the polarizability in the case of grating should be compared to those given in Figure 3 for the same quantity. For the continuous sheet the enhancement of the imaginary part of $\alpha$ is about that found in the present case. This is attributed to the approximate description of the reflection coefficients of the grating. Indeed, we would expect the renormalization to be larger in the case of the grating as the latter supports excitation of plasmons by far-field radiation, whereas in the case of the continuous graphene sheet the excitation of plasmons is due to near-field excitation only. We also note that the resonance peak in the imaginary part of the polarizability is not of broad band when compared to the same quantity in the continuous case.

On other hand, the frequency where the maximum of the resonance is located is larger in the present case. This happens since we can tune the position of the resonance in the grating by varying both the Fermi energy and the geometric parameters of the grating. Therefore, the grating system has a versatility that cannot be found in the continuous sheet case. Indeed, using gratings with smaller periods, the resonance can be tuned across the electromagnetic spectrum, from the THz spectra to the infrared. We also note that the renormalization of ${\alpha }_{zz}$ component (black dotted line) is substantially larger than the ${\alpha }_{xx}$ component (red solid line) and the ${\alpha }_{yy}$ one (brown dashed line). This happens because the $zz-$component of the Green’s function is about twice as large compared to the $xx-$component.

As noted above, tuning the grating parameters allows for an additional degree of freedom to control the impact of graphene on the polarizability of a nanoparticle. The tuning can be twofold: (1) changing the period of the grating and keeping the aspect ratio to 1/2 (constant filling factor); and (2) changing the aspect ratio, keeping the period fixed. In the first case, the procedure allows for large changes in the spectral position of the resonance; this is because the momentum of the SPP is essentially $q=\pi /L$. Since the SPP dispersion is proportional to $\sqrt{q}$, L fixes the position of the resonance. In the second case, a fine tuning of the position of the resonance is achieved. However, the effect of the renormalization of the polarizability is greater for the half-filled case. This is because in this regime the plasmonic resonance has the maximum intensity [16,62]

Finally, we have verified that when $w\to L$, we recover the results of a continuous graphene sheet.

4. Extension of the Formalism When the Quantum Emitter Has Both an Electric and a Magnetic Dipole

A current density ${\mathbf{j}}_f(\mathbf{r},\omega )$ of a particle can be described in terms of its moments in a multipole expansion [63]. A small particle, however, can often be described using only the multipole moments of the lowest orders. In the case of a metallic nanoparticle, its response is dominated by the electric dipole moment. Nevertheless, it is known that in some cases it is necessary to go beyond the electric dipole approximation and consider higher-order moments [6]. In particular, it has been shown that silicon nanoparticles with size between the tens and hundreds of nanometers can have strong responses in the infrared and visible spectra due to higher order moments [4,5,6,23,64], with the magnetic dipole moment contributing the most, even though the particles are not magnetic by themselves. This motivates us to generalize the formalism of the previous sections to the case of a point-like nanoparticle (or quantum emitter) with both electric and magnetic dipole moments. Although the Green’s functions technique has been used before in this problem [6,15,38], some details regarding the behavior of the Green’s functions at coincidence, that is, when ${\mathbf{r}}^{\prime }=\mathbf{r}$, have been overlooked. Therefore, we carefully present the full formalism, that is, accounting for both electric and magnetic dipole contributions, below.

4.1. Free-Space Electric, Magnetic and Mixed Green’s Functions

Our starting points are the inhomogeneous Helmholtz equations for the electric and the magnetic fields (in fact the magnetic field induction $\mathbf{B}(\mathbf{r},\omega )$) in the presence of a source current density (see Appendix A for the derivation):

$$\begin{array}{cc}\hfill -{\nabla }^2\mathbf{E}(\mathbf{r},\omega )-k_n^2\mathbf{E}(\mathbf{r},\omega )& =i\omega {\mu }_n{\mu }_0\left[{\mathbf{j}}_f(\mathbf{r},\omega )+\frac{1}{k_n^2}\nabla \left(\nabla ·{\mathbf{j}}_f(\mathbf{r},\omega )\right)\right]\hfill \end{array}$$

(76)

$$\begin{array}{cc}\hfill -{\nabla }^2\mathbf{B}(\mathbf{r},\omega )-k_n^2\mathbf{B}(\mathbf{r},\omega )& ={\mu }_n{\mu }_0\nabla \times {\mathbf{j}}_f(\mathbf{r},\omega ).\hfill \end{array}$$

(77)

As before, we can write the solution for the inhomogeneous Helmholtz equations as

$$\begin{array}{cc}\hfill \mathbf{E}(\mathbf{r},\omega )& ={\mathbf{E}}_0(\mathbf{r},\omega )+i\omega {\mu }_n{\mu }_0{\int }_{\backslash V_{\delta }\left(\mathbf{r}\right)}{d}^3{\mathbf{r}}^{\prime }{g}_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)\left(\overleftrightarrow{I}+\frac{1}{k_n^2}{\nabla }^{\prime }{\nabla }^{\prime }\right){\mathbf{j}}_f({\mathbf{r}}^{\prime },\omega )\hfill \end{array}$$

(78)

$$\begin{array}{cc}\hfill \mathbf{B}(\mathbf{r},\omega )& ={\mathbf{B}}_0(\mathbf{r},\omega )+{\mu }_n{\mu }_0{\int }_{\backslash V_{\delta }\left(\mathbf{r}\right)}{d}^3{\mathbf{r}}^{\prime }{g}_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right){\nabla }^{\prime }\times {\mathbf{j}}_f({\mathbf{r}}^{\prime },\omega ).\hfill \end{array}$$

(79)

In the same spirit of Equation (38), we write the current in terms of polarization, ${\mathbf{P}}_f$, and magnetization, ${\mathbf{M}}_f$, densities as

$${\mathbf{j}}_t(\mathbf{r},\omega )=-i\omega {\mathbf{P}}_f(\mathbf{r},\omega )+\nabla \times {\mathbf{M}}_f(\mathbf{r},\omega ).$$

(80)

Inserting the latter result into Equations (78) and (79) we obtain

$$\begin{array}{cc}\hfill \mathbf{E}(\mathbf{r},\omega )& ={\mathbf{E}}_0(\mathbf{r},\omega )+{\omega }^2{\mu }_n{\mu }_0{\int }_{\backslash V_{\delta }\left(\mathbf{r}\right)}{d}^3{\mathbf{r}}^{\prime }{g}_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)\left(\overleftrightarrow{I}+\frac{1}{k_n^2}{\nabla }^{\prime }{\nabla }^{\prime }\right){\mathbf{P}}_f({\mathbf{r}}^{\prime },\omega )\hfill \\ & +i\omega {\mu }_n{\mu }_0{\int }_{\backslash V_{\delta }\left(\mathbf{r}\right)}{d}^3{\mathbf{r}}^{\prime }{g}_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right){\nabla }^{\prime }\times {\mathbf{M}}_f({\mathbf{r}}^{\prime },\omega ),\hfill \end{array}$$

(81)

$$\begin{array}{cc}\hfill \mathbf{B}(\mathbf{r},\omega )& ={\mathbf{B}}_0(\mathbf{r},\omega )-i\omega {\mu }_n{\mu }_0{\int }_{\backslash V_{\delta }\left(\mathbf{r}\right)}{d}^3{\mathbf{r}}^{\prime }{g}_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)\nabla \times {\mathbf{P}}_f({\mathbf{r}}^{\prime },\omega )\hfill \\ & +{\mu }_n{\mu }_0{\int }_{\backslash V_{\delta }\left(\mathbf{r}\right)}{d}^3{\mathbf{r}}^{\prime }{g}_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)\left(-{\nabla }^{\prime 2}+{\nabla }^{\prime }{\nabla }^{\prime }\right){\mathbf{M}}_f({\mathbf{r}}^{\prime },\omega ),\hfill \end{array}$$

(82)

where we have used the fact that ${\nabla }^{\prime }·\left({\nabla }^{\prime }\times {\mathbf{M}}_f({\mathbf{r}}^{\prime },\omega )\right)=0$ and ${\nabla }^{\prime }\times {\nabla }^{\prime }\times {\mathbf{M}}_f({\mathbf{r}}^{\prime },\omega )={\nabla }^{\prime }\left({\nabla }^{\prime }·{\mathbf{M}}_f({\mathbf{r}}^{\prime },\omega )\right)-{\nabla }^{\prime 2}{\mathbf{M}}_f({\mathbf{r}}^{\prime },\omega )$. We now proceed as in Section 2.1, using integration by parts, while taking into account the boundary terms due to the excluded volume $V_{\delta }$ enclosing the point ${\mathbf{r}}^{\prime }=\mathbf{r}$, in the same form that we have already dealt with for the electric field Green’s function. The crossed terms relating the magnetization to the electric field and the polarization to the magnetic field only involve one derivative of the Helmholtz Green’s function and therefore the generated boundary term vanishes in the limit of infinitesimal excluded volume. Therefore, we may simply write

$${\int }_{\backslash V_{\delta }\left(\mathbf{r}\right)}{d}^3{\mathbf{r}}^{\prime }{g}_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right){\nabla }^{\prime }\times {\mathbf{M}}_f({\mathbf{r}}^{\prime },\omega )={\int }_{\backslash V_{\delta }\left(\mathbf{r}\right)}{d}^3{\mathbf{r}}^{\prime }\nabla g_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)\times {\mathbf{M}}_f({\mathbf{r}}^{\prime },\omega ),$$

(83)

where we have used the fact that in a translation-invariant system ${\nabla }^{\prime }{g}_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)=-\nabla g_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)$. Finally, the term that relates the magnetization to the magnetic field (magnetic field induction) can be treated in a similar way as the one for the electric field Green’s function. The only difference is that we also have to use integration by parts for the Laplacian term. The steps to treat this term are exactly the same as the ones to treat the ${\nabla }^{\prime }{\nabla }^{\prime }$ term in Section 2.1 and we obtain

$$\begin{array}{c}{\int }_{\backslash V_{\delta }\left(\mathbf{r}\right)}{d}^3{\mathbf{r}}^{\prime }{g}_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)\left(-{\nabla }^{\prime 2}+{\nabla }^{\prime }{\nabla }^{\prime }\right){\mathbf{M}}_f({\mathbf{r}}^{\prime },\omega )=\hfill \\ ={\int }_{\backslash V_{\delta }\left(\mathbf{r}\right)}{d}^3{\mathbf{r}}^{\prime }\left(-{\nabla }^{\prime 2}+{\nabla }^{\prime }{\nabla }^{\prime }\right)g_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right){\mathbf{M}}_f({\mathbf{r}}^{\prime },\omega )+L_{V_{\delta }}{\mathbf{M}}_f(\mathbf{r},\omega )-{\overleftrightarrow{L}}_{V_{\delta }}·{\mathbf{M}}_f(\mathbf{r},\omega ),\hfill \end{array}$$

(84)

where ${\overleftrightarrow{L}}_{V_{\delta }}$ is given by Equation (9) and $L_{V_{\delta }}=\mathrm{Tr}\left({\overleftrightarrow{L}}_{V_{\delta }}\right)$, (see Equation (A28)). This quantity is just the solid angle of excluded volume $V_{\delta }$ centered at ${\mathbf{r}}^{\prime }=\mathbf{r}$ divided by $4\pi$, which equals 1 for any surface (see Appendix B). We also point out that for ${\mathbf{r}}^{\prime }\ne \mathbf{r}$ we have $-{\nabla }^{\prime 2}{g}_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)=k_n^2{g}_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)$. These results allow us to write

$$\begin{array}{cc}\hfill \mathbf{E}(\mathbf{r},\omega )& ={\mathbf{E}}_0(\mathbf{r},\omega )+{\omega }^2{\mu }_n{\mu }_0\int d^3{\mathbf{r}}^{\prime }{\overleftrightarrow{G_0}}^{EE}\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)·{\mathbf{P}}_f({\mathbf{r}}^{\prime },\omega )\hfill \\ & +\omega {\mu }_n{\mu }_0{k}_n\int d^3{\mathbf{r}}^{\prime }{\overleftrightarrow{G_0}}^{EM}\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)·{\mathbf{M}}_f({\mathbf{r}}^{\prime },\omega ),\hfill \end{array}$$

(85)

$$\begin{array}{cc}\hfill \mathbf{B}(\mathbf{r},\omega )& ={\mathbf{B}}_0(\mathbf{r},\omega )-\omega {\mu }_n{\mu }_0{k}_n\int d^3{\mathbf{r}}^{\prime }{\overleftrightarrow{G_0}}^{ME}\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)·{\mathbf{P}}_f({\mathbf{r}}^{\prime },\omega )\hfill \\ & +{\mu }_n{\mu }_0{k}_n^2\int d^3{\mathbf{r}}^{\prime }{\overleftrightarrow{G_0}}^{MM}\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)·{\mathbf{M}}_f({\mathbf{r}}^{\prime },\omega ),\hfill \end{array}$$

(86)

where we have the electric field and magnetic field Green’s functions

$$\begin{array}{cc}\hfill {\overleftrightarrow{G_0}}^{EE}\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)& =\mathrm{P}.\mathrm{V}{.}_{V_{\delta }}\left[\overleftrightarrow{I}+\frac{1}{k_n^2}\nabla \nabla \right]g_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)-\frac{1}{k_n^2}{\overleftrightarrow{L}}_{V_{\delta }}\delta \left(\mathbf{r}-{\mathbf{r}}^{\prime }\right),\hfill \end{array}$$

(87)

$$\begin{array}{cc}\hfill {\overleftrightarrow{G_0}}^{MM}\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)& =\mathrm{P}.\mathrm{V}{.}_{V_{\delta }}\left[\overleftrightarrow{I}+\frac{1}{k_n^2}\nabla \nabla \right]g_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)+\frac{1}{k_n^2}\left(\overleftrightarrow{I}-{\overleftrightarrow{L}}_{V_{\delta }}\right)\delta \left(\mathbf{r}-{\mathbf{r}}^{\prime }\right),\hfill \end{array}$$

(88)

and we have the mixed Green’s functions defined as

$${\overleftrightarrow{G_0}}^{EM}\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)={\overleftrightarrow{G_0}}^{ME}\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)={\mathrm{P}.\mathrm{V}.}_{V_{\delta }}\left[\begin{array}{ccc}0& -{\partial }_z& {\partial }_y\\ {\partial }_z& 0& -{\partial }_x\\ -{\partial }_y& {\partial }_x& 0\end{array}\right]\frac{i}{k_n}{g}_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right).$$

(89)

These describe magnetoelectric effects, which can be important when the nanoparticle sits on a substrate [6]. The dyadic $\left(\overleftrightarrow{I}-{\overleftrightarrow{L}}_{V_{\delta }}\right)\delta \left(\mathbf{r}-{\mathbf{r}}^{\prime }\right)$ in Equation (88) can be interpreted as a demagnetization term. For the case for a spherically-symmetric excluded volume, we have $\left(\overleftrightarrow{I}-{\overleftrightarrow{L}}_{V_{\delta }}\right)\delta \left(\mathbf{r}-{\mathbf{r}}^{\prime }\right)=\overleftrightarrow{I}2/3\delta \left(\mathbf{r}-{\mathbf{r}}^{\prime }\right)$. The factor of $2/3$ is well-known as being the demagnetization factor of a spherical particle [63], however, to the best of our knowledge, this term has not been discussed in the literature before in the context of the application of Green’s functions to electromagnetic problems. Correctly taking this term into account is essentially to describe self-field effects in the magnetization of a particle (analogous to the self-field effects in the polarization in the (electric-only) case considered before).

For the case of nanoparticle characterized by a permittivity ${ϵ}_{\mathrm{np}}$ and permeability ${\mu }_{\mathrm{np}}$, the free polarization and magnetization densities inside the nanoparticle volume read

$$\begin{array}{cc}\hfill {\mathbf{P}}_f(\mathbf{r},\omega )& ={\mathbf{P}}_{\mathrm{np}}(\mathbf{r},\omega )-{\mathbf{P}}_n\left(\omega \right)={ϵ}_0\left({ϵ}_{\mathrm{np}}-{ϵ}_n\right)\mathbf{E}(\mathbf{r},\omega ),\hfill \end{array}$$

(90)

$$\begin{array}{cc}\hfill {\mathbf{M}}_f(\mathbf{r},\omega )& ={\mathbf{M}}_{\mathrm{np}}(\mathbf{r},\omega )-{\mathbf{M}}_n\left(\omega \right)={\mu }_0^{-1}\left({\mu }_n^{-1}-{\mu }_{\mathrm{np}}^{-1}\right)\mathbf{B}(\mathbf{r},\omega ),\hfill \end{array}$$

(91)

where ${\mathbf{P}}_{\mathrm{np}}(\mathbf{r},\omega )$ and ${\mathbf{P}}_n\left(\omega \right)$ are the polarization densities of the nanoparticle and host medium, and ${\mathbf{M}}_{\mathrm{np}}(\mathbf{r},\omega )$ and ${\mathbf{M}}_n\left(\omega \right)$ are their densities, respectively. We used the linear constitutive relations ${\mathbf{P}}_n(\mathbf{r},\omega )={ϵ}_0\left({ϵ}_n-1\right)\mathbf{E}(\mathbf{r},\omega )$ and ${\mathbf{M}}_n(\mathbf{r},\omega )={\mu }_0^{-1}\left(1-{\mu }_n^{-1}\right)\mathbf{B}(\mathbf{r},\omega )$. The previous relation between the magnetic field induction and the magnetization follows from the equations $\mathbf{H}=\mathbf{B}/{\mu }_0-\mathbf{M}$ and $\mathbf{M}=\chi \mathbf{H}=\chi (\mathbf{B}/{\mu }_0-\mathbf{M})$, where $\chi$ is the magnetic susceptibility, then $\mathbf{M}(1+\chi )=\chi \mathbf{B}/{\mu }_0\iff \mathbf{M}=\frac{\chi }{1+\chi }\mathbf{B}/{\mu }_0\iff \mathbf{M}=\frac{\chi +1-1}{1+\chi }\mathbf{B}/{\mu }_0\iff \mathbf{M}=(1-{\mu }^{-1})\mathbf{B}/{\mu }_0$. The same reasoning provides the relation between the polarization and the electric field. Inserting the two previous equations in Equations (85) and (86), we obtain

$$\begin{array}{cc}\hfill \mathbf{E}(\mathbf{r},\omega )& ={\mathbf{E}}_0(\mathbf{r},\omega )+{\omega }^2{\mu }_n{\mu }_0{ϵ}_0\left({ϵ}_{\mathrm{np}}-{ϵ}_n\right){\int }_V{d}^3{\mathbf{r}}^{\prime }{\overleftrightarrow{G}}_0^{EE}\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)·\mathbf{E}({\mathbf{r}}^{\prime },\omega )\hfill \\ & +\omega {\mu }_n{k}_n\left({\mu }_n^{-1}-{\mu }_{\mathrm{np}}^{-1}\right){\int }_V{d}^3{\mathbf{r}}^{\prime }{\overleftrightarrow{G}}_0^{EM}\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)·\mathbf{B}({\mathbf{r}}^{\prime },\omega ),\hfill \end{array}$$

(92)

$$\begin{array}{cc}\hfill \mathbf{B}(\mathbf{r},\omega )& ={\mathbf{B}}_0(\mathbf{r},\omega )+{\mu }_n{k}_n^2\left({\mu }_n^{-1}-{\mu }_{\mathrm{np}}^{-1}\right){\int }_V{d}^3{\mathbf{r}}^{\prime }{\overleftrightarrow{G}}_0^{MM}\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)·\mathbf{B}({\mathbf{r}}^{\prime },\omega )\hfill \\ & -\omega {\mu }_n{\mu }_0{k}_n{ϵ}_0\left({ϵ}_{\mathrm{np}}-{ϵ}_n\right){\int }_V{d}^3{\mathbf{r}}^{\prime }{\overleftrightarrow{G}}_0^{ME}\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)·\mathbf{E}({\mathbf{r}}^{\prime },\omega ).\hfill \end{array}$$

(93)

The set of coupled Equations (92) and (93) are the Lippmann–Schwinger equations for electromagnetic scattering. Solving them, we can obtain the electric and magnetic fields scattered by the nanoparticle.

4.2. Weyl’s or Angular Spectrum Representation of Magnetic and Mixed Green’s Functions

Now we will see what is the Weyl’s (or angular spectrum) representation of the magnetic and mixed Green’s functions. The magnetic Green’s function is almost the same as the electric Green’s function, the only difference being the different additional $\overleftrightarrow{I}\delta \left(\mathbf{r}-{\mathbf{r}}^{\prime }\right)/k_n^2$ self-field term, which is isotropic and independent of the chosen excluded volume. Therefore, we can write

$${\overleftrightarrow{G_0}}^{MM}\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right)={\widehat{e}}_s{\widehat{e}}_s\frac{i}{2{\beta }_n}{e}^{i{\beta }_n\left|z-z^{\prime }\right|}+{\widehat{e}}_{p,n}^{\pm }{\widehat{e}}_{p,n}^{\pm }\frac{i}{2{\beta }_n}{e}^{i{\beta }_n\left|z-z^{\prime }\right|}+\frac{1}{k_n^2}\left(\overleftrightarrow{I}-{\widehat{e}}_z{\widehat{e}}_z\right)\delta \left(z-z^{\prime }\right).$$

(94)

We point out that the demagnetization term $\left(\overleftrightarrow{I}-{\widehat{e}}_z{\widehat{e}}_z\right)$ was previously obtained in [38]. As for the mixed Green’s function, their Weyl’s representation can be obtained by making the replacements: $g_0\left(\mathbf{r},{\mathbf{r}}^{\prime },\omega \right)\to g_0\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right)$, $\left({\partial }_x,{\partial }_y\right)\to i{\mathbf{p}}_{\parallel }$ and ${\partial }_z\to \pm i{\beta }_n$ for $z\gtrless z^{\prime }$. Therefore, we obtain

$${\overleftrightarrow{G_0}}^{EM}\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right)={\overleftrightarrow{G}}_0^{ME}\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right)=\frac{1}{k_n}\left[\begin{array}{ccc}0& \sigma {\beta }_n& -p_y\\ -\sigma {\beta }_n& 0& p_x\\ p_y& -p_x& 0\end{array}\right]\frac{i}{2{\beta }_n}{e}^{i{\beta }_n\left|z-z^{\prime }\right|},$$

(95)

where $\sigma =\pm 1$ for $z\gtrless z^{\prime }$. As for the electric and the magnetic Green’s functions, the mixed Green’s functions in the Weyl representation can also be written in terms of the s- and p-polarization vectors. It is straightforward to verify that

$$\frac{1}{k_n}\left[\begin{array}{ccc}0& \sigma {\beta }_n& -p_y\\ -\sigma {\beta }_n& 0& p_x\\ p_y& -p_x& 0\end{array}\right]={\widehat{e}}_{p,n}^{\pm }{\widehat{e}}_s-{\widehat{e}}_s{\widehat{e}}_{p,n}^{\pm },$$

(96)

which allows us to write

$${\overleftrightarrow{G_0}}^{EM}\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right)={\overleftrightarrow{G}}_0^{ME}\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right)=\left[{\widehat{e}}_{p,n}^{\pm }{\widehat{e}}_s-{\widehat{e}}_s{\widehat{e}}_{p,n}^{\pm }\right]\frac{i}{2{\beta }_n}{e}^{i{\beta }_n\left|z-z^{\prime }\right|}.$$

(97)

This representation is useful, as it allows for a simple interpretation of the emitted fields generated by the electric and magnetic dipoles in terms of s- and p–polarized electromagnetic waves.

If we are interested in the problem of scattering at a planar interface between two dielectric media with ${ϵ}_1$ for $z>0$, and ${ϵ}_2$ for $z<0$, we can construct reflected and transmitted Green’s functions expressed in terms of reflection and transmission coefficients, as done previously for ${\overleftrightarrow{G_0}}^{EE}\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right)$. However, some care must be taken in what the polarization vectors mean in the Green’s function, considering that the polarization of an electromagnetic field is usually defined by the polarization of the $\mathbf{E}$ field. The quantity ${\overleftrightarrow{G_0}}^{EM}\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right)$ gives us the electric field generated by a point magnetic dipole located at $z^{\prime }$. Therefore, the reflected and transmitted Green’s functions are constructed in the same way as for ${\overleftrightarrow{G_0}}^{EE}\left({\mathbf{p}}_{\parallel },z,z_0,\omega \right)$, and for $z_0>0$ we obtain

$$\begin{array}{cc}\hfill {\overleftrightarrow{G_r}}^{EM}\left({\mathbf{p}}_{\parallel },z,z_0,\omega \right)& =r_p\frac{i}{2{\beta }_1}{\widehat{e}}_{p,1}^{+}{\widehat{e}}_s{e}^{i{\beta }_1(z+z_0)}-r_s\frac{i}{2{\beta }_1}{\widehat{e}}_s{\widehat{e}}_{p,1}^{-}{e}^{i{\beta }_1(z+z_0)},\hfill \end{array}$$

(98)

$$\begin{array}{cc}\hfill {\overleftrightarrow{G_t}}^{EM}\left({\mathbf{p}}_{\parallel },z,z_0,\omega \right)& =t_p\frac{i}{2{\beta }_1}{\widehat{e}}_{p,2}^{-}{\widehat{e}}_s{e}^{-i{\beta }_{2}z}{e}^{i{\beta }_1{z}_0}-t_s\frac{i}{2{\beta }_1}{\widehat{e}}_s{\widehat{e}}_{p,1}^{-}{e}^{-i{\beta }_{2}z}{e}^{i{\beta }_1{z}_0}.\hfill \end{array}$$

(99)

For the magnetic Green’s functions, ${\overleftrightarrow{G}}^{MM}\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right)$ and ${\overleftrightarrow{G}}^{ME}\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right)$, we must take into account that these describe a field $\mathbf{B}$ generated by, respectively, a point magnetic and electric dipole. For electric and magnetic dipoles, ${\mathbf{d}}_0$ and ${\mathbf{m}}_0$, located at $z_0$, the primary magnetic field emitted for $z_0>z>0$ is given by

$$\begin{array}{cc}\hfill {\mathbf{B}}_0\left({\mathbf{p}}_{\parallel },z,\omega \right)& ={\mu }_1{\mu }_0{k}_n^2{\overleftrightarrow{G}}_0^{MM}\left({\mathbf{p}}_{\parallel }z,z_0,\omega \right)·{\mathbf{m}}_0-\omega {\mu }_1{\mu }_0{k}_n{\overleftrightarrow{G}}_0^{ME}\left({\mathbf{p}}_{\parallel }z,z_0,\omega \right)·{\mathbf{d}}_0\hfill \\ & =B_{0,s}{e}^{i{\beta }_1\left|z-z_0\right|}{\widehat{e}}_s+B_{0,p}{e}^{i{\beta }_1\left|z-z_0\right|}{\widehat{e}}_{p,1}^{-},\hfill \end{array}$$

(100)

with

$$\begin{array}{cc}\hfill B_{0,s}& ={\mu }_n{\mu }_0{k}_1^2\frac{i}{2{\beta }_1}\left({\widehat{e}}_s·{\mathbf{m}}_0\right)+\omega {\mu }_1{\mu }_0{k}_1\frac{i}{2{\beta }_n}\left({\widehat{e}}_{p,1}^{-}·{\mathbf{d}}_0\right),\hfill \end{array}$$

(101)

$$\begin{array}{cc}\hfill B_{0,p}& ={\mu }_n{\mu }_0{k}_1^2\frac{i}{2{\beta }_1}\left({\widehat{e}}_{p,1}^{-}·{\mathbf{m}}_0\right)-\omega {\mu }_1{\mu }_0{k}_1\frac{i}{2{\beta }_n}\left({\widehat{e}}_s·{\mathbf{d}}_0\right).\hfill \end{array}$$

(102)

The corresponding electric field can be obtained from Maxwell’s equations as ${\mathbf{E}}_0({\mathbf{p}}_{\parallel },z,\omega )=-{\omega }^{-1}{\mathbf{p}}_n^{\pm }\times {\mathbf{B}}_0({\mathbf{p}}_{\parallel },z,\omega )$. More explicitly (for $z>0$) we have for the primary field

$${\mathbf{E}}_0\left({\mathbf{p}}_{\parallel },z,\omega \right)=v_1{B}_{0,s}{e}^{i{\beta }_1\left|z-z_0\right|}{\widehat{e}}_{p,1}^{-}-v_1{B}_{0,p}{e}^{i{\beta }_1\left|z-z_0\right|}{\widehat{e}}_s.$$

(103)

This primary electric field is scattered by the interface at $z=0$, giving origin to a reflected field for $z>0$, which reads

$${\mathbf{E}}_r\left({\mathbf{p}}_{\parallel },z>0,\omega \right)=r_p{v}_1{B}_{0,s}{e}^{i{\beta }_1\left(z+z_0\right)}{\widehat{e}}_{p,1}^{+}-r_s{v}_1{B}_{0,p}{e}^{i{\beta }_1\left(z+z_0\right)}{\widehat{e}}_s,$$

(104)

and to a transmitted field for $z<0$

$${\mathbf{E}}_t\left({\mathbf{p}}_{\parallel },z<0,\omega \right)=t_p{v}_1{B}_{0,s}{e}^{-i{\beta }_{2}z}{e}^{i{\beta }_1{z}_0}{\widehat{e}}_{p,2}^{-}-t_s{v}_1{B}_{0,p}{e}^{-i{\beta }_{2}z}{e}^{i{\beta }_1{z}_0}{\widehat{e}}_s.$$

(105)

The corresponding magnetic fields can be obtaining using Faraday’s law applied to Equations (104) and (105). For example: if $\mathbf{E}=E_0{\widehat{e}}_{p,n}$ then $i\omega \mathbf{B}=\nabla \times \mathbf{E}=i{\mathbf{p}}_n^{-}\times {\widehat{e}}_{p,n}^{-}{E}_0=i{k}_n{\widehat{e}}_s{E}_0$, where $k_n=\left|p_n^{-}\right|$ and $E_0$ is the amplitude of the s-component of the field. The obtained magnetic fields are given by

$$\begin{array}{cc}\hfill {\mathbf{B}}_r\left({\mathbf{p}}_{\parallel },z>0,\omega \right)& =r_p{B}_{0,s}{e}^{i{\beta }_1\left(z+z_0\right)}{\widehat{e}}_s+r_s{B}_{0,p}{e}^{i{\beta }_1\left(z+z_0\right)}{\widehat{e}}_{p,1}^{+},\hfill \end{array}$$

(106)

$$\begin{array}{cc}\hfill {\mathbf{B}}_t\left({\mathbf{p}}_{\parallel },z<0,\omega \right)& =t_p\frac{v_1}{v_2}{B}_{0,s}{e}^{-i{\beta }_{2}z}{e}^{i{\beta }_1{z}_0}{\widehat{e}}_s+t_s\frac{v_1}{v_2}{B}_{0,p}{e}^{-i{\beta }_{2}z}{e}^{i{\beta }_1{z}_0}{\widehat{e}}_{p,2}^{-}.\hfill \end{array}$$

(107)

From the above Equations (106) and (107), we can, after replacing Equations (101) and (102) in Equations (106) and (107), obtain the reflected and transmitted magnetic Green’s functions, which are given by

$$\begin{array}{cc}\hfill {\overleftrightarrow{G}}_r^{MM}\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right)& =\frac{i}{2{\beta }_1}{e}^{i{\beta }_1\left(z+z_0\right)}\left[r_p{\widehat{e}}_s{\widehat{e}}_s+r_s{\widehat{e}}_{p,1}^{+}{\widehat{e}}_{p,1}^{-}\right],\hfill \end{array}$$

(108)

$$\begin{array}{cc}\hfill {\overleftrightarrow{G}}_r^{ME}\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right)& =\frac{i}{2{\beta }_1}{e}^{i{\beta }_1\left(z+z_0\right)}\left[r_p{\widehat{e}}_s{\widehat{e}}_{p,1}^{-}-r_s{\widehat{e}}_{p,1}^{+}{\widehat{e}}_s\right],\hfill \end{array}$$

(109)

$$\begin{array}{cc}\hfill {\overleftrightarrow{G}}_t^{MM}\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right)& =\frac{i}{2{\beta }_1}{e}^{-i{\beta }_{2}z}{e}^{i{\beta }_1{z}_0}\frac{v_1}{v_2}\left[t_p{\widehat{e}}_s{\widehat{e}}_s+t_s{\widehat{e}}_{p,2}^{-}{\widehat{e}}_{p,1}^{-}\right],\hfill \end{array}$$

(110)

$$\begin{array}{cc}\hfill {\overleftrightarrow{G}}_t^{ME}\left({\mathbf{p}}_{\parallel },z,z^{\prime },\omega \right)& =\frac{i}{2{\beta }_1}{e}^{-i{\beta }_{2}z}{e}^{i{\beta }_1{z}_0}\frac{v_1}{v_2}\left[t_p{\widehat{e}}_s{\widehat{e}}_{p,1}^{-}-t_s{\widehat{e}}_{p,2}^{-}{\widehat{e}}_s\right].\hfill \end{array}$$

(111)

Notice that the Fresnel reflection and transmission coefficients are defined for the electric field. With the last four equations we conclude the development of the formalism for the calculation of the renormalized polarizability [6] of a nanoparticle possessing both electric and magnetic dipoles.

5. Conclusions

In this paper we have studied the influence of two plasmonic structures on the effective polarizability of a nanoparticle made of either a metal (with a nearly dispersionless bare polarizability) or a polar dielectric or semiconductor (with a resonant polarizability due to polar optical phonons). The two studied structures are a continuous graphene sheet and a plasmonic graphene-based grating. In both cases a significant enhancement of the imaginary part of the polarizability has been observed. The two media possess plasmonic resonances which, however, occur at different frequencies. In the particular case of the grating, the resonance is tunable in two different ways: by adjusting the gate voltage and by changing the geometric parameters of the grating. In this case, it is possible to scan the resonance from the THz to the mid-IR range, whereas for the continuous graphene sheet the resonance is always in the THz range for the currently achieved values of electronic doping using a gate. The approach pursued here was to model the nanoparticle by a point-like dipole. The main motivation for this approach lies in its ability to make analytic progress. However, in real systems, one has a finite-size particle which can be modeled as an assembly of many point-like dipoles. These are determined by the coupled dipole equations [49]. In this case, the particle, even a spherical one, has other multipole resonances that can couple to the incoming radiation and contribute to the extinction cross-section. The two lowest multipoles, besides the electric dipole, are the magnetic dipole and the electric quadrupole. It can be shown numerically that for semiconductor nanoparticles such as spheres, cubes, pyramids, disks, and cylinders, the extinction cross-section has a strong magnetic-dipole resonance [4,5,6,23]. We note, however, that for semiconductor nanoparticles, if we consider interband transitions, that is, exciton resonances that are characteristic of semiconductors, the relevance of higher multipole resonances depends much more on the underlying band structure than on the shape. The formalism used in this paper to describe the renormalization of the electric dipole resonances can be extended to include the problem of magnetic dipole resonance [6], as we have seen in the previous section. The contribution to the extinction cross section of the magnetic dipole is given by ${\sigma }_{\mathrm{ext}}^{\mathrm{m}}=\frac{\omega \mu }{2S_{\mathrm{inc}}}\Im ([{\mathbf{H}}_0^{*}\left({\mathbf{r}}_{\mathbf{0}}\right)·\mathbf{m}\left({\mathbf{r}}_{\mathbf{0}}\right)]$, where $S_{\mathrm{inc}}$ is the power per unit area of the incoming radiation and $\mathbf{m}\left({\mathbf{r}}_{\mathbf{0}}\right)={\overleftrightarrow{\alpha }}_{MM}{\mathbf{H}}_0\left({\mathbf{r}}_{\mathbf{0}}\right)$, with ${\overleftrightarrow{\alpha }}_{\mathrm{MM}}$ being the effective magnetic polarizability of the particle and ${\mathbf{H}}_0$ the incoming magnetic field. The effective magnetic polarizability can be derived as done before for the electric dipole case. To that end, we will need the dyadic magnetic Green’s function which can be obtained by writing the wave equation for the magnetic field using the procedure outlined in Section 4. This study will be pursued in a forthcoming paper.