Boundary element method for 2D materials and thin films

M. Hrtoň; V. Křápek; T. Šikola

doi:10.1364/OE.25.023709

1. Introduction

The critical point in the research of 2D materials was without doubt the successful fabrication of graphene by Geim and Novoselov in 2004 [1]. Since then this field underwent a rapid growth and 2D materials have become a major research area. The efforts were driven mostly by graphene and its superb mechanical stability and durability, excellent transport properties, and exceptionally strong interaction with light [2,3].

Among other things, graphene is also a key material for nanophotonics, as it supports plasmonic oscillations with very long lifetimes in a wide range of energies [4]. The key element of graphene plasmonics is the possibility to dynamically tune the optical properties of graphene by changing its Fermi level via electrostatic gating. Alongside the relatively low losses and broad tunability, plasmons in graphene exhibit high degree of spatial confinement accompanied with huge electric field enhancement [5]. These unique features have already found their use in mid-infrared optoelectronics [6], bio-sensing [7, 8], fingerprinting [9] and beam shaping [10] and there are many other applications, in which graphene shows a lot of promise [11].

Besides graphene, there is a new emerging class of materials such as topological insulators [12], black phosphorus [13–15], hexagonal boron nitride (hBN) [16] and transition metal dichalcogenides [17]. These materials exhibit 2D character and also support polaritonic excitations.

The ongoing research of interaction between light and all these 2D materials is inseparably connected with the need for an efficient and robust simulation tool for solving Maxwell’s equations of electromagnetism. There are several techniques suited for this task such as the finite-difference time domain method (FDTD) [18], the discontinuous Galerkin time-domain method (DGTD) [19], the finite element method (FEM) [20, 21] or the boundary element method (BEM) [22–25], each of them with its own advantages and limitations [26,27]. When dealing with graphene and other 2D materials, we face the problem of modelling an atomically thick material within 3D simulations. One solution is to treat the 2D material as a very thin 3D film with a dielectric function constructed from its 2D dynamical conductivity σ_2D(ω). This approach proved to be successful in predicting the optical response of many graphene structures ranging from ribbons and discs [4] to plasmonic dimers [28] and graphene-loaded gold antennas [29, 30]. The other solutution is based on incorporating the 2D conductivity directly into the boundary conditions between individual discretization elements and this concept has been already implemented in FDTD [31], DGTD [32] and FEM [33].

The boundary element method, on the other hand, still lacks such an extension, even though it could be very beneficial: The time-complexity of BEM scales with N³ (N being the number of discretization elements) and the reduction in N associated with the transition from a very thin 3D film to a single 2D interface is substantial. Therefore, to provide a more efficient approach for modelling of materials such as graphene within boundary element method, we decided to develop a modification that would inherently incorporate 2D materials into the boundary conditions for the electromagnetic fields at the object interfaces. We succeeded in formulating a new set of integral Maxwell’s equations that correctly describe the interaction of light with both extended and structered graphene of an arbitrary shape.

The following sections briefly introduce the apparatus of BEM and outline our solution to the aforementioned issue. Subsequently, we test it using our own implementation in Matlab on problems with known analytical solutions. In order to manifest the ability of the modified BEM algorithm to model light scattering involving inhomogeneously doped graphene [34], we employ it to study highly localized plasmons in side-gated graphene ribbons predicted by Thongrattanasiri et al. [35].

2. Principles of BEM

This section comprises the most basic definitions and equations of electromagnetism involved in the boundary element method. In the course of the following derivation we shall mostly adopt the notation introduced by Abajo and Howie [36], with differences arising mainly due to giving preference to SI units. The starting point are the Maxwell’s equations of electromagnetism in the absence of magnetic materials and written in the frequency domain ω

\begin{array}{l} \nabla \cdot (ε E) & = & \frac{ρ}{ε_{0}} & \nabla \times B & = & μ_{0} j - i \frac{ω}{c^{2}} ε E, \\ \nabla \cdot B & = & 0, & \nabla \times E & = & i ω B, \end{array}

where the electric field E and the magnetic field B are related to the charge ρ and current j densities and to each other (c denotes the speed of light). BEM allows to solve the Maxwell’s equations for a finite system of discrete homogeneous objects (each described with a local and spatially invariant dielectric function ε (ω)). Within the method, Maxwell’s equations are transformed to a set of integro-differential equations for the auxiliary charges and currents at the object boundaries. In this way, the dimensionality of the problem is reduced to two at a price of increased number of independent variables and limited applicability to non-discrete systems.

The last two Maxwell’s equations allow us to express E and B in terms of scalar and vector potentials Φ and A

E = - \nabla Φ + i k A,

B = \frac{1}{c} \nabla \times A,

i k ε Φ = \nabla \cdot A,

with

k = \frac{ω}{c}

. Note the factor 1/c in Eq. (2) which was introduced to preserve the notation encountered in Ref. [36].

Assuming that the investigated area is composed of a finite number of discrete objects, one can write the solution to the above equations in terms of auxiliary surface charge and current densities σ and h defined at the object interfaces

Φ_{j} (r, ω) = Φ_{j}^{ext} (r, ω) + \int_{\partial Ω_{j}} d s G_{j} (r, s, ω) σ_{j} (s, ω),

A_{j} (r, ω) = A_{j}^{ext} (r, ω) + \int_{\partial Ω_{j}} d s G_{j} (r, s, ω) h_{j} (s, ω),

where the subscript j is used to discriminate between the distinct objects and the terms

Φ_{j}^{ext}

,

A_{j}^{ext}

account for the external fields, as well as for possible charges and currents embedded within the volume of the objects. The integration is performed over the entire surface of the particular object ∂Ω_j. The Green’s function G_j (r, r′, ω) is solution to the wave equation

(\nabla^{2} + k_{j}^{2}) G_{j} (r, r^{'}, ω) = - δ (r - r^{'}),

with

k_{j} = k \sqrt{ε_{j}}

. The values of the auxiliary surface charge and current densities are determined via the requirement that the electromagnetic fields and potentials satisfy the customary boundary conditions at the interfaces between two media

Φ_{2} - Φ_{1} = 0,

A_{2} - A_{1} = 0,

n_{s} \cdot (ε_{2} E_{2} - ε_{1} E_{1}) = 0,

n_{s} \times (B_{2} - B_{1}) = 0,

where the orientation of the surface normal vector n_s is chosen to be directed from the medium 1 towards the medium 2. Even though one could apply the above conditions as they are, it proves to be beneficial to transform the last one into a more suitable form.

Starting with the identity ∇ (n_s · A) = (n_s · ∇) A + n_s × (∇ × A) + (A · ∇) n_s + A × (∇ × n_s) and exploiting the definition of the vector potential [Eq. (2)] and its continuity across any interface [Eq. (8)] one obtains

\begin{array}{r} \nabla [n_{s} \cdot (A_{2} - A_{1})] = (n_{s} \cdot \nabla) (A_{2} - A_{1}) + \\ + c n_{s} \times (B_{2} - B_{1}) . \end{array}

Let us now take projections of the above equation along the surface normal n_s and the two independent tangential vectors

t_{s}^{1}

and

t_{s}^{2}

n_{s} \cdot \nabla [n_{s} \cdot (A_{2} - A_{1})] = n_{s} \cdot (n_{s} \cdot \nabla) (A_{2} - A_{1}),

\begin{array}{r} t_{s}^{ν} \cdot \nabla [n_{s} \cdot (A_{2} - A_{1})] = t_{s}^{ν} \cdot (n_{s} \cdot \nabla) (A_{2} - A_{1}) + \\ + c t_{s}^{ν} \cdot [n_{s} \times (B_{2} - B_{1})] . \end{array}

Finally, by employing the identity ∇ · [n_s (n_s · A)] = (n_s · A) (∇ · n_s) + n_s · ∇ (n_s · A), the continuity of the vector potential and its tangential derivatives, the Lorentz gauge condition [Eq. (3)] and the continuity of the tangential magnetic field [Eq. (10)] one arrives at

n_{s} \cdot (n_{s} \cdot \nabla) (A_{2} - A_{1}) - i k (ε_{2} Φ_{2} - ε_{1} Φ_{1}) = 0,

\begin{array}{l} t_{s}^{ν} \cdot \nabla [n_{s} \cdot (A_{2} - A_{1})] - \\ - t_{s}^{ν} \cdot (n_{s} \cdot \nabla) (A_{2} - A_{1}) = 0 . \end{array}

In the standard implementation of BEM, Eqs. (7), (8), (9), (14) and (15) constitute together a system of eight equations for the eight unknown auxiliary quantities σ₁, σ₂, h₁ and h₂.

3. η-BEM: Additional surface current density approach

The presence of graphene or any other 2D material can be established in two different manners. The first, more general approach requires introduction of an additional surface current density, thus increasing both the number of unknown quantities and the number of equations by three. The alternative approach resides in incorporating the conductivity of the 2D material directly into the boundary conditions, which leaves the number of unknowns at eight, but proves to be inapplicable in certain types of problems.

Starting with the first of the two considered approaches, the current density η induced within graphene is proportional to the local electric field via the 2D optical conductivity σ_2D

η = μ_{0} c σ_{2 D} E_{‖},

with the symbol ‖ indicating that only the component parallel to the 2D sheet is considered. The electric field can be split into three distinct contributions, i.e. the external electric field E^ext, the electric field due to auxiliary surface charge and current densities, and the electric field produced by 2D material itself

η = μ_{0} c σ_{2 D} {[E^{ext} + E (σ, h) + i k \int_{S_{G}} d s \overset{\leftrightarrow}{G} (r, s), η (s)]}_{‖} .

where G⃡ denotes the free space dyadic Green’s function [37] and the integration is performed over the area of the 2D material sheet S_G. Importantly, the presence of the 2D material is fully contained within the new auxiliary current density, while the other quantities σ and h only serve to adjust the potentials to the discontinuities in the dielectric properties of space as in the original BEM. Notice the ambiguity in the calculation of the tangential electric field. Since it is continuous across the interface, one can evaluate it on either side of the interface with the same final result. We have chosen to take the average of the two fields, as it gives the system of equations certain sense of symmetry. Taking this into consideration the previous equation reads

\begin{array}{l} η & = & \frac{1}{2} μ_{0} c σ_{2 D} [E_{1}^{ext} + E_{2}^{ext} - (\nabla G_{1} σ_{1} + \nabla G_{2} σ_{2}) + \\ {+ i k (G_{1} h_{1} + G_{2} h_{2}) + i k ({\overset{\leftrightarrow}{G}}_{1} + {\overset{\leftrightarrow}{G}}_{2}) η]}_{‖}, \end{array}

where the matrix notation has been adopted and for the sake of brevity shall be used from now on. In this representation matrix and vector indices play the role of spatial coordinates and matrix-vector products are equivalent to integration over surface (for a more detailed description see Ref. [36]).

Apart from Eq. (18), the current density η also enters the expressions for the scalar and vector potentials

Φ_{j} = Φ_{j}^{ext} + G_{j} σ_{j} + \frac{1}{i k ε_{j}} \nabla G_{j} \cdot η,

A_{j} = A_{j}^{ext} + G_{j} h_{j} + G_{j} η,

where the last term in Eq. (19) follows from the Lorentz gauge condition. Moreover, the tangential magnetic field now exhibits a discontinuity at the interface due to the non-vanishing surface current η

n_{s} \times (B_{2} - B_{1}) = \frac{1}{c} η .

Adopting these changes, the original system of boundary conditions becomes

(\frac{1}{i k ε_{1}} \nabla G_{1} - \frac{1}{i k ε_{2}} \nabla G_{2}) \cdot η + G_{1} σ_{1} - G_{2} σ_{2} = Φ_{2}^{ext} - Φ_{1}^{ext},

G_{1} h_{1} - G_{2} h_{2} + (G_{1} - G_{2}) η = A_{2}^{ext} - A_{1}^{ext},

\begin{array}{l} H_{2} ε_{2} σ_{2} - H_{1} ε_{1} σ_{1} + i k n_{s} \cdot (G_{1} ε_{1} h_{1} G_{2} ε_{2} h_{2}) + \\ + i k n_{s} \cdot [({\overset{\leftrightarrow}{G}}_{1} ε_{1} - {\overset{\leftrightarrow}{G}}_{2} ε_{2}) η] = D_{2}^{ext} - D_{1}^{ext}, \end{array}

\begin{array}{l} n_{s} \cdot (H_{1} h_{1} - H_{2} h_{2}) + i k (G_{2} ε_{2} σ_{2} - G_{1} ε_{1} σ_{1}) + \\ + [n_{s} (H_{1} - H_{2}) + \nabla G_{2} - \nabla G_{1}] \cdot η = α_{2} - α_{1}, \end{array}

\begin{array}{l} (n_{s} T_{1}^{ν} - t_{s}^{ν} H_{1}) \cdot h_{1} - (n_{s} T_{2}^{ν} - t_{s}^{ν} H_{2}) \cdot h_{2} + \\ + [n_{s} (T_{1}^{ν} - T_{2}^{ν}) - t_{s}^{ν} (H_{1} - H_{2} - 𝕀)] \cdot η = β_{2}^{ν} - β_{1}^{ν}, \end{array}

where the following quantities have been introduced

𝕀 = identity matrix,

H_{j} = n_{s} \cdot \nabla G_{j},

T_{j}^{ν} = t_{s}^{ν} \cdot \nabla G_{j},

D_{j}^{ext} = ε_{j} n_{s} \cdot (- \nabla Φ_{j}^{ext} + i k A_{j}^{ext}),

α_{j} = n_{s} \cdot (n_{s} \cdot \nabla) A_{j}^{ext} - i k ε_{j} Φ_{j}^{ext},

β_{j}^{ν} = t_{s}^{ν} \cdot \nabla (n_{s} \cdot A_{j}^{ext}) - t_{s}^{ν} \cdot (n_{s} \cdot \nabla) A_{j}^{ext} .

The major setback of the presented approach is apparently the coupling of previously uncoupled equations, which substantially increases the number of algebraic operations needed to solve the problem.

Before proceeding with the alternative approach we should address the generalization of the method to an arbitrary number of media. First of all, the Green’s function G_j and its derivatives are nonzero only if the two points it connects belong to the same homogenous region of space filled with the medium j. Denoting j₂ and j₁ the materials in the directions parallel and anti-parallel to the surface normal, the expressions such as G₁ε₁σ₁ − G₂ε₂σ₂ then read

\begin{array}{l} \int d s^{'} {G_{j_{1} (s)} ε_{j_{1} (s)} [σ_{1} δ_{j_{1} (s) j_{1} (s^{'})} + σ_{2} δ_{j_{1} (s) j_{2} (s^{'})}] - \\ - G_{j_{2} (s)} ε_{j_{2} (s)} [σ_{1} δ_{j_{2} (s) j_{1} (s^{'})} + σ_{2} δ_{j_{2} (s) j_{2} (s^{'})}]} = \\ = \int d s^{'} {G_{j_{1} (s^{'})} ε_{j_{1} (s^{'})} [δ_{j_{1} (s) j_{1} (s^{'})} - δ_{j_{2} (s) j_{1} (s^{'})}] σ_{1} - \\ - G_{j_{2} (s^{'})} ε_{j_{2} (s^{'})} [δ_{j_{2} (s) j_{2} (s^{'})} - δ_{j_{1} (s) j_{2} (s^{'})}] σ_{2}}, \end{array}

where the Kronecker delta δ_ij = 1 if i = j and zero otherwise and argument in the subscript specifies the location at which the dielectric properties are taken. Similarly, the terms from Eq. (18) such as ∇G₁σ₁ + ∇G₂σ₂ become

\begin{array}{l} \int d s^{'} {\nabla G_{j_{1} (s^{'})} [δ_{j_{1} (s) j_{1} (s^{'})} + δ_{j_{2} (s) j_{1} (s^{'})}] σ_{1} + \\ + \nabla G_{j_{2} (s^{'})} [δ_{j_{2} (s) j_{2} (s^{'})} + δ_{j_{1} (s) j_{2} (s^{'})}] σ_{2}} . \end{array}

Apparently, one can retain the simple form of Eqs. (18), (22) – (26) by replacing the Green’s function and its derivatives with

\begin{array}{l} G_{1} \to {{}^{\pm}G}_{1} = G_{j_{1} (s^{'})} [δ_{j_{1} (s) j_{1} (s^{'})} \pm δ_{j_{2} (s) j_{1} (s^{'})}] \\ G_{2} \to {{}^{\pm}G}_{2} = G_{j_{2} (s^{'})} [δ_{j_{2} (s) j_{2} (s^{'})} \pm δ_{j_{1} (s) j_{2} (s^{'})}] . \end{array}

so that the system of equations becomes

\begin{array}{l} η = \frac{1}{2} μ_{0} c σ_{2 D} [E_{1}^{ext} + E_{2}^{ext} - (\nabla {{}^{+}G}_{1} σ_{1} + \nabla {{}^{+}G}_{2} σ_{2}) + \\ {+ i k ({{}^{+}G}_{1} h_{1} + {{}^{+}G}_{2} h_{2}) + i k ({}^{+}{\overset{\leftrightarrow}{G}}_{1} + {}^{+}{\overset{\leftrightarrow}{G}}_{2}) η]}_{‖}, \end{array}

(\frac{1}{i k ε_{1}} \nabla {{}^{-}G}_{1} - \frac{1}{i k ε_{2}} \nabla {{}^{-}G}_{2}) \cdot η + {{}^{-}G}_{1} σ_{1} - {{}^{-}G}_{2} σ_{2} = Φ_{2}^{ext} - Φ_{1}^{ext},

{{}^{-}G}_{1} h_{1} - {{}^{-}G}_{2} h_{2} + ({{}^{-}G}_{1} - {{}^{-}G}_{2}) η = A_{2}^{ext} - A_{1}^{ext},

\begin{array}{l} {{}^{-}H}_{2} ε_{2} σ_{2} - {{}^{-}H}_{1} ε_{1} σ_{1} + i k n_{s} \cdot ({{}^{-}G}_{1} ε_{1} h_{1} - {{}^{-}G}_{2} ε_{2} h_{2}) + \\ + i k n_{s} \cdot [({}^{-}{\overset{\leftrightarrow}{G}}_{1} ε_{1} - {}^{-}{\overset{\leftrightarrow}{G}}_{2} ε_{2}) η] = D_{2}^{ext} - D_{1}^{ext}, \end{array}

\begin{array}{l} n_{s} \cdot ({{}^{-}H}_{1} h_{1} - {{}^{-}H}_{2} h_{2}) + i k ({{}^{-}G}_{2} ε_{2} σ_{2} - {{}^{-}G}_{1} ε_{1} σ_{1}) + \\ + [n_{s} ({{}^{-}H}_{1} - {{}^{-}H}_{2}) + \nabla {{}^{-}G}_{2} - \nabla {{}^{-}G}_{1}] \cdot η = α_{2} - α_{1}, \end{array}

\begin{array}{l} (n_{s} {}^{-}{T_{1}^{ν}} - t_{s}^{ν} {{}^{-}H}_{1}) \cdot h_{1} - (n_{s} {}^{-}{T_{2}^{ν}} - t_{s}^{ν} {{}^{-}H}_{2}) \cdot h_{2} + \\ + [n_{s} ({}^{-}{T_{1}^{ν}} - {}^{-}{T_{2}^{ν}}) - t_{s}^{ν} ({{}^{-}H}_{1} - {{}^{-}H}_{2} - 𝕀)] \cdot η = β_{2}^{ν} - β_{1}^{ν} . \end{array}

4. BC-BEM: 2D conductivity incorporated into the boundary conditions

The system of equations presented in the previous section allows for the presence of a 2D material but at the cost of an increase in the number of unknowns and their mutual coupling. The alternative is to adjust the original customary boundary conditions so that they contain the 2D material inherently, i.e. let the current induced within the 2D material act as a discontinuity in the tangential magnetic field and the accumulated 2D charge density as a discontinuity in the normal displacement (the latter being the consequence of the continuity equation)

n_{s} \times (B_{2} - B_{1}) = μ_{0} σ_{2 D} E_{‖},

n_{s} \cdot (ε_{2} E_{2} - ε_{1} E_{1}) = - \frac{i}{ω ε_{0}} \nabla \cdot (σ_{2 D} E_{‖}) .

It is instructive to expand the right hand side of the last equation in terms of the surface normal vector

\begin{array}{l} - \frac{i}{ω ε_{0}} \nabla \cdot (σ_{2 D} E_{‖}) & = & - \frac{i}{ω ε_{0}} \nabla σ_{2 D} \cdot [E - n_{s} (n_{s} \cdot E)] - \\ - \frac{i σ_{2 D}}{ω ε_{0}} \nabla \cdot [E - n_{s} (n_{s} \cdot E)] . \end{array}

The first term containing the gradient of the 2D optical conductivity is largely responsible for the limited applicability of the approach based on adjustment of the boundary conditions. In case of structured planar 2D materials such as ribbons and discs, the abrupt changes in the conductivity at its edges give rise to singularities which are difficult to treat. Moreover, in case of structures with curved surface, e.g. cylinders or spheres, the spatial derivatives of surface normal vector are no longer zero and new terms specific to each geometry have to be included. Nevertheless, this method can still find its use in geometries containing for example extended graphene sheets with constant or smoothly varying Fermi level.

Restricting ourselves to infinite planar 2D materials and adopting the above boundary conditions, Eqs. (7), (8), (9), (14) and (15) can be recast as

{{}^{-}G}_{1} σ_{1} - {{}^{-}G}_{2} σ_{2} = Φ_{2}^{ext} - Φ_{1}^{ext},

{{}^{-}G}_{1} h_{1} - {{}^{-}G}_{2} h_{2} = A_{2}^{ext} - A_{1}^{ext},

\begin{array}{l} {{}^{-}H}_{2} ε_{2} σ_{2} - {{}^{-}H}_{1} ε_{1} σ_{1} + i k n_{s} ({{}^{-}G}_{1} ε_{1} h_{1} - {{}^{-}G}_{2} ε_{2} h_{2}) - \\ - \frac{1}{2} \frac{i σ_{2 D}}{ω ε_{0}} [δ (s - s^{'}) (σ_{1} + σ_{2}) + {(n_{s} \cdot \nabla)}^{2} ({{}^{+}G}_{1} σ_{1} + {{}^{-}G}_{2} σ_{2})] - \\ - \frac{1}{2} \frac{σ_{2 D}}{c ε_{0}} n_{s} \cdot ({{}^{+}H}_{1} h_{1} + {{}^{+}H}_{2} h_{2}) = D_{2}^{ext} - D_{1}^{ext} + \frac{1}{2} \frac{i σ_{2 D}}{ω ε_{0}} \nabla \cdot {(E_{1}^{ext} + E_{2}^{ext})}_{‖}, \end{array}

n_{s} \cdot ({{}^{-}H}_{1} h_{1} - {{}^{-}H}_{2} h_{2}) + i k ({{}^{-}G}_{2} ε_{2} σ_{2} - {{}^{-}G}_{1} ε_{1} σ_{1}) = α_{2} - α_{1},

\begin{array}{l} (n_{s} {}^{-}{T_{1}^{ν}} - t_{s}^{ν} {{}^{-}H}_{1}) \cdot h_{1} - (n_{s} {}^{-}{T_{2}^{ν}} - t_{s}^{ν} {{}^{-}H}_{2}) \cdot h_{2} - \frac{1}{2} μ_{0} c σ_{2 D} ({}^{+}{T_{1}^{ν}} σ_{1} + {}^{+}{T_{2}^{ν}} σ_{2}) + \\ + \frac{1}{2} i ω μ_{0} σ_{2 D} t_{s}^{ν} \cdot ({{}^{+}G}_{1} h_{1} + {{}^{+}G}_{2} h_{2}) = β_{2}^{ν} - β_{1}^{ν} - \frac{1}{2} μ_{0} c σ_{2 D} t_{s}^{ν} \cdot (E_{1}^{ext} + E_{2}^{ext}) . \end{array}

5. Numerical procedure

The numerical procedure resides in discretization of the surface integrals by dividing the interfaces into N small elements of area S_n within which the auxiliary surface charges and currents exhibit only small variations and can be considered to be nearly constant. As a consequence, the auxiliary quantites are represented by complex vectors (N × 1 for σ and N × 3 for h and η), while the Green’s function and its derivatives are approximated by N × N matrices so that the integrations over surface become true matrix-vector products. Importantly, both the Green’s function and its derivatives diverge in the s′ → s limit and generally can exhibit substantial variations at small distances. This can be resolved by averaging those quantities over the small interface elements S_n such as

{{}^{\pm}G}_{j, m n} = \int_{S_{n}} d s^{'} {{}^{\pm}G}_{j} (s_{m}, s^{'}),

where m, n serve as the matrix indices and s_m is the position of the m-th interface element.

The time complexity of the BEM algorithm is largely derived from the number of matrix multiplications and inversions necessary for the calculation of auxiliary charges and currents. Since the system of equations in the proposed modification with an additional surface current density η (η-BEM) is larger by three and the equations themselves are fully coupled to each other, the number of required algebraic operations is rather high compared to the original BEM algorithm. On the other hand, the equations in the modification with revised boundary conditions (BC-BEM) are coupled only partially and the time complexity is close to that of the traditional BEM. Table 1 then provides quantitative comparison between the three approaches both in 2D and 3D.

Table 1. Comparison of the three different approaches for modelling of 2D materials. Based on our implementation of BEM.

View Table

Apparently, the proposed modifications will be of benefit when the number of discretization elements needed to achieve convergence is reduced accordingly. Assuming O (N³) complexity for both the matrix inversion and multiplication, the rough estimate of the required reduction factor for the η-BEM is 2.5. Another benefit offered by the proposed modifications is the possibility to assign to each element a different dynamical conductivity, thus allowing for simulations including inhomogeneously doped 2D materials without the need of dividing those materials into many separate small segments.

6. Benchmark of the methods

In order to demonstrate the superiority of the presented modifications to the original version of BEM, we compared their performance on several problems with known analytical solutions, namely a self-standing graphene ribbon, and an infinite dielectric rod and a dielectric sphere, both coated with a very thin metallic film.

Figure 1 shows a resonance peak in the extinction cross-section spectra of a self-standing graphene ribbon calculated by three different methods. In case of the ordinary BEM (red line), the ribbon was modelled as a t = 1 nm thick film with an effective dielectric function

ε_{3 D} = 1 + \frac{i σ_{2 D}}{ε_{0} ω t},

whereas the η-BEM (green line) treated it as a conductive interface between two dielectric half-spaces. The numerically calculated spectra were tested against the semi-analytical solution (blue line) obtained by projecting Eq. (17) into the Fourier space, a method introduced by Koppens et al. [4].

Fig. 1 Comparison between extinction cross-section spectra of a self-standing graphene ribbon calculated by the ordinary BEM (red), η-BEM (green) and semi-analytical Fourier expansion method (blue). The ribbon has an infinite length, a width of 100 nm and in case of BEM a thickness of 1 nm (infinitesimal thickness is assumed for η-BEM). The doping level within the ribbon was set to E_F = 0.2 eV. The inset shows the shifts in the spectral position of the observed plasmonic peak as we increase the number of boundary elements. The 2D conductivity of graphene was modelled using Eq. 51 from Appendix A.

Download Full Size | PDF

The numerical solutions were considered to be converged once the relative change between two consecutive iterations was less than 1%. In this particular case, the number of discretization elements was doubled with each subsequent iteration. The inset of Fig. 1 shows the spectral position of the observed resonance peak as a function of the number of boundary elements N. Both the BEM-based approaches ultimately converged to values very close to the analytical solution (with a relative error of 0.6% for η-BEM and 1.3% for BEM), although the number of discretization elements needed by η-BEM to satisfy the convergence criterion (essentially the measure of how fast the method converges) is roughly 4 times smaller, a high enough factor to compensate for its increased time-complexity due to the large number of algebraic operations.

Next, we study the gold-coated dielectric particles (2D cylindrical rods and 3D spheres). The particles are composed of a dielectric core with the relative permittivity ε_Core = 2 and the diameter D_in = 72 nm, and a thin gold layer with the thickness t = 4 nm. The extinction cross-section spectra calculated by BEM, BC-BEM, and η-BEM are plotted in Fig. 2 and compared to the exact analytical solution obtained from the Mie theory [38]. Clearly, the use of BC-BEM and η-BEM implies that we need to convert the 3D dielectric function of gold into 2D surface conductivity, which is easily achieved by rearranging Eq. (49):

σ_{2 D} = i ε_{0} ω t (1 - ε_{3 D}) .

To assess the validity of treating a thin layer as a conductive interface, two analytical solutions were derived. One solution with the gold-shell represented by a thin film (Mie-TF), the other with the gold shell replaced by a single conductive interface of matching optical properties (Mie-CI). The analytical approaches coincide rather well, with the difference in the dipole peak wavelength of ≈ 7.5 nm for the coated rod and ≈ 0.5 nm for the coated sphere. This indicates that a thin 3D layer can be very well approximated by a conductive interface and its optical response can be therefore studied using the BC-BEM and η-BEM variants. Provided that in each subsequent simulation N is doubled (2D rods), resp. increased by 50% (3D spheres), the number of the boundary elements in particular simulations ensures the convergence of the peak energies with a relative error of 1%.

Fig. 2 (a) Comparison between extinction cross-section spectra of a gold-coated dielectric rod and sphere calculated by the ordinary BEM (blue), BC-BEM (green), η-BEM (red) and analytical Mie theory (black solid line for the thin gold film, black dashed line for the 2D conductive interface). The diameter of the dielectric core (ε = 2) was set to D_in = 72 nm and the thickness of the gold shell to t = 4 nm. The insets shows the spectral position of the dipolar plasmonic peak as a function of the number of discretization elements for the coated sphere (b) and cylinder (c). Note that the absence of data points for smaller N in the case of ordinary BEM in (b) is due to its failure to produce the dipolar plasmonic peak inside the simulated frequency range. This further underlines the poor convergence of the original BEM for this type of problems.

Download Full Size | PDF

In the case of the infinitely long rod, the ordinary BEM and η-BEM are equally capable to correctly calculate both the spectral position and strength of the dipolar plasmonic peak, using roughly the same number of discretization elements. Meanwhile, the BC-BEM seems to be numerically unstable and the solution tends to deteriorate with increasing N, likely due to the limited numerical accuracy of the term (n_s · ∇)² ⁺G_1,2 in Eq. (45). Turning our attention towards the coated sphere, the η-BEM gains the upper hand as both the traditional BEM and BC-BEM require substantially larger number of elements to achieve convergence. As one might expect, the convergence of the methods is quite sensitive to the thickness of the shells; generally, for thicker films the situation becomes more favourable for the ordinary BEM, whereas for thinner ones it is the other way around.

The proposed modifications to the current BEM algorithm proved to be of benefit in a number of problems involving 2D materials and thin coatings, since the considerable increase in time complexity is redeemed by a sufficient reduction in the number of discretization elements. Moreover, the η-BEM approach permits us to assign to each interface element a different dynamical conductivity, giving us a viable tool to simulate light scattering on a whole new class of structures including inhomogeneously doped graphene and corrugated thin films.

7. Inhomogeneous doping

The group of García de Abajo has recently shown that even for the simplest gating configuration, i.e. back-gated graphene structures such as ribbons [35] and discs [39], the distribution of induced charge can be strongly inhomogeneous. Interestingly, this has only a little effect on the charge distribution and spectral position of plasmonic modes.

Among other configurations, Abajo et al. also investigated plasmonic properties of side-gated graphene structures in the quasistatic limit. When an electrically isolated graphene ribbon is exposed to a strong static electric field E_SG with orientation parallel to the ribbon, its edges become oppositely charged to counteract the external field. This redistribution of charge within the ribbon is naturally accompanied by a shift in the Fermi level, highest at the edges and gradually decreasing as we move towards the center of the ribbon. Consequently, the central part, where the induced charge density changes the sign and E_F approaches zero, acts as a very narrow spacer between two optically conductive regions. In analogy to optical dimer antennas, such structure should support plasmonic oscillations with very large field enhancement in the “gap”.

Here, we study the side-gated graphene using our modified BEM, with emphasis on the role of interband transitions which were omitted in the previous work by the group of García de Abajo. The procedure is rather straightforward: Once the spatial dependence of the Fermi level is known (for the detailed calculation see Ref. [35]), we can insert it into the model for the optical conductivity of extended graphene set down in Eq. 51 in Appendix A and subsequently calculate the optical response of the side-gated ribbon using the η-BEM modification.

Figure 3 shows extinction spectra of a side-gated ribbon for several values of the static electric field E_SG and two different models of the optical conductivity - with and without the interband transitions. The dimensions of the ribbon and its orientation towards the static electric field and the incident electromagnetic wave are depicted in Figure 3(b)) and spatial dependence of the Fermi level for the three values of E_SG in Figure 3(c). Clearly, the built-up of charge at the edges of the ribbon is accompanied by a prominent drop in the Fermi level right at its center, virtually dividing the ribbon into two separate, optically conductive regions.

Fig. 3 (a) Extinction cross-section spectra of a side-gated self-standing ribbon for three intensities of the static electric field E_SG. The solid and dashed lines discriminate between two different models for dynamical conductivity, with the dashed one corresponding to a model with enabled interband transitions. Based on the maps of local electric field and distribution of induced charge density n_pl we identified the two marked peaks as “single antenna” and “dimer antenna” plasmonic resonances. (b) Schematical drawing of the simulated structure - a single self-standing graphene ribbon in a parallelly oriented static electric field illuminated from the top by infrared radiation. The ribbon is infinite in the direction perpendicular to the figure plane. (c) The spatial distribution of the absolute value of the Fermi level over the graphene ribbon for the three intensities of the static electric field E_SG. Note that in the region 0 < x < 50 nm, there is depletion of electrons and the Fermi level is therefore negative.

Download Full Size | PDF

The differences between the two sets of spectra (solid and dashed) demonstrate the importance of the interband transitions. While in the case of homogeneously doped or even back-gated ribbon the neglection of interband processes is justifiable, for the strongly inhomogeneous doping profile of the side-gated ribbon the situation dramatically changes, as the shift in the Fermi energy is not universally high enough to ensure Pauli blocking [40] in the entire ribbon. Even though this effect is largely restricted to the central part of the ribbon, where the Fermi level is close to zero, the extinction spectra plotted in Fig. 3 indicate that the interband transitions are indeed a strong loss channel. Namely, the spectra calculated using the conductivity model without interband transitions (solid lines) display prominent plasmonic peaks for electric field intensities E_SG = 30 MV· m⁻¹ and higher, whereas the spectra with interband transitions turned on (dashed lines) remain flat regardless of the strength of the applied electric field.

To verify the above findings and gain more insight into the nature of the observed plasmonic resonances, we further calculated the plasmon charge distributions and the maps of local electric field intensity at the energies corresponding to those resonances. Taking into account the fact that in resonance the motion of electrons should be delayed by 90 degrees with respect to the driving electromagnetic field, we associate the spatial profile of a particular plasmon mode with the imaginary part (green) of the plasmon charge distribution n_pl, while the real part (blue) contains the off-resonance contributions from other plasmon modes. Figure 4 shows the charge distributions of the first four plasmon modes excited in a self-standing graphene ribbon subjected to a static electric field of strength E_SG = 30 MV· m⁻¹ (green in Fig. 3). The mode with the lowest energy ħω = 0.059 eV clearly corresponds to a simple dipolar antenna. In the second one (ħω = 0.113 eV) we recognize the predicted ”dimer antenna” mode with a substantial amount of charge piled up in the middle of the ribbon. The remaining two resonances at ħω = 0.155 eV and ħω = 0.183 eV then belong to higher order modes. For completeness we have also included the spatial profiles that are obtained when the conductivity model with enabled interband transitions (dashed lines) is used. Apparently, the newly introduced decay mechanism virtually quenches all plasmonic oscillations in the parts with a sufficiently low Fermi level and as we move towards the higher energies, larger and larger portion of the ribbon becomes affected. This further deepens the impression given by the absence of plasmonic resonances in the extinction spectra, that the side-gating configuration is not well-suited for applications where sharp resonance features and high field-enhancement factors are demanded. The suppression of the “dimer antenna” plasmonic mode when the interband transitions are turned on is also clearly visible on the maps of the local electric field presented in Fig. 5. While the field distribution for the lowest mode with the highest intensity at the ribbon edges remains mostly unaffected, there is a clear difference between those two conductivity models as we turn our attention towards the maps taken at the photon energy ħω = 0.113 eV corresponding to the resonant excitation of the “dimer antenna” plasmon. Due to the quenching of plasmonic oscillations, the hot spot at the center of the ribbon completely vanishes.

Fig. 4 The spatial profile of the lowest four plasmon modes observed in Fig. 3 for E_SG = 30 MV · m⁻¹. The imaginary part of the plasmon charge density n_pl corresponds to the mode in resonance, while the real part represents the off-resonant contributions from the rest of the modes. The type of the line reflects the conductivity model used in the calculation.

Download Full Size | PDF

Fig. 5 Maps of the local electric field intensity for the two plasmon modes with the lowest energy from Fig. 4. The calculations with interband transitions turned off (on) are on the left (right). The hot spot associated with the build-up of the charge in the “gap” separating the two optically conductive regions completely vanishes when the decay via interband transitions is enabled.

Download Full Size | PDF

8. Conclusion

In this work we have addressed the issue of modeling graphene and other 2D materials in light scattering simulations using the boundary element method. We introduced a modification of the traditional BEM that allows for a substantial reduction in the number of required discretization elements. As another important advantage of the presented extension to the current BEM formulation we see the possibility to assign to each discretization element a different dynamical conductivity, thus enabling, for example, simulations of inhomogeneously doped graphene structures. In order to demonstrate this ability we carried out a study focusing on a side-gated graphene ribbon which was predicted to support highly localized plasmon oscillations. Our calculations show that these oscillations are strongly damped via interband processes and the expected tunable hot-spot is therefore entirely missing. Nonetheless, our modified boundary element method represents a valuable tool that facilitates more efficient modeling of 2D materials and opens a door to the study of systems and devices in which the spatial tailoring of optical properties can lead to new interesting effects and functional properties.

Appendix A: 2D conductivity of graphene

In its natural state, graphene behaves as a zero-gap semiconductor, with the Fermi level positioned directly at the crossing between the valence and conduction band, i.e. the Dirac point [41]. Via doping or interaction with a substrate, the Fermi level can be shifted and graphene acquires metallic character. This naturally also affects its optical properties. In terms of the Fermi energy E_F (measured from the Dirac point), the energy of the absorption edge reads 2 |E_F|.

The optical properties of extended graphene are described by its non-local dynamical conductivity σ_2D(q, ω), with q and ω denoting spatial and temporal frequency, but Thongrattanasiri et al. [42] have recently shown that the majority of graphene structures can be adequately treated within the local limit q → 0. In this work we use for σ_2D(ω) the following frequently encountered formula derived within random phase approximation [43,44]

σ (ω) = \frac{e^{2} | E_{F} |}{π ℏ^{2}} \frac{i}{(ω + i τ^{- 1})} + \frac{e^{2}}{4 ℏ} [θ (ℏ ω - 2 | E_{F} |) + \frac{i}{π} ln | \frac{ℏ ω - 2 | E_{F} |}{ℏ ω + 2 | E_{F} |} |],

where τ stands for the relaxation time of intraband scattering processes and θ is the Heaviside step function. The i/ (ω + iτ⁻¹) dependence indicates that at low frequencies graphene behaves as a Drude metal, while for photon energies larger than 2 |E_F|, the second term corresponding to the onset of interband transitions becomes dominant.

Funding

Grant Agency of the Czech Republic (grant No. 15-21581S), MEYS CR (project No. LQ1601 CEITEC 2020), and Technology Agency of the Czech Republic (grant No. TE01020233).

References and links

1. K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, “Electric field effect in atomically thin carbon films,” Science 306, 666–669 (2004). [CrossRef] [PubMed]

2. T. Low and P. Avouris, “Graphene plasmonics for terahertz to mid-infrared applications,” ACS Nano 8, 1086–1101 (2014). [CrossRef] [PubMed]

3. A. C. Ferrari, F. Bonaccorso, V. Fal’ko, K. S. Novoselov, S. Roche, P. Boggild, S. Borini, F. H. L. Koppens, V. Palermo, N. Pugno, J. A. Garrido, R. Sordan, A. Bianco, L. Ballerini, M. Prato, E. Lidorikis, J. Kivioja, C. Marinelli, T. Ryhanen, A. Morpurgo, J. N. Coleman, V. Nicolosi, L. Colombo, A. Fert, M. Garcia-Hernandez, A. Bachtold, G. F. Schneider, F. Guinea, C. Dekker, M. Barbone, Z. Sun, C. Galiotis, A. N. Grigorenko, G. Konstantatos, A. Kis, M. Katsnelson, L. Vandersypen, A. Loiseau, V. Morandi, D. Neumaier, E. Treossi, V. Pellegrini, M. Polini, A. Tredicucci, G. M. Williams, B. Hee Hong, J.-H. Ahn, J. Min Kim, H. Zirath, B. J. van Wees, H. van der Zant, L. Occhipinti, A. Di Matteo, I. A. Kinloch, T. Seyller, E. Quesnel, X. Feng, K. Teo, N. Rupesinghe, P. Hakonen, S. R. T. Neil, Q. Tannock, T. Lofwander, and J. Kinaret, “Science and technology roadmap for graphene, related two-dimensional crystals, and hybrid systems,” Nanoscale 7, 4598–4810 (2015). [CrossRef] [PubMed]

4. F. H. L. Koppens, D. E. Chang, and F. J. García de Abajo, “Graphene plasmonics: A platform for strong light-matter interactions,” Nano Letters 11, 3370–3377 (2011). [CrossRef] [PubMed]

5. S. Thongrattanasiri and F. J. García de Abajo, “Optical field enhancement by strong plasmon interaction in graphene nanostructures,” Physical Review Letters 110, 187401 (2013). [CrossRef] [PubMed]

6. M. Freitag, T. Low, W. Zhu, H. Yan, F. Xia, and P. Avouris, “Photocurrent in graphene harnessed by tunable intrinsic plasmons,” Nature Communications 42951 (2013). [CrossRef] [PubMed]

7. D. Rodrigo, O. Limaj, D. Janner, D. Etezadi, F. J. G. de Abajo, V. Pruneri, and H. Altug, “Mid-infrared plasmonic biosensing with graphene,” Science 349, 165–168 (2015). [CrossRef] [PubMed]

8. T. Wu and L. Wei, “Tunable resonant graphene plasmons for mid-infrared biosensing,” Optics Express 24, 26241–26248 (2016). [CrossRef] [PubMed]

9. H. Hu, X. Yang, F. Zhai, D. Hu, R. Liu, K. Liu, Z. Sun, and Q. Dai, “Far-field nanoscale infrared spectroscopy of vibrational fingerprints of molecules with graphene plasmons,” Nature Communications 7, 12334 (2016). [CrossRef] [PubMed]

10. E. Carrasco, M. Tamagnone, J. R. Mosig, T. Low, and J. Perruisseau-Carrier, “Gate-controlled mid-infrared light bending with aperiodic graphene nanoribbons array,” Nanotechnology 26, 134002 (2015). [CrossRef] [PubMed]

11. F. J. García de Abajo, “Graphene plasmonics: Challenges and opportunities,” ACS Photonics 1, 135–152 (2014). [CrossRef]

12. J. Qi, H. Liu, and X. C. Xie, “Surface plasmon polaritons in topological insulators,” Physical Review B 89155420 (2014). [CrossRef]

13. T. Low, R. Roldán, H. Wang, F. Xia, P. Avouris, L. M. Moreno, and F. Guinea, “Plasmons and screening in monolayer and multilayer black phosphorus,” Physical Review Letters 113106802 (2014). [CrossRef]

14. X. Wang and S. Lan, “Optical properties of black phosphorus,” Advances in Optics and Photonics 8, 618–655 (2016). [CrossRef]

15. J. Wang and Y. Jiang, “Infrared absorber based on sandwiched two-dimensional black phosphorus metamaterials,” Optics Express 25, 5206–5216 (2017). [CrossRef] [PubMed]

16. T. Low, A. Chaves, J. D. Caldwell, A. Kumar, N. X. Fang, P. Avouris, T. F. Heinz, F. Guinea, L. Martin-Moreno, and F. Koppens, “Polaritons in layered two-dimensional materials,” Nature Materials 10384792 (2016). [CrossRef]

17. J. B. Khurgin, “Two-dimensional exciton-polariton-light guiding by transition metal dichalcogenide monolayers,” Optica 2, 740–742 (2015). [CrossRef]

18. A. F. Oskooi, D. Roundy, M. Ibanescu, P. Bermel, J. Joannopoulos, and S. G. Johnson, “Meep: A flexible free-software package for electromagnetic simulations by the FDTD method,” Computer Physics Communications 181, 687–702 (2010). [CrossRef]

19. J. Niegemann, M. Knig, K. Stannigel, and K. Busch, “Higher-order time-domain methods for the analysis of nano-photonic systems,” Photonics and Nanostructures - Fundamentals and Applications 7, 2–11 (2009). [CrossRef]

20. V. W. Brar, M. S. Jang, M. Sherrott, J. J. Lopez, and H. A. Atwater, “Highly confined tunable mid-infrared plasmonics in graphene nanoresonators,” Nano Letters 13, 2541–2547 (2013). [CrossRef] [PubMed]

21. J. Yang, J. Yang, W. Deng, F. Mao, and M. Huang, “Transmission properties and molecular sensing application of cgpw,” Optics Express 23, 32289–32299 (2015). [CrossRef] [PubMed]

22. V. Myroshnychenko, E. Carbó-Argibay, I. Pastoriza-Santos, J. Pérez-Juste, L. M. Liz-Marzán, and F. J. G. de Abajo, “Modeling the optical response of highly faceted metal nanoparticles with a fully 3D boundary element method,” Advanced Materials 20, 4288–4293 (2008). [CrossRef]

23. U. Hohenester and A. Trügler, “MNPBEM – a matlab toolbox for the simulation of plasmonic nanoparticles,” Computer Physics Communications 183, 370–381 (2012). [CrossRef]

24. U. Hohenester and J. Krenn, “Surface plasmon resonances of single and coupled metallic nanoparticles: A boundary integral method approach,” Physical Review B 72195429 (2005). [CrossRef]

25. J. Waxenegger, A. Trügler, and U. Hohenester, “Plasmonics simulations with the MNPBEM toolbox: Consideration of substrates and layer structures,” Computer Physics Communications 193, 138–150 (2015). [CrossRef]

26. V. Myroshnychenko, J. Rodriguez-Fernandez, I. Pastoriza-Santos, A. M. Funston, C. Novo, P. Mulvaney, L. M. Liz-Marzan, and F. J. Garcia de Abajo, “Modelling the optical response of gold nanoparticles,” Chemical Society Reviews 37, 1792–1805 (2008). [CrossRef] [PubMed]

27. B. Gallinet, J. Butet, and O. J. F. Martin, “Numerical methods for nanophotonics: standard problems and future challenges,” Laser & Photonics Reviews 9, 577–603 (2015). [CrossRef]

28. J. Mertens, A. L. Eiden, D. O. Sigle, F. Huang, A. Lombardo, Z. Sun, R. S. Sundaram, A. Colli, C. Tserkezis, J. Aizpurua, S. Milana, A. C. Ferrari, and J. J. Baumberg, “Controlling subnanometer gaps in plasmonic dimers using graphene,” Nano Letters 13, 5033–5038 (2013). [CrossRef] [PubMed]

29. Y. Yao, M. A. Kats, P. Genevet, N. Yu, Y. Song, J. Kong, and F. Capasso, “Broad electrical tuning of graphene-loaded plasmonic antennas,” Nano Letters 13, 1257–1264 (2013). [CrossRef] [PubMed]

30. Y. Yao, M. A. Kats, R. Shankar, Y. Song, J. Kong, M. Loncar, and F. Capasso, “Wide wavelength tuning of optical antennas on graphene with nanosecond response time,” Nano Letters 14, 214–219 (2014). [CrossRef]

31. V. Nayyeri, M. Soleimani, and O. M. Ramahi, “Modeling graphene in the finite-difference time-domain method using a surface boundary condition,” IEEE Transactions on Antennas and Propagation 61, 4176–4182 (2013). [CrossRef]

32. P. Li, L. J. Jiang, and H. Bagci, “A resistive boundary condition enhanced DGTD scheme for the transient analysis of graphene,” IEEE Transactions on Antennas and Propagation 63, 3065–3076 (2015). [CrossRef]

33. Z. Li, K. Yao, F. Xia, S. Shen, J. Tian, and Y. Liu, “Graphene plasmonic metasurfaces to steer infrared light,” Scientific Reports 5, 12423 (2015). [CrossRef] [PubMed]

34. L. Ye, Y. Chen, G. Cai, N. Liu, J. Zhu, Z. Song, and Q. H. Liu, “Broadband absorber with periodically sinusoidally-patterned graphene layer in terahertz range,” Optics Express 25, 11223–11232 (2017). [CrossRef]

35. S. Thongrattanasiri, I. Silveiro, and F. J. García de Abajo, “Plasmons in electrostatically doped graphene,” Applied Physics Letters 100, 201105 (2012). [CrossRef]

36. F. J. García de Abajo and A. Howie, “Retarded field calculation of electron energy loss in inhomogeneous dielectrics,” Physical Review B 65, 115418 (2002). [CrossRef]

37. L. Novotny and B. Hecht, Principles of Nano-Optics (Cambridge University, 2006). [CrossRef]

38. C. F. Bohren and D. R. Huffman, Absorption and Scattering of Light by Small Particles (Wiley-VCH, 1998). [CrossRef]

39. I. Silveiro and F. J. García de Abajo, “Plasmons in inhomogeneously doped neutral and charged graphene nanodisks,” Applied Physics Letters 104, 131103 (2014). [CrossRef]

40. Q. Bao, H. Zhang, Y. Wang, Z. Ni, Y. Yan, Z. X. Shen, K. P. Loh, and D. Y. Tang, “Atomic-layer graphene as a saturable absorber for ultrafast pulsed lasers,” Advanced Functional Materials 19, 3077–3083 (2009). [CrossRef]

41. T. Ando, “Theory of electronic states and transport in carbon nanotubes,” Journal of the Physical Society of Japan 74, 777–817 (2005). [CrossRef]

42. S. Thongrattanasiri, A. Manjavacas, and F. J. García de Abajo, “Quantum finite-size effects in graphene plasmons,” ACS Nano 6, 1766–1775 (2012). [CrossRef] [PubMed]

43. B. Wunsch, T. Stauber, F. Sols, and F. Guinea, “Dynamical polarization of graphene at finite doping,” New Journal of Physics 8, 318 (2006). [CrossRef]

44. E. H. Hwang and S. Das Sarma, “Dielectric function, screening, and plasmons in two-dimensional graphene,” Physical Review B 75, 205418 (2007). [CrossRef]

Boundary element method for 2D materials and thin films

Abstract

1. Introduction

2. Principles of BEM

3. η-BEM: Additional surface current density approach

4. BC-BEM: 2D conductivity incorporated into the boundary conditions

5. Numerical procedure

6. Benchmark of the methods

7. Inhomogeneous doping

8. Conclusion

Appendix A: 2D conductivity of graphene

Funding

References and links

Cited By

Figures (5)

Tables (1)

Equations (54)

Optics Express

2D	Multiplications	Inversions	Total
BEM	9	12	21
BC-BEM	10	22	32
η-BEM	95	24	119

3D	Multiplications	Inversions	Total
BEM	11	14	25
BC-BEM	22	30	52
η-BEM	246	48	292