FDTD modelling of nanostructured OLEDs: analysis of simulation parameters for accurate radiation patterns

Nanostructures in OLEDs allow to manipulate the radiation pattern and direct light of a specific wavelength into a specific angle in the far field by resonant light outcoupling. We present a theoretical study investigating simulation conditions for accurate radiation patterns employing the FDTD method. Radiation patterns of single dipole emitters placed in a sufficiently large simulation domain are superimposed to model the incoherent emission of a continuous emission layer. The necessary simulation domain size is derived by considering propagation length and field profile of the guided mode of a corresponding unstructured OLED. A 5 nm mesh resolution is determined as a good compromise between accuracy and simulation time. The position of an emitter with respect to the grating has a large influence on the radiation pattern. Nevertheless, we demonstrate for nine different OLED structures that the radiation pattern of a homogeneous emission layer is approximated with an amplitude error of less than 5% by superimposing only four equally spaced dipoles at positions ±(1/8,3/8,5/8,7/8)·Λ/2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pm \,(1/8, 3/8, 5/8, 7/8) \cdot \varLambda /2$$\end{document}, which avoid points of high symmetry.


Introduction
Organic light-emitting diodes (OLEDs) are light sources with a broad range of applications, including general lighting, displays, automotive, and compact optical sensors. One of their major advantages is the prospective possibility to fabricate OLEDs using large scale, cost-efficient roll-to-roll printing technologies. This enables new application areas such as 1 3 139 Page 2 of 20 single-use sensor devices for point-of-care diagnostics using organic optoelectronics (Søndergaard et al. 2013;Smith et al. 2016;Jahns et al. 2017).
OLEDs are thin-film devices made of multiple layers of organic and inorganic (semi-) conducting materials. The layers have thicknesses in the range of 10 nm to 200 nm and typical refractive indices of around 1.8, thereby forming a slab waveguide (Lu and Sturm 2002). From a typical OLED stack, only 20% of the generated light is extracted to the outside due to trapping in the waveguide or in the substrate and absorption inside the OLED (Meerheim et al. 2010). The radiation pattern of a planar OLED is approximately Lambertian (Kim et al. 2000). A photograph of a yellow-emitting OLED is shown in Fig. 1b. Nanostructures in photonic devices are periodic or aperiodic, deterministic or stochastic geometric structures with dimensions in the range of the light's wavelength. They manipulate the propagation of light and can be used for optical filters, sensors, or-in the case of OLEDs-to control the light outcoupling (Geyer et al. 2008;Schwab et al. 2014). Tailoring the broadband Lambertian emission radiation pattern such that the operating wavelength is mainly directed into one specific emission angle is particularly interesting for compact optical sensors. This tailoring is achieved by including a periodic modulation of the refractive index (i.e. a grating) into the OLED. Figure 1a shows a simplified view of an OLED stack with an 1D grating. The spontaneous emission is continuously distributed inside the emission layer (EML). Light travelling along the grating is outcoupled by the grating, resulting in the characteristic 'resonance cross' in the wavelength-dependent radiation pattern shown in Fig. 1c.
Integrating a nanostructure into an OLED offers many new degrees of freedom, which are difficult to explore purely experimentally. It is therefore reasonable to perform numerical simulations in order to study the influence of different parameters on the device performance. In this work we discuss the simulation of such nanostructured OLEDs using the FDTD method (Taflove and Hagness 2005) and give guidelines for a correct simulation setup. The direct simulation is especially challenging due to small feature sizes in the 10 nm regime, long propagation lengths up to hundreds of micrometres and incoherent and randomly distributed emitters, which break the periodicity of the structure.
The challenge of incoherent emitters has been addressed by using many equally (Lee et al. 2003;Kim et al. 2014) or randomly (Lee et al. 2005) distributed dipoles with random initial phases or by multiple runs with single dipoles at different positions which are incoherently superimposed afterwards (Callens et al. 2014). The long propagation lengths require long domain sizes (Callens et al. 2014). Due to memory limitations, these are not realizable in 3D simulations and the domain is truncated by perfect mirrors after a few micrometers (Lee et al. 2003(Lee et al. , 2005Kim et al. 2014). These memory limitation can be overcome by using the inverse Floquet transformation (Zschiedrich et al. 2013) or the reciprocity theorem (Janssen et al. 2010;Zhang et al. 2015). Apart from FDTD, the finite element method (FEM) and the rigorous coupled-wave analysis (RCWA) are commonly used to analyse nanostructured OLEDs. Paulsen et al. (2017) conducted a comparison of these three methods for photonic structures, showing comparable simulation results, but differences in the computation time. Figure 1e shows an exemplary grating OLED unit cell with a continuous emission layer, which has to be approximated by simulating the radiation of multiple dipoles (Fig. 1f).
As described in Lumerical's documentation, the FDTD method is inherently coherent (Lumerical Inc. 2018b). Therefore, exciting the structure with multiple dipoles simultaneously gives incorrect results. It is recommended to run the simulation for each dipole individually and to obtain the incoherent result by summing the intensities afterwards (Lumerical Inc. 2018e). It is also stated that if temporal incoherence is achieved by running a long simulation with multiple dipoles with random phase shifts, it is not possible to use the results for a nearto-far-field transformation (Lumerical Inc. 2018f).
The radiation pattern of each individual dipole depends on the dipole's position with respect to the grating (an example is shown in Fig. 1g). We therefore address the question of how to place the dipoles to emulate the continuous emission layer.
This work is structured as follows: In Sect. 2 the general simulation setup is presented. Then, in Sect. 3 we determine how much structural detail has to be included in the model by analysing two structures with different levels of detail. In Sect. 4 we discuss the mesh resolution and simulation domain size. The placement of dipoles is discussed in Sect. 5.

Simulation method and structure under investigation
We use the commercially available software FDTD Solutions by Lumerical Inc. for all near-field simulations (Lumerical Inc. 2018a). This work focuses on the outcoupling of light propagating inside the stack perpendicular to the grating lines. Therefore, the simulations are restricted to the 2D case. In this study, OLED stacks with a Super Yellow (PDY-132, Merck) emission layer are considered. The two main Super Yellow OLED stacks considered here are shown in Fig. 2, both as a planar structure and as a unit cell of a nanostructured version. Since Super Yellow is a long-chain molecule, its orientation will be mainly parallel to the layers in a real grating-OLED (McBranch et al. 1995). Therefore and because TE modes have less interaction with the absorbing metal electrode, only TE polarized dipoles are considered in this study. For all simulations, the emission zone is assumed to be a thin layer with a vertical distance of 20 nm from the upper Super Yellow boundary (see Fig. 1e).
The simulation region is surrounded by perfectly matched layers (PMLs) to suppress reflections from the domain boundaries. One side of the OLED stack is terminated by a thick ( ≥ 70 nm ) metallic electrode. The small transmission through this layer is neglected and the boundary condition is set to perfect electric conductor (PEC). On the opposite side, where the light is extracted through air or a glass substrate, a 16-layer PML with a predefined 'steep angle' profile, which enhances the absorption of light travelling nearly parallel to the boundary, is used. The mesh is set to a fixed size in the grating area, which is discussed later, and a gradual increase of the vertical mesh size is allowed in the relatively thick outermost layer.
The structure is excited by a broadband dipole source placed at different locations in the emission layer. This dipole source mimics the oscillating dipole moment that forms during a radiative state transition (Henderson 1979;Lumerical Inc. 2018d). Hence FDTD is a time-domain method, the full spectral response can be calculated from a single short pulse (Schneider 2010) and any desired spectral emitter characteristic may be applied afterwards.
All fields are collected in the near field using a horizontal field monitor close to the PML. From that data, the far-field radiation patterns are calculated using an in-house python implementation of a standard near-to-far-field transformation algorithm (Schneider 2010). The farfield is calculated in the outermost layer (encapsulation or substrate layer, , an electron transport layer (ETL) and a silver cathode (Ag). The substrate is transparent (e.g. glass). b The same stack with a nanostructure. The primary refractive index modulation occurs between the substrate and the ITO. c Stack of a planar Super Yellow consisting of a reflective cathode (silver, Ag), the emitter Super Yellow (SY), a hole transport layer (HTL), a transparent low-index polymer anode (PEDOT:PSS, n = 1.4 ) and an encapsulation layer (enc). d The same stack with a nanostructure. The primary refractive index modulation occurs between the metal and the Super Yellow and between the HTL and the anode respectively) and projected into the free space ( n = 1 ) using Snell's law and the Fresnel equations.

Required degree of detail
Fabricated nanostructured OLEDs may have sub-wavelength feature sizes ranging down to 10 nm or even lower. Furthermore, when fabricating such nanostructured OLEDs, the original structure geometry is not perfectly replicated into the OLED because corners are rounded and steep edges are flattened. In order to define the simulated structure with the necessary degree of detail and to obtain results that are comparable to results obtained from measuring fabricated devices, it is important to know whether small details in the geometry will influence the final results. It has been shown that the general effect of resonant light-outcoupling is mainly described be the grating's first-order Fourier coefficient (Kluge et al. 2014). Yet higher-order details also have a relevant influence on the results. We performed simulations of two similar structures, which are shown in Fig. 3, to investigate the influence of small details on the radiation pattern. The differences between the structures are that structure (a) represents a perfect grating replication with steep edges, constant layer thicknesses and grating heights through the layers, and structure (b) resembles more closely a fabricated device where the layer thicknesses are not uniform and the edge steepness decreases through the layers. During fabrication, a perfectly rectangular grating in the first layer may even deform into a sinusoidal shape in subsequent layers.
The results are shown in Fig. 4. Although the change in the structures is rather small, the change in the far-field spectrum, especially around 0 • , is significant. We conclude that it is necessary to include detailed information about the nanostructure geometry into the model to accurately predict radiation patterns of fabricated devices. Furthermore, it is necessary to use a sufficiently small mesh size which is capable of resolving such details.

Mesh and simulation domain size
The smallest and largest length scales in the model determine the mesh size and the simulation domain size. Both have direct influence on the calculation time, memory requirement, and the result's accuracy. In this section we present a convergence analysis to find a compromise between computation cost and accuracy.

Mesh resolution
With the FDTD method, it is not possible to freely assign different mesh resolutions to parts of the simulated structure. Using a finer mesh, e.g. for the grating structures, influences the mesh globally because the mesh elements have to maintain their rectangular shape. The rectangular mesh also makes it difficult to simulate sloped or curved boundaries. The FDTD solutions software provides a proprietary mesh refinement method (Lumerical Inc. 2018c) which is an extension of the Yu-Mittra method (Yu and Mittra 2001).
In order to find a compromise between accuracy and calculation time, we conducted a simple convergence study by varying the mesh size and calculating the radiation patterns for a single dipole. We varied the mesh size from 100 nm to 2 nm and used (1) no mesh refinement ('staircase'), (2) mesh refinement only for dielectric interfaces ('conformal variant 0'), (3) mesh refinement for all interfaces including metals ('conformal variant 1') (cf. Lumerical Inc. 2018c). The electric field's convergence behaviour at different positions is shown in Fig. 5.
Good convergence is reached for mesh resolutions around 10 nm and below if Lumerical's conformal variant 1 (method (3)) is used. The other methods do not provide fully converged results down to mesh sizes of 2 nm . As higher resolution would require unreasonably long simulation times (see Table 1), we did not further investigate the other methods' convergence behaviour systematically.
We conclude that a mesh resolution of 5 nm combined with the mesh refinement option 'conformal variant 1' provides good accuracy while maintaining reasonable calculation times.

Simulation domain size
The simulation domain size is determined by two properties of the quasi-guided mode: its propagation length and its evanescent field. Because we are analysing resonant outcoupling through the grating, the full interaction of the propagating mode with the grating must be simulated. The propagation length is determined by absorption losses and radiation losses. The radiation losses are not known beforehand, but the absorption may be estimated from the corresponding planar structure. Modes are calculated using e.g. the transfer-matrix method (Kwon 2009). Here, we choose a lateral simulation domain size such that the mode is attenuated by absorption to 1% of its initial intensity at the edge of the simulation domain, when excited at the center. This leads to the following requirement for the minimum domain size L: Additional scattering losses result in shorter propagation lengths and thus this is a worstcase scenario.
To illustrate the requirement for such long domain sizes, we simulated the OLED structure shown in Fig. 2d with different lateral domain sizes. The results are shown in Fig. 6. In order to simulate the resonant outcoupling correctly, the quasi-guided mode needs a sufficiently long interaction length with the grating.  From that L/2 can be calculated to 17 μm , 33 μm and 41 μm , respectively. Figure 5 shows that the relative errors fall below 1% around L∕2 ≈ 20 μm , 30 μm , and 40 μm , respectively, which corresponds well with the estimated propagation lengths (Due to the oscillations that can be seen in Fig. 6c, the relative errors do not decrease smoothly and it is not meaningful to determine exact crossing points with the dashed 1% line). We conclude that long lateral domain sizes are indeed necessary and that they can be estimated using the given equation. In a 3D simulation it is usually not possible to use such large domain sizes. In this case, short domain sizes of the order of 5-10 μm in combination with perfect mirror boundaries have been used (Lee et al. 2005;Kim et al. 2014). However, special care must be taken to ensure that cavity effects or coherent interference effects do not influence the final result, as Fig. 6e, f shows that naively replacing the PML boundaries with perfect mirrors leads to incorrect results.
The vertical extension of the simulation domain, i.e. the thickness of the most outer layer (substrate or cladding), must be large enough that the guided mode's evanescent field does not interfere with the horizontal near-field monitor or the PML to avoid numerical errors in the near-to-far-field transformation and distortions of the mode structure by the PML. The required thickness is estimated from the corresponding planar structure's guided mode's field profile, which is also calculated with the transfer-matrix method. The thickness is typically in the range of 500 nm to 1 μm . With these mesh and domain sizes and using 12 cores of a Intel Xeon E5 v3 or v4 processor equipped with 128 GB of RAM or more leads to simulation times on the order of tens of seconds to minutes per simulation.

Dipole positions
The light in an OLED is generated by spontaneous emission inside a thin film in the emissive material. The emitting molecules are continuously distributed and may-depending on the molecule-be randomly oriented or have a preferred orientation. In the simulation this emission is modelled by incoherent dipole sources. The emission zone position in the emission layer depends on the OLED's electrical characteristics. With the incorporation of a nanostructure, the OLED geometry no longer has a continuous translational invariance. To mimic the continuous emission film of a real OLED, the radiation of many uniformly distributed incoherent dipoles is superimposed. To achieve incoherency, the simulation is carried out for each dipole individually and the final incoherent result is obtained by where E 1 , … , E N are the individual fields calculated for N dipole positions (Lumerical Inc. 2018b). Figure 7 shows radiation patterns for dipoles located at different positions with respect to the grating. The radiation patterns are asymmetric and depend strongly on the dipole position. It is desirable to simulate as few dipoles as possible in order to reduce the simulation time.

Radiation patterns for dipoles at different position
Placing dipoles explicitly into positions of high symmetry, e.g. above the centre of the grating ridge or groove, leads to an overestimation of resonance effects compared with the background emission. This is seen in Fig. 8, where the radiation patterns of dipoles at highly symmetric positions are compared to the emission of the complete continuous layer, which is modelled by dipoles equidistantly placed with 1 nm lateral distance along the assumed emission zone, as shown in Fig. 1e.

Statistical analysis
We perform a statistical analysis to investigate whether a small number of randomly placed dipoles gives a reliable result. Small sample sizes of N = 1, 2, 4, … , 64 dipoles are randomly placed and the results superimposed for 1000 iterations. The results are shown for arbitrarily chosen wavelength and outcoupling angle ( = 550 nm , = 5 • ) in Fig. 9a. Using four or less dipoles, it is unlikely to obtain a reliable result close to the correct value.
With 32 or more dipoles the value of twice the standard deviation of the mean is so small that the error will be less than 5% with a probability of 95% . Therefore, the approach to superimpose randomly placed dipoles gives unreliable and often wrong results, unless one uses an unreasonably high number of dipoles.

The four-dipole approximation
Next, we show that it is possible to achieve accurate results by specific placement of only four dipoles. We analysed three typical OLED stacks, whose characteristics are given in Table 2. The distinct difference between these stacks is that for the PEDOT OLEDs, the emission layer itself has the highest index and forms the waveguide. For the ITO OLEDs, the ITO has the highest index and forms most of the waveguide. For all stacks and parameter variations, we calculate the radiation patterns of a continuous emission layer by superimposing the results from dipoles regularly spaced with 1 nm separation, forming a quasi-continuous layer as in Fig. 1e, and take these results as the "correct" results (this neglects the electrical effects of the nanostructure on the local current and charge density and a possible inhomogeneous emitter distribution Fujita et al. 2005). Then, we calculate the superposition of only four dipoles at different positions and compare it to the correct result. We optimized the placement of four dipoles to get the best approximation (here defined as minimizing the maximal error in the radiation pattern) of the correct result using the algorithm described in Appendix 1. Figure 10 shows exemplary comparisons between the correct result and the optimized four-dipole approximation. The difference between the results is so small that it is practically invisible. A maximum relative error below 5% is obtained for all structures. Because all the structures under investigation have the same periodicity, we can look at the distribution of optimal dipole positions, which is shown in Fig. 11. It is striking that although the analysed structures have different duty cycles and thereby different grating edge positions, the optimal positions agglomerate in four different regions (with the exception of a few outliers). These four regions are subsumed quite well by the equally spaced positions (a) (b) Fig. 8 Emission spectra at = 5 • for a single dipole centred above the grating ridge and centred in the grating groove compared and the correct result for a continuous emission layer for a the structure shown in Fig. 7a and b the same structure, but with a much wider ridge (duty cycle 0.7 instead of 0.4). In both cases, using the dipoles at the symmetry points overestimates the resonance intensity and its ratio to the background intensity ± (1∕8, 3∕8, 5∕8, 7∕8) ⋅ ∕2 , which avoid points of high symmetry. Figure 12 shows the comparison between the correct results and the approximations using these four positions for chosen structures. The difference between the results is again small. The full comparison figures are shown in Appendix 2. Even for the structures to which the outliers belong (A1 4th dipole; A2 2nd, 3rd and 4th dipole; C2 1st dipole), the approximation is good. It is therefore possible to obtain accurate radiation patterns for all investigated structures by superimposing only four dipoles placed at positions which are independent of the actual grating structure.
Of course, this four-dipole approximation cannot hold for arbitrary structures and we presume that no rigorous justification and no generally valid error estimation exist. Nevertheless, the numeric experiments show that the approximation works for the investigated structures, which have the following properties in common: period length in the range of the wavelength, not more than one structure (grating ridge) per period, no strong refractive index contrasts close to the emitters, not too sharp resonances, TE excitation. To understand better the range of applicability, we give a heuristic explanation why four dipoles yield such good approximation and show four cases to which the afore-mentioned properties to not apply. Figure 13 shows the individual farfield spectra for all 185 individual dipoles that emulate the quasi-continuous layer in structure A1 at the central resonance ( = 0 • , = 574 nm) . The dipoles' positions with respect to the grating (above ridge, groove or edge) are colour-coded. One can see that only too main curve shapes exist, one associated with dipoles above the ridge and the other associated with dipoles in the groove. Therefore, we can expect that only a few individual dipoles of each kind are needed to reconstruct all features of the overall farfield radiation pattern. In general, each dipole's radiation is determined by the available electromagnetic modes and their profiles. Thus the reasoning given above is only valid if the modes' field distributions are approximately aligned with the grating geometry and have no sharp features. We investigated structures with sharp edges and stronger refractive index contrasts close to the emission layer (structure D1, D2, E1, E2 cf. Fig. 14) in order to test the limits of the four-dipole approximation. In structure E the emission layer is interrupted. There we used only the two dipole position (5∕8, 7∕8) ⋅ ∕2 , that lie inside the emitter material (SY, yellow colour). Details of results are shown in Fig. 15 and the full results are shown in Appendix 2. In comparison with the previous results (Fig. 12), the agreement between the four dipole approximation and the continuous EML is slightly reduced. For nanostructure  Table 2. The red points are equally spaced postions at 23 nm , 69 nm , 115 nm , and 161 nm . They lie close to the centres of the individual agglomerations designs typically used in OLEDs, the test structures shown in Fig. 14 represent borderline cases. We therefore conclude that the four-dipole approximation can be used for typical nanostructured OLED design that are similar to the structures analysed in this work.

Summary
For calculating accurate radiation patterns of nanostructured OLEDs it is important to include details of the nanostructure into the simulation. It follows that a mesh resolution much smaller than the typical value of ∕10 is necessary for converged results. On the large scale, the required simulation domain size is determined by the quasi-guided mode's field profile (vertical extension) and its propagation length (lateral extension). Because this leads to a ratio of domain size to mesh resolution of around 10 4 … 10 5 , simulation times are long and an extension to full 3D simulations is not feasible. Nevertheless, the direct approach is useful for obtaining first results when exploring new structures and getting reference results to compare more advanced simulation methods. Approximating the continuous emission layer with few dipole emitters is possible for reducing simulation times. We demonstrated that four dipoles are sufficient to obtain accurate results for all investigated structures. Although this heuristic result is obtained with quite different OLED structures, a more rigorous analysis is necessary for a general statement. The fact that the same dipole positions give accurate results for all structures regardless of the grating geometry suggests a more theoretical study on this phenomenon. From a theoretical point of view, it is also interesting if the phenomenon extends to 3D simulations, although the large simulation domain size renders this approach impractical for the 3D case.  In E2 and E2*, the relative error is rather high, but the absolute error is still small due to the small amplitudes. The amplitudes are small because no resonant outcoupling occurs for = 5 • . The structures E1 and E2 have been calculated with only the two of the four equally spaced dipoles that actually lie in the emission area. Additionally, E1* and E2* show the results for four dipoles, equidistantly placed just inside the emitter material, which results again in a good approximation