Dipole radiation within one-dimensional anisotropic microcavities: a simulation method.

We present a simulation method for light emitted in uniaxially anisotropic light-emitting thin film devices. The simulation is based on the radiation of dipole antennas inside a one-dimensional microcavity. Any layer in the microcaviy can be uniaxially anisotropic with an arbitrary orientation of the optical axis. A plane wave expansion for the field of an elementary dipole inside an anisotropic medium is derived from Maxwell's equations. We employ the scattering matrix method to calculate the emission by dipoles inside an anisotropic microcavity. The simulation method is applied to calculate the emission of dipole antennas in a number of cases: a dipole antenna in an infinite medium, emission into anisotropic slab waveguides and waveguides in liquid crystals. The dependency of the intensity and the polarization on the direction of emission is illustrated for a number of anisotropic microcavities.


Introduction
Optically anisotropic materials exhibit different material properties depending on the orientation of the electric field. Traditionally they are used in polarizing beam splitters [1] and most importantly in liquid crystal displays (LCDs) which were responsible for the breakthrough of flat-panel displays in a large number of applications: TV's, mobile phones, all types of electronic devices,... The anisotropic properties of liquid crystals can be controlled by applying a moderate voltage over the device or by other external means. In the last years more and more liquid crystal based optical devices have been reported outside of their use in displays, such as: liquid crystal based lasers [2], polymerized liquid crystal devices [3], tunable waveguide devices [4], switchable/tunable filters [5],... Many other applications have been demonstrated or suggested, a more complete overview can be found in a review paper [6].
Anisotropic materials are also used in other types of photonic devices such as organic light emitting diodes (OLEDs) which are attracting a lot of attention for high quality light sources [7] and displays [8]. OLEDs feature some unique properties like: low power consumption, potential for cheap manufacturing, excellent colors, new design possibilities, ... Anisotropy has been reported in OLEDs using polymer molecules [9], whose long chains are naturally anisotropic, but recently also in small molecule OLEDs [10]. OLEDs made with polymerized liquid crystals have also been fabricated [11].
With the growing number and increasing diversity of applications for optically anisotropic materials, there is a need for accurate simulation and design tools for such applications. In this article we present a simulation method to calculate the radiation of an electrical dipole antenna inside a thin anisotropic film or in a stack of uniaxially anisotropic layers. This method simulates light emission, taking into the account the anisotropy of the materials, inside liquid crystal devices, anisotropic OLEDs and many other thin film devices. Thin layer architectures are often encountered in devices made by spincoating, deposition from vapor phase or liquid filling of cells.
Our method is based upon the decomposition of the dipole field into a superposition of plane and evanescent waves. The radiation pattern of the dipole inside the layer structure is modified by reflection, transmission and interferences of these waves by interfaces inside the layer stack. Simulation of one-dimensional thin film stacks using plane wave decomposition was pioneered by Lukosz [12] and has been applied in many fields [13] for stacks containing only isotropic materials. Wasey describes an extension of this model to take in account anisotropic layers but only allows anisotropic materials with optical axes perpendicular to the layer interfaces [14]. The method we present here can simulate the radiation of an electric dipole antenna inside a stack of one-dimensional anisotropic layers with arbitrary orientation of the optical axes. Earlier methods are restricted to stacks of isotropic materials or to stacks of uniaxial anisotropic materials with specific orientations of the optic axis [14]. Other methods such as the finite-difference-time-domain (FDTD) method or rigorous coupled wave analysis (RCWA) can be used to model the radiation from non-planar structures at the expense of more computation resources [15,16].
In section 2 the plane wave decomposition of the dipole field in an infinite medium and the alterations of the dipole field caused by layered stacks are described. The method described in section 2 is then applied to examples of anisotropic devices in section 3. First the radiation of dipole in a homogeneous anisotropic medium is compared to that of a dipole in an isotropic layer enclosed by an anisotropic medium. The radiation of dipoles in an anisotropic waveguiding layer bound by a metal electrode is discussed with particular attention to the polarization and propagation direction of the waveguided modes. The effects of coupling between the ordinary and extra-ordinary waves is also illustrated. Finally we calculate the radiation of a dipole into the modes of a waveguide in a liquid crystal. We summarize our results in section 4.

Physical background & notations
We consider the problem of a radiating elementary electrical dipole inside a stack of uniaxial anisotropic one-dimensional layers. The problem of a classical elementary electrical dipole antenna is equivalent to the emission of photons by an ensemble of excited states decaying through an electric dipole transition [17]. The coordinate system is chosen so that the x-and yaxis are parallel to the layer stack and the z-axis is normal to the stack. The orientation of the dipole moment p is defined by an inclination angle ν (with respect to the z-axis) and an azimuth angle ζ in the xy-plane. We choose the origin of the xyz-system to coincide with the location of the dipole. Each uniaxial medium i inside the layer stack is characterized by the orientation of the optical axis c i (extra-ordinary polarization) and the two eigenvalues of the dielectric tensor ε ⊥ and ε , where ε ⊥ and ε are the eigenvalues respectively in a plane perpendicular (ordinary polarization) and parallel to c i . The orientation of c i is also determined by an inclination angle α i and azimuth angle β i . A sketch of the coordinate system and layer stack is presented in Fig. 1. Any monochromatic electric field E can be written as a superposition of plane and evanescent waves. The electric field of a single plane or evanescent wave is written as: and there is a similar expression for the magnetic field H. ω is the angular frequency of the dipole antenna. Each plane or evanescent wave is characterized by its wave-vector k: Because of Snell's law waves with given k x and k y only couple to waves with the same k x and k y in the other media. This transverse part of k can be grouped to In an anisotropic medium the value for k z depends on the medium and the polarization of the plane wave. In uniaxial materials the electric field of the ordinary polarization is perpendicular to c, the extra-ordinary polarization has an electric field component parallel to c. The wave-vector can also be written as a sum of a part perpendicular k ⊥ and a part parallel k to c. The amplitude of the wave-vector in medium i is given by . λ is the wavelength of the light in vacuum and n i is the refractive index. For uniaxial anisotropic media there is a different refractive index for the ordinary waves n i,o and extra ordinary waves n i,e . The value of n i,o is the same for every direction in the medium but n i,e depends on the direction of propagation. Throughout this paper the subscripts , ⊥ and z are used to denote components respectively parallel to c, perpendicular to c and parallel to z.
A plane wave with κ < i k travels in a certain direction, the inclination angle of propagation is found with When κ > i k , z k is an imaginary number and the wave is an evanescent wave that decays exponentially with z. In non-absorbing media no power is dissipated by evanescent waves. No direction of propagation can be associated with evanescent waves since these only have a non-zero field in a small area. Evanescent waves represent the near-field of the dipole.

Radiation of elementary dipoles in a homogeneous anisotropic medium
In the following paragraphs the plane wave decomposition of the dipole field is derived. Reference [18] rewrites Maxwell's equations as uncoupled differential equations for 4 scalar Hertz potentials Θ , Ψ , Π and Φ . An uncoupled differential equation is obtained for Θ and Ψ (Eq. (74) and 75 in [18]): Where ⊥ ∇ denotes the components of the ∇ operator perpendicular to c and / ∂ ∂c is the partial derivative along c. In these equations the following source terms appear: J is the electric current, J is the component of J along c, u and v are auxiliary functions. In this case the only source is an oscillating elementary dipole exp( ) ω j t p and J , u and v are given by: A solution for Θ and Ψ can be found by performing a spatial and temporal Fourier transform on Eq. (3) to (7). Using the relations: / an expression for Θ and Ψ in the Fourier domain can be easily derived: From the scalar Hertz potentials Ψ and Θ the two remaining scalar Hertz potentials Φ and Π can be found (by performing the Fourier transform of Eq. (76) and 77 in [18]): Taking the Fourier transform of Eq. (80) in [18] then allows to calculate E in the Fourier domain from the scalar Hertz potentials Θ and Ψ and source term u: An explicit expression for E(k) is obtained by inserting Eq.
The electric field in xyz-coordinates is then found by the inverse Fourier transform of Eq. (13) along k x , k y and k z .
Only the inverse Fourier transform along k z is explicitly calculated to obtain the plane wave decomposition E(k x , k y ). The integral of the field in Eq. (13) over k z can be split in two terms and separately solved using contour integration. The first term of Eq. (13) has two poles: ε ε ε ⊥ ∆ = − and · = t c κ c 1 . These poles can either be purely real ( κ κ ≤ crit ) or contain an imaginary part ( κ κ > crit ).
The poles of the second term of Eq. (13) are given by: 2 2 µε ω κ ⊥ = ± − z k (17) and can also be either real ( 2 ). The contour integration of the two terms of Eq. (14) is performed by choosing a contour so that exp( , when z>0 this is the lower half of the complex k z plane. In case the poles are purely real a small loss term ε ′ − j is added to ε ⊥ and ε and after integration the limit 0 ε ′ → is taken. For each term only one pole lies inside the contour and only this pole is evaluated when performing the integration. For z>0 and purely real poles, the poles inside the contour are: In this way expressions are found for E (for z>0): (21) and a similar expression for z<0. This approach has the advantage that the field is decomposed in plane wave eigenmodes which can be readily used in multilayer stack algorithms, as explained in section 2.3. Another approach is to start from the formulas for the radiation of a dipole in vacuum and to anisotropically transform the problem using the method of Clemmow [19]. The results from both approaches are in agreement. For the field in (20), is independent of the direction of c, this is called the ordinary polarization. The power flux through a plane with constant z radiated by a dipole is then found by integrating the z-component of the Poynting vector S (unit W/m 2 ) over that plane: The * denotes complex conjugation. E and H in Eq. (22) are two separate double integrals over dk x dk y . When integrating over a plane of constant z we can bring both E and H under the same double integral. Because of orthogonality an extra factor 2 4π is needed in the new integrand K: The integrand K (unit W.m 2 ) can be split in an ordinary K o and an extra ordinary part K e . K is calculated for the field in a plane with z>0 (K + ) or a plane with z<0 (K -). The total power radiated by a dipole F (in W) is: From Eq. (20) and (21) it is also possible to derive the unit vectors of the ordinary and extra-ordinary eigenmodes of a uniaxial medium. The electric field of a mode is found by multiplying the unit vector with a complex amplitude. Looking at Eq. (20) and (21) it is clear that the ordinary field cannot have a component along c but the extra-ordinary field can. In an anisotropic medium the electric fields of the eigenmodes are no longer perpendicular to each other (unlike in an isotropic medium). (26)

Radiation of dipole inside an anisotropic microcavity
The radiation pattern of a dipole antenna can be significantly altered by placing the dipole inside an optical microcavity [20]. Interferences between reflections at the different layer interfaces cause variations in the local density of states at the location of the dipole antenna. As a result the angular emission pattern of the dipole is modified, this is also called the Purcell effect. Such microcavities are often made by depositing a series of thin films. The lateral dimensions of the films are much larger than their thicknesses and so the microcavities are in good approximation one-dimensional.

A dipole in an isotropic microcavity
For isotropic microcavities a model for dipole radiation based on plane wave decomposition was developed by Lukosz [12] and has been applied to simulate a wide range of applications [13]. The starting point of the Lukosz method is the plane wave decomposition of the field of a dipole antenna in an infinite medium [21]. The electric field of the dipole in an infinite medium is then altered by interference from the reflections of the various interfaces of the microcavity. A schematic can be seen in Fig. 2. The combined reflections lead to the following formula for the amplitudes , / cav TE TM E of the TE and TM waves in an isotropic microcavity based on the amplitude ∞ E of the field in an infinite medium (found with Eq. (20) and (21)) and the complex field reflection coefficients of the top and bottom parts of the cavity + A and − A :  27) and (28) the numerator expresses the effect of wide angle interference, i.e. interference of the upward emitted wave with the reflection of the downward emitted wave (or vice versa). The denominator embodies multiple beam interference: it sums all reflections between the top and bottom parts of the cavity. + A and − A may be calculated as described in [13]. Since ± A and ± ∞ E are different for the transverse magnetic (TM) and transverse electric (TE) polarization, a separate calculation is done for each polarization.

A dipole in an anisotropic microcavity
For an anisotropic microcavity a similar method can be applied to calculate the radiation of a dipole in the cavity. The eigenmodes of anisotropic materials are two linearly polarized waves, the ordinary (o) and the extraordinary (e) wave (instead of the TE and TM polarization). The normalized fields of the eigenmodes are given by Eq. amplitudes of the ordinary and extra-ordinary waves and their difference in phase. One must be careful to simulate all changes in polarization that occur during propagation and reflection or transmission in an anisotropic cavity. In anisotropic media e-and o-waves are coupled when reflection or transmission at an interface takes place, therefore the reflection coefficients have to be replaced by reflection matrices. Equation (27) and (28) should be replaced by: The reflection matrices ± A are: where lm A is the amplitude reflection coefficient for an incoming wave with polarization l into a reflected wave with polarization m.  (27) and (28) respectively.

The scattering matrix method
To calculate the reflection matrix we employ the scattering matrix method introduced by Ko [22]. Other methods, such as the Berreman 4x4 matrix method [23], can also be applied but we have chosen the scattering matrix method because it is numerically more stable when dealing with evanescent waves and total internal reflection.
. iN S is a block matrix of four 2x2 matrices. A sketch of the input and output waves in the scattering matrix method is shown in Fig. 3.
, , , , ,21 An analogous procedure is used to determine − A and − T . The fields emitted to the outside layer (N or 0) out E can be calculated from cav E : , , , , of a stack with an additional layer when ij S and the matrix j I are known [22]. j I relates the fields above the j + 1/j interface to the fields above the interface j/j-1.   We construct j I in the following way (a sketch is shown in Fig. 4). First the complex amplitudes of the input waves , , ,  At the interface j/j+1 all transverse fields have to be constant, this is expressed by the following boundary conditions: The matrix elements of the first two rows of In summary (see Fig. 4), j D relates the amplitudes of the eigenmodes in medium j at the j-1/j interface to the amplitude of the eigenmodes in medium j at the j+1/j interface, j B (or 1 + j B ) relates the amplitudes of the four eigenmodes to the four transverse field components in medium j (or j+1) and j I relates the amplitudes of the eigenmodes at the j-1/j interface in medium j to the amplitudes of the eigenmodes at the j+1/j interface in medium j+1. The boundary condition at the j+1 /j interface can then be expressed as:   , ,  , 1  ,  , ,   , ,  , 1  ,  , ,  1 , , , , ,

Calculation examples
In this section we demonstrate the abilities of the method described in section 2 by applying it to three problems: the radiation of a dipole in an infinite medium, the emission inside a multilayer structure containing an anisotropic layer and one-dimensional waveguides in liquid crystal. For these calculations the method is implemented as a computer program.

Emission from a thin isotropic film into an anisotropic medium
As a reference we calculate the total power F radiated by a dipole emitting with a wavelength 0.53 λ = µm in an infinite anisotropic medium for various orientations of the dipole moment. We consider a dipole oriented along the x, y or z-axis (respectively called p x , p y and p z ) inside an anisotropic medium ( n 1 = o , n 2 = e ) with c x . We calculate F by evaluating the integral in Eq. (24). We then also calculate F for a dipole situated in the middle of an isotropic layer (n=1) with thickness d which is bound on top and bottom side by a semi-infinite anisotropic medium with the same properties as described above. F is then calculated as a function of d (a schematic is shown in Fig. 5). The dipole moment is chosen so that this dipole radiates a power of 1 Watt in vacuum. In Fig. 6 F is shown as a function of the thickness d of the isotropic layer. For a dipole in a homogeneous anisotropic medium an analytical formula for F (Watt) can be found that depends on the refractive indices of the material and the angle ρ between p and c:  Figure 6 shows F in an anisotropic medium for dipoles p x , p y and p z . The solid lines represent F calculated with the anisotropic plane wave decomposition. The dashed lines show F for a dipole in an isotropic layer of different thicknesses. The results of the simulation for a dipole in a homogeneous anisotropic medium are the same as the results of the analytical formula (Eq. (45)). The calculation with a dipole in an isotropic layer between anisotropic media agrees with the analytical solution for 0 → d . The reflections at the layer interfaces lead to interference effects, which are a function of the thickness d of the layer and of the projections of the wave vector k x and k y . The constructive and destructive interferences in the emission pattern, lead, after integration, to an oscillating behavior of the emitted power F . As the layer thickness increases from 100nm to the µm range the number of minima and maxima in the integration increases which causes a damping of the oscillation. For very thick isotropic layers, the anisotropic material no longer influences F.

Emission and waveguiding in a thin anisotropic layer
As an example of an anisotropic multilayered structure we calculate the radiation pattern of a dipole into a film of anisotropic material deposited on a thick glass substrate and covered with a metal mirror. The configuration (depicted in Fig. 7) is an optically thick glass substrate covered with a thin isotropic emitting layer and an anisotropic layer with higher refractive index than the substrate (the parameters are for the liquid crystal 5CB) and an Aluminum mirror. The following material parameters are used: 1.5196 =  The radiation of a dipole is investigated by considering the power density ( , ) x y K k k per interval x y dk dk . The unit of K is W.µm 2 , so that after integration over x y dk dk (unit µm −2 , see Eq. (24)), we obtain the total power emitted by the dipole. Again we have chosen the dipole so that 1 Watt would be emitted by the dipole in vacuum.  In Fig. 9 + K is shown for an anisotropic waveguiding layer with c z (a total overview is shown in Fig. 9a), Fig. 9b)  κ > e CB k n no more modes are observed, because both ordinary and extra-ordinary waves are now evanescent. Since the layer stack and the orientation of the optical axis are invariant for rotation around the z-axis, the angular distribution is the same for every azimuth angle of propagation φ . In a structure which is no longer invariant to rotation around the z-axis, the radiation of the dipole will also depend on the azimuth angle φ . Figure 10 and Fig. 11 show + K for a waveguiding layer with c x . Figure 10a)

Coupling between ordinary and extra-ordinary waves
In an anisotropic multilayered structure the o-and e-waves are coupled by reflection and transmission at the interfaces in the cavity. This coupling influences the polarization and direction in which a dipole emits. We have simulated the emission from a dipole in the following layer structure (see Fig. 12): an optically thick glass substrate with a 1000nm thick anisotropic layer, a 5nm anisotropic emitting layer and another 1000nm anisotropic layer, above that there is a thick layer of air. The optical axes of the anisotropic layers are parallel to the xy-plane. The optical axes of the anisotropic layers on the glass and air side have an azimuth angle respectively 45 β°= − and 45 β°= + , for the emitting layer 0 β°= . The refractive indices of the materials and the wavelength are the same as in the previous example. A small absorption term 0.0001j is added to the refractive indices of 5CB, n e,5CB and n o,5CB , so that waveguided modes have a non-zero width.

Anisotropic waveguides
Waveguides in liquid crystals are another example of anisotropic microcavities. A onedimensional waveguide in a liquid crystal can be created by applying a voltage to a liquid crystal cell [24]. The orientation of the liquid crystal is locally altered by the applied electric field and in this way a refractive index profile is created for e-waves inside the cell. A dipole antenna inside the liquid crystal can emit light to the outside world and into the waveguide. We calculate the radiation pattern of randomly oriented dipoles inside a one dimensional waveguide, with the configuration sketched in Fig. 14. A 7.1µm thick layer of 5CB is sandwiched between two optically thick substrates. The inclination angle α of 5CB is assumed to follow a Gaussian profile around the center (z=0) of the LC layer. The profile of the inclination angle becomes 2 2 90 90 exp( / 2 ) α σ°°= − −z . The full width at half maximum of this profile is determined by σ ( 2.355σ = FWHM ). The Gaussian profile is built up as a stack of homogeneous 100nm thick layers with the inclination angle of each layer given by the Gaussian distribution for the middle of that layer. The refractive indices of 5CB and wavelength ( 0.53 λ = µm ) are the same as in the previous example. Again a small absorption term 0.0001j is added to the refractive indices of 5CB , 5 n e CB and ,5 n o CB so that waveguided modes have a non-zero width.  Figure  15b) shows the effective indices of refraction for the waveguided modes as a function of the FWHM of the Gaussian profile. For a small FWHM the waveguide is single-mode, as the FWHM increases, higher order modes arise.

Conclusions
In this paper we presented a simulation method for light emission inside uniaxial anisotropic thin film devices. In this method the electric field of an elementary dipole inside an anisotropic layer is written as a superposition of plane and evanescent waves. A scattering matrix formalism for plane and evanescent waves is used to calculate the emission inside onedimensional microcavities. The dipole emission is calculated for a number of interesting examples. A comparison between the emission from an anisotropic material and the emission from an isotropic layer embedded in anisotropic material was made. We illustrated the dependence of waveguided modes in anisotropic layers on the polarization and the direction of propagation with respect to the optical axis of the material. As a third example we calculated the modes of a one dimensional waveguides in an inhomogeneous liquid crystal and the effective indices of refraction for such modes.