We develop accurate finite-difference time-domain (FDTD) modeling of polymer bulk heterojunction solar cells containing Ag nanoparticles between the hole-transporting layer and the transparent conducting oxide-coated glass substrate in the wavelength range of 300 nm to 800 nm. The Drude dispersion modeling technique is used to model the frequency dispersion behavior of Ag nanoparticles, the hole-transporting layer, and indium tin oxide. The perfectly matched layer boundary condition is used for the top and bottom regions of the computational domain, and the periodic boundary condition is used for the lateral regions of the same domain. The developed FDTD modeling is employed to investigate the effect of geometrical parameters of Ag nanospheres on electromagnetic fields in devices. Although negative plasmonic effects are observed in the considered device, absorption enhancement can be achieved when favorable geometrical parameters are obtained. Plasmonics [1] have received increasing attention for a variety of applications related to biosensing [2] and compact nanophotonics [3] – [5] due to their ability to produce large field enhancement and subwavelength field confinement. Recently, plasmonic light-trapping geometries using metal nanoparticles have been employed to improve the efficiency of solar cells [6] . Metal nanoparticles have been widely utilized for inorganic thin-films [7] , organic thin-films [6] , [8] , [9] , and dye-sensitized solar cells [10] . Among various solar cell technologies, organic photovoltaics (OPVs) have been of particular interest for their use in the production of large-area flexible modules [11] due to high-throughput, low temperature processes [12] for low-cost roll-to-roll manufacturing with an improved environmental stability [13] . Plasmonic OPVs have been investigated both numerically and experimentally. Among the numerical techniques used, the finite-difference time-domain (FDTD) method [14] – [17] has been widely employed because of its accuracy, robustness, and matrix-free characteristics. Moreover, a single FDTD simulation can compute a wideband response by using a Fourier transform, since it is a time-domain method. In FDTD, the frequency-dependent permittivity of materials in plasmonic OPVs should be incorporated by an appropriate dispersion model. However, the previous FDTD analyses have not considered dispersive properties of OPV materials due to the difficulty involved in doing so. In fact, because the real part of the relative permittivity of Ag is negative, a dispersive FDTD algorithm should be applied to Ag so that the resulting FDTD algorithm does not suffer from instability. Without a proper FDTD dispersion model, one should perform many simulations over the wavelengths of interest by way of adispersive FDTD for Ag and a non-dispersive FDTD for other OPV materials (with the corresponding permittivity and conductivity at a specific wavelength), which leads to overwhelming computational costs. Therefore, it is of great interest to develop accurate FDTD dispersive modeling for the optical analysis of plasmonic OPVs. In this work, we develop — based on the Drude dispersion model — FDTD dispersive modeling for plasmonic OPVs. The perfectly matched layer (PML) [18] – [19] and the periodic boundary condition (PBC) [14] are used for the termination and lateral regions of the computational domain, respectively. We also employ the proposed FDTD algorithm to investigate the effect of the geometrical parameters of Ag nanospheres on electromagnetic fields in the photoactive layer. It is worth noting that the purpose of this paper is to develop FDTD dispersive modeling suitable for plasmonic OPVs, not to optimize plasmonic OPVs for improved performance. We consider polymer:fullerene bulk heterojunction (BHJ) solar cells. For the plasmonic structure, the self-assembled Ag nanospheres are formed between the indium tin oxide (ITO)-coated glass substrate and the poly(3, 4-ethylene dioxythiophene) polystyrene sulfonate (PEDOT:PSS) — the latter acting as the hole-transporting layer. The photoactive layer is poly(3-hexylthiophene) (P3HT) and [6,6]-phenyl-C ₆₁ -butyric acid methyl ester (PCBM) blend film. The final structure considered here is ITO/Ag nanospheres/PEDOT:PSS/ P3HT:PCBM. The schematics of the side and top views of the considered structure are shown in Fig. 1. Here, we assume the hexagonal configuration of Ag nanospheres. As previously alluded to, FDTD has been widely used to study plasmonic structures [4] – [6] , [20] – [22] . In this work, we develop broadband-accurate FDTD modeling of plasmonic polymer BHJ solar cells. For FDTD dispersion modeling, we apply the Drude dispersion model for Ag, PEDOT:PSS, and ITO as follows: (1)$${\epsilon }_r\left(\omega \right)={\epsilon }_{\infty }+\frac{{\omega }_p^2}{j\omega \left(\Gamma +j\omega \right)},$$ where ε_∞ is the relative permittivity at the high-frequency limit, ω_p is a plasma frequency, and Γ is a damping coefficient. These Drude parameters can be obtained by using a nonlinear optimization technique. The extracted Drude parameters and the corresponding root-mean-square errors (RMSEs) are listed in Table 1 . It should be noted that the Drude dispersion model has been widely employed for Ag [5] , [20] – [24] . More complex dispersion models (for example, multispecies Drude-Lorentz dispersion models) can improve computational accuracy [25] but tend to lead to heavy computational burdens. Therefore, the Drude dispersion model is employed in this work. To the best of our knowledge, the Drude dispersion model is successfully applied to PEDOT:PSS [26] and ITO [27] for the first time. Figure 2 shows the relative permittivity of PEDOT:PSS and ITO — both of which were under our Drude dispersion model and the corresponding experimental data. The Drude dispersion model agrees well with the experimental data for PEDOT:PSS, however, deviates a little from the experimental data for ITO. Albeit with this error, the Drude dispersion model for ITO is superior to that of the already widely-used Drude dispersion model for Ag (see Table 1 ). Note that the relative permittivity of P3HT:PCBM is set as 3.4, which is the same as in [28] . As previously alluded to, without our Drude dispersion model, many simulations should be performed over the wavelengths of interest under a sine-wave excitation, which leads to overwhelming computational costs. It is also noted that more complex dispersion models should be employed when a wider range of wavelengths is considered. Now, let us derive the FDTD update equations for the Drude dispersion model. Toward this purpose, we first consider Maxwell’s Ampere’s law in the frequency domain as follows: (2)$$\nabla \times H\left(\omega \right)=j\omega {\epsilon }_0{\epsilon }_r\left(\omega \right)E\left(\omega \right).$$ Inserting (1) into (2), we have (3) $$\nabla \times H\left(\omega \right)=j\omega {\epsilon }_0{\epsilon }_{\infty }E\left(\omega \right)+\frac{{\epsilon }_0{\omega }_p^2}{j\omega +\Gamma }E\left(\omega \right).$$ Introducing the equivalent current J ( ω ) [20] and then rearranging the resulting equation, we have (4)$$\nabla \times H\left(\omega \right)=j\omega {\epsilon }_0{\epsilon }_{\infty }E\left(\omega \right)+J\left(\omega \right),$$ (5) $$j\omega J\left(\omega \right)+\Gamma J\left(\omega \right)={\epsilon }_0{\omega }_p^{2}E\left(\omega \right).$$ Applying the inverse Fourier transform to the above equations, we obtain (6) $${\epsilon }_0{\epsilon }_{\infty }\frac{\partial }{\partial t}E(t)+J(t)=\nabla \times H(t),$$ (7) $$\frac{\partial }{\partial t}J(t)+\Gamma J(t)={\epsilon }_0{\omega }_p^{2}E(t).$$ Applying the central difference scheme (CDS) to the temporal derivatives, we can obtain the FDTD update equations for E ⁿ⁺¹ and J ⁿ⁺¹ below (8) $$E^{n+1}=E^n+\frac{\Delta t}{{\epsilon }_0{\epsilon }_{\infty }}\nabla \times H^{n+1/2}-\frac{\Delta t}{2{\epsilon }_0{\epsilon }_{\infty }}\left(J^{n+1}+J^n\right),$$ (9)$$J^{n+1}=C_{Ja}{J}^n+C_{Jb}\left(E^{n+1}+E^n\right),$$ where (10)$$C_{Ja}=\frac{1-0.5\Delta t\Gamma }{1+0.5\Delta t\Gamma },$$ (11)$$C_{Jb}=\frac{0.5\Delta t{\epsilon }_0{\omega }_p^2}{1+0.5\Delta t\Gamma }.$$ Note that A ⁿ indicates a vector field component A at the time n Δ t , where Δ t is the FDTD step size [14] . Note that we cannot update directly E ⁿ⁺¹ and J ⁿ⁺¹ from (8) and (9), since field values at simultaneous times are involved. Therefore, we insert (9) into (8) and then manipulate the resulting equation, which leads to the final FDTD update equation for E ⁿ⁺¹ , which is given as (12)$$E^{n+1}=C_{Ea}{E}^n+C_{Eb}\nabla \times H^{n+1/2}+C_{Ec}{J}^n,$$ where (13)$$C_{Ea}=\frac{C_{\alpha }-C_{Jb}}{C_{\alpha }+C_{Jb}},$$ (14)$$C_{Eb}=\frac{2}{C_{\alpha }+C_{Jb}},$$ (15)$$C_{Ec}=-\frac{1+C_{Ja}}{C_{\alpha }+C_{Jb}},$$ with C_α = 2 ε ₀ ε _∞ /Δ t . In the above equations, spatial discretization should be performed by applying the CDS to the curl operator. The FDTD update equation for H ^n+1/2 can be simply obtained by using the standard difference scheme in Maxwell’s Faraday’s law. Note that in every time step we update H ^n+1/2 , E ⁿ⁺¹ , and J ⁿ⁺¹ sequentially. In FDTD, it is necessary to apply appropriate boundary conditions to truncate the computational domains. In this work, the PML is used for the top and bottom regions to avoid spurious reflections from the external grid boundaries by employing the complex-coordinate stretching technique [18] – [19] . In lateral regions, the PBC [14] is used to reduce computational burdens by using the periodicity of the structure. In this section, we apply our homemade FDTD modeling to analyze the optical responses of plasmonic OPVs. The step size is Δ x = Δ y = Δ z = 0.2 nm to 1.6 nm, depending on the radius ( r ) of Ag nanospheres used to accurately model various nanospheres (that is, 10 cells per the radius of a sphere). The time step size is given by Δ t = 0.95Δ x / c ₀ /√3 to satisfy the stability condition, where c ₀ is the vacuum light speed in m/s. An x-polarized uniform plane wave is uniform in the space domain, and a pulse is modulated by a sine-wave in the time domain. To obtain the spectral results ( E_x ( λ ), E_y ( λ ), and E_z ( λ )), a discrete Fourier transform is used for the temporal results ( E_x ( t ), E_y ( t ), and E_z ( t )) in the photoactive layer. Plasmonic effect is quantitatively estimated using a field intensity ratio (FIR) — that is, the ratio of the square of the electric fields found (| E ( λ )| ² = | E_x ( λ )| ² +| E_y ( λ )| ² +| E_z ( λ )| ² ) in the photoactive layer for the plasmonic device (with Ag nanospheres) and the control (without Ag nanospheres) structure. Before proceeding with plasmonic OPVs that are based on the hexagonal periodicity of Ag nanospheres in three-layered media, we analyze a plasmonic OPV based on a single Ag nanosphere in three-layered media. In this case, the plasmonic OPV is simulated by replacing the PBC by the two-stage PML [21] in the lateral regions. Figure 3 shows FIRs for various radii. For comparison, we also simulate the structure for r = 10 nm using a Drude dispersive FDTD for Ag and a non-dispersive FDTD for ITO and PEDDOT:PSS (with the corresponding permittivity and conductivity at five different wavelengths), the results of which are indicated by diamond-shaped symbols in Fig. 3 . Please note that only a single simulation is performed for our proposed FDTD, however, five different such simulations should be performed for the comparative study. As shown in Fig. 3 , the resulting graph plots agree with each other. As the radii of Ag nanospheres increase, the peaks and bases of the graph plots are clearly observed. For example, at r = 10 nm, a single peak (base) occurred at λ ≈ 422 nm ( λ ≈ 401 nm), which is consistent with the existence of forward-direction enhancement at λ > λ _sp and backward-direction enhancement at λ > λ _sp [29] . In addition, red shifts are observed when r is increased. Now, let us consider a plasmonic OPV based on the hexagonal periodicity of Ag nanospheres in three-layered media, as shown in Fig. 1 . First, we illustrate the effect of r on FIR at d / r = 2, where d is the particle-to-particle spacing of Ag nanospheres (see Fig. 1 ). Figure 4 shows FIR versus wavelength for various sizes of Ag nanospheres, revealing that FIR is highly dependent on r . As the radii of Ag nanospheres increase, FIR decreases. Differently from the single Ag nanosphere case, no red shifts are observed when r is increased, which is due to complicated coupling between neighboring Ag nanospheres. Note that the single Ag nanosphere case is shown in Fig. 3 . Also note that ripples are observed and that they are more distinct for smaller r . Next, Fig. 5 illustrates the effect of d / r on FIR with a fixed radius of r = 10 nm. Similar to the size of Ag nanospheres, FIR depends on the particle-to-particle spacing of Ag nanospheres. As the particle-to-particle spacing decreases, smaller FIR is observed due to strong coupling between neighboring Ag nanospheres. As shown in Fig. 5 , FIR ≈ 0 (that is, | E | ² ≈ 0) is observed. In other words, not much light can penetrate into the photoactive layer. This phenomenon may result from complicated coupling between neighboring Ag nanospheres. We also examine FIR for various simultaneous r and d / r . We consider r = 2 nm to 16 nm and d / r = 0.55 nm to 4nm. The maximum r is set as 16 nm because we assume that metal nanoparticles are embedded in PEDOT:PSS with some geometrical margins. Note that the thickness of PEDOT:PSS is 40 nm for the device. Figure 6 shows FIR versus r and d / r . In this case, we integrate FIR over the considered spectrum (300 nm to 800 nm). As shown in Fig. 6 , it is noted that larger r and smaller d / r lead to smaller FIR. It is worth noting that the absorbance of the photoactive layer is proportional to | E | ² ; thus, the short circuit current is also proportional to | E | ² for the ideal carrier transport condition (that is, the perfect internal quantum efficiency) [30] . Figure 7 shows the time-averaged | E | ² distribution in the considered plasmonic OPV, with r = 10 nm and d / r = 2 at the wavelength of 413 nm. Note that this time-averaged | E | ² distribution is normalized by the time-averaged | E | ² distribution of the control structure (without Ag nanospheres). Three snapshots are depicted on the xy-plane at the center of the Ag nanosphere ( z = 90 nm), the yz -plane at x = 0 nm, and the xz -plane at y = 0 nm. Strong field intensity is observed near and between Ag nanospheres, but weak field intensity is shown inside the photoactive layer (120 nm ≤ z ≤ 200 nm). This explains why negative plasmonic effects are observed in the considered device. It is well known that a single simulation cannot exactly emulate the experiment. In this section, we discuss the differences between our simulation and the experiment. These differences are summarized in Table 2 . In our developed FDTD, we assume a perfect geometry. However, this is impossible in reality due to manufacturing errors, especially for Ag nanospheres. In addition, modeling errors are inherent to FDTD algorithms because of staircase approximations, but such errors can be alleviated by using a conformal-path model [14] . In our simulation, we have not considered the cathode of plasmonic OPVs, because of computational costs (memory and CPU time). In fact, the wave-guiding mode between the cathode and the P3HT:PCBM/PEDOT:PSS interface can increase the field intensity inside the photoactive layer. To fully consider a plasmonic OPV and optimize it, hardware-acceleration techniques, such as GPU-FDTD [14] or MPI-FDTD [31] , can be employed. As explained previously, we have used the constant permittivity of P3HT:PCBM, as done in [28] . However, in reality the real part of the permittivity of the P3HT:PCBM is not constant and also the imaginary part of the permittivity of the P3HT:PCBM does in fact exist. To accurately consider dispersive characteristics of P3HT:PCBM, new complex dispersion models, such as the complex conjugate dispersion model [32] or the quadratic complex rational function dispersion model [17] , [33] , can be employed. The improvement of our dispersive FDTD algorithm is currently under investigation to overcome the above-mentioned limitations. As a part of the development of a systematic design tool for optimizing plasmonic OPVs, broadband-accurate FDTD modeling has been developed for plasmonic polymer:fullerene BHJ solar cells. The Drude dispersion model has been applied for Ag, PEDOT:PSS, and ITO. To apply proper boundary conditions, the PML and the PBC are employed for the termination regions and the lateral regions, respectively. We have examined the effects of the size of Ag nanospheres and their inter-particle spacing on field intensity in the photoactive layer. Although negative plasmonic effects are observed in the considered device, the absorption enhancement in a photoactive layer can be achieved by optimizing the shape, size, and inter-particle spacing of Ag nanoparticles within the wide range of geometrical parameters available and also by changing the location of Ag nanoparticles [34] – [35] . We hope that our FDTD dispersive modeling can be a foundation stone of a design tool for optimizing plasmonic-enhanced OPVs. Kyung-Young Jung received his BS and MS degrees in electrical engineering from Hanyang University, Seoul, Rep. of Korea, in 1996 and 1998, respectively and his PhD degree in electrical and computer engineering from The Ohio State University, Columbus, USA, in 2008. From 2008 to 2009, he was a postdoctoral researcher with The Ohio State University, and from 2009 to 2010, he was an assistant professor with the Department of Electrical and Computer Engineering, Ajou University, Suwon, Rep. of Korea. Since 2011 he has been with Hanyang University, where he is now an assistant professor with the Department of Electronic Engineering. He is a recipient of the graduate study abroad scholarship from the National Research Foundation of Korea, the presidential fellowship from The Ohio State University, and the distinguished lecturer from Hanyang University. His current research interests include computational electromagnetics, bio electromagnetics, and plasmonics. Woo-Jun Yoon received his BS and MS degrees in metallurgical engineering from Korea University, Seoul, Rep. of Korea, in 1999 and 2001, respectively and his MS and PhD degrees in electrical engineering from The Ohio State University, Columbus, USA, in 2006 and 2009, respectively. His research experience encompasses modeling, design, processing, and characterization of optoelectronic devices for photovoltaic applications, with emphasis on organic photovoltaics. He is currently working as a research associate at the Naval Research Laboratory, Washington, DC, USA. He has coauthored over 30 journal and conference publications. Yong Bae Park received his BS (summa cum laude), MS, and PhD degrees in electrical engineering from the Korea Advanced Institute of Science and Technology, Daejeon, Rep. of Korea, in 1998, 2000, and 2003, respectively. From 2003 to 2006, he was with the Korea Telecom Laboratory, Seoul, Rep. of Korea. In 2006, he joined the School of Electrical and Computer Engineering, Ajou University, Suwon, Rep. of Korea, where he is now an associate professor. His research interests include electromagnetic field analysis and electromagnetic interference and compatibility. Paul R. Berger is a professor of electrical & computer engineering at Ohio State University, Columbus, USA. He is the founder of the Nanoscale Patterning Laboratory. He received his BSE in engineering physics and his MSE and PhD (1990) degrees in electrical engineering all from the University of Michigan, Ann Arbor, USA. Currently, he is actively working on conjugated polymer-based optoelectronic and electronic devices; molecular electronics; Si/SiGe nanoelectronic devices and fabrication processes; Si-based resonant interband tunneling diodes and quantum functional circuitry; bioelectronics; and the fabrication and growth of semiconductor materials. Formerly, he worked at Bell Laboratories, Murray Hill, NJ (1990–1992) and taught electrical and computer engineering at the University of Delaware (1992–2000). In 1999, he took sabbatical leave and went on to work firstly at the Max-Planck Institute for Polymer Research, Mainz, Germany, while supported by Prof. Dr. Gerhard Wegner, and then at Cambridge Display Technology, Ltd., Cambridge, United Kingdom, working under Dr. Jeremy Burroughes. In 2008, he spent an extended sabbatical leave at the Interuniversity Microelectronics Center in Leuven, Belgium, while appointed as a visiting professor in the Department of Metallurgy and Materials Engineering, Katholieke Universiteit Leuven, Belgium. He has authored in excess of 100 articles and five book sections and has been issued 17 patents, with five more pending from 50+ disclosures. Some notable recognitions for him were an NSF CAREER Award (1996), a DARPA ULTRA Sustained Excellence Award (1998), a Lumley Research Award (2006, 2011), and a Faculty Diversity Excellence Award (2009). He has been on the program and advisory committees of numerous conferences, including the IEDM and ISDRS meetings. He is currently the chair of the Columbus IEEE EDS/LEOS Chapter and Faculty Advisor to Ohio State’s IEEE Student Chapters. He is a fellow and distinguished lecturer of IEEE EDS and a senior member of the Optical Society of America. Fernando L. Teixeira received his BS and MS degrees in electrical engineering from the Pontifical Catholic University of Rio de Janeiro, Brazil and his PhD degree in electrical engineering from the University of Illinois, Urbana-Champaign, USA. He was a postdoctoral associate with the Massachusetts Institute of Technology, and since 2000 he has been with The Ohio State University, Columbus, USA, where he is now a professor with the Department of Electrical and Computer Engineering and affiliated with the ElectroScience Laboratory. He is a recipient of the NSF CAREER Award, the triennial Booker Fellowship from the International Union of Radio Science, and the outstanding young engineer award from the IEEE Microwave Society. He currently serves as an associate editor for the IEEE Antennas and Wireless Propagation Letters. His current research interests include computational electromagnetics; plasmonics and metamaterials; electromagnetic sensing; and inverse scattering.