Dipole emission in stratified media with multiple spherical scatterers: Enhanced outcoupling from OLEDs

Abstract Scattering particles find application in organic light emitting diodes (OLEDs) for an enhanced outcoupling of the generated light. This paper presents a computational scheme to exactly model the electromagnetic fields and the power outcoupling efficiency of a typical OLED geometry, comprising a thin film system with spherical scattering particles inside. The model is based on the expansion of the fields in plane and spherical vector wave functions, as well as the scattering matrix formalism for the layer system reflections. In a numerical application example, the effect of 1000 spherical high index scattering particles on the internal outcoupling from a realistic OLED structure is discussed.


Introduction
White organic light emitting diodes (WOLEDs) are an appealing technology for general lighting applications. They are characterized by a broad emission spectrum and large active areas. In principle, it is possible to manufacture OLEDs in a roll to roll printing process on flexible substrates, which can in future lead to reduced production costs. Regarding device efficiency, one substantial challenge is the limited light outcoupling: The excitation of waveguide and surface plasmon polariton modes in the OLED thin film system confines a significant part of the generated light to the layer stack. Only a minority of the generated photons, typically 20-30%, can escape and contribute to usable light [1]. Since state of the art OLEDs already show internal quantum efficiencies near unity [2], poor light extraction in fact turns out to be the most severe efficiency limiting factor.
The trapped light can be partly outcoupled by substrate modification techniques like micro lenses or surface roughening. However, if one wants to access the photons bound to the OLED thin film system, it is necessary to include scattering geometries right inside the thin films. One possibility is scattering layers based on nanoparticles, which partly scatter the bound modes to the far field, increasing the light output of the device [1,3].
A quantitative description of the dipole emission in a thin film system with scattering particles is of great importance for the design of improved light emitting devices. Approximate methods like ray tracing have been suggested [4], but in case of dense scattering centers, geometrical optics cannot fully account for coherent multiple scattering. Purely numerical 3D approaches, including FEM and FDTD, can only handle volumes of limited size [5]. Therefore, only a limited number of scattering particles can be considered. Semi analytical methods like the volume integral method were also applied to the problem [6]. This avoids discretizing the whole layer system, but a discretization scheme is used to model the scatterers. If the scattering particles represent spheres, a purely analytical approach by means of an expansion of the scattered field in spherical vector wave functions becomes feasible. In the context of scattering layers in OLEDs, this method was recently suggested by Tishchenko [7]. The objective of this paper is to present a comprehensive and exact computational procedure to model the electromagnetic dipole radiation, including the full scattering interaction of the geometry of layer interfaces and spherical particles. The method is built on three cornerstones: the dyadic Green's function for the light emission process [8], the scattering matrix for the response of the layer system and the transition matrices (T-matrix) for the spherical scatterers.
Researchers from different fields, including astronomy, atmospheric physics, material characterization and biophysics, have studied the scattering by clusters of particles by means of T-matrix methods. A comprehensive database of papers dealing with related problems can be found in [9]. Its application to electromagnetic scattering by a single particle near an interface was elaborated in 1980 by Kristensson [10]. In 2008, Mackowski published the solution to the problem of multiple spheres near a single interface [11]. The general method, applicable also for non-spherical scatterers above films with several interfaces, was discussed by Doicu et al. [12]. Very recently, Kristensson presented a treatment of scattering by a single perfectly conducting particle embedded in a parallel plate waveguide, including dipole radiation as the incident field [13].
The present study is similar to the approach presented in [11]. The direct scattering interaction of the spheres is considered by employing an addition theorem for the spherical vector wave functions, whereas the interface mediated scattering interaction requires a transformation of the electric field representation from spherical waves to plane waves. Each partial plane wave is then scattered by the layer interfaces and retransformed into a series of spherical wave functions again. In contrast to [11], the plane interfaces are located on either side of the scattering centers. For this configuration, the layer system transition matrix (see Section 5) plays the role of a generalized reflection coefficient. In addition, the incident field of the present analysis is given by a point dipole source, also located inside the layer system, whereas in [11], the incident field is given by a plane wave.
The paper is organized as follows: In Section 2, an overview on the problem and the computational method is given. We introduce the vector wave functions and the relevant transformation relations in Section 3. Sections 4 and 5 treat the scattering by a spherical particle and by the layer system, respectively. Section 6 is dedicated to the calculation of the initial field coefficients. In Section 7, we present the multiple scattering relations, before in Section 8 the linear system for the scattered field coefficients is compiled. Section 9 treats the evaluation of the far field pattern and the power flux, as well as the dipole life time. Finally, in Section 10, we use two application examples to validate the model and to demonstrate the use of the method for a realistic OLED structure before we draw conclusions in Section 11.

Overview
The investigated configuration is depicted in the lefthand side of Fig. 1. It consists of a stratified medium with interfaces normal to the z-direction. The outermost layers are confined by only one interface and are thus semiinfinitely extended. A single Hertzian dipole is included in one of the layers and an arbitrary number of scattering spheres are located in the same or another layer than the dipole. The task is to calculate the total electric field for a given dipole strength and orientation.
The here proposed method relies on the expansion of the incident and the scattered electric field of each sphere in spherical vector wave functions (SVWFs). The key step of the analysis is to evaluate the incident field for each scattering sphere. It is given by the sum of the following three contributions (see the right-hand part of Fig. 1): the initial dipole field (including the layer system response), the scattered field of all other spheres and finally the layer system reflection of the scattered field. When the incident field on a sphere is known, the scattered field can be evaluated using the transition matrix of that sphere. The above implies that the incident field on each sphere is a function of the scattered field coefficients of all spheres. On the other hand, the scattered field coefficients naturally depend on the incident field. Therefore, the incident and scattered field coefficients cannot be determined one after the other but they are calculated simultaneously by solving a self-consistent linear system [12]. This way, multiple scattering is automatically accounted for exactly and no scattering order approximation is necessary. Finally, when the scattered field coefficients of the spheres are known, all quantities of interest can be calculated from these coefficients.
This approach is subject to the following limitations: We restrict ourselves to isotropic and non-magnetic media. The layers containing the spheres and the dipole must be lossless. The infinite series expansion of the fields in spherical waves is truncated in order to get a finite set of equations. However, if the spherical waves are cut-off at a high enough multipole order, the error introduced by this truncation is very small [14]. Further, the assembly of the linear system includes the numerical evaluation of Sommerfeld integrals.
The geometry of the problem includes both spherical and planar structures. Therefore, the electrical field will be expanded in SVWFs and plane vector wave functions (PVWFs) and the conversion of the one representation to the other plays a key role in the approach.

The vector wave functions
We seek the solution of the time harmonic Helmholtz equation, with a time dependency expðÀiωtÞ implicitly understood for all fields. The boundary conditions are given by the radiation condition for z-71 and continuous parallel components of the electric and the magnetic field at all interfaces. The fields will be expanded in basis solutions of the homogeneous Helmholtz equations, the so-called vector wave functions (VWFs). The task is then to find the proper coefficients for this expansion. The vector wave functions and their transformation relations are taken from [15] and we follow the notation and normalization conventions of that reference.

Definition of the VWFs
The spherical vector wave functions (SVWFs) are defined as where the radial wave function z ðνÞ l stands either for the spherical Bessel function of order l, z ð1Þ l ¼ j l , or the spherical Hankel function of first kind, z ð3Þ l ¼ h indices of Ψ ðνÞ τσml stand for τ the polarization (TE ¼ 1, TM ¼ 2), σ the parity (1 ¼even, 2 ¼odd) of the ϕ dependence, m ¼ 0…l the angular index with respect to ϕ and l ¼ 1; 2; … the angular index with respect to θ. A multiindex n is introduced to subsume them all, ðτσmlÞ-n. The number ðνÞ indicates if the SVWF is of regular kind ðν ¼ 1Þ or represents an outgoing wave ðν ¼ 3Þ.
A second set of basis solutions is given by the plane vector wave functions (PVWFs) [15]: Here, ðk; α; βÞ are the spherical coordinates of the wave vector k. The index j of Φ j indicates the polarization, with j¼1 corresponding to TE waves and j¼2 corresponding to TM waves.

Transformations of the VWFs
The spherical wave functions can be transformed into plane waves and vice versa. We will need the expansions Ψ ð3Þ The transformation matrix B is defined in the appendix. A shorthand notation for the double integration over the angles α and β was introduced: k . The contour C used in (2) runs from 0 to π=2 and then parallel to the imaginary axis to π=2Ài1. Note that in (2) for z o 0 the contour C À from [15] was replaced by C þ ¼ C by means of the variable substitution α-α À ¼ π À α, dα-Àdα À , C À -À C þ . In the following, we will use the shorthand notation α À for π À α. It marks the waves propagating in the À z-direction.
Translations of PVWFs are just a trivial phase shift: In case of SVWFs, the following addition theorem accounts for translations of the coordinate origin [15]: The translation operator P nn 0 ðdÞ is defined in the appendix.

Scattering by a sphere
Consider a spherical scatterer of radius R and refractive index n S , embedded in a medium with refractive index n M .
The corresponding wave numbers are k S ¼ n S ω=c and k M ¼ n M ω=c. Let r S be the sphere's center position. The incident field can be expanded in terms of regular SVWFs relative to r S : n ðr Àr S Þ If we write the scattered field in terms of outgoing SVWFs relative to r S , the coefficients of the scattered field are related to the coefficients of the incoming field by the transition matrix For spherical scatterers, the matrix T S is diagonal in all indices τ; σ; m and l. It can be determined by considering the continuity conditions of the electric and the magnetic field at the boundary of the sphere. An explicit expression is given in the appendix, see (C.1).

Response of the layer system
We call z 1 o ⋯ o z N the positions of the layer interfaces and L 0 ; …; L N the layers, where z i bounds L i from below. The layers are characterized by their refractive indices n 0 ; …; n N and thicknesses d 1 ; …; d N À 1 (the outermost layers are semi infinite in size).
Reflections at the planar layer interfaces couple the forward propagating waves (plane waves with αrπ=2) to their backward propagating counterpart, running in the α À ¼ π À α direction. This coupling is most conveniently represented by a matrix formalism. In the following, we therefore introduce a 2-vector notation for the coefficients of forward and backward propagating plane waves. The layer system response to any excitation can then be written in terms of a generalized 2 Â 2 transition matrix. Expanding the electrical field in layer i in PVWFs, we write The coefficients g þ i and g À i belong to the forward and backward propagating plane waves of layer i, respectively. They refer to the field strength at position R dα sin α. The integration variable α need not coincide with the propagation angle α i : if the field is generated in another layer i exc with refractive index n exc and then propagated into layer i, the integration variable α will be the angle of propagation in the emitting layer. α i and α À i ¼ π Àα i are then understood to be functions of α, as the refraction of waves through different layers obeys Snell's law: sin α i ¼ n exc =n i sin α. We can set α i ¼ arcsinðn exc =n i sin αÞ such that 0 rRα i rπ=2 and Iα i r0.

Layer system scattering matrix
In the scattering matrix scheme the coefficients of the outgoing plane waves are expressed as a linear function of the incoming amplitudes. In contrast to the also used transfer matrix scheme, this approach leads to a numerically more stable algorithm [7]. The incoming and outgoing wave amplitudes of a layer system are grouped into vectors. The scattering matrix is then defined by the relation It can be calculated in a recursive algorithm following the steps described in [16].

Waves generated inside the layer system
Consider an electromagnetic radiation source or a scattering center, located at r exc in layer i exc . The excitation can be a dipole source D or the scattering by a sphere S. The excited field in the layer containing the excitation is then of the form In other words, at positions above the excitation origin the field is expressed in terms of plane waves propagating in positive z-direction (with propagation angle α þ ¼ α), whereas at positions under the excitation, it is represented by waves propagating in negative z-direction (with angle α À ¼ π À α). This discontinuous definition corresponds to the nature of the excitation as an outgoing field. However, if the evanescent fields are correctly accounted for in the expansion, the total field is again smooth at all points except r exc , where it is singular.
Both, the upward and downward propagating waves experience refraction and reflection, such that the total response of the layer system is a combined reflection of the upward and downward propagating partial waves.
To model this coupling of up and down propagating partial waves, the layer system is split into two subsystems by adding a virtual interface at the position of the excitation, z exc (see Fig. 2, where in the left, exc ¼ D representing dipole excitation and in the right, exc ¼ S for the scattering case). The bottom subsystem consists of all layers below the excitation origin, i.e. the interfaces z 1 ; z 2 ; …; z iexc ; z exc . Its top layer is a hypothetical semi-infinite layer with refractive index n iexc , bounded from below by z exc . On the other hand, the top subsystem starts with an interface z exc and contains all interfaces with z 4 z exc , i.e. z exc ; z iexc þ 1 ; …; z N . Let S 0;exc j and S exc;N j denote the scattering matrices of the bottom and the top subsystem, respectively.
We write the layer system's response to E exc , evaluated in the layer containing the excitation, as The subscript R indicates the layer system response. In the following, we will suppress the polarization index j and the angular arguments α and β for a clearer notation. First, we consider the scattering matrix equation for the lower subsystem. The field propagating into the lower subsystem from z exc is given by the exciting field plus the reflected field coming from the top layer subsystem, represented by g À R;exc . The reflected field, coming from the bottom subsystem, is again given by g þ R;exc : The zeros on the right-hand side of Eq. (9) are due to the fact that no light is incident on the lower subsystem from below. The analog scatter equation for the top subsystem reads Eqs. (9) and (10) can be used to determine the linear relation between the excitation and the response of the layer system: with The subscript exc indicates that the transition matrix represents the response of the layer system, evaluated at the excitation position. When the g 7 R;exc are known, we can also evaluate the electric field in other layers, caused by the excitation. It is represented by the coefficients g 7 R;i : with again α i as a function of α ¼ α iexc in the sense of Snell's law. For layer i with z i 4 z exc , we employ the scattering matrix S exc;i for the layer subsystem with the interfaces z exc ; …; z i : For z i oz exc the scattering matrix equation reads with S i;exc the scattering matrix for the layer system z i ; …; z exc . The linear dependency of the coefficients g 7 R;i in layer i on the excitation coefficients g 7 exc can again be written in a transition matrix equation: The transition matrix T L exc;i is a generalization of the layer system transition matrix T L exc , incorporating the layer system response to an excitation, evaluated in layer i. It can be calculated by solving Eqs. (13) and (14) The field coefficients g and the scatter and layer system transition matrices S and T L have to be calculated separately for the two polarization states j¼1,2 and for each propagation angle α. Although the here used scattering matrix algorithm is more stable than the so-called transfer matrix scheme, the exact calculation of the transition matrix T L exc;i causes trouble if the z-component of the wave number has a large imaginary part (i.e. in the limit α-π=2Ài1). In particular the propagation of the layer system response to another layer turns out to introduce instability. However, the response of the layer system is damped by expði cos α kΔzÞ-expðÀ sin α kΔzÞ with Δz the scattering path length. Thus for large sin α it is a good approximation to ignore multiple scattering by using Eq. (16) where we set the reflection coefficient of the layer interfaces to zero in the calculation of the matrices S.

Dipole field
Using the formalism developed in the last section, it is a straightforward matter to evaluate the dipole field at the sphere positions. The dipole is characterized by its current density jðrÞ ¼ ð1=iωÞδðr À r D Þd, located at r D ¼ ð0; 0; z D Þ in layer i D . The exciting field of this dipole is then given by the dyadic Green function: with [15] Gðr; for z≷z D : The functions Φ † have the explicit i in (1) changed to À i.
Eq. (17) has the form of equation (8) with exc ¼ D and One can use the layer system transition matrix (15) to get the propagated initial field coefficients g 7 R;iS;j for layer i S , i.e. the layer containing the spheres, such that the dipole excited field reads k ½Φ j ðα iS ; β; r Àr iS Þ; Φ j ðα À iS ; β; r À r iS Þ Using (4), the PVWFs are translated to the origin of sphere S: ðβ Àϕ S Þ, where ðρ S ; ϕ S ; z S Þ represent the cylindrical coordinates of r S , we can analytically evaluate the β-integration. ! id z sin αI 1 m ðp n ; q n ; k i S ρ S sin α i S ; ϕ S Þ þ À 1 1 ! i cos αI 2 m1 ðp n ; q n ; d x ; d y ; k i S ρ S sin α i S ; ϕ S Þ:

It results in
The integral expressions abbreviated by I 1 and I 2 are given in the appendix, as well as the Boolean expressions p n and q n , see (A.3)-(A.5). Finally, we can write the incoming field coefficients at sphere S for the layer system mediated dipole excitation: Eq. (19) contains a Sommerfeld type integral, which has to be carried out numerically. Note that α is the angle in the dipole layer and α i S is a function of α via Snell's law.

Multiple scattering
In this section, we will write the scattered field from sphere S 0 , as well as its reflection from the layer system, in the form of an incident field on sphere S.
The scattered field reads In order to express it in terms of regular spherical waves around the center of sphere S, we make use of the translation formula (5) with d ¼ r S À r S 0 : In the second step, the matrix W S;S 0 was introduced. It represents the linear relation between the scattered field coefficients from sphere S 0 and the incoming field coefficients at sphere S.
Up to now, only the direct scattered field from sphere S 0 was considered as an incoming field at S. To account for the full interaction of the spheres, also the reflections of the layer system have to be taken into account. The first step is to expand E S 0 in plane waves, using (2): n 0 ∬ C d 2ê k B n 0 ;j ðα 7 ; βÞΦ j ðα 7 ; β; r Àr S 0 Þ; z≷z S 0 : ð21Þ As discussed in Section 5.2, the layer system's response to this field is calculated by means of the layer system transition matrix, with exc ¼ S 0 : To translate the plane waves to the origin of sphere S, use (4). Afterwards, the PVWFs Φ are replaced by an expansion in regular SVWFs around S, using (3).
where the coupling matrix W S;RS 0 is given by dα sin α ½B pol † n;j ðαÞe ik cos αz S 0 S ; B pol † n;j ðα À Þe À ik cos αz S 0 S T L S 0 ;j ðαÞ B pol n 0 ;j ðαÞ B pol n 0 ;j ðα À Þ 2 4 3 5 I 2 mm 0 ðp n ; q n ; p n 0 ; q n 0 ; kρ S 0 S sin α; ϕ S 0 S Þ: ð24Þ I 2 , p and q are defined in the appendix and ðρ S 0 S ; ϕ S 0 S ; z S 0 S Þ represent the cylindrical coordinates of r S À r S 0 .

Self-consistent equation
Combining the results from the preceding sections, we can write down the linear system that has to be solved in order to find the electromagnetic field coefficients. For sphere S we have, from (6), If we subsume the coefficients b n S into an overall coefficient vector b, we can formally rewrite (25) in a matrix-vector Once the coefficient vector b is determined, all quantities of interest can be derived from it, including the scattered near and far field, as well as the dipole life time and the power flux through any surface.
9. Far field, power flux and dipole life time The layer system mediated electric field at a position r is calculated as with i the layer containing r. Note that in the layers i D and i S , the direct radiated or scattered field has to be accounted for, too. The coefficients g 7 i;j in Eq. (27) are a superposition of the contributions from the dipole and the scattered fields from the spheres: The Jacobi factors are necessary, because the fields propagated through the layer system are now calculated by an integral over the angle in the layer containing the field point, in contrast to Eq. (12). Using Snell's law, we can write d 2ê The magnetic field is then (we set μ ¼ 1) The complex (conjugated) Poynting vector is defined as It can be used to calculate the time averaged power flux in z-direction: The r J dependency of the integrand is given by expðiðk J À k 0 J Þ Á r J Þ and thus the integral yields Further, the integrand vanishes for j a j 0 , because k but it is parallel to Φ 1 and therefore it has no z-component. Carrying out theê k 0 integral yields If the power flux through the bottom or the top layer is calculated, either g þ or g À is zero. Then, the near field carries no power and the integral can be carried out over the real angles 0 rα r π=2: The integrand of Eq. (30) is the power emitted into top/ bottom layer per solid angle. It is thus the far field pattern of the solution. The β integral can be carried out analytically. The resulting expressions are quite lengthy and can be found in a supplementary file.

Dipole lifetime and outcoupling efficiency
The decay rate of the emission dipole is proportional to the work that it does against its own field. It can be determined either by the power radiated through two surfaces, above and below the dipole, or by the electrical power dissipated by the source current against the total electrical field at the dipole position [8,17]: In an infinite medium, the total radiated power that a dipole emits is calculated as In the presence of interfaces and scattering centers, the contributions from the reflected and the scattered field need to be added: They can be evaluated by means of Eq.
An important figure of merit for the optical design of light emitting diodes is the outcoupling efficiency. It is defined by the power that the dipole and the scatterers feed into the far field modes of the system, divided by the total electromagnetic power of the dipole:

Validation and application example
The computational method presented in Sections 2-9 was implemented in a MATLAB s routine. For the evaluation of the Sommerfeld integrals (19) and (24), the values of the most time consuming functions along a fixed grid in α are stored in a lookup table. The trapezoidal rule is used for integration. For the evaluation of the Bessel functions, a two dimensional lookup table is employed to allow for an interpolation in both α and ρ S 0 S . To reduce the necessary sampling rate in α, the integration path is deformed towards 0-1-π=2 À 0:1i-π=2 À 1:5i in order to avoid narrow peaks in the integrand which stem from resonances corresponding to guided modes of the layer system. This deformation of the integration path is justified, because the integrand is an analytic function of α in the domain 0 r Rα r π=2, Iα r 0. The linear system (26) is solved by means of LU decomposition.

Validation
The first test geometry for our implementation of the here presented method is a simple lossless layer structure that supports no thin film waveguide modes. This way, we can verify the conservation of energy in a straightforward manner. The structure consists of a high index substrate (n¼1.7), a scattering layer (thickness d ¼300 nm, n M ¼ 1:5), a luminescent layer (d ¼200 nm, n¼1.6) and an almost perfectly conducting mirror material (n ¼ 0:0001 þ 1000i). In the scattering layer, nine scattering spheres with radius 120 nm are located in a square grid pattern, at positions x; y ¼ À300; 0; 300 nm, in the middle of the scattering layer. The field is excited by an electrical point dipole, emitting at a vacuum wavelength of 520 nm, located at ðx; yÞ ¼ ð0; 0Þ in the middle of the luminescent layer. For the refractive index n S of the spheres, we consider the values 2.2, 2.6 and 3.
The scattering configuration is depicted in the left-hand side of Fig. 3, and the right-hand side shows the calculated far field pattern for a parallel dipole orientation and n S ¼ 3. With regard to the conservation of energy, we check that the total far field power equals the dipole dissipated power, 〈P far 〉 ¼ 〈P D 〉, using Eqs. (30) and (31). For all six cases (vertical/parallel dipole, n ¼ 2:2; 2:6; 3), we verify the energy conservation with a relative precision better than 10 À 3 . Fig. 4 shows the electric near field, evaluated along a probing line that runs diagonally across the sphere array, inside the substrate with a distance of 100 nm to the substrate/scattering layer interface, see Fig. 3. The case of a vertically oriented dipole and a parallel dipole is shown on the left and the right side, respectively. For comparison, we compute the near field also by means of the FEM, using a commercially available software (COMSOL Multiphysics s ).
The results are in perfect agreement. So we have shown that our implementation respects the conservation of energy and is consistent with near field calculations from established software packages.

Application example: OLED outcoupling
To demonstrate the applicability of the presented method, we use it to model a generic bottom emitting OLED structure with a scattering film for enhanced internal outcoupling, see the left-hand side of Fig. 5. The OLED consists of the following layers: glass substrate (n ¼1.5)/ scattering layer (d¼500 nm, n M ¼ 1:55)/transparent anode (d ¼100 nm, n ¼ 1:95 þ 0:0035i)/organic (d ¼ 125 nm, n ¼ 1.75) /metal cathode (n ¼ 0:68 þ 5:32i). In the scattering layer, a set of high refractive index spheres with radius 120 nm is distributed randomly with a particle density of 20=μm 3 . This configuration is identical to the green bottom emitting OLED described in [1], except that we assume a scattering layer thickness of 0:5 μm instead of 2 μm. As the scattering effect on the internal waveguide modes is exponentially damped with increasing distance to the waveguide core layers, the outcoupling due to a given number of particles can be better studied with a lower scattering film thickness.
The figure of merit of our analysis is the substrate coupling efficiency, i.e. the power radiated into the substrate divided by the total dissipated electromagnetic power. For the special case of a vertically oriented dipole, located in the middle of the organic layer and emitting at 520 nm vacuum wavelength, the substrate coupling efficiency without scattering is very low: only $ 2% of the dissipated power goes into the substrate far field, whereas 490% of the power is coupled into the surface plasmon polariton mode (SPP). This is in contrast to the case of a parallel dipole, which shows a substrate coupling efficiency of η sub 4 70% already without scattering particles.
The authors of [1] compare the decay length of the respective waveguide modes to the calculated bulk mean free path in the scattering layer. For the fundamental TE waveguide mode, the energy decay length is $ 13 μm and larger than the scattering length such that an extensive scattering of the energy out of this waveguide mode can be expected. The substrate coupling efficiency is then mainly limited by the reflection and absorption losses of the direct  dipole radiation and of the scattered radiation from the particles. Much more difficult is the extraction of the SPP mode: first, the mode overlap with the scattering layer is considerably smaller and, second, the energy decay length of the SPP mode is only $ 1:6 μm. Coupling only to TM waves, the vertical dipole is therefore an ideal testing ground for the presented method: as most of the relevant scattering happens inside of a circular domain of 1-2 SPP decay lengths around the dipole position, we can expect the estimated substrate coupling efficiency to converge for a moderate number of model scattering particles.
The left side of Fig. 5 shows the layer system, together with the electric field profile of the SPP mode. The solid line refers to the mode profile in the original layer system, whereas the dashed line refers to the case of a scattering matrix material with an increased refractive index of n M ¼ 1:8. To define the sphere ensemble (right-hand side of Fig. 5), we successively generate evenly distributed random sphere center positions in a cylindrical domain of radius 5:6 μm, aiming at a final sphere density of 20=μm 3 . The minimal distance of the sphere centers to the scattering layer interfaces is set to 120 nm in order to prevent intersection. In addition, the algorithm checks for each newly added sphere if it intersects with one of the preceding spheres, in which case the new sphere is discarded. The algorithm terminates, when the ensemble includes the desired number of 1000 spheres. To study the influence of the sphere number on the outcoupling (see Fig. 6, left), smaller sphere samples are realized by considering only the N nearest spheres to the origin.
In principle, one has to average over the random realization for the sphere positions. Instead, we keep the sphere ensemble fixed and average over nine dipole positions in a quadratic region of size 1 μm 2 , assuming that the local dipole/sphere configuration is the dominant influencing factor for the substrate coupling efficiency of a specific random realization. Fig. 6 shows the calculated substrate coupling efficiency for a vertical dipole emitting at a vacuum wavelength of 520 nm, as a function of the number of spheres for a fixed sphere density (left). Hence, the size of the area containing scattering spheres increases accordingly. Due to the short decay length of the SPP modes, the value for 1000 spheres is near to the asymptotic value for an infinite number of spheres. The other TM waveguide modes decay more slowly, but are only weakly excited by the dipole. Together  with the TE waveguide modes that were excited by scattering, they are responsible for the residual slope in the substrate coupling efficiency at 1000 spheres. The right-hand side of Fig. 6 shows the dependence of the substrate coupling efficiency as a function of the sphere refractive index for 1000 spheres. The error bars stem from the averaging over the nine different dipole positions. For comparison, the case of a high index scattering layer (n M ¼ 1:8 instead of 1.55) is shown, too. Due to the more gentle decay of the evanescent mode field, a larger mode overlap with the scattering spheres leads to a more efficient outcoupling, given that the particles are made of a material with a sufficiently high refractive index in order to ensure a considerable index contrast to the host material. For very large refractive indices, a substrate coupling efficiency of 4 20% can be expected, which is more than twice the maximal value for the case n M ¼ 1:55. Due to the resonant behavior of the Mie scattering efficiency, the substrate coupling efficiency has a maximum at n S % 2:9.
An experimental external quantum efficiency enhancement by a factor of approximately 1.5 was reported by the authors of [1]. This enhancement can be attributed to the positive effect of the scattering particles on both, the outcoupling of the thin film bound modes and on the extraction of the light from the substrate into air (which is beyond the scope of this study). In addition, the statistical distribution of the dipole moment orientation has to be considered. However, for the case of a parallel dipole orientation, the largest loss channel is the coupling to the fundamental TE waveguide mode with a larger decay length than the SPP mode. Therefore, a larger scattering sample would be required to reliably estimate the substrate coupling efficiency.

Conclusions
We have developed a computational scheme to model the propagation of light in a layered structure with spherical scattering particles embedded inside the layer system. This approach can be used to describe light outcoupling in OLEDs containing scattering layers. In comparison to FEM based calculations of the near field and by checking the conservation of energy, we have demonstrated the validity of the approach. As the fields are represented by a series in outgoing spherical vector wave functions (36 parameters per sphere for l max ¼ 3), the exact calculation of fields in OLED structures including a significant number of particles becomes feasible with moderate memory cost. This is an advantage over purely numerical approaches like FEM or FDTD, which are limited to a small computational domain. Furthermore, if the dipole position or orientation changes, only the initial field coefficients a n RD have to be recalculated, whereas the sphere interaction matrices W and W R remain unchanged. Therefore, an average over dipole position and orientation is possible with very little extra computational cost. The same holds for a variation of the sphere refractive index or radius (which only affect the transition matrix T). In fact, the scattering spheres can be exchanged by any species of scattering particles for which the transition matrix T is known. The limiting factor of the presented approach with respect to both, computation time and memory, is the large dimension of the interaction matrices W and W R , which is proportional to the number of considered spheres. The most time consuming task is the assembly of the layer system mediated interaction matrix W R which involves the numerical evaluation of an integral for each matrix element. For an improved computational performance, a sophisticated algorithm for the efficient evaluation of these Sommerfeld integrals is essential.

Acknowledgments
We wish to thank Christian Yorck, Jan Mescher, Carola Moosmann and Siegfried Kettlitz for inspiring discussions. A.E. gratefully acknowledges financial support from the Karlsruhe School of Optics & Photonics (Grant no. GSC 21).

Appendix A. Transformation matrix B
The transformation matrix B for the transformation between the SVWFs and the PVWFs [15] can be split into two parts depending only on α or β: Here, the Boolean coefficients p and q select the cos or the sin βÀdependency. They read p τσ;j ¼ δ τj δ σ1 þð1 À δ τj Þδ σ2 ; q τσ;j ¼ δ τj δ σ2 Àð1 À δ τj Þδ σ1 ðA:3Þ The matrix B † has all explicit i in (A.2) changed to À i.