Novel semi-analytical optoelectronic modeling based on homogenization theory for realistic plasmonic polymer solar cells

Numerical-based simulations of plasmonic polymer solar cells (PSCs) incorporating a disordered array of non-uniform sized plasmonic nanoparticles (NPs) impose a prohibitively long-time and complex computational demand. To surmount this limitation, we present a novel semi-analytical modeling, which dramatically reduces computational time and resource consumption and yet is acceptably accurate. For this purpose, the optical modeling of active layer-incorporated plasmonic metal NPs, which is described by a homogenization theory based on a modified Maxwell–Garnett-Mie theory, is inputted in the electrical modeling based on the coupled equations of Poisson, continuity, and drift–diffusion. Besides, our modeling considers the effects of absorption in the non-active layers, interference induced by electrodes, and scattered light escaping from the PSC. The modeling results satisfactorily reproduce a series of experimental data for photovoltaic parameters of plasmonic PSCs, demonstrating the validity of our modeling approach. According to this, we implement the semi-analytical modeling to propose a new high-efficiency plasmonic PSC based on the PM6:Y6 PSC, having the highest reported power conversion efficiency (PCE) to date. The results show that the incorporation of plasmonic NPs into PM6:Y6 active layer leads to the PCE over 18%.

Apart from the researches mentioned above, there are many other experimental studies on the plasmonic PSCs [59][60][61][62][63][64] , while there are a few theoretical studies for the simulation of plasmonic PSCs that almost all of them cannot be successfully applied to realistic plasmonic PSCs. To investigate the effects of certain properties of incorporating NPs, some simplifications have considered in the theoretical studies in the literature [65][66][67][68][69][70][71][72][73] that should be addressed for the modeling of realistic plasmonic PSCs. Two of the simplifications include the incorporation of plasmonic NPs with the same size in an ordered array, while the experimental reports have documented the randomly blending of NPs with various sizes in the PSCs. Both of these factors cause significant changes in optical properties and consequently the PCE of plasmonic PSCs. Furthermore, the effects of interference and reflection introduced by different layers on the absorption of NPs and the effects of increased trap-assisted recombination due to NPs on the electrical properties have to be considered for realistic modeling. The implementation of these conditions to the modeling of realistic plasmonic PSCs through numerical methods leads to prohibitively long-time complex computation process due to the consideration an enormous collection of NPs, which is far beyond the reach of even modern computers. To overcome this limitation, in this study, we present a novel semi-analytical modeling for predicting the performance of the realistic structure of plasmonic PSCs to dramatically reduce computational time and resource consumptions and yet is acceptably accurate. Hence, this paper is structured as follows. The modeling is explained in Section II, where the geometrical parameters of plasmonic BHJ PSC and some assumptions we make to implement modeling are expressed in section II-A, the effect of embedded plasmonic NPs into PSCs on optical properties is described by an analytical modeling based on homogenization theory (HT) in section II-B, and the electrical properties of plasmonic PSCs are obtained in section II-C with solving the coupled equations of Poisson, continuity, and drift-diffusion. In section III, to evaluate the applicability of the semi-analytical modeling, its results are compared with the experimental results of a fabricated plasmonic PSC. Besides, the influence of the NPs parameters, including the amount of size dispersion and concentration of the NPs, on the performance of plasmonic PSCs is discussed. In section IV, based on the reported PSC with the best performance so far, a new high efficiency plasmonic PSC is proposed and investigated. Optical modeling for the effect of incorporated plasmonic NPs. Optical properties of an ordered array of metal NPs incorporated into the active layer can be accurately obtained by considering a few NPs and defining their geometry exactly and then solving the Maxwell equations using numerical techniques such as finite difference time domain (FDTD) 86 , boundary element method (BEM) 87 , discrete dipole approximation (DDA) 88 , etc. In the case of irregularly incorporated NPs with various sizes, a large number of NPs must be considered for defining the geometry of active layer, and the aforementioned numerical methods are much less tractable due to the high dependence of calculation time on the size of the system. To surmount the shortcomings of numerical methods in obtaining the optical properties of a vast collection of plasmonic NPs, analytical approaches based on the HTs are methodical [89][90][91] . In the HTs, a complex medium formed by the inclusion of NPs into a material is replaced with a homogeneous medium that has the same optical properties as the complex medium [92][93][94][95] . Therefore, the optical properties of the active layer:NPs composite can be expressed through the HTs with the complex dielectric function of this homogeneous medium, simply, ε HM .
A conventional HT is Maxwell-Garnett theory (MGT), derived from the Lorentz-Lorenz relations or the Clausius-Mossotti formula 96 . It averages over the induced electric dipole moments of individual NPs without considering the interaction between NPs to calculate the LSPR band of NPs with the same sizes 97 , so it fails to apply for NPs with size dispersion or with low interparticle distance (small distance among neighboring NPs leads to the high volume fraction of NPs). Furthermore, the MGT is based on the quasistatic limit and ignores retardation effects, so it produces significant errors for NPs with diameters that the electrostatic limit is no longer valid and the retardation effects become predominant, namely, it is restricted to the spherical NPs with the diameter (d) much smaller than the wavelength of incident photons (λ), i.e., d < < λ 98-100 . For example, d < 6 nm for Au NPs and d < 3 nm for Ag NPs are acceptable sizes for applying the MGT 101 . Moreover, the dependence of intrinsic confinement effects, induced for very small NPs, smaller than the mean free path of conduction electrons (typically d < 4 nm) 102,103 , on the NPs size is not taken into account in the MGT. Indeed, the size of NPs does not explicitly appear in the MGT.
To remove the restrictions mentioned above, different extensions of the MGT have been proposed 97,[104][105][106][107][108][109][110][111] . For example, a corrected version of MGT which accounts for a dipolar interaction between NPs is obtained by Markel et al. 97 , or an extension of MGT considering both intrinsic confinement and retardation effects, called Maxwell-Garnett-Mie theory (MGMT), is achieved by replacing quasistatic electric dipole polarizability with that of obtained by the Mie theory [109][110][111] .
In this paper, to model the optical properties of active layer:NPs composite, a modified MGMT is developed by considering the size dispersion of NPs, size-dependent intrinsic confinement for very small NPs, and retardation effects for large NPs. Therefore, it can predict the LSPR band of NPs distributed over a wide range of sizes, from very small to relatively large. Besides, light scattering within the active layer by NPs is considered by an additional contribution to the modified MGMT.
The MGMT gives the complex dielectric function of active layer:NPs composite in terms of volume fraction of NPs in the active layer (f), mean radius of NPs (R̅ ), and frequency-and size-dependent Mie polarizability (α Mie (ω,R̅ )) as 90 where ψ 1 and ζ 1 are the first order of Riccati-Bessel functions of the first and second kind, respectively, and m ω and x ω are defined as: where ε NP (ω,R̅ ) is the size-dependent dielectric function of NPs. It should be noted that the effect of intrinsic confinement is considered in the MGMT through the implementation of ε NP , instead of using the dielectric function of bulk metal ( ε bm ), in the M e (1) (ω,R̅ ). By assuming that the only effect of NPs size is on the free electrons, ε NP can be derived from the Matthiessen rule by modifying ε bm , described by Lorentz-Drude Model, as 112 : where ω p is the plasma frequency, Г 0 is the damping constant, A s is the parameter depending on the scattering process of the electrons of NP surface 113 , and v f is the Fermi velocity of free electrons.
The effect of size dispersion of NPs can be included in the MGMT, Eq. (1a), by considering the mean-field theory. Size dispersion leads to various electric dipole moments for NPs with different radii. Therefore, the average polarizability is calculated by weighting the polarizabilities over the relative abundance of each NP. Hence, by considering Gaussian distribution with the standard deviation of γ for the NPs radii, M e (1) (ω,R̅ ) in Eq. (1b) would be replaced with the following expression: where R min and R max stand for the smallest and largest radius of NPs in the size dispersion.
The homogenized dielectric function ( ε HM ), describing the optical properties of the active layer:NPs composite, must address all the mechanisms resulted from inserting NPs. Therefore, in addition to the impact of plasmonic near-field due to LSPR excitation, considered through ε MGMT , the effect of light scattering by the embedded NPs within the active layer must be reflected in the ε HM . Because of the size dispersion of NPs, the absorption mechanism is partly attributed to enhanced LSPR near-field around the small size NPs and partly attributed to light scattering from the large size NPs 114-117 that disperse the electromagnetic waves of the incident light. The reemitting of incident light in different directions inside the active layer leads to an increase in the optical path length 20 . The effect of enhanced optical path length by the specific angular spread of scattered light can be expressed by the Percus-Yevick correction term 89,90 . This term is added to the ε MGMT , Eq. (1a), to obtain ε HM for the active layer:NPs composite as: Electrical modeling of plasmonic BHJ PSCs. The mechanism of generating electron-hole pairs and their transport in plasmonic BHJ PSCs, like pristine BHJ PSCs (without plasmonic NPs), is that the absorbed photons by the active layer cause the transition of electrons from the highest occupied molecular orbital (HOMO) of electron-donating material to its lowest unoccupied molecular orbital (LUMO) and creating neutral Frenkel excitons (FEs) with the generation rate of G F . Generated FEs diffuse to the donor:acceptor interface (10-20 nm) and then dissociate into electrons on the LUMO of electron-accepting material and holes on the HOMO of electron-donating material on either side of the interface with the Coulomb interaction between them, called charge transfer excitons (CTEs). CTEs will undergo recombination to FEs after a finite time unless induce to separate 118,119 . The motion of electrons causes the dissociation of CTEs into free electrons and holes moving towards the corresponding electrodes by incoherent hopping between localized states randomly distributed in www.nature.com/scientificreports/ space due to the field arising from the difference of the energy levels of intermediate layers or the work functions of electrodes. Therefore, for electrical modeling, following mechanisms have to be taken into account: (1) the generation, dissociation, and recombination of CTEs, (2) generation and recombination of free charges, (3) drift and diffusion of charges, and (4) the extraction of charges at the electrodes. To consider these mechanisms, several one-dimensional electrical models differing in the choice of their components, the definition of boundary conditions, and the method of solving drift-diffusion equations have been developed in the literature 80,83,[120][121][122][123][124][125][126] . In the following, while expressing the coupled equations of continuity, drift-diffusion, and Poisson for obtaining the density of electrons and holes (n and p) and electric potential (φ) in the plasmonic BHJ PSCs, the effect of plasmonic NPs on the aforementioned mechanisms will be clarified.
where Eqs. (6) and (7) are the continuity equations for electrons and holes, respectively, and Eq. (8) is the Poisson equation, j n and j p , comprised of drift and diffusion components, are the electron and hole current densities, respectively, q is the elementary charge, μ n and μ p are the mobility of electrons and holes, respectively, N A , N D , and n RC are the densities of ionized acceptors, donors, and trapped charges in recombination centers, respectively, and ε HM = Real[ε HM (ω = 0,f,γ,R̅ ,R min ,R max )] is the homogenized dielectric constant of active layer:NPs composite, which shows that the plasmonic NPs affect clearly the electrical modeling of plasmonic BHJ PSCs through ε HM . The right-hand side of Eqs. (6) and (7) describes generation and recombination processes, where G CT is the amount of CTEs generated in the active layer which is considered as equal to the G F , i.e., the conversion efficiency of FEs to CTEs is considered to be unity. G F is equivalent to useful absorption, i.e., the portion of incident photons of sunlight absorbed by the active layer. To calculate G F , the portion of parasitic photons, including absorbed photons by the non-active layers and scattered photons escaping from the PSC in all directions, is subtracted from the total number of incident photons to the plasmonic PSC. For this purpose, transfer matrix formalism that considers all optical interference effects is implemented to calculate the attenuation in each layer and the transmission and reflection at each interface layer of the plasmonic PSC, with the inputs of thickness and complex refractive index (n͂ (ω) = η(ω) + iκ(ω)) of each layer [127][128][129] . Therefore, the effect of NPs on the G CT is through the complex refractive index of the active layer:NPs composite, n͂ HM , defined as: It is noted that the portion of absorbed photon by NPs does not contribute to creating FEs and, therefore, is parasitic, but it is not taken into account in the semi-analytical optoelectronic modeling because Morawiec et al. 130 have shown that it is insignificant in the visible part of the AM1.5G spectrum. P CT→e-h in Eqs. (6) and (7) is the probability of dissociation from CTEs to free electrons and holes defined by: where k F is the rate constant of the decaying of CTEs to FEs, and k D is the rate constant of CTEs separation to free electrons and holes. The analytical expression of k D reported by Mihailetchi et al. 131 is implemented in our optoelectronic modeling, which is defined as: where a and E b are the separation distance and binding energy of bound electron-hole pairs; respectively, k B is the Boltzmann constant, T is temperature, J 1 is the first order of Bessel function, and ∇φ is the electric field strength in the active layer:NPs composite. As seen in Eq. (10b), the impact of embedded NPs on k D and consequently on P CT→e-h is through the homogenized dielectric constant and electric potential of the active layer:NPs composite. The second term on the right-hand side of Eqs. (6) and (7) describes the recombination process where R Lan and R trap stand for the Langevin bimolecular and trap-assisted monomolecular recombination, respectively 125,132 . The recombination of two free opposite charges created from different CTEs refers to Langevin bimolecular, and the recombination of a free charge with an immobilized charge at a trap state refers to monomolecular, where the first is defined as 133 : www.nature.com/scientificreports/ where C Lan = q(μ n + μ p )/ε HM , stands for the recombination coefficient predicted by Langevin model, ξ is an additional reduced factor taking into account experimentally derived reduced Langevin factor [134][135][136][137] , and n 1 and p 1 are the characteristic electron and hole concentrations, respectively, the product of which is equal to the square of intrinsic carrier density, i.e., n 1 p 1 = n i 2138 . Trap-assisted monomolecular recombination is described by a modification of Shockley-Read-Hall rate equation (r SRH (E)) 139 in which Gaussian density of state is considered for recombination centers (DOS RC (E)), as follows 140 : where the first curly bracket in the integral refers to DOS RC (E), E RC is the center of Gaussian distribution of recombination centers considered in the middle of the band gap, δ RC is the width of Gaussian distribution, the second curly bracket in the integral refers to r SRH (E), N RC is the total density of recombination centers including defects, impurities, and NPs, τ n (τ p ) is electron (hole) lifetime, and C trap stands for trap-assisted recombination coefficient. As reported by Wu et al. 132 , incorporating NPs into the active layer causes the increase in the recombination centers at the interfacial region of donor:acceptor. Therefore, in addition to defects and impurities density of donor:acceptor blend, the density of NPs is included in N RC 89 . Consequently, embedded NPs also affect τ n and τ p because these are inversely proportional to N RC 141 . It is to be noted that the incorporation of NPs in the PSCs changes the number of photo-generated carriers. As a result, both recombination processes in plasmonic PSCs differ from NPs-free counterpart PSCs. In addition to the number of photo-generated carriers, Eqs. (11) and (12) show that R Lan and R trap are respectively influenced by embedded NPs through ε HM and N RC .
Embedded NPs also influence on the Poisson equation, Eq. (8), through trapped charges in recombination centers, n RC , calculated as: where the curly bracket in the integral refers to the possibility of a recombination center being occupied by one electron.
To obtain a unique solution to Eqs. (6) to (8), it is necessary to specify appropriate boundary conditions 126,142 . They are defined for n(z), p(z), and φ(z), by assuming that the contact of active layer with ABL is hole ohmic and with CBL or cathode is electron ohmic, at z = 0 and z = t active as: where V is the applied bias, WF C and WF A are the work functions of cathode and anode, respectively, E gap is the effective energy band gap of active layer defined by the difference of the HOMO energy level of the donor polymer (E HOMO-do ) and the LUMO energy level of the acceptor molecule (E LUMO-ac ), and N c(v) is the effective density of states for electrons (holes).

Evaluation of the validity of semi-analytical optoelectronic modeling
Choosing a fabricated plasmonic BHJ PSCs as a test case. Plasmonic BHJ PSCs based on blending of poly-(3-hexylthiphene) and phenyl-C 61 -butyric acid methyl ester, P3HT:PCBM, with the incorporation of metallic NPs with different volume fractions, sizes, and shapes have been fabricated and extensively studied 20,40,147,148 .
In the following, we will focus on the conventional structure of P3HT:PCBM PSC, reported by Stratakis and Kymakis groups 22,28,149 , and P3HT:PCBM:NPs PSC with following geometrical parameters and materials: the weight ratio of P3HT:PCBM is 1:1, t active = 100 nm, the ABL is poly (3,4- Comparison to experimental data and discussion of modeling results. The validation of our semi-analytical optoelectronic modeling is investigated by Fig. 2. At the first stage, the pristine PSC, the ITO/ (11) R Lan = ξ C Lan np − n 1 p 1 = ξ q µ n + µ p ε HM np − n 1 p 1 www.nature.com/scientificreports/ PEDOT:PSS/P3HT:PCBM/Al structure without metal NPs, is simulated with the parameters indicated in Table 1. As can be found in the experimental literature, the aspects of the fabrication process of PSCs lead to a range of possible values for a parameter. From these ranges, the values of parameters in Table 1 are chosen from literature by comparing the experimental J-V curve with the modeling results to reach the best fitting. At the second stage, the experimental data for J-V characteristics of P3HT:PCBM:NPs PSC with 5% Au NPs concentration 22 , reported by Paci et al. 22 and shown with pink solid circles in Fig. 2, are compared with the simulated J-V characteristics obtained by our semi-analytical modeling. Since ref. 22 has reported that the NPs diameters with an average of 10 nm are distributed in the range of 1.5 to 20 nm, R̅ = 5 nm, R min = 0.75 nm, R max = 10 nm, and γ = 5 nm are considered for the optical modeling of P3HT:PCBM:NPs layer. The good agreement of simulated photovoltaic parameters of P3HT:PCBM:NPs PSC with the experimental ones shown in Table 2 confirms the validity of the semi-analytical optoelectronic modeling to examine the performance of plasmonic PSCs. It can be found from Fig. 2 and Table 2 that incorporating Au NPs into the active layer induces an improvement of PCE by 39%, which is the result of increase in the J sc and FF, while the V oc remains unchanged. To understand the origin of improved J sc , the absorption coefficient of the composite active layer of P3HT:PCBM:NPs, calculated by the HT, is shown in Fig. 3a. For comparison, the figure also shows the absorption coefficient of P3HT:PCBM. It can be seen that incorporating NPs improves the optical absorption of the active layer, resulting from the excitation of LSPR modes and multiple light scattering by the NPs simultaneously. Besides, the wide range of NPs diameter (from 1.5 to 20 nm) leads to broadband optical absorption enhancement, originating from the red shift of resonance peak of LSPR for large size NPs in size dispersion due to retardation effect and the blue shift of resonance peak for small size NPs due to intrinsic confinement effect. Another reason for improving J sc is the increase in the absorbed photons percentage in the active layer (useful absorption) of plasmonic PSC. To illustrate it, the percentage of absorbed photons in each layer and the percentage of reflection, including the interference induced by electrodes and scattered light escaping from the PSC, for the P3HT:PCBM and P3HT:PCBM:NPs PSCs, calculated by transfer matrix method, are shown in Fig. 4. It is found that after incorporating Au NPs in the active layer, the percentage of reflection decreases, and the percentage of absorbed photon in  www.nature.com/scientificreports/ the active layer increases for wavelengths above 650 nm. It should be noted that the incorporation of Au NPs also affects recombination processes. The presence of Au NPs in the active layer leads to increased trap density and consequently higher trap-assisted recombination, and strong local field around NPs results in a higher density of photo-generated excitons, G F , thereby leading to increased bimolecular recombination, which both recombination processes deplete the photo-generated carriers. However, the increase in the photo-generated carriers    In the P3HT:PCBM:NPs PSC, the V oc increase caused by the increase in the G F due to the LSPR effect of embedded NPs are counteracted with the V oc decrease caused by the increase in the C trap due to the addition of embedded NPs to the recombination centers.
FF is pronouncedly influenced by the electrical resistivity of active layer and electrodes 152 . The inverse of electrical resistivity, electrical conductivity (σ), can be calculated by the real part of complex refractive index (η) and absorption coefficient (α) as follows 153 : where c 0 is the velocity of light in vacuum and ε 0 is the vacuum permittivity.
The refractive index of the composite active layer of P3HT:PCBM:NPs calculated by the HT (η HM (ω)) is shown in Fig. 3b. It is found from Eq. (16) and Fig. 3a,b that the conductivity of P3HT:PCBM:NPs is more than the P3HT:PCBM. Increasing the conductivity is beneficial to carrier transport, leading to the increased FF of the P3HT:PCBM:NPs PSC.
The influence of the concentration (volume fraction) of Au NPs on the J-V characteristics of the P3HT:PCBM:NPs PSCs is simulated and depicted in Fig. 5, and the calculated photovoltaic parameters are listed in Table 2. Low volume fractions of NPs, varying from 0.04 to 0.07, are considered, because incorporating NPs with high volume fraction may deteriorate the morphology of donor:acceptor blend of active layer, and may create a short circuit leading to the degradation of the electrical properties of P3HT:PCBM:NPs PSCs.
With increasing the volume fraction (f) of NPs in the P3HT:PCBM:NPs active layer, no variation in V oc is observed from Fig. 5, but J sc increases. On the one hand, increasing f typically increases the number of NPs in the active layer and, consequently, increases recombination centers and, as a result, decreases the number of photo-generated carriers. On the other hand, increasing f improves the effective absorption coefficient of the active layer, as shown in Fig. 6 (right axis), and decreases the percentage of parasitic absorption (the sum of reflection and non-active layers absorption), as shown in Fig. 6 (left axis), leading to an increase in the number of photo-generated carriers. Photo-generated carrier enhancement due to the last two reasons outweighs their loss through recombination, which leads to an increase in J sc with increasing f. Figure 5 and Table 2 also show the experimental data of J-V characteristics and respective photovoltaic parameters obtained by Spyropoulos et al. 149 . It is seen that for the P3HT:PCBM:4%NPs and P3HT:PCBM:5%NPs PSCs, there is good agreement between calculated and experimental data, but not for the P3HT:PCBM:6%NPs PSC. The poor performance of P3HT:PCBM:6%NPs PSC reported by Spyropoulos et al. 149 is due to poor properties of NPs dispersion into the active layer resulting in NPs aggregation, affecting the plasmonic effects.
To investigate the effect of size dispersion on the performance of plasmonic PSCs, the J-V characteristics of P3HT:PCBM:5%NPs PSCs for Au NPs of 10 nm in mean diameter (2R̅ = 10 nm) and various radius dispersions (γ) are shown in Fig. 7. The blue, black, and red solid J-V curves are respectively for the NPs with equal diameters of 10 nm (γ = 0), diameters in the 5 to 15 nm range (2γ = 5 nm), and in the 1.5-20 nm range (2γ = 10 nm). With increasing γ, V oc is almost constant, and FF slightly decreases from 64.57 to 64.52%. An improvement in the www.nature.com/scientificreports/ J sc and PCE by 2.8% and 2.6%, respectively, for 2γ = 5 nm and by 12.9% and 12.7%, respectively, for 2γ = 10 nm compared to the γ = 0 is achieved. J sc enhancement is due to the increase in the absorption coefficient of the active layer, calculated by HT and shown in the right axis of Fig. 8a, and the decrease in the percentage of parasitic photons, shown in the left axis of Fig. 8a. Figure 7 also shows the J-V characteristics of P3HT:PCBM:5%Ag NPs PSC (the pink solid circles). The parameters of incorporated Ag NPs are chosen to be the same as incorporated Au NPs reported in refs. 22,149 (f = 0.05, 2γ = 10 nm, 2R̅ = 10 nm, 2R min = 1.5 nm, and 2R max = 20 nm). The P3HT:PCBM PSC incorporated with Ag NPs exhibits higher PCE of 4.03% compared to the P3HT:PCBM:Au NPs PSC, which is 3.81%. PCE enhancement is due to improved J sc , from 9.87 mA/cm 2 for Au NPs to 10.41 mA/cm 2 for Ag NPs, implying that the improved photocurrent results from enhanced absorption coefficient of the active layer in the wavelength range of 350-550 nm and the decrease in the percentage of parasitic photons below 480 nm for the P3HT:PCBM:Ag NPs PSC compared to the P3HT:PCBM:Au NPs PSC (see Fig. 8 (b)).

Proposing a high efficiency plasmonic BHJ PSC
As can be found from the modeling results and the experimental data of previous section, the efficiency of the plasmonic P3HT:PCBM PSCs is very low which is due to the much poorer performance of the reference P3HT:PCBM PSC, mainly resulted from intrinsic shortcomings of fullerene PCBM acceptor. The best reported PCE of single-junction PSCs with fullerene derivative acceptors, certified by National Renewable Energy www.nature.com/scientificreports/ Laboratory (NREL), is 11.5% 154 . Therefore, the designing of high-performance PSCs based on non-fullerene acceptors has attracted tremendous efforts in the recent years [155][156][157][158][159][160][161][162][163][164][165] . These efforts along with developing polymer, optimizing several aspects of BHJ morphology, and interface engineering have not only promoted the PCE of the non-fullerene PSCs to a high level of 17.23% (NREL-certified value 16.77%), but have also improved stability compared to fullerene PSCs 165 . The BHJ PSC with the highest efficiency so far, reported by Li group 165 , is a conventional structure with the layers of ITO/PEDOT:PSS/PM6:Y6/PDINN/Ag, where conjugated polymer PM6 and Y6 molecule of PM6:Y6 blend are applied as p-type donor and acceptor, respectively, and aliphatic amine-functionalized perylene-diimide (PDINN)/Ag as a bilayer cathode. Based on experimental reports in the literature and our modeling results in the previous section, we foresee that, by means of incorporating plasmonic NPs into PM6:Y6 active layer, it is possible to achieve even higher efficiency than 17.23%. Therefore, our semi-analytical modeling, providing a realistic prediction, is employed to investigate the performance of PM6:Y6:NPs PSCs. First, our modeling results have been fitted in Fig. 9 with the experimental J-V characteristics of reference PM6:Y6 PSC, reported by Li group 165 . In the electrical modeling, we set WF C = 3.72 eV 165 , which is the Ag WF modified by PDINN. Electron and hole mobilities of PM6:Y6 active layer with the weight ratio of 1:1.2 and thickness of 150 nm, determined by the space charge limited current method and reported in ref. 164 , are 5.90 × 10 -4 cm 2 V −1 s −1 and 2.00 × 10 -4 cm 2 V −1 s −1 , respectively. Second, through the HM calculation of the absorption spectrum of PM6:Y6 film embedding Ag, Au, Al, and Cu NPs along with the size optimization of the NPs, the most appropriate metal with optimized size for incorporating into PM6:Y6 blend is found. Third, the J-V characteristics of PM6:Y6:NPs PSC with the optimized conditions of NPs, which is Ag NPs with mean size of 20 nm ranged from 5 to 35 nm, is calculated www.nature.com/scientificreports/ with our semi-analytical modeling and shown in Fig. 9. From this figure and the photovoltaic parameters summarized in  Fig. 10a) and the increased percentage of absorbed photons in the active layer (see Fig. 10b), prevail the negative influence, which is increased trap-assisted recombination due to increase in trapping states emanating from the embedded Ag NPs.  www.nature.com/scientificreports/ The increased FF of PM6:Y6:NPs PSC most likely arises from the improved conductivity of the PM6:Y6:NPs compared to the PM6:Y6, as can be found from Eq. (16) and Fig. 10a, leading to the better electron transport within the active layer in the presence of Ag NPs 114,166 .

Conclusion
In conclusion, a semi-analytical optoelectronic modeling that could predict the performance of plasmonic BHJ PSCs where spherical NPs were incorporated into the active layer was demonstrated. Firstly, the effect of incorporation of NPs into the active layer on the optical properties was analytically modeled by the homogenization theory, considering a disordered array of NPs with size dispersion. Secondly, the percentage of useful absorption by the active layer was calculated by transfer matrix method, in which the number of photons related to absorption in the non-active layers, interference induced by electrodes, and scattered light escaping from the PSC in all directions were subtracted from the total absorbed photons in the PSC. Finally, J-V characteristics of plasmonic PSCs were modeled by coupled Poisson, continuity, and drift-diffusion equations. Then, by comparing the results obtained by the semi-analytical modeling with the experimental data reported in the literature for the photovoltaic parameters of P3HT:PCBM:NPs PSCs, which showed a good agreement, our modeling approach was verified. Therefore, for realistic prediction of the photovoltaic parameters of a new high efficiency plasmonic PM6:Y6 PSC, the modeling was applied and yielded that incorporating Ag NPs into PM6:Y6 active layer led to 10.9% improvement in the PCE, from 17 to 18.86%.

Data availability
The source code of this work will be made available from the corresponding author upon reasonable request.