Towards fully self-consistent optoelectronic simulation of planar devices

Photon recycling plays an important role in various optoelectronic devices, such as thin-film solar cells and intracavity double-diode structures presently studied for electroluminescent cooling. However, a complete description of photon recycling requires developing fully self-consistent simulation tools of photon and electronic charge transport in realistic device structures. In this paper, we describe a possible route towards such simulation tools in planar devices by combining the radiative transfer (RT) equation of photon transport with the drift-diffusion (DD) equations of electron and hole dynamics. We investigate the feasibility of the approach by selected proof-of-principle device studies. The results indicate that combining the RT and DD models not only yields the expected device characteristics, but also produces new insight into the pertinent local emission and absorption properties within the device structures. In particular, the model allows to accurately investigate how emission directivity and saturation affect the coupling coefficient in the double diode structures, showing that the main contribution to the coupling coefficient arises from the reflectivity of the top contact. For completeness, we also discuss how interference affects photon transport in resonant devices by an example calculation using the quantized fluctuational electrodynamics framework, paving way for self-consistent modeling tools that also account for wave-optical effects.


Introduction
Photonics and optoelectronics are among the key enabling technologies of the modern society allowing displays, high-efficiency light-emitting diodes (LEDs), solar cells, global optical fibre networks and multiple other technologies.However, continuous research is needed even for these basic technologies to tap their full potential and to find new application areas.One of the topics in need of further research is photon recycling, which has been acknowledged as an important process in optoelectronic devices like solar cells and LEDs already through several decades (Stern and Woodall 1974;Pazos-Outón et al. 2016;De Neve et al. 1997;Benisty et al. 1998;Saliba et al. 2015;Dupont et al. 2000;Shi et al. 2016) but that still presents a challenge for numerical modeling.A more recent example of a structure where photon recycling needs particular attention is the double-diode structure (DDS), which is presently being studied for the possibility of electroluminescent cooling (Olsson et al. 2016).The operation of the DDS is fundamentally based on emission and reabsorption of photons at different parts of the same integrated device structure.For a complete description of such devices, new modeling tools need to be developed where the optical and electrical processes are coupled fully self-consistently.This would enable e.g.optimizing the detailed layer structures of electroluminescent cooling devices, thin-film solar cells, and numerous other devices.
Several works have combined the simulation of optical and electrical properties to yield useful insight into different devices.However, typically the optical properties are included in full-device simulations without complete self-consistency with the carrier dynamics.Earliest works to address the problem estimated the effects of photon recycling using the Beer-Lambert law [see e.g.Durbin and Gray (1994)].More recently in Wilkins et al. (2016), Wilkins et al. accounted for photon recycling and luminescent coupling by calculating coupling matrices from transfer matrices and adding new terms in the radiative recombination rates.In Wang et al. (2013), Wang et al. treated photon recycling in GaAs solar cells using a ray tracing model.Particularly in the case of nanostructures, specialized and computationally heavy frameworks have been developed relying e.g. on the finite-difference time-domain method or the scattering matrix method (Walker et al. 2015;Wang et al. 2014;Kivisaari et al. 2018a).In this work, we explore the possibility to perform fully self-consistent simulations of photon and carrier transport in planar devices.We also illustrate how such simulations can provide useful insight into their emission and absorption processes.
To perform self-consistent simulations of optical and electrical transport, we make use of the drift-diffusion (DD) model of charge transport and the radiative transfer (RT) model of photon transport.Such models enable full self-consistency, if the absorption and emission in RT are coupled with the carrier distributions, and the total generation and recombination rates in DD are calculated by integrating over the photon numbers from RT.In this work we consider only planar devices, which enable writing the RT model for a single position coordinate by taking advantage of the translational symmetry in the lateral plane.The feasibility of combining RT and DD models is investigated by proof-of-principle studies of selected device structures.As an additional test, we simulate photon transport using the quantized fluctuational electrodynamics framework to expand the self-consistent simulations towards accounting for interference effects.

Theory
We start constructing the self-consistent models for the optics and electronics of planar devices by combining the DD and RT models.The conventional RT model does not account for interference effects, and it dictates how the spectral photon number ± (z, E, K) depends on the properties and excitation level of the structure.The conven- tional RT equation can be written for the present case as (Chandrasekhar 1960) where ± refers to upward and downward propagating modes (respectively), z is the position perpendicular to the layer interfaces, E is the photon energy, K is the lateral component of the full wavevector parallel to the layer interfaces, 0 is the absorption coefficient of the material in equilibrium, is the propagation angle given by sin −1 (K∕| |) , and f c , f v are the local Fermi-Dirac functions for conduction band electrons and valence band holes calculated at energy E. The Fermi-Dirac functions are given by where E Fn , E Fp are the quasi-Fermi levels and E e , E h are the electron and hole energies corresponding to photon energy E, assuming direct transitions.The electron and hole energies are calculated as where E c , E v are the conduction and valence band edge energies, E g is the bandgap energy, and m e , m h are the electron and hole effective masses.In steady state conditions, Eq. (1) would then also directly lead to the modified Bose-Einstein distribution with a chemical potential E Fn − E Fp (Würfel 1982).
The RT model is coupled with the DD model through f c and f v , which are thereby determined by the quasi-Fermi levels and the conduction and valence band edge energies.The DD model is written along the z direction as where is the static permittivity, U is the electrostatic potential, e is the elementary charge, n, p are the electron and hole densities, N d , N a are the ionized donor and acceptor densities, J n , J p are the electron and hole current densities, n , p are the electron and hole mobili- ties, and R is the net generation-recombination rate including radiative and nonradiative (1) processes.Details of the DD model used here are described in more detailed in our previous works [see e.g.Sadi et al. (2018) and references therein].
Equations ( 1) and ( 4) are all solved on a 1D domain, but for Eq. ( 4) we use the 3D densities of states, and Eq. ( 1) is solved sequentially for each E and K (and for as determined by K).Therefore the model describes the physics of wide planar structures where no changes occur in the lateral direction.The DD model in particular is a highly nonlinear equation system, whose solution schemes have attracted considerable interest for a long time (Ghione and Benvenuti 1997).Here, we are using an in-house finite-element method (FEM) implementation with the equation-based tools of Comsol Multiphysics, where the linear part is solved with the MUMPS solver and the nonlinear iteration is performed with the Newton method using an automatically adjusted damping factor.The implementation results in a satisfactory convergence for all the cases studied in this paper, but it could possibly be further improved e.g. by using the well-known Scharfetter-Gummel solution scheme (Scharfetter and Gummel 1969).
The coupled RT-DD model becomes fully self-consistent by calculating the net radiative recombination-generation as a double integral of the photon numbers as where v = c∕n r is the speed of light in the medium and g(E, ) = E 2 ∕( 2 ℏ 3 v 3 ) sin ∕2 is the energy-and angle-resolved optical density of states, accounting for both independent polarizations.In other words, Eq. (1) has to be solved self-consistently with the DD equations for all the K and E values of interest to get the total generation/recombination rate.Equation ( 5) is obtained by integrating R rad,E, = 1 E d dz S E, , where S E, is the spectral Poynting vector obtained from the photon numbers as S E, = Ev cos g(E, )( + − − ).
While the combined RT-DD model can give useful insight to the complete operation of devices whose layers are thick with respect the wavelength, it becomes less accurate with thinner layers that exhibit resonance effects.In such devices, the photon numbers ± can be replaced by the photon number expectation values ⟨n ±, ⟩ calculated as a suitably weighted average of the emission throughout the structure using the (quantized) fluctuational electrodynamics (Q)FED approach.Both FED and QFED have been developed to calculate the emission, propagation and absorption of photons by making use of the dyadic Green's function and the corresponding optical densities of states.In our case, we give a short summary of QFED, because it has been used to formulate the interference-exact radiative transfer model described shortly below.The photon number expectation values from QFED are given by Partanen et al. (2017) where ∈ {TE, TM} refers to the polarization of light, (z, K, ) is the optical local density of states (LDOS), NL±, (z, K, , z � ) are the optical nonlocal densities of states (NLDOS) of the right and left propagating fields, and ⟨ η ⟩ is the expectation value of the so-called source-field photon number, which accounts for the emission at the source location z ′ .The NLDOS of the right and left propagating fields depend on the full NLDOS as NL±, (z, K, , z � ) = NL, (z, K, , z � ) ± IF, (z, K, , z � ) , where IF, (z, K, , z � ) is the optical interference density of states (IFDOS).All the different densities of states can be expressed with the dyadic Green's functions as specified in Partanen et al. (2017) and ( 5) references therein.While there is no direct correspondence between from the conventional RT model above and ⟨n⟩ due to their different origin, they can be compared e.g. by relating the Poynting vector expressed with above with the one resulting from ⟨n⟩ .In this case the spectral Poynting vector is calculated as S K,, = Ev   2 (⟨n +, ⟩ − ⟨n −, ⟩) , and finally to obtain the full Poynting vector, S K, , has to be integrated over d and dK 2 .
The photon number in Eq. ( 6) does not lend itself easily to full device simulations, because it requires performing the integral for all possible source points and e.g.constructing nonlocal coupling terms for the emission calculation.Therefore we are also investigating the use of the so-called interference-exact radiative transfer (IFRT) model, transforming the equations for the photon number to first order differential equations (Partanen et al. 2017) where ±, (z, K, ) and ±, (z, K, ) are position-dependent damping and scattering coefficients, accounting for emission and absorption as well as photon exchange between the right and left propagating modes due to interference and internal reflections.Essentially, the IFRT model is constructed by comparing Eq. ( 7) and the photon numbers from Eq. ( 6) to establish the ±, and ±, coefficients from the dyadic Green's functions.
As mentioned above, QFED and the models derived from it require calculating the dyadic Green's function and the corresponding optical densities of states.Typically these are calculated for planar structures for example by using transfer matrices.However, as a potentially more straightforward alternative for transfer matrices, in our recent work we formulated the QFED using optical admittances in order to enable a general and straightforward framework to obtain the dyadic Green's functions and other quantities for QFED and IFRT (Kivisaari et al. 2018b).Also in this work, we use the optical admittances to calculate results where interference is accounted for using the equations presented in Kivisaari et al. (2018b).
Before we proceed to the results, we comment shortly on the expected time requirements of the electro-optical simulation framework described in this paper.Both the electrical and optical simulations are done in a 1D domain.Considering a single bias point, the the simulation is started by an electrical simulation that is carried out only for one set of parameters and is nearly instantaneous.The optical 1D simulation then follows for all propagation angles and photon energies, but as they are independent of each other, those 1D simulations can be run in parallel.Then depending on the quality of initial conditions, typically only a few iterations between the optical and electrical simulations are required for convergence.Overall, with full parallelization, the simulation of a single bias point should therefore take less than a minute depending on the size of the geometry and the spectral and angular discretization of the optical problem.On a standard multi-core desktop computer the speedup available from parallelization is limited, but even without parallelization, the simulation for a single bias point can be presently performed within minutes.These considerations should hold for both RT-DD modeling and the corresponding model constructed later with the IFRT, as calculating the and coefficients required for Eq. ( 7) is faster than performing the 1D electrical and optical simulations. (7) 1 3 95 Page 6 of 13

Results and discussion
Here we investigate the fully self-consistent RT-DD simulation tool described above by carrying out selected proof-of-principle calculations.We start by performing RT-DD simulations of a GaAs solar cell device, where photon energies spanning the full solar energy spectrum above the GaAs bandgap energy are accounted for.Then, we proceed to RT-DD simulations of double-diode structures (DDSs), where the RT is solved only for one nearbandedge photon energy for simplicity and the energy integration is replaced by multiplication by k B T .Here the model is used to study the dependence of the optical coupling coef- ficient on the applied bias.For completeness, in the end we demonstrate the effects from interference by solving the right and left propagating photon numbers in a simplified DDS from Eq. ( 6).

Thin-film solar cells
As the first example, the RT-DD method is applied to study an example shallow-junction GaAs solar cell reported in Bauhuis et al. (2016).The structure has a 3.05 μm thick p-GaAs base layer (hole density 3 × 10 22 m −3 ) in the bottom and a 0.15 μm n-GaAs emitter layer (electron density 2 × 10 24 m −3 ) on top.Below the base layer, in our simulations there is a bottom contact which also reflects 95 % of the incoming light at all angles.On top of the structure we have the top contact, which in this proof-of-principle study is assumed to be perfectly transparent for angles smaller than the total internal reflection (TIR) angle.At angles larger than the TIR angle, the top contact reflects specularly all the light coming from inside the sample.To fully account for the energy spectrum of solar radiation, here we also perform the integration over all photon energies of interest, i.e., from 1.42 to 4.43 eV.In the calculations, we assume normally incident AM1.5 solar spectrum in the RT equation.The purpose of this study is to see if such RT-DD calculations with the full energy integration are numerically feasible and result in the expected qualitative solar cell behavior.The absorption coefficient as a function of photon energy is directly obtained from Rakić and Majewski (1996) and Polyanskiy (2018).Other parameters are the standard ones and we chose not to list them fully in this proof-of-principle study, as the purpose is not to report detailed device characteristics but rather to illustrate the modeling framework in more general terms.Most importantly, however, the effective masses are m c = 0.067m 0 , m h = 0.606m 0 ( m 0 being the free-electron mass), and the mobilities are n = 8500 cm 2 ∕Vs , h = 400 cm 2 ∕Vs (Vurgaftman and Meyer 2001; Levinshtein et al. 1996).
The current-voltage characteristics as simulated using the RT-DD model are shown in Fig. 1 both under illumination ("Photocurrent") and under dark current conditions.The curve under illumination exhibits typical solar cell behaviour with a short-circuit current The RT-DD model allows also studying the internal recombination-generation processes in the structure for each bias voltage.As an example, in Fig. 2 we show the internally emitted (a) upward and (b) downward propagating photon numbers as a function of propagation angle and position at a bias voltage of 1.1 V and for photon energy 1.42 eV.The incoming solar intensity at 1.42 eV and normal incidence is also plotted in Fig. 2b in arbitrary units (starting from its maximum at the top contact and decreasing towards zero).The applied bias in Fig. 2 has been chosen to illustrate the contrast between the internally emitted and incident fields and thereby highlight the effects that can be studied with the RT-DD method.Starting from Fig. 2a, the upward propagating photon number starts to increase from zero already in the p-GaAs layer due to the quasi-Fermi level difference there.However, it increases most at the pn junction before arriving at the top contact where it is extracted or reflected back into downward propagating modes.In Fig. 2b, the downward photon number is first zero at the top contact at angles up to the TIR angle due to the boundary conditions specified above.The downward photon number then increases at the pn junction due to the large electron/hole concentration, and is almost fully absorbed within the underlying p-type GaAs layer due to photon recycling before arriving at the bottom contact.Note that without the self-consistent simulation of optical and electrical properties considered in this paper, these optical fields would only be present in the solar modes (here only at the 0 angle).In that case, emission inside the solar cell would be accounted for by a radiative recombination term that only depends on the quasi-Fermi level separation but not on the optical fields.In this paper, the self-consistent formulation combines the effects of photon recycling and the applied bias, and therefore also the upward modes exist at all angles in Fig. 2a.In this way, optical effects such as emission saturation (i.e., complete reabsorption of emitted photons) are also directly included in the calculation of the radiative recombination.This can be expected to provide more realistic radiative recombination profiles than e.g. using only a material-specific radiative recombination coefficient.

Double-diode structures
In this Subsection, we apply the RT-DD model to the DDS introduced by Olsson et al. (2016).The layer compositions, thicknesses and doping levels are chosen following Sadi et al. (2018): excluding some of the extremely thin etch-stop layers from this list to make it more compact, the photodiode on the bottom consists of 1 μm of p-doped GaAs ( 10 18 cm −3 ), 3 μm of lightly p-doped GaAs ( 3.3 × 10 17 cm −3 ), and 700 nm of n-doped GaAs ( 10 18 cm −3 ).In the LED side integrated on top of the photodiode, there is 1 μm of lightly n-doped Al 0.3 Ga 0.7 As ( 1.3 × 10 17 cm −3 ), a 300 nm thick intrinsic GaAs active layer, 300 nm of lightly p-doped Al 0.3 Ga 0.7 As ( 3 × 10 17 cm −3 ), 100 nm of p-doped Al 0.3 Ga 0.7 As ( 10 19 cm −3 ), and finally 20 nm of p-doped GaAs ( 10 19 cm −3 ).In accordance with the exper- imental setup, we have three contacts in the device, one at the bottom, one at the border between the photodiode and the LED, and one on top, all described by Dirichlet boundary conditions.The LED is biased with a forward bias, while the photodiode is short-circuited.
Here we simulate carrier transport along a straight vertical line through the middle of the structure and again solve the RT equation for all propagation angles from Eq. ( 1).This is a reasonable approximation, as the DDSs are typically 0.1-1 mm wide but the total thickness of the LED and photodiode parts is not more than a couple of μm .The main objective of this Subsection is to study how accounting for the optical propagation in the structure affects the optical coupling coefficient of the DDS.To obtain a first-order estimate of the coupling coefficient, we solve the RT equation only for the photon energy 1.42 eV corresponding to the bandgap of GaAs for all angles, and replace the energy integration in Eq. ( 5) by multiplying by k B T .We assume an absorption coefficient of 10 6 1∕m for GaAs and 10 3 1∕m for the other layers with the chosen photon energy, roughly following the val- ues reported in Rakić and Majewski (1996) and Polyanskiy (2018).
Figure 3 shows (a) the upward propagating photon number and (b) the downward propagating photon number as a function of propagation angle and position from the RT-DD simulation.The layer structure is juxtaposed with the colormaps to help interpreting them.The bias voltage over the LED is 1.3 V in Fig. 3, while the photodiode is short-circuited.Starting from Fig. 3a, the upward propagating photon number increases at the uppermost GaAs layer, which has both a large absorption coefficient and a large quasi-Fermi level separation resulting in a net emission of photons.The photon number increases more at larger propagation angles, because there the photons travel over a larger distance in the active GaAs layer.
In Fig. 3b, the downward photon number is reflected from the top contact, where we assume a constant reflectivity of 95% for all propagation angles for simplicity.The photon number increases once more at the uppermost GaAs layer and then travels with very small losses through the larger-bandgap AlGaAs layers.In the lower GaAs layers of the photodiode, the photon number decreases to zero due to absorption.In the chosen DDS structure, absorption is weaker in the n-doped GaAs layer of the photodiode due to the n-doping and the fairly large value of f c that decreases the absorption as determined by Eq. (1).Absorption becomes much stronger in the intrinsic GaAs layers below the n-GaAs, since the values of both f c and f v are much smaller there.
To evaluate the optical coupling coefficient, Fig. 4 shows the net radiative recombination rate in the LED GaAs layer, the net generation rate in all GaAs layers of the photodiode, the absorption rate in the uppermost GaAs contact layer, and the coupling coefficient defined as the ratio between the photodiode net generation and the LED net recombination.It can be seen that with the structure and parameters chosen here, the coupling coefficient is roughly 0.94 up to bias voltages around 1.4 V, and therefore it is mostly determined by the reflectivity of the top contact.For voltages larger than 1.4 V, the relative share of emission towards larger angles increases due to the saturation of the effective absorption coefficient caused by band filling.These larger-angle modes are more strongly absorbed in the thin GaAs top contact layer, and this causes the coupling coefficient to decrease slightly at the larger applied biases.This is an example of a photon recycling effect directly included in the RT-DD model, which would be challenging to account for with a simpler model.

Calculations including interference
As a last proof-of-principle example of how wave optics affect photon transport in optoelectronic devices of interest, here we study a simplified DDS made of InGaAs/InP instead of GaAs/AlGaAs.The purpose here is to explore the basic field properties that the IFRT model predicts in the optical transport and pave way for self-consistent simulations that also account for wave-optical effects.The structure to be simulated is shown schematically above both columns of Fig. 5.It contains two In 0.47 Ga 0.53 As layers separated by 2 μm of InP, and the considered photon energy is 0.81 eV, corresponding to the bandgap of InGaAs.The simulations in this Subsection are still done without full optoelectronic coupling so that electron-hole transport is not simulated.In the biased InGaAs layer, we assume a quasi-Fermi level separation of 0.8 eV, and in other parts of the structure we assume no electrical excitation.As boundary conditions, the InP layers at both ends of the structure are assumed to continue to infinity, and we assume no incoming photon fluxes.The relative permittivity for the chosen photon energy is 9.95 + 0.03i for InP and 12.46 + 0.54i for InGaAs (Polyanskiy 2018;Adachi 1989).The left and right propagating photon numbers are solved from Eq. ( 6), using the densities of states given in Partanen et al. (2017) and the Green's functions calculated from optical admittances as in Kivisaari et al. (2018b).6) for the (a) and (b) TM modes as a function of position and the lateral K vector.It can be seen that as in the case without interference, the largest increase in the photon numbers takes place at the biased active layer (InGaAs in this case).However, especially at larger K values, the right propagating photon numbers have nonzero values even before the biased InGaAs layer, and this happens mainly due to coupling with the left propagating photon numbers through internal reflections.In Fig. 5c, d, we show the left-propagating modes calculated in a similar way for TE and TM.Also here it is clear that the photon numbers increase at the biased InGaAs layer.We also see a slight decrease in the left propagating photon numbers in the unbiased InGaAs layer due to absorption.In addition, the leftpropagating photon numbers likewise show an interference pattern between the InGaAs layers especially at larger K values due to interference with the right-propagating photon numbers.
To study the energy transport in the example DDS, in Fig. 5e, f we show the magnitude of the Poynting vector for TE and TM, calculated from the right-and left-propagating photon numbers following (Partanen et al. 2015).It can be seen that the interference patterns are not anymore visible in the Poynting vectors in Fig. 5e, f, and their values are negative (positive) on the left (right) side of the biased InGaAs layer.This is the expected behavior, because the biased layer should in this case emit net radiation into both directions away from it.Furthermore, the absolute value of the Poynting vector increases towards larger K values, as light again propagates a larger distance in the emitting GaAs layer.It can also be seen that the absolute value decreases in the unbiased InGaAs layer due to absorption.Minor resonant effects can be seen in the Poynting vector for TE in Fig. 5e, as its value is not a monotonic function of K∕k 0 .We expect such effects to become more important when the layer thicknesses decrease more towards the wavelength of light.

Conclusions
In this paper we described our efforts towards the development of fully self-consistent modeling tools for optical and electrical properties of planar devices.Combining the radiative transfer model of photon transfer with the drift-diffusion model of electron-hole dynamics allows directly connecting the predicted device characteristics with the pertinent local emission and absorption properties within the device.We believe that this will allow useful new insight for optimizing various optoelectronic device structures for different existing and emerging applications.Completing the models with wave-optical effects such as interference and models for scattering surfaces such as plasmonic or other gratings will further allow studying the performance limits of very thin devices exhibiting strong resonance effects.
Acknowledgements Open access funding provided by Aalto University.We acknowledge financial support from the Academy of Finland and the European Research Council (Horizon 2020 research and innovation programme, Grant Agreement No. 638173).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creat iveco mmons .org/licenses/by/4.0/),which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Fig. 1
Fig.1Current-voltage characteristics of the GaAs solar cell under illumination and in dark current conditions simulated with the RT-DD method, accounting for all photon energies and propagation angles

Fig. 2
Fig.2Photon numbers from the RT-DD simulation for a upward and b downward propagation as a function of propagation angle and position.The colormaps include only photons emitted inside the structure, and the incoming solar intensity at 1.42 eV and normal incidence is plotted in arbitrary units in (b), starting at its maximum at the top contact and decreasing towards zero.Bias voltage is 1.1 V and photon energy is 1.42 eV.The pn-junction is marked with a dashed horizontal line

Fig. 3 aFig. 4
Fig. 3 a Upward propagating photon number and b downward propagating photon number in the DDS from the RT-DD simulation as a function of the propagation angle and position.The simplified layer structure is juxtaposed to make it easier to interpret the figures.The bias voltage is 1.3 V over the LED, and the photodiode is short-circuited

Fig. 5
Fig. 5 Photon numbers and Poynting vectors from QFED for the structure illustrated above the two panels as a function of position and lateral K vector: a right propagating photon numbers for TE, b right propagating photon numbers for TM, c left propagating photon numbers for TE, d left propagating photon numbers for TM, e magnitude of the Poynting vector for TE, and f magnitude of the Poynting vector for TM