Analytical expressions for the Electromagnetic Dyadic Green’s Function in Graphene and thin layers

Analytical expressions for the Electromagnetic Dyadic Green’s Function in Graphene and thin layers

A. Yu. Nikitin, F. J. Garcia-Vidal, and  L. Martin-Moreno A. Yu. Nikitin and L. Martin-Moreno are with the Instituto de Ciencia de Materiales de Aragón and Departamento de Física de la Materia Condensada, CSIC-Universidad de Zaragoza, E-50009, Zaragoza, Spain (e-mail:; J. Garcia-Vidal is with the Departamento de Física Teórica de la Materia Condensada, Universidad Autónoma de Madrid, E-28049, Madrid, Spain (e-mail:

An analytical general analysis of the electromagnetic Dyadic Green’s Function for two-dimensional sheet (or a very thin film) is presented, with an emphasis on on the case of graphene. A modified steepest descent treatment of the fields from a point dipole given in the form of Sommerfeld integrals is performed. We sequentially derive the expressions for both out-of-plane and in-plane fields of both polarizations. It is shown that the analytical approximation provided is very precise in a wide range of distances from a point source, down to a deep subwavelength region ( of wavelength). We separate the contribution from the pole, the branch point and discuss their interference. The asymptotic expressions for the fields are composed of the plasmon, Norton wave and the components corresponding to free space.

Graphene, thin films, plasmon, Dyadic Green’s Function.

I Introduction

Electromagnetic properties of graphene have recently received a lot of attention due to a variety of application in photonics[1]. One of the attractive properties of graphene is its capacity to support highly-localized (nanometric) surface modes, i.e. graphene surface plasmons (GSPs) in terahertz (THz) and micro-wave frequency ranges [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]. GSPs and their gate tunability open interesting possibilities for construction of tunable meta-materials [4], and merging photonics and electronics. In particular, for applications related to strong light-matter interactions and biosensing, the interaction of a graphene sheet with a point emitter presents a special interest [5, 6, 7, 8, 9, 10, 11, 12]. Localized excitation has allowed experimental demonstration of GSPs [11, 12]. In comparison with other experimental techniques, local excitation of GSPs is more favorable due to very high values of GSP momentums.

The computation of patterns of the electromagnetic fields in graphene created by point source (as well as spontaneous emission rates), requires the knowledge of the Dyadic Green’s Function (DGF) [2, 3, 8, 9, 10]. Its calculation involves notoriously difficult Sommerfeld-type integrals, with the integrands containing quickly oscillating functions, poles and branch cuts [14, 13, 15].

In this paper we will perform an analytical analysis of DGF, providing an asymptotic series expansion with the modified steepest-decent method [13, 15]. We will show that this expansion, while being exact for long distances, in practice is very precise outside of its formal validity range. We will explicitly provide contributions for GSP, out-of-plane propagating waves and field components decaying algebraically along the graphene sheet. Notice that while all the examples correspond to the conductivity of graphene, the analytical expressions are applicable to any two-dimensional (2D) sheet. The analytical expression can be useful for treatment of more complicated problems, related for instance to Lippmann-Schwinger integral equation [14].

Fig. 1: The schematic of the studied system. A point source (dipole) is placed onto a graphene sheet, characterized by a two-dimensional conductivity . The colorplot represents an example of the spatial distribution of the electric field modulus for shown orientation of the dipole moment. Representative parameters for the graphene sheet corresponding to the formation of GSP are (see Appendix B for the definitions): , K, ps, eV.

Ii Formulation of the problem

Let us consider a dipole with an arbitrary oriented dipole moment placed at the point , with being the distance from a free-standing graphene sheet which covers the plane , see Fig.1. Without loss of generality, we will suppose that . The time dependency is supposed to be , where is the angular frequency. Throughout the article, we express coordinates in the in-plane () and normal () directions to the graphene sheet in dimensionless units as and , with being the free-space wavevector.

Graphene is represented by its in-plane complex conductivity . Throughout this paper, all examples are performed for based on the random-phase-approximation [16, 17, 18], see Appendix B.

The electric field emitted by our electric dipole, is given through the DGF by the following relation


In the next section we will provide the exact general expressions for .

Iii The general form of the Green’s dyadic

We have performed the analysis for a graphene placed onto the boundary of two different dielectrics. However, we have found that there is no much qualitative difference between this general case and free-standing (suspended) graphene. Since the analytical formula in the general case are quite lengthy, in this paper we will consider the DGF for a suspended graphene (taking also into account that free-standing samples are widely used in the experiments).

DGF satisfy the following differential equation


where is a diagonal unit matrix and is the Dirac delta function. The DGF must be complemented by the boundary conditions at the graphene sheet that we present below.

Iii-a Angular representation of DGF in Cartesian coordinates

The solution for can be can be expressed (see e.g. [14, 13]) in terms of plane waves in free-space , characterized by their in-plane momentum and polarization ). The unitary vectors characterizing the polarization of each mode are:


where “” (“”) upward (downward) propagation along respectively. is the normalized z-component of the wavevectors.

The DGF reads


with being the DGF in free space (FS),


and , being the contributions due to the reflection and transmission in our 2D system


In these expressions the superscript “T” stands for transposition.

In the above expressions and are the reflection and transmission coefficients for graphene. These coefficients can be found by matching the magnetic and electric fields through the boundary conditions:


where () and () stay for the electric (magnetic) fields in the regions of negative and positive , respectively, and is the unitary vector along the direction. As a result of the matching we have


with , being the dimensionless 2D conductivity.

Explicitly, we have for


Analogously, reads


and finally, for the transmission part we have


Eqs. (9)-(11) present the angular representation of DGF in Cartesian coordinates. By transforming a cylindrical coordinate system, these expression can be greatly simplified.

Fig. 2: Transformation of the coordinates.

Iii-B DGF in cylindrical coordinates

As was previously shown, the Purcell factor (the total decay rate normalized to the free space decay rate) of the point emitter placed directly over the graphene monolayer diverges due to losses through the evanescent waves with large -components [6, 7, 10]. This means that the dipole placed directly on the graphene surface is quenched. However, the parameter that accounts for the efficiency of the coupling to GSP, (-factor, defined as the ratio of the emitter’s decay rate through GSP to its total decay rate) has an optimum value in the region of very small distances from the dipole to graphene: (see [10]). Then, taking into account that the problems related to the high values of -factor are relevant, the range of small distances presents a special interest. Another important point is that the field patterns at the distances do not differ essentially upon the field patterns created by a dipole lying directly on the monolayer. Therefore, for the above two reasons, in this paper we will consider the dipole placed directly onto the graphene sheet.

We would like to notice, that some physical systems can be reduced to a problem of a dipole lying directly on the graphene sheet. For instance, a subwavelength aperture in graphene sheet can be represented by an effective dipole placed directly onto the sheet, (for comparison with the case of metal films see [21, 22]).

Additionally, since the analytical treatment of both reflection () and transmission () parts of the DGF is similar, we will derive the expressions for the transmission part, .

So, supposing that in the previous expressions, the Green’s dyadic can be simplified. The symmetry of the problem makes it convenient to work in cylindrical coordinates , see Fig. 2 (a):


In this system of coordinates the Green’s dyadic can be obtained from the one in cartesian coordinates through




Having performed this transformation, the DGF in cylindrical coordinates is expressed in the form of Sommerfeld integrals as :


where the subscripts “s” and “p” correspond to TE and TM polarizations respectively. In this expressions, and are Bessel functions of th order. Equations (15) can be treated numerically. We show some recipes for the integration in the complex plane in Appendix A. Let us now proceed with the asymptotic expansion of the DGF.

Iv Asymptotic expansion of the DGF

In this section we will derive explicit asymptotic expressions for the elements of DGF following the steepest-decent method, modified in order to take into account the presence of both poles and branch points close to the integration path [13, 21].

First, using the identity


where are Hankel functions, we can extend the limit of integration to the whole real -axis in (15)


where according to Eqs. (15) function is odd/even for even/odd values of . In (17) we have used the symmetry of the Hankel functions .

Second, we use the asymptotic form of the Hankel functions for large arguments (). Notice that the region that provides the major contribution to the integral corresponds to . Then the formal condition of the asymptotic expansion validity for the lower value of the contributing reads as . However, as we will show below, the true region of distances where the asymptotic approximation is valid is much less restricted. Retaining the first two terms in , the expansion reads


For convenience, let us normalize the DGF as follows


Then using (16)-(17), we find from (15) the expressions for :


with and the superscripts of indicate distinct dependencies upon . The denominators in the integral are


and nominators


where , , represent the unit vectors.

Fig. 3: Integration contours in the complex plane of the integration variables [panel (a)] and [panel (b)] for different values of angle . In panels (a) and (b) the initial integration contour corresponds to and respectively. The pole position is shown by a circular symbol. In the complex plane of the initial path is the same and the steepest-decent path depends upon , while in the complex plane the steepest-decent path is the same, , and the initial integration path changes with . The position of the pole in the complex plane is dependent upon .

As we see, the integrands contain branch points corresponding to , and branch cuts corresponding to . In order to remove these problematic peculiarities, we perform the following standard change of integration variable and coordinates (see Fig. 2 (b)):


The integrals (20) then transforms as


where the integration contour pases through the complex plane , see Fig. 3 (a) and corresponds to the real axis in the complex -plane. At this stage we can perform a steepest-decent integration. Then the integration path must be transformed to , see Fig. 3 (a), with the saddle point . However, the steepest-decent integration is much easier in another complex variable plane. This variable is given as follows


so that . With this change of variable becomes


where the integration contour is shown in see Fig. 3 (b). The steepest-decent integration path is now simply given by , and the saddle point is located in the origin, . Notice that the branch points (corresponding to ) are given by in the complex plane , while in the -plane they are located at . For the branch point and saddle point in -plane coincide.

An important point here is that the elements of the integrant dyadic in are singular due to the presence of the poles [ and ] in the denominators. These poles are located at and respectively. They correspond to TM-surface wave (GSP) and TE surface wave [2]


The positions of the poles depend essentially upon the value of the normalized conductivity . TM waves correspond to , while TE ones correspond to , so that TM and TE surface waves cannot exist at the same frequency see [2]. High values of correspond to large while low values of correspond to large . In principle, taking into account a wide range of metamaterials that are available at present, a wide range of is also accessible. A mathematical treatment of the problem corresponding to a three-dimensional (3D) layer of a very thin thickness with the dielectric permittivity can be reduced to the case of a 2D sheet with an effective 2D normalized conductivity . The relation between 2D effective conductivity and 3D permittivity is established from the comparison of the Fresnel coefficients and reads as .

Returning to the case of graphene, due to small values of in graphene, TE surface waves are very weakly bounded () and therefore they virtually do not couple to the point emitter. This means that, in practice, the contribution from TE pole can be neglected in graphene [8]. Nevertheless, taking into account a wide range of possible 2D sheets, in our analysis we retain the contribution from both poles.

In order to proceed with the series expansion, the singular terms in the integrands must be separated. The separation for the dyadics into a pole and a smooth part is as follows:


where are dyadics with the elements corresponding to the residues of , which after some algebra can be computed from Eq. (26) as:


Let us explicitly write out the expressions for the residue dyadics


Now we can deform the integration contour in the complex plane, into the real axis. Then the singular terms in Eq. (28) for that enter to the integral (28) can be integrated analytically. The result for can be presented in the form of a sum


where the term with the complementary error function is due the singularity. This function includes the contribution from the pole that is automatically taken into account when the pole is crossed by the transformation of the integration contour. The second term in Eq. (31), presents a nonsingular contribution


The integral appearing in Eq. (32) is of the Gauss type, so the functions can be expanded in Tailor series close to and integration of every term is easily performed with the following result


where is Gamma function. Retaining in this expression the terms up to order (which is enough for the most of the cases) yields


Recalling that the saddle point corresponds to and therefore to , we can explicitly write out the dyadics


where we have taken into account the identity . The explicit expression for the third term in (34) is more cumbersome for arbitrary , therefore we give its formal expression, involving derivatives of previously defined functions


The expressions (31), with given by (34), present the analytical approximation of the DGF. Let us now analyze different terms in this expression and the applicability of the approximation.

V Analysis of the analytical solution

In its general form, the analytical solution presents a non-trivial combination of the contribution from the pole and saddle point. The “interaction” between these contributions depends both upon the distance between the saddle point and the pole in the complex plane and upon the physical distance responsible for the oscillations of the integrand. The parameter that measures the interaction between the saddle point and the pole is called by Sommerfeld “numerical distance” and its square presents the argument of the complementary error function, .

In the -plane, the saddle point corresponds to the condition of the extremum of the exponential phase in (20), i.e. to , or . This can be considered as an equation for the rays in “Ray Optics” (RO). If the contribution of the pole is neglected, then the same result can be derived following the standard stationary phase evaluation. In the far field, the leading RO contribution corresponds to the first term in (34) with replaced by . It reads


Returning to the DGF via Eq. (19), we have explicitly [using Eqs. (21), (22) and the relation ]


When the graphene sheet disappears, , we recover the spherical wave term () of DGF corresponding to a dipole in FS. Notice that at (for the fields along ), the elements , , in vanish independently upon whether the graphene sheet is present or not. In contrast, the element in at is very sensitive to the presence of graphene. It takes non-zero values for free space, and vanishes for . This property of the element is similar to the diffraction shadow effect in metals due to the presence of surface modes and formation of the Norton waves[19]. In case of metals, however, the diffraction shadow appears for the TM part of the dyadic, while in thin films this takes place for the TE part.

The fact that some elements of the dyadic vanish indicates that other terms in the general solution for must be considered in order to provide the correct far-field representation of the DGF and electric fields. Let us consider in details the case of .

V-a Asymptotic expansion of DGF in graphene plane,  ()

Here we present simplified expressions for at . Recall that the general expression is given by the Eq. (31), with given by Eq. (34). The functions at read


Introducing for a brevity of notations the following dyadic


The derivatives (36) simplify as


Now we can explicitly write the full expression for the TE dyadic


and for the TM one


Recall that in Eqs. (42)-(43) the residue dyadics are given by Eq. (30), and have the form of Eq. (41). The locations of the poles in the complex -plane at read


We have checked that the Eqs. (42), (43) transform to the DGF of FS in the limit . To perform this limit, one should carefully expand all the coefficients taking into account that for small the poles become , . In particular, for TM part of the DGF the expansion for large arguments of the complementary error function must be taken.

V-B Numerical check on the validity of the analytical approximation

Let us present some illustrative examples that demonstrate the validity of the analytical approximation. For this we directly compare numerical and analytical calculations for two components of DGF in the most unfavorable situation, i.e. for (), when the “RO” contribution disappears. We perform the precise (converged) numerical calculations according to Appendix A. The analytical approximation is given by the expressions (42), (43).

The conductivity of graphene is a function of frequency, , chemical potential , temperature and scattering time , see Appendix B. For the illustration, we consider the room temperature K and eV (48THz), which are typical experimental values.

Fig. 4: Comparison between the numeric and analytic calculations of and DGF elements at . The main figures show both real and imaginary parts of and as a function of distance. The values of these elements have been normalized to the maximal values of their modules in the shown range of distances. The parameters for graphene in both panels are: K, eV, ps, THz. The insets show the dependencies of the relative error upon the frequency for different distances. In the insets and are the same as in main figures.

The comparison is shown in Fig. 4. The range of the distances corresponds to the subwavelength region (outside of this region the the numerical and analytical curves are virtually undistinguishable). In order to characterize the difference between numerical, , and analytical, , results, in the insets we have represented the relative error in THz frequency range. According to the insets to Fig. 4, the error is a non-monotonous function of frequency. However, it has a decaying tendency with frequency increase. To understand such behavior, let us notice that for higher frequencies decreases so that both the real and imaginary parts of increase. In the lower limit ( THz) so that , while in the upper limit ( THz) so that . Thus, the propagation length of the GSP, , decreases due to increase of . This leads to a strong spacial decay of the terms in the solution (at the deep sub-wavelength distance), related to the GSP field components. Since the analytical solution recovers the FS DGF (up to ), the coincidence between the analytical and exact solutions improves for higher frequencies.

Vi Long distance limit. Surface modes and algebraically-decaying components

In the region of parameters, where the argument of the complementary error function is a large number (the numerical distance is large, ), this function can be substituted by a few terms from its asymptotic expansion