Evolutionary optimization of light-matter coupling in open plasmonic cavities

The recent advent of polaritonic devices [optoelectronic devices based on exciton-polaritons (EPs) arising from the strong coupling of excitons in semiconductors and photons in optical cavities] is attracting significant interest due to the remarkable properties of EP states. The large binding energy and the associated high transition dipole moments of Frenkel excitons in organic semiconductors lead to many intriguing phenomena emerging from strong light–matter coupling at room temperature. The hybrid character of organic EPs makes them exhibit photon and exciton properties that result in the modification of intrinsic semiconductor properties,^1–3 such as the enhancement of the exciton diffusion length,^4–9 increase of the long-range energy transfer,^10–12 or reduction of photobleaching.^13,14 Enhanced non-linearities in EPs are also responsible for polariton lasing^15–18 and quantum condensation due to the bosonic character of EPs,^19–23 which can lead to low threshold coherent emission.

Organic photovoltaic (OPV) devices can profit from strong light–matter coupling and EPs. OPV cells are easy to fabricate, lightweight, and mechanically flexible, and can be partially transparent, which makes them very appealing for large-area energy conversion. However, the power conversion efficiency of OPV cells is low and limited by the thickness of the active layer, which is determined by exciton diffusion lengths and carrier recombination losses.^24,25 Other limiting factors are the charge transfer across the interface separating the donor–acceptor materials and the photo-degradation due to intersystem crossing.²⁶ The modified energies of EPs, compared to excitons, could be used to reduce the triplet density or enhance reversed intersystem crossing,^11,27–29 modify the excited-state life time,^13,30–32 or enhance exciton and charge transport.^{4–9,33–35} These phenomena could be exploited to solve the problems limiting a more extended utilization of OPV, and the first steps in this direction have been made.^36–38 However, most of the works on strong light–matter coupling are based on closed Fabry–Perot microcavities with the active semiconductor enclosed between two mirrors.³⁹ Such a closed geometry limits the application of strong light–matter coupling in OPV, where incident light needs to be absorbed as efficiently as possible in the active material.

In this manuscript, we investigate the coupling of excitons in regioregular poly(3-hexylthiophene-2,5-diyl) (P3HT) semiconducting polymer with optical modes in resonant cavities formed by square arrays of plasmonic nanodisks made of aluminum (Al). P3HT is an electron donor polymer, widely employed in OPV. Nanoparticle arrays define open optical cavities with optical modes that result from the radiative coupling of localized surface plasmon resonances to diffraction orders in the plane of the array. These optical modes are known as Surface Lattice Resonances (SLRs).⁴⁰ SLRs can have remarkably narrow linewidths (high Q-factors) due to low radiation losses and high field enhancements over large volumes.^41–44 When the nanoparticle array is covered by a P3HT layer, SLRs can couple to excitons in the organic semiconductor as evidenced by the formation of lower and upper EP states in the absorption spectrum with an associated Rabi splitting.⁴⁵ We have used a Particle Swarm Optimization (PSO) algorithm to find the global best sample dimensions for light–matter coupling.^46–49 We run the PSO algorithm together with Finite-Difference in Time-Domain (FDTD) simulations using 30 stochastic populations with randomly chosen sample parameters that evolve during 100 generations.⁵⁰ Samples have been fabricated with the optimized dimensions, and the optical extinction dispersion has been measured, confirming the results obtained by the PSO algorithm.

A schematic representation of the investigated system is depicted in Fig. 1(a). This system consists of a square array of aluminum (Al) nanoparticles with a layer of P3HT on top of a glass substrate (refractive index n = 1.51). The choice of Al is motivated by the plasmonic response of this metal at short wavelengths in the visible spectrum and its stability. The thickness of the P3HT layer is 100 nm, comparable to the typical thickness of the active layer in organic solar cells. The Al nanoparticles, with diameter d and height h, are arranged in an infinite periodic square array with a lattice constant a. The nanoparticle array is designed to support SLRs at the P3HT exciton energy for normal incidence, i.e., light with a wavelength close to 500 nm and that incident along the normal direction to the surface of the arrays is diffracted in the plane of the array. The optical properties of P3HT, defined by the refractive index and extinction coefficient, are obtained from ellipsometry measurements and fitted to a Lorentz oscillator model that describes the excitonic response, as shown in Fig. 1(b). (The values of the fit can found in Sec. S1 of the supplementary material.) From this fit, we retrieve an exciton energy of 2.38 eV, corresponding to a wavelength of 521 nm and a background refractive index (refractive index far from the exciton energy) of 1.72.

Two questions are addressed in this manuscript: (i) What is the maximum Rabi splitting that can be achieved by coupling of P3HT excitons and SLRs supported by a square array of Al nanoparticles? (ii) Which parameters of the array lead to this maximum Rabi splitting? We seek the answers to these questions using a PSO algorithm and FDTD. PSO algorithms are designed for numerical optimization problems that look for the solution by defining a random population of candidate solutions (“particles”) that can move (or evolve) over the search space by changing position and velocity. Each “particle” movement is influenced by its own local best position and the best known positions of the other “particles”. Subsequent generations converge to a global best solution. A more detailed description of the PSO method can be found in Sec. S2 of the supplementary material.

The properties of SLRs in the Al nanoparticle arrays are controlled by the size of the nanoparticles and the lattice constant. We tune the geometrical parameters, including the height h, diameter d of nanodisk, and the lattice constant of the square array a, independently, to find the lattice with a SLR absorption that matches spectrally and spatially the absorption spectrum and the location of the P3HT layer. The fitness value that we maximize with the PSO algorithm is defined by

where P_abs is given below and corresponds to the overlap integral of the absorption spectra of the bare (uncoupled cavity and exciton states) leading to the coupled system and γ_c and γ_p are the decay rates of these bare states. The fitness value is inversely proportional to the modulus of the difference between the loss rates of the bare states, which is justified by the fact that in a coupled oscillator model, the Rabi splitting is maximum when γ_c = γ_p (see Sec. S3 of the supplementary material).

The systems supporting the bare states are schematically represented in Fig. 2(a) and correspond to (i) SLRs in an empty cavity, i.e., a nanoparticle array on glass with a 100 nm thick layer of a material with a refractive index of n = 1.72 + 0.0001i on top, similar to the refractive index of P3HT at energies far from the exciton energy and (ii) excitons in a P3HT layer with a thickness of 100 nm on top of glass. These bare states have been simulated using FDTD by determining the absorbed power spectrum in the volume occupied by the layers and are given by $Pac(ω)$ and $Pap(ω)$ for the bare cavity and bare P3HT layer, respectively. The normalized absorbed power spectra are displayed in Fig. 2(b), where we use h = 49 nm, d = 130 nm, and a = 280 nm. As shown in Sec. III, these values correspond to the array with the largest coupling strength. The absorbed powers are defined by integrating over the volume occupied by the layers in a unit cell of the array: $Pac,p(x,y,z,ω)=(1/2)ϵ0ω|Ec,p(x,y,z,ω)|2εc,p′′(ω)$⁠, where E_c,p(x, y, z, ω) is the electric field at frequency ω and position (x, y, z) in the layer on top of the array and in P3HT and $εc,p′′(ω)$ is the imaginary component of the dielectric function of the layers. Note that we allow for a sufficiently small imaginary component of the refractive index in the layer on top of the Al nanoparticle array to calculate an absorbed optical power in this layer without introducing a significant change to the SLR mode profile. The absorbed power peak at 533 nm corresponds to the SLR emerging from the (1,0) and (−1,0) diffraction orders, while the narrower peak around 428 nm is due to the (0,±1) diffraction orders. P_abs in Eq. (1) is given by

where V is the volume of the polymer layer in a unit cell of the array. The fitness value given by Eq. (1) is thus proportional to the overlap integral of the absorbed power of the bare systems calculated in the volume occupied by the polymer. Correspondingly, we consider the regions in the sample where the SLR mode overlaps with the excitons in P3HT. The choice of this fitness value ensures that both the absorption of the bare states and their overlap is maximum, which should also lead to a maximum coupling strength and Rabi splitting.

We also show with the dashed curve in Fig. 2(b) the normalized scattering efficiency of a single Al nanoparticle with the same dimensions as the nanoparticles in the array to illustrate the significant narrowing of the SLR as a result of the collective response of the array. The tunability of the linewidth of the SLR, which is controlled by the localized plasmonic resonances in the individual nanoparticles and the diffraction orders in the plane of the array,^51–53 turns out to be critical to achieve the largest possible coupling strength and Rabi splitting.

Figure 3 shows simulations of the electric field intensity enhancement (local field intensity normalized by the incident field intensity) in the bare cavity at 533 nm, i.e., the SLR peak wavelength (simulations of the field intensity enhancement at 428 nm corresponding to the other peak in the absorption spectrum are shown in Sec. S4 of the supplementary material). These simulations are performed by considering a plane wave illumination at normal incidence and polarization along the x-direction. Figure 3(a) corresponds to the x–y plane at the top of the nanoparticles, while (b) and (c) are the cross sections in the z–x and z–y planes indicated by the dashed lines in (a). We can observe in these figures a high field intensity enhancement close to the nanoparticles due to the localized surface plasmon resonance (note that to highlight the SLRs, the color scale is saturated at an intensity enhancement of 10) and a field intensity enhancement extending over the array and in the thin layer on top as a result of the (1,0) and (−1,0) in-plane diffraction orders.

It is worth mentioning that it is not possible to define the fitness value as the energy difference between the lower polariton (LP) and the upper polariton (UP) in the coupled system. The Rabi splitting is given by the minimum splitting between the LP and the UP, which is defined at wave vectors in the dispersion relation of zero detuning for the bare states. A PSO algorithm based on a fitness value defined for the coupled system would not know a priori the energy and dispersion of the bare states, converging to a wrong set of parameters. This point is discussed in more detail in Sec. III, and it is the reason why we chose a fitness value that is defined for the bare states.

Some constraints are imposed to the geometrical values and the fitness value calculation by the PSO algorithm: (a) a > d + 50 nm. This condition is imposed to ensure that the nanoparticles do not overlap and are sufficiently separated to be possible to fabricate them using electron-beam lithography. (b) F = 0 when the wavelength difference between the SLR and the P3HT center exciton energy is Δλ > 50, 40, 30, and 20 nm in the 2–10, 11–20, 21–30, and 31–100 generations, respectively. This condition is imposed to ensure that the populations evolve to a solution in which the fitness value is maximum for overlapping resonances. (c) $F=2Pabs0.001$ if |γ_c − γ_p| < 0.001. This condition is imposed to avoid the divergence of F when γ_c = γ_p.

The PSO algorithm was run using 30 initial random populations of Al nanodisk arrays. The values of h, d, and a were bounded to the ranges 10–100, 50–200, and 200–400 nm, respectively. These ranges of values could be eventually realized experimentally.

The largest fitness value in the first generation was 3.73 × 10², and after 100 iterations, the fitness value increased about 45 times to obtain a global best F = 1.69 × 10⁴, as can be seen in Fig. 4(a) with a black line. In less than 50 generations, the PSO algorithm was able to find a geometry with dimensions that increase F to values close to convergence, showing the effectiveness of this method. The mean value of F of all populations is shown in Fig. 4(a) with a red line. This mean value approaches the global best value as the number of generations increases, as it could be expected when this global best attracts the other populations.^46,50 Figure 4(b) displays the fitness values (in a log-color scale) of each population and for all the generations, showing the increase in F for all the populations as they converge toward the best geometry.

Figure 5(a) shows the normalized spectra of $Pac(ω)$ and $Pap(ω)$ for the 11 best populations during the full range of generations (the number gives the generation at which F increases, as it is also indicated in Fig. 4), while Fig. 5(b) displays the simulations of the optical power absorption enhancement calculated for the coupled systems with the same structural parameters as these best populations. The power absorption enhancement is defined as the power absorption in the P3HT layer of the coupled system, normalized by the power absorption of the bare P3HT layer to suppress the contribution of dark and uncoupled excitons. The best populations for the first nine generations correspond to detuned SLR–exciton systems, as can be appreciated by the large wavelength difference between the P3HT absorption maximum [black dotted curves in Fig. 5(a)] and the maxima in absorption spectra of the bare particle array corresponding to SLRs [red solid curves in Fig. 5(a)]. As a result of this detuning, indicated by δ in Fig. 5(a), the difference between the LP and the UP shown in Fig. 5(b) does not correspond to the Rabi splitting. Only when the detuning between the bare states becomes sufficiently small, we can identify the maxima in the absorption spectra of the coupled system with the LP and UP states for zero detuning, defining the Rabi splitting. The Rabi splitting reaches a maximum of 0.7996 eV after 85 generations. This value is larger than the decay rates of the bare states obtained from the FWHM of the power absorption spectra, γ_c = 0.5874 eV and γ_p = 0.5881 eV, and satisfies the criteria for strong coupling of coupled oscillators given in Sec. S3 of the supplementary material. The Rabi splitting corresponds to 34% of the exciton energy, which places the system at the onset of the ultra-strong coupling regime.⁵⁴ 

In PSO algorithms, the value of F depends on several parameters that gradually converge toward those defining the global best solution as F increases. To illustrate the convergence of the parameters in the optimization process, we show in Figs. 6(a)–6(c) the mean (black curve) and standard deviation (SD) (red curve) of h (particle height), d (particle diameter), and a (lattice constant) for each generation, respectively. After 100 generations, we obtain the maximum value of F corresponding to an Al nanoparticle array with h = 49 ± 8, d = 130 ± 16, and a = 279 ± 10 nm. The large particle diameter that is retrieved after the optimization leads to a broad SLR and an optimum overlap with the absorption spectrum of P3HT. From these values, we can see that the SD of the lattice constant is less than 4% of the mean, being the relative values of the SD of the diameter and height of the nanoparticles 12% and 16% of their respective means. This low relative SD suggests that the lattice constant is the most sensitive parameter for the optimization of the coupling strength and the one converging the fastest, which can be understood by the strong dependence of the linewidth of SLRs with the lattice constant in nanoparticle arrays.⁵⁵ 

The relation between the Rabi energy and the power absorption integral defined in Eq. (2) is illustrated in Fig. 7(a) for several populations in different generations, where we can see that Ω_R increases first linearly with P_abs, converging for large values of $ΩR$ to ∼0.8 eV. As expected, the last generation corresponds to the largest Rabi splittings as marked by the red ellipse in Fig. 7(a). The global best population is indicated with the red data point.

Based on the results of the PSO algorithm, we have fabricated a sample with structural dimensions as close as possible to the optimal dimensions for coupling with excitons in P3HT. The sample was fabricated using electron-beam lithography (EBL) and reactive-ion etching (RIE). A 50 nm thick polycrystalline Al layer was deposited on a glass substrate (Eagle XG, Corning) by an electron-beam deposition at room temperature. A resist (nega-type, NEB22A2, Sumitomo, Japan) was spin-coated on the Al film, and the array patterns were written by EBL (F7000S-KYT01, Advantest, Japan). After the development of the patterns by a developer solution (SD-1), the Al film with a patterned resist was etched using RIE. After that, the remaining resist on top of the array was removed by an ashing process where oxygen plasma selectively etches the resist away. There is a 3 nm thick natural oxide layer (Al₂O₃) formed on top of the nanoparticles that gives the long-term stability needed for applications. An optical microscope image of the 200 × 200 μm² array is shown in Fig. 8(a), while a scanning electron microscope image of the nanoparticles with a = 280 ± 10 and d = 130 ± 10 nm is displayed in Fig. 8(b).

The sample has been characterized by measuring the dispersion of the optical extinction using Fourier plane imaging spectroscopy.⁵⁶ For these measurements, a low-power white-light beam from a tungsten–halogen lamp is focused onto the sample with a 40× objective (NA = 0.6). The transmitted light through the sample is collected by another 60× objective (NA = 0.7), and the back focal plane of this objective is imaged at the entrance slit of a spectrometer. The imaging spectrometer retrieves the transmission spectra as a function of angle. These measurements can be plotted in a dispersion diagram, representing the optical extinction as a function of the photon energy and the wave vector parallel to the surface, i.e., k_‖ = (2π/λ)sin θ. The related extinction is defined as 1 − T/T₀, where T is the transmission through the sample and T₀ is a reference transmission measured through bare P3HT without the sample.

The extinction measurements are shown in Fig. 8(c) and compared with the simulations obtained with FDTD in (d). These measurements and simulations correspond to the bare optical cavity, which has been fabricated by spin-coating a film of P3HT with a thickness of 100 nm on top of the Al nanoparticle array that has been photobleached using the expanded beam from a continuous wave (CW) laser of λ = 445 nm. The photobleaching is measured in a pump–probe setup where a thin film of P3HT is pumped non-resonantly with a 445 nm CW laser and probed with a white-light source. The measured spectra over time are used to characterize the photostability of the sample. During the photobleaching, P3HT is photo-excited into singlet states that can relax in two different ways to the ground state: direct relaxation or triplet states. The triplet states can react with ambient triplet oxygen, creating products such as singlet oxygen and hydrogen peroxide, which bleach the polymer and reduce the absorption. By photobleaching, we remove the excitonic resonance in the P3HT layer while maintaining a thin dielectric layer on top of the Al nanoparticle array to reproduce experimentally the bare cavity shown in Fig. 3(a). We note, however, that the refractive index of the bleached P3HT layer is ∼1.5 (measured with ellipsometry), thus lower than the background refractive index of P3HT. The SLRs of the array in Figs. 8(c) and 8(d) correspond to the bands of high extinction. We have also marked in these figures the in-plane non-degenerate (1,0) and (−1,0) and the degenerate (0,±1) diffraction orders satisfying the grating equation.⁵⁷ The in-plane diffraction of the (1,0) and (−1,0) orders is responsible of the low-energy SLR that can couple to the P3HT excitons at 2.35 eV. This SLR is indicated in the figures and calculated using a damped coupled oscillator model with one oscillator representing the localized surface plasmon polaritons of the individual nanoparticles and other two representing the (1,0) and (−1,0) orders,⁵⁸ neglecting the mutual coupling between the diffraction orders.

The coupling of the SLR to excitons in P3HT is demonstrated by the measurements and simulations shown in Figs. 8(e) and 8(f). In these measurements and simulations, which correspond to the nanoparticle array with the P3HT layer on top before photobleaching, we can see the formation of the LP band at lower energies than the energy of the SLR of the bare particle array. Note that for the simulations, we have used the complex refractive index of P3HT retrieved from ellipsometry and shown in Fig. 1 with the dashed curves. The SLR here is slightly red-shifted with respect to the SLR in the photobleached sample due to the higher refractive index of P3HT. The LP band has been fitted to the coupled oscillator model, described in Sec. S3 of the supplementary material and assuming a constant loss rate of the SLR equal to the value at k_‖ = 0, to retrieve a Rabi energy of Ω_R = 0.74 eV. This value of the Rabi splitting is smaller than the maximum value of 0.8 eV obtained from the PSO algorithm and the simulations, which can be explained by the fact that extinction measurements give the contribution of absorption and scattering to the spectrum, while the PSO algorithm was applied only to absorption. As discussed in the literature, absorption measurements provide a more reliable value of the Rabi splitting than reflection, transmission, or extinction measurements.^59–61 However, accurate absorption measurements need to be done using an integrating sphere to collect all the reflected and transmitted light intensities by the sample, which typically require a size of the sample larger than the one used in our experiments. We also note in Figs. 8(e) and 8(f) that the UP band calculated with the coupled oscillator model does not coincide with the band of large extinction over the whole range of wave vectors. This discrepancy can be attributed to the coupling to the (0,±1) SLRs that cannot be neglected at high energies. For completeness, we have also indicated with white dashed curves and lines in Figs. 8(e) and 8(f) the dispersion of the SLR in the bare cavity and the P3HT central exciton energy. The nearly flat band at 2.1 eV in Fig. 8(f) corresponds to 0–1 intrachain excitons in P3HT, as can be also appreciated in the experimentally determined values of the optical constants of Fig. 1(b). These excitons are not considered in the PSO calculations.

Using a particle swarm optimization algorithm and finite-difference in time-domain simulations, we have determined the optimal size and lattice constant of square arrays of Al nanoparticles for strong coupling of surface lattice resonances with excitons in P3HT films. Through 100 generations of 30 populations, we retrieve a maximum Rabi splitting of 0.8 eV in these open optical cavities. We have also fabricated the Al nanoparticle array with optimized geometrical parameters and measured the optical extinction using back focal plane imaging spectroscopy. The experimental results show the formation of exciton polaritons in the strong coupling limit, in agreement with the simulated extinction. This work illustrates an efficient approach to optimize resonant structures for light–matter interaction in open systems that can be used in polaritonic devices.

The supplementary material showing the fit to the optical constants of P3HT, the description of the PSO algorithm, the coupled oscillator model for strong light–matter coupling, and the electric field distribution of the high energy lattice resonance is available.

This work was financially supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) (Vici Grant No. 680-47-628). P.B. is sponsored by the China Scholarship Council. S.M. acknowledges financial support from MEXT, Japan (Kakenhi, Grant No 17KK0133).

The data that support the findings of this study are available from the corresponding author upon reasonable request.