1. Introduction

Numerical methods for the computation of electromagnetic (EM) fields, such as the finite-difference time-domain (FDTD) method [1], finite-element method (FEM) [2], and boundary element method [3], have been developed to analyze the performance and determine the optimal design of various optical devices. The quality factor of a photonic crystal slab cavity, for example, can be significantly improved by designing the cavity structure based on the results of FDTD simulations [4]. Furthermore, these EM-field simulations are of immense benefit in the fascinating research fields of plasmonics and metamaterials.

The interesting optical properties that are found in plasmonics originate from highly intense and strongly localized EM fields mediated by a surface plasmon, and similar optical properties would be expected from exciton resonance. However, established numerical methods are limited to dielectrics and metals, and an EM simulation that takes into account an exciton confined in arbitrary geometry has yet to be developed in a general form. One exception is the extended version of the discrete dipole approximation (DDA) [5], which has been used to study the optical selection-rule breakdown [6] and finite exciton momentum [7] in a carbon nanotube placed in a nanogap between gold nanoparticles. The extended DDA requires advance knowledge of a large number of exciton wavefunctions that account for the electron–hole exchange interaction [8]. In this paper, we present an FEM-based EM simulation method that can be applied to exciton resonance in an arbitrary confinement geometry without needing to calculate exciton wavefunctions.

Excitons are bound states of an optically excited electron–hole pair in a semiconductor and play an important role in the optical properties of semiconductor nanostructures. In fact, the spectral peaks of a nanostructure depend on its size and geometry because of the quantum confinement effect, which arises from the strong dispersion relation between the exciton energy and wave vector.

In contrast to an exciton, the dispersion relation of a plasmon in a metal is negligibly small, i.e., the plasmon energy is almost independent of its wave vector. As a result, there is at most only one plasmon polariton mode in a metal at a given frequency. This means that EM simulation methods can be directly applied to metallic nanostructures by simply using a frequency-dependent dielectric function. For an exciton, however, which has a strong dispersion relation, there can be at most two transverse exciton polariton modes propagating in a semiconductor at a given frequency. In addition, there is one longitudinal exciton mode. The multimode excitations in semiconductor nanostructures mean that the number of boundary conditions for the associated Maxwell’s equations is less than the number of variables (the field amplitudes). To resolve this problem, two methods have been developed: the introduction of an additional boundary condition (ABC) for the exciton polarization [9–14], and a microscopic nonlocal theory [15,16] in which the exciton polarization and EM field are determined self-consistently. We adopt the ABC in our FEM-based EM simulations.

The EM simulation method proposed here for a confined exciton has been developed for one-and two-dimensional (1D and 2D) systems with the ABC [17]. The method can be applied to the limiting case in which the longitudinal component of the exciton is not excited by incident light polarized in a direction parallel to the semiconductor surface (S polarization), and cannot be applied for P polarized incident light even for 1D and 2D systems. For a three-dimensional arbitrary confinement geometry, the longitudinal component of the exciton arises in general; for the exciton in a spherical quantum dot (QD), for example, the longitudinal and transverse mixed mode is strongly coupled to a plane-wave EM field when the size of the sphere is much smaller than the wavelength of the EM field [18]. The present FEM-based EM simulation method can also be applied to the general case in which the longitudinal component of the exciton is excited.

2. Variational formulation

Maxwell’s equations lead to vector wave equations (Maxwell’s wave equations) for the electric and magnetic fields and the boundary conditions at the interface of materials with different permittivities or permeabilities. This boundary-value problem can be formulated using a functional ℱ, and the stationary condition of the functional δℱ = 0 provides the governing Maxwell’s wave equations and Neumann boundary conditions.

2.1. Electromagnetic field modes

We consider a situation in which an incident field E^(inc) irradiates a semiconductor nanostructure (region Ω₁) surrounded by a dielectric material with a dielectric constant of ε, as illustrated schematically in Fig. 1 (a). The surface of region Ω₁ is denoted by Γ₁. In an open-region scattering problem, the scattered field propagates to infinity, but to restrict the computational region, we introduce an artificial interface Γ₀ within which the EM fields are numerically calculated by subdividing the region into small volume elements. The dielectric area between Γ₀ and Γ₁, denoted by Ω₀, and the exterior (far-field) region, denoted by Ω_∞, have dielectric constants of ε, and thus there is no reflection or refraction at Γ₀.

 

Fig. 1 (a) Schematic illustration of a semiconductor nanostructure (region Ω₁) surrounded by a dielectric (regions Ω₀ and Ω_∞). The interface between the semiconductor nanostructure and dielectric is denoted by Γ₁, and an artificial spherical surface, denoted by Γ₀, is introduced within which the region is subdivided into small volume elements for the numerical computation. Various wave modes are also depicted. (b) The dispersion relations of two transverse exciton polaritons and one longitudinal exciton in bulk CuCl (the CuCl parameters for the calculation are summarized in Table 1).

Download Full Size | PDF

In the semiconductor region Ω₁, the dielectric function including the exciton effect is expressed as

The total electric field E⁽¹⁾ in region Ω₁ is obtained as

We consider a plane-wave incident field given by E^(inc)(r) = I₀ exp(ikz) with I₀ being polarized in the x-direction. The incident field can be expanded as

2.2. Boundary conditions

At the interface Γ₁, the following Maxwell’s boundary conditions (MBCs) are satisfied:

At the artificial interface Γ₀, the MBCs are given by

2.3. Functionals

The electric field E⁽⁰⁾ in region Ω₀ obeys the following Maxwell’s wave equation:

In region Ω₁,the electric-field components E^(t₁) and E^(t₂) of the two exciton polaritons are governed by the Maxwell’s wave equation

3. Finite element analysis

To numerically calculate the solutions of the boundary-value problem, the computational region is subdivided into small volume elements as shown in Fig. 2. As we choose a tetrahedral geometry for the volume elements, the interface is subdivided into triangular elements. All the volume elements are labeled by a set of integers e = 1, 2,···, and the nodes of each volume element are labeled as i = 1, 2, 3, 4. The ith node in the eth volume element is denoted by a local label (e;i). Any node can be indicated by multiple local labels because the node belongs to multiple volume elements. For one-to-one mapping between the nodes and integers, we create a global label s(e;i) = 1, 2,··· for all nodes. The triangular element at the interface belonging to the eth volume element is denoted by f_e.

 

Fig. 2 Example of the subdivision of the mesh by tetrahedral elements.

Download Full Size | PDF

In the conventional FEM formalism, scalar or vector fields are assigned at the nodes of the volume elements, referred to as nodal elements. For transverse vector fields, however, a higher accuracy can be obtained if the fields are assigned at the edges of the volume elements (edge elements) [22, 23]. Here, we develop a hybrid method, i.e., the scalar field Φ is represented by nodal elements and the transverse vector fields E^(α) with α = {0, t₁, t₂} are represented by edge elements. In more detail, Φ is expanded by the nodal-element basis L_(e;i)(r):

Maxwell’s wave equations can be derived from the stationary condition of the functional ℱ, whose explicit expression is given in the Appendix B.2, with respect to $E_{s(e;i)s(e;j)}^{(\alpha )}$ and Φ_s(e;i). To include the MBC [Eq. (10)], we find a stationary condition

The other stationary conditions are given by

The MBC [Eq. (10)], ABC [Eq. (12)], and the MBC [Eq. (15)] are included in ℱ̃ (see the Appendix B.2). The MBC [Eq. (14)] is automatically satisfied for edge elements because the tangential electric field at the interface expresses common fields that belong to the different regions of the interface. However, this is not the case for the MBC [Eq. (9)] at the interface Γ₁ because the electric field in a semiconductor area consists of E^(t₁), E^(t₂), and (−∇Φ). Therefore, we set the three variables $E_{s(e;i)s(e;j)}^{(t_1)}$, $E_{s(e;i)s(e;j)}^{(t_2)}$, and $E_{s(e;i)s(e;j)}^{(0)}$ at each edge element on the interface Γ₁ and explicitly impose the MBC [Eq. (9)] as follows:

Let us check that the number of variables is equal to the number of equations. The numbers of equations defined by Eqs. (30) and (31) are equal to the numbers of variables a_mn and b_mn, respectively, and the numbers of equations defined by Eqs. (26) and (27) are equals to the numbers of variables $E_{s(e;i)s(e;j)}^{(\alpha )}$ [s(e; i), s(e; j) ∉ Γ₁] and Φ_s(e;i), respectively. The remaining variables are $E_{s(e;i)s(e;j)}^{(\alpha )}$ [s(e; i), s(e; j) ∈ Γ₁], and the number of these variables is 3M, where M is the number of edges on Γ₁. The remaining equations are explicit boundary conditions [Eqs. (28) and (29)] and the stationary condition [Eq. (25)], each of which provides M equations, and thus there are a total of 3M equations. Therefore, $E_{s(e;i)s(e;j)}^{(\alpha )}$ and Φ_s(n;i) can be uniquely determined.

4. Numerical results

We apply the present FEM to light scattering problems in semiconductor nanostructures. The scattering cross section σ_s is calculated using the expansion coefficients a_mn and b_mn of ${\mathbf{E}}_{\infty }^{(\text{scat})}(\mathbf{r})$ as follows:

To begin, we calculate the scattering cross section of a spherical QD. The scattering cross section for this geometry has also been calculated using Mie theory, in which the EM fields inside and outside the sphere are expanded by vector spherical harmonics and the expansion coefficients are determined from the MBCs and Pekar-type ABC [24]. The validity of the present FEM is confirmed by a comparison of the cross section with that calculated by Mie theory. We then study the scattering cross section and field distributions of a hexagonal-disk QD, which cannot be calculated using established methods.

4.1. Spherical QD

We consider a spherical QD of CuCl with a 10-nm diameter. Spherical QDs can be synthesized in a glass matrix [25], and recently, spherical CuCl QDs have been fabricated from a bulk sample using laser ablation in superfluid helium; size selective transportation of the QDs has also been reported [26]. The lowest exciton of CuCl, known as the Z₃ exciton, has a simple electronic structure and a Bohr radius of 0.7 nm. Since the exciton Bohr radius is sufficiently smaller than the QD diameter, we can use the bulk exciton parameters. We take the Z₃ exciton into account in the following calculations using the parameters listed in Table 1.

Table 1. Parameters of the Z₃ exciton of CuCl and Γ₅(B) exciton of ZnO; m_e denotes the free electron mass.

View Table

Figure 3 (a) shows the volume element mesh in the QD region. A total of 6609 nodes and 35429 elements were used in the calculation. A sufficient cut-off angular momentum is N_c = 3. Figure 3 (b) shows the calculated scattering cross sections for an incident plane wave polarized in the x-direction using the present FEM (crosses) and Mie theory (solid line); the exciton damping energy was set to h̄γ = 1 meV. The FEM results exhibit excellent agreement with those of Mie theory, which confirms the validity of the present FEM. Figures 3 (c) and (d) show the distributions of the electric field and exciton polarization under the resonant photon energy of 3.206 eV, respectively. A strong enhanced field appears around the edge of the QD, as shown in Fig. 3 (c). This enhanced-field distribution appears in metallic nanostructures because of the surface-plasmon resonant scattering, whereas in semiconductor nanostructures, it originates from the exciton resonance. The profile of the polarization field in Fig. 3 (d) indicates that the excited state is a 1s-like nodeless exciton.

 

Fig. 3 (a) Volume element mesh for a spherical QD. (b) Scattering cross section of a spherical QD of CuCl with a 10-nm diameter calculated by the present FEM (crosses) and by Mie theory (solid line). (c) Vector plots of the calculated electric field at the resonant photon energy indicated by the arrow in (b). (d) Vector plots of the calculated exciton polarization at the resonant photon energy indicated by the arrow in (b).

Download Full Size | PDF

4.2. Hexagonal-disk QD

To demonstrate the advantages of the present FEM, we consider a hexagonal-disk QD. In contrast to the spherical geometry, there are no suitable vector basis functions. Furthermore, other established numerical techniques cannot include the exciton effect for nanostructures with an arbitrary geometry. However, we show here that the present FEM is capable of calculating the light scattering from a hexagonal-disk QD of ZnO. The exciton of ZnO has a large binding energy of 60 meV, and by virtue of this fact, room-temperature laser emission from excitons in ZnO microcrystallite thin films has been realized; the hexagonal ZnO QDs were grown by laser molecular beam epitaxy [27].

Figure 4 (a) shows the volume element mesh in the QD region, and in this case, 12947 nodes, 72169 elements, and N_c = 3 were used. Figure 4 (b) shows the calculated scattering cross section for an incident plane wave propagating in the z-direction and polarized in the x-direction; the exciton damping energy was set to h̄γ = 1 meV. The electric field distributions at the resonant photon energies, indicated by the left and right arrows, are mapped in Fig. 4 (c) and Fig. 4 (f), respectively. The magnitudes of the field have been normalized to the magnitude of the incident field. Both scattered fields are enhanced by a factor of 10 at the edges of the hexagonal-disk QD near the x-axis. Although these field distributions are similar, the exciton polarization patterns are quite different, as can be seen in Fig. 4 (d) and (g) (x-component) and Fig. 4 (e) and (h) (y-component). These plots indicate that 1s-like and 2p-like excitons are excited at resonant photon energies of 3.39063 eV and 3.39227 eV, respectively. Thus, the present FEM is capable of resolving the excited modes of the excitons.

 

Fig. 4 (a) Volume element mesh for a hexagonal-disk QD. (b) Scattering cross section calculated by the present FEM of a hexagonal-disk QD of ZnO with a length of 30 nm and a height of 10 nm. At the resonant photon energy of 3.39063 eV indicated by the left arrow in (b), (c) the magnitudes of the electric fields and the (d) x- and (e) y-components of the exciton polarization are plotted in the x–y plane of the top surface of the QD. At the resonant photon energy of 3.39227 eV indicated by the right arrow in (b), (f) the magnitudes of the electric fields and the (g) x- and (h) y-components of the exciton polarization are plotted in the same plane.

Download Full Size | PDF

5. Conclusion

We have developed an FEM-based EM simulation method for semiconductor nanostructures with arbitrary geometries. The EM field in a semiconductor propagates as an exciton polariton with a transverse character. At a given frequency, there are at most two exciton-polariton and one longitudinal-exciton modes in the semiconductor, and thus, an ABC is necessary in addition to the MBCs for determining the EM fields. To represent the weak-form functional in terms of variables defined in subdivided volume elements, we have developed a hybrid edge–node-element method to improve the accuracy. The transverse fields and the scalar field providing the longitudinal field are represented by edge and node elements, respectively. We have found the stationary condition [Eq. (25)] from which the MBC [Eq. (10)] can be added to the derived linear equations. EM fields have been obtained by solving the equations derived from the stationary conditions [Eqs. (25)–(27)] for ℱ̃ and the equations of the explicit boundary conditions [Eqs. (28)–(31)].

The present method has been confirmed to be valid through a comparison with the scattering cross section of a spherical semiconductor calculated by Mie theory. Furthermore, we have applied the present FEM to the calculation of the electric-field and exciton-polarization distributions of a hexagonal-disk QD of ZnO. At each resonant photon energy, the characteristic polarization patterns reflecting the size-quantized exciton are well resolved by the calculations. We expect that the developed FEM will be a powerful tool for designing and analyzing semiconductor optical devices utilizing exciton resonance.

A. Vector spherical harmonics

An EM field in region Ω_∞ can be expanded by the basis vector functions for an outgoing wave satisfying the Maxwell’s wave equations. Since the artificial boundary Γ₀ is a sphere, it is useful to represent the vector functions in spherical polar coordinates. The incident plane wave can also be expanded by the basis vector functions for the outgoing wave. These vector functions are given by

(33)$${\mathbf{M}}_{mn}(\mathbf{r})=\nabla \times [\mathbf{r}{\mathrm{\Psi }}_{mn}(\mathbf{r})],$$

(34)$${\mathbf{N}}_{mn}(\mathbf{r})=\frac{1}{k}\nabla \times \nabla \times [\mathbf{r}{\mathrm{\Psi }}_{mn}(\mathbf{r})],$$

with

(35)$${\mathrm{\Psi }}_{mn}(\mathbf{r})=h_n^{(1)}(kr)P_n^m(\text{cos}\theta )e^{im\varphi },$$

where $h_n^{(1)}$ is the nth spherical Hankel function of the first kind, and $P_n^m$ is the associated Legendre function of order n and degree m. These vector functions have a transverse character of

(36)$$\nabla \cdot {\mathbf{M}}_{mn}(\mathbf{r})=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\nabla \cdot {\mathbf{N}}_{mn}(\mathbf{r})=0,$$

and relations of

(37)$$\nabla \times {\mathbf{N}}_{mn}(\mathbf{r})=k{\mathbf{M}}_{mn}(\mathbf{r}),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\nabla \times {\mathbf{M}}_{mn}(\mathbf{r})=k{\mathbf{N}}_{mn}(\mathbf{r}).$$

Explicit expressions for the vector spherical harmonics are given by

(38)$${\mathbf{M}}_{mn}(\mathbf{r})=i\frac{m}{\text{sin}\theta }{h}_n^{(1)}(kr)P_n^m(\text{cos}\theta )e^{im\varphi }\widehat{\theta }-h_n^{(1)}(kr)\frac{d{P}_n^m(\text{cos}\theta )}{d\theta }{e}^{im\varphi }\widehat{\varphi },$$

(39)$$\begin{array}{r}\hfill {\mathbf{N}}_{mn}(\mathbf{r})=\frac{n(n+1)}{kr}{h}_n^{(1)}(kr)P_n^m(\text{cos}\theta )e^{im\varphi }\widehat{r}+\frac{1}{kr}{\xi }_n^{\prime }(kr)\frac{d{P}_n^m(\text{cos}\theta )}{d\theta }{e}^{im\varphi }\widehat{\theta }\\ \hfill +i\frac{1}{kr}\frac{m}{\text{sin}\theta }{\xi }_n^{\prime }(kr)P_n^m(\text{cos}\theta )e^{im\varphi }\widehat{\varphi },\end{array}$$

where ξ_n(x) = xh_n(x) is the Ricatti–Bessel function and ξ′_n(x) = dξ(x)/dx.

The integrals of the inner product of the vector spherical harmonics on the spherical surface Γ₀ of radius R are calculated as

(40)$${\iint }_{{\mathrm{\Gamma }}_0}{d}^2\mathbf{r}{\mathbf{M}}_{mn}^{*}\cdot {\mathbf{N}}_{m^{\prime }{n}^{\prime }}=0,$$

(41)$${\iint }_{{\mathrm{\Gamma }}_0}{d}^2\mathbf{r}{\mathbf{M}}_{mn}^{*}\cdot {\mathbf{M}}_{m^{\prime }{n}^{\prime }}={\delta }_{m{m}^{\prime }}{\delta }_{n{n}^{\prime }}{\left|h_n^{(1)}(kR)\right|}^2{R}^2{D}_{mn}\equiv S_{mn}^{(M)}{\delta }_{m{m}^{\prime }}{\delta }_{n{n}^{\prime }},$$

(42)$$\begin{array}{ll}{\iint }_{{\mathrm{\Gamma }}_0}{d}^2\mathbf{r}{\mathbf{N}}_{mn}^{*}\cdot {\mathbf{N}}_{m^{\prime }{n}^{\prime }}\hfill & ={\delta }_{m{m}^{\prime }}{\delta }_{n{n}^{\prime }}{R}^2{D}_{mn}\left[{\left|\frac{{\xi }^{\prime }(kR)}{kR}\right|}^2+n(n+1){\left|h_n^{(1)}(kR)\right|}^2\right]\hfill \\ \hfill & \equiv S_{mn}^{(N)}{\delta }_{m{m}^{\prime }}{\delta }_{n{n}^{\prime }},\hfill \end{array}$$

with

(43)$$D_{mn}=4\pi \frac{n(n+1)}{2n+1}\frac{(n+m)!}{(n-m)!}.$$

B. Hybrid-type finite element method

B.1. Basis functions

Scalar fields in the tetrahedral element are expressed by the nodal-element basis L_(e;i)(r) defined as

(44)$$L_{(e;i)}(\mathbf{r})=a_{(e;i)}^x+a_{(e;i)}^{y}y+a_{(e;i)}^{z}z+b_{(e;i)},$$

where the coefficients $a_{(e;i)}^x$, $a_{(e;i)}^y$, $a_{(e;i)}^z$, and b_(e;i) are obtained from the condition L_(e;i)[r_(e;j)] = δ_ij, with r_(e;i) being the position of node (e; i).

In contrast, the vector fields are expressed by the edge-element basis W_(e;ij)(r) defined as

(45)$${\mathbf{W}}_{(e;ij)}(\mathbf{r})=\frac{1}{l_{ij}}\left[L_{(e;i)}(\mathbf{r})\nabla L_{(e;j)}(\mathbf{r})-L_{(e;j)}(\mathbf{r})\nabla L_{(e;i)}(\mathbf{r})\right],$$

where l_ij is the length of the edge between nodes (e; i) and (e; j). Figure 5 shows the fields of the nodal- and edge-element bases.

 

Fig. 5 Schematic illustrations of the (a) nodal-element basis L_(e;i)(r) and (b) edge-element basis W_(e;ij)(r).

Download Full Size | PDF

The expansion coefficients Φ_s(e;i) in Eq. (22) and $E_{s(e;i)s(e;j)}^{(\alpha )}$ in Eqs. (23) and (24) represent field amplitudes at specific points:

(46)$${\mathrm{\Phi }}_{s(e;i)}=\mathrm{\Phi }[{\mathbf{r}}_{(e;i)}],$$

(47)$$E_{s(e;i)s(e;j)}^{(\alpha )}={\widehat{t}}_{(e;ij)}\cdot {\mathbf{E}}^{(\alpha )}[{\mathbf{r}}_{(e;ij)}],$$

where t̂_(e;ij) denotes a unit vector parallel to the edge between nodes (e; i) and (e; j), and r_(e;ij) = [r_(e;i) + r_(e;j)]/2. Note that the edge-element basis W_(e;ij)(r) automatically satisfies the divergence-free condition [∇ · W_(e;ij) = 0] of transverse fields [23], and thus, the transverse character of E^(α) is guaranteed. This is why edge elements avoid spurious solutions and provide more accurate transverse fields than those calculated using nodal elements.

B.2. Functionals

The functional ℱ₀ in the FEM is obtained by substituting Eq. (15) into Eq. (18) and using Eq. (23) as follows:

(48)$$\begin{array}{ll}{ℱ}_0=\hfill & -\frac{1}{2}\sum _{e\in {\mathrm{\Omega }}_0}\sum _{ijkl=1}^4{K}_{ijkl}^{(e)}(k)E_{s(e;i)s(e;j)}^{(0)*}{E}_{s(e;k)s(e;l)}^{(0)}\hfill \\ \hfill & \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{k}{2}\sum _{e;f_e\in {\mathrm{\Gamma }}_0}\sum _{ij=1}^4\sum _{n=1}^{N_c}\sum _{m=-n}^n{E}_{s(e;i)s(e;j)}^{(0)*}\hfill \\ \hfill & \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\times \{(I_0{p}_{mn}{\delta }_{m,\pm 1}+a_{mn}){\iint }_{{\mathrm{\Gamma }}^{(f_e)}}{\text{d}}^2\mathbf{r}\hspace{0.17em}\hspace{0.17em}{\mathbf{W}}_{(e;ij)}\cdot {[{\widehat{n}}_0\times {\mathbf{N}}_{mn}(\mathbf{r})]}_{f_e}\hfill \\ \hfill & \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+(I_0{q}_{mn}{\delta }_{m,\pm 1}+b_{mn}){\iint }_{{\mathrm{\Gamma }}^{(f_e)}}{\text{d}}^2\mathbf{r}\hspace{0.17em}\hspace{0.17em}{\mathbf{W}}_{(e;ij)}\cdot {[{\widehat{n}}_0\times {\mathbf{M}}_{mn}(\mathbf{r})]}_{f_e}\}\hfill \\ \hfill & +\frac{1}{2}\sum _{e;f_e\in {\mathrm{\Gamma }}_1}\sum _{ij=1}^4{E}_{s(e;i)s(e;j)}^{(0)*}{\iint }_{{\mathrm{\Gamma }}^{(f_e)}}{\text{d}}^2\mathbf{r}\hspace{0.17em}\hspace{0.17em}{\mathbf{W}}_{(e;ij)}\cdot {[{\widehat{n}}_0\times \nabla \times {\mathbf{E}}^{(0)}]}_{f_e},\hfill \end{array}$$

with

(49)$$K_{ijkl}^{(e)}(k)={\iiint }_{{\mathrm{\Omega }}^{(e)}}{\text{d}}^3\mathbf{r}\left\{\left[\nabla \times {\mathbf{W}}_{(e;ij)}\right]\cdot \left[\nabla \times {\mathbf{W}}_{(e;kl)}\right]-k^2{\mathbf{W}}_{(e;ij)}\cdot {\mathbf{W}}_{(e;kl)}\right\},$$

where Ω^(e) is the region of the eth volume element, Γ^(f_e) is the area of the interface f_e, and ∑_{e;fe∈Γ0(Γ1)} denotes the summation over the volume elements e that have one triangular surface f_e belonging to Γ₀ (Γ₁).

The other functional ℱ₁ is obtained by substituting Eq. (12) into Eq. (21) and using Eqs. (22) and (24) as follows:

(50)$$\begin{array}{l}{ℱ}_1=-\frac{1}{2}\sum _{\mu =1}^2\sum _{e\in {\mathrm{\Omega }}_1}\sum _{ijkl=1}^4{K}_{ijkl}^{(e)}(k_{t_{\mu }})E_{s(e;i)s(e;j)}^{(t_{\mu })*}{E}_{s(e;k)s(e;l)}^{(t_{\mu })}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{1}{2}\sum _{\mu =1}^2\sum _{e;f_e\in {\mathrm{\Gamma }}_1}\sum _{ij=1}^4{E}_{s(e;i)s(e;j)}^{(t_{\mu })*}{\iint }_{{\mathrm{\Gamma }}^{(f_e)}}{\text{d}}^2\mathbf{r}\hspace{0.17em}\hspace{0.17em}{\mathbf{W}}_{(e;ij)}\cdot {[{\widehat{n}}_1\times \nabla \times {\mathbf{E}}^{(t_{\mu })}]}_{f_e}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-\frac{1}{2}\sum _{e\in {\mathrm{\Omega }}_1}\sum _{ij=1}^4{J}_{ij}^{(e)}(k_l){\mathrm{\Phi }}_{s(e;i)}^{*}{\mathrm{\Phi }}_{s(e;j)}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{1}{2}\sum _{\mu =1}^2\sum _{e;f_e\in {\mathrm{\Gamma }}_1}\sum _{ij=1}^4\frac{\chi (k_{t_{\mu }},\omega )}{\chi (k_l,\omega )}{\mathrm{\Phi }}_{s(e;i)}^{*}{E}_{s(e;i)s(e;j)}^{(t_{\mu })}{\iint }_{{\mathrm{\Gamma }}^{(f_e)}}{\text{d}}^2\mathbf{r}\hspace{0.17em}\hspace{0.17em}{L}_{(e;i)}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{[{\widehat{n}}_1\cdot {\mathbf{W}}_{(e;ij)}]}_{f_e},\end{array}$$

with

(51)$$J_{ij}^{(e)}(k_l)={\iiint }_{{\mathrm{\Omega }}^{(e)}}{\text{d}}^3\mathbf{r}\hspace{0.17em}\hspace{0.17em}\left\{\left[\nabla L_{(e;i)}\right]\cdot \left[\nabla L_{(e;j)}\right]-k_l^2{L}_{(e;i)}{L}_{(e;j)}\right\}.$$

In the stationary condition [Eq. (25)], the terms coming from the functional derivative of the last term in Eq. (48) and the second term in Eq. (50) can be eliminated using the MBC [Eq. (10)]. Therefore, we can include the MBC [Eq. (10)] by removing these terms, i.e., we use the following functional ℱ̃₀ and ℱ̃₁:

(52)$$\begin{array}{ll}{\widetilde{ℱ}}_0=\hfill & -\frac{1}{2}\sum _{e\in {\mathrm{\Omega }}_0}\sum _{ijkl=1}^4{K}_{ijkl}^{(e)}(k)E_{s(e;i)s(e;j)}^{(0)*}{E}_{s(e;k)s(e;l)}^{(0)}\hfill \\ \hfill & +\frac{k}{2}\sum _{e;f_e\in {\mathrm{\Gamma }}_0}\sum _{ij=1}^4\sum _{n=1}^{N_c}\sum _{m=-n}^n{E}_{s(e;i)s(e;j)}^{(0)*}\hfill \\ \hfill & \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\times \{(I_0{p}_{mn}{\delta }_{m,\pm 1}+a_{mn}){\iint }_{{\mathrm{\Gamma }}^{(f_e)}}{\text{d}}^2\mathbf{r}\hspace{0.17em}\hspace{0.17em}{\mathbf{W}}_{(e;ij)}\cdot {[{\widehat{n}}_0\times {\mathbf{N}}_{mn}(\mathbf{r})]}_{f_e}\hfill \\ \hfill & \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+(I_0{q}_{mn}{\delta }_{m,\pm 1}+b_{mn}){\iint }_{{\mathrm{\Gamma }}^{(f_e)}}{\text{d}}^2\mathbf{r}\hspace{0.17em}\hspace{0.17em}{\mathbf{W}}_{(e;ij)}\cdot {[{\widehat{n}}_0\times {\mathbf{M}}_{mn}(\mathbf{r})]}_{f_e}\},\hfill \end{array}$$

and

(53)$$\begin{array}{l}{\widetilde{ℱ}}_1=-\frac{1}{2}\sum _{\mu =1}^2\sum _{e\in {\mathrm{\Omega }}_1}\sum _{ijkl=1}^4{K}_{ijkl}^{(e)}(k_{t_{\mu }})E_{s(e;i)s(e;j)}^{(t_{\mu })*}{E}_{s(e;k)s(e;l)}^{(t_{\mu })}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}-\frac{1}{2}\sum _{e\in {\mathrm{\Omega }}_1}\sum _{ij=1}^4{J}_{ij}^{(e)}(k_l){\mathrm{\Phi }}_{s(e;i)}^{*}{\mathrm{\Phi }}_{s(e;j)}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{1}{2}\sum _{\mu =1}^2\sum _{e;f_e\in {\mathrm{\Gamma }}_1}\sum _{ij=1}^4\frac{\chi (k_{t_{\mu }},\omega )}{\chi (k_l,\omega )}{\mathrm{\Phi }}_{s(e;i)}^{*}{E}_{s(e;i)s(e;j)}^{(t_{\mu })}{\iint }_{{\mathrm{\Gamma }}^{(f_e)}}{\text{d}}^2\mathbf{r}\hspace{0.17em}\hspace{0.17em}{L}_{(e;i)}{[{\widehat{n}}_1\cdot {\mathbf{W}}_{(e;ij)}]}_{f_e}.\end{array}$$

The total functional ℱ̃ is then given by ℱ̃ = ℱ̃₀ + ℱ̃₁.

Acknowledgments

The authors would like to gratefully acknowledge Gang Bao (MSU, USA) for creating opportunity to start this study and for stimulating discussions. This work was supported by Grant-in-Aid for Scientific Research (C), No. 25400325.

References and links

1. K. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag. 14, 302–307 (1966). [CrossRef]  

2. J.-M. Jin, The Finite Element Method in Electromagnetics (John Wiley, 2002).

3. P. K. Banerjee and R. Butterfield, Boundary Element Methods in Engineering Science (McGraw-Hill, 1981).

4. Y. Akahane, T. Asano, B.-S. Song, and S. Noda, “High-Q photonic nanocavity in a two-dimensional photonic crystal,” Nature 425, 944–947 (2003). [CrossRef]   [PubMed]  

5. Y. Mizumoto and H. Ishihara, “Theory of resonant radiation force exerted on single organic molecules near metallic nanogap,” Proc. SPIE 8097, 80971C (2011). [CrossRef]  

6. M. Takase, H. Ajiki, Y. Mizumoto, K. Komeda, M. Nara, H. Nabika, S. Yasuda, H. Ishihara, and K. Murakoshi, “Selection-rule breakdown in plasmon-induced electronic excitation of an isolated single-walled carbon nanotube,” Nat. Photonics 7, 550–554 (2013). [CrossRef]  

7. S. Uryu, H. Ajiki, and H. Ishihara, “Model of finite-momentum excitons driven by surface plasmons in photoexcited carbon nanotubes covered by gold metal films,” Phys. Rev. Lett. 110, 257401 (2013). [CrossRef]   [PubMed]  

8. H. Ajiki, “Calculation method for exciton wavefunctions with electron–hole exchange interaction: Application to carbon nanotubes,” J. Phys. Soc. Jpn. 82, 4701 (2013). [CrossRef]  

9. S. I. Pekar, “The theory of electromagnetic waves in a crystal in which excitons are produced,” Sov. Phys. JETP 6, 785–796 (1958).

10. J. Birman, E. Rashba, and M. Sturge, Excitons (North-Holland, 1982).

11. V. M. Agranovich and V. L. Ginzburg, Crystal Optics with Spatial Dispersion, and Excitons (Springer, 1984). [CrossRef]  

12. A. d’Andrea and R. Del Sole, “Wannier-mott excitons in semi-infinite crystals: Wave functions and normal-incidence reflectivity,” Phys. Rev. B 25, 3714–3730 (1982). [CrossRef]  

13. P. Halevi and R. Fuchs, “Generalized additional boundary condition for non-local dielectrics. I. Reflectivity,” J. Phys. C Sol. Stat. Phys. 17, 3869–3888 (1984). [CrossRef]  

14. K. Cho and M. Kawata, “Theoretical analysis of polariton interference in a thin platelet of CuCl. I. Additional boundary condition,” J. Phys. Soc. Jpn. 54, 4431–4443 (1985). [CrossRef]  

15. K. Cho, “Nonlocal theory of radiation-matter interaction: Boundary-condition-less treatment of Maxwell equations,” Prog. Theor. Phys. Suppl. 106, 225–233 (1991). [CrossRef]  

16. K. Cho, Optical Response of Nanostructures: Microscopic Nonlocal Theory (Springer, 2003). [CrossRef]  

17. Y. Sun, H. Ajiki, and G. Bao, “Computational modeling of optical response from excitons in a nano optical medium,” Commun. Comput. Phys. 4, 1051–1068 (2008).

18. H. Ajiki and K. Cho, “Longitudinal and transverse components of excitons in a spherical quantum dot,” Phys. Rev. B 62, 7402–7412 (2000). [CrossRef]  

19. K. Cho and H. Ishihara, “ABC-theory of an exciton-polariton in a slab,” J. Phys. Soc. Jpn. 59, 754–764 (1990). [CrossRef]  

20. J. Lagois, “Depth-dependent eigenenergies and damping of excitonic polaritons near a semiconductor surface,” Phys. Rev. B 23, 5511 (1981). [CrossRef]  

21. M. Zamfirescu, A. Kavokin, B. Gil, G. Malpuech, and M. Kaliteevski, “ZnO as a material mostly adapted for the realization of room-temperature polariton lasers,” Phys. Rev. B 65, 161205 (2002). [CrossRef]  

22. A. Bossavit, “A rationale for ‘edge-elements’ in 3-D fields computations,” IEEE Trans. Mag. 24, 74–79 (1988). [CrossRef]  

23. A. Bossavit, “Solving Maxwell equations in a closed cavity, and the question of ‘spurious modes’,” IEEE Trans. Mag. 26, 702–705 (1990). [CrossRef]  

24. R. Ruppin, “Optical absorption by excitons in microcrystals,” J. Phys. Chem. Solids 50, 877–882 (1989). [CrossRef]  

25. A. Ekimov, A. L. Efros, and A. Onushchenko, “Quantum size effect in semiconductor microcrystals,” Solid State Commun. 56, 921–924 (1985). [CrossRef]  

26. K. Inaba, K. Imaizumi, K. Katayama, M. Ichimiya, M. Ashida, T. Iida, H. Ishihara, and T. Itoh, “Optical manipulation of CuCl nanoparticles under an excitonic resonance condition in superfluid helium,” Phys. Stat. Solidi B 243, 3829–3833 (2006). [CrossRef]  

27. Z. Tang, G. K. Wong, P. Yu, M. Kawasaki, A. Ohtomo, H. Koinuma, and Y. Segawa, “Room-temperature ultraviolet laser emission from self-assembled ZnO microcrystallite thin films,” Appl. Phys. Lett. 72, 3270–3272 (1998). [CrossRef]