In this paper a fast, yet accurate method to estimate the spectral and angular distribution of the scattered radiation of a fluorescent material is described. The proposed method is an extension of the adding-doubling algorithm for non-fluorescent samples. The method is validated by comparing the spectral and angular transmittance and reflectance characteristics obtained with the extended algorithm with the results obtained using Monte Carlo simulations. The agreement using both methods is within 2%. However, the adding-doubling method achieves a reduction of the calculation time by a factor of 400. Due to the short calculation time, the extended adding-doubling method is very useful when fluorescent layers have to be optimized in an iterative process.
© 2012 OSA
The behavior of light propagating through a fluorescent material is the subject of investigation in a number of fields. In the bio-medical field, fluorescence methods are valuable tools in diagnostics [1–3]. In lighting applications, a phosphor mixture is used to create white light from blue LEDs or from uv emission in a gas discharge [4–6]. In photovoltaics, fluorescent coatings are used to obtain a better match between the solar spectrum and the spectral response of the solar cell [7–9].
The optimization of the fluorescent material for lighting or photovoltaic applications is usually done either empirically [4, 10, 11], using Monte Carlo simulations [12–14] or by applying a simplified one dimensional model [15, 16]. The empirical method requires a lot of effort and time and allows usually only one parameter at the time to be investigated. Monte Carlo simulations take up a lot of computation time, since many rays must be traced in order to reduce the stochastic noise within acceptable margins . The 1-D models usually depart from the well-known Lambert-Beer law which is not well suited to take scattering into account.
In this paper, a fast and accurate method to estimate the wavelength dependent angular distribution of light scattered and emitted by a fluorescent material is described. The proposed method is an extension of the existing adding-doubling algorithm for non-fluorescent samples .
Any macroscopic material can be considered as the conjunction of a number of plane parallel layers of the same material. If the reflection and transmission characteristics of one layer are known, the reflection and transmission characteristics of a combination of such layers can be calculated with the adding-doubling method. This method has been applied already in 1862 by Stokes . The transmission and reflection characteristics of a plane parallel “single-scatter” layer can be calculated by solving the radiative transfer equation.
The adding-doubling method is commonly used to determine the angular light distribution of materials with volume scattering [20–22]. Since fluorescence can be interpreted as inelastic volume scattering, it should be possible to extend this method to determine the characteristics of materials showing both elastic and inelastic scattering.
An extended Kubelka Munk theory for fluorescence has been used to determine the reflection and transmission characteristics for the individual layers. However, the Kubelka Munk theory has some restrictions: an ideal diffuse irradiation is required, the flux propagating in the material must be diffuse, Fresnel reflections at the surface are not taken into account and the correlation between the concentration of the diffusing or fluorescent particles and the absorption and scatter constants is not straightforward [23–25]. Rosema et al. explored the possibility of the adding-doubling method for describing fluorescent layers using a two flux approach . This approach returns total reflection and transmission but no angular information about the scattered radiation.
In this paper we exploit the full potential of the adding-doubling algorithm by using an n-flux model and starting from the radiative transfer equation for light propagation in a fluorescent layer to determine the reflection and transmission characteristics of a “single-scatter” layer. The n-flux model subdivides the flux propagating through a material into n fluxes each propagating within a conical segment with polar angle θ. The optical properties defined in the radiative transfer equation, in particularly absorption and scattering coefficient, are linearly dependent on the concentration . This is advantageous if concentration is an optimization parameter. Furthermore, the use of the n-flux model provides angular information on the scattered and converted light.
2. General approach
For the adding-doubling method the incident flux is divided into a number of channels. The channels are conical segments, each labeled with the polar angle θ. The conical segments are represented in Fig. 1 .
Channel 1 contains the light travelling in directions lying between the normal on the sample and polar angle θ1, channel 2 contains the light travelling in directions lying between θ1 and θ2 and so forth. The flux in the positive or the negative Z direction can be represented by a vector, where each element represents the amount of light within the conical segment. The vector contains n elements, equal to the number of channels the flux is divided in.
The reflected flux in each conical segment originating from the incident flux in each other conical segment can be written in matrix notation as given by Eq. (1):
Herein is I+(θj) the incident flux in channel j propagating in the positive Z direction, I-(θi) the reflected flux in channel i, propagating in the negative Z direction and R(θi, θj) gives the reflection from channel j into channel i. In a similar way a transmission matrix can be composed. The angle θ can take values from 0 to 90°.
3.1 Adding-doubling for non-fluorescent slabs
The adding-doubling method for slabs where only scattering and absorption occurs is extensively described in literature [20–22]. If the reflection (R) and transmission (T) matrices of two given slabs at a particular wavelength are known, the new R and T matrices of a slab consisting of the two slabs joined together can be calculated with the adding method. The doubling method is similar as the adding method for the combination of two identical slabs. For one slab the following equations can be written for each wavelength with Eq. (2) and Eq. (3):
Here, Ix± is the flux within the channels in vector notation, the plus and minus sign denote the propagation direction along the positive and negative Z-axis respectively. The subindex represents the location of an interface starting from 0. Rxy respectively Txy is the reflection matrix and transmission matrix for a flux incident on interface x of the slab defined by interface x and y. This is schematically shown in Fig. 2 .
3.2 Adding-doubling for fluorescent slabs
At excitation wavelengths, the regular adding-doubling can be applied to obtain the reflection and transmission characteristics at each individual excitation wavelength interval. At an emission wavelength, the adding-doubling method described above must be adjusted for the contribution of the fluorescence. Indeed, the flux in forward and backward direction at an emission wavelength is also dependent on the incident flux at the excitation wavelengths.
For fluorescent slabs, the relevant wavelength ranges of excitation and emission are subdivided into respectively N wavelength intervals ΔλXi and L wavelength intervals ΔλMj, with i going from 1 to N and j going from 1 to L. Each wavelength region ΔλXi can have a contribution to the flux at emission wavelength interval ΔλMj, as schematically shown in Fig. 3 .
Rxy and Txy are the reflection and transmission matrices at the selected emission wavelength, representing scattering and absorption. The dependence on wavelength has been omitted from the equations for readability reasons. The matrices Rcxy(ΔλXi, ΔλMj) and Tcxy(ΔλXi, ΔλMj) represent the conversion of light from an excitation wavelength interval ΔλXi towards a particular emission wavelength interval ΔλMj for each incident and scattered conical segment.
Similar equations can be written for a second layer 1-2. The reflection and transmission matrices for the joined layer 0-2 are given in Eqs. (10-13). Details about the calculations can be found in appendix A.
4. Initial matrices
4.1 Initial matrices for non-fluorescent slabs
In the previous section, it has been described how the transmission and reflection matrices for a stack of layers are obtained when the initial R and T matrices are known for each layer. In this section we explain how the initial R and T matrices can be determined. For non-fluorescent applications this is already extensively described in literature . A short summary will be given. The initial matrices are calculated for thin “single-scatter” layers by solving the radiative transfer equation.
It is assumed that the sample is a uniform infinite plane parallel slab of finite thickness. The incident beam and the scattering behavior of the material are assumed to be azimuthally symmetric. The calculations start from the time-independent azimuthally-averaged radiative transfer equation in the case of unpolarized radiation in plane-parallel slab geometry without sources given by Eq. (14) :
With I the flux within a selected channel, ν the cosine between the average propagation direction in the selected channel and the normal on the slab (ν can take values from −1 to 1). P(ν,ν’) is the azimuthally averaged phase function (normalized to 1), µa and µs are respectively the absorption and scattering coefficient. The phase function models the spatial redistribution of scattered light.
A commonly used phase function is the Henyey-Greenstein phase function . This phase function was originally developed to model diffuse radiation in the galaxy, but it is also commonly used for medical  and biological  applications. The Henyey-Greenstein phase function PHG(Υ) is given in Eq. (15), where Υ is the cosine of the angle between incident and scattered direction.
Herein g represents the anisotropy factor. Using the Henyey-Greenstein phase function as a fixed phase function for all materials, a non-fluorescent material can be described with three optical parameters: µa, µs, g.
Wiscombe et al. extensively described the solution of the radiative transfer equation given in Eq. (14) . After discretization and integration over a slab 0-1 with thickness Δd, Eq. (14) can be rewritten as Eq. (16) and Eq. (17):
With νi the value of ν in conical segment i. The weights in matrix c are determined by the used quadrature scheme. Prahl et al. described the use of a Radau and Gaussian quadrature scheme for the discretization . Further manipulation of Eqs. (16) and (17) and identification with Eqs. (2) and (3) results in the initial matrices given by Eq. (23) and Eq. (24) :
With Γ defined by Eq. (25):
4.2 Initial matrices for fluorescent slabs
At the excitation wavelength and other non-emission wavelengths, the initial R and T matrices for a fluorescent slab are calculated as described above for non-fluorescent applications. For the calculation of the initial R and T matrices at an emission wavelength interval, an additional contribution of the fluorescence has to be taken into account .
Equation (26) gives the radiative transfer equation for light within the emission wavelength interval ΔλMj with central wavelength λMj.
In Eq. (26) I(ν) is the channel flux at the selected emission wavelength interval, I(ν’, ΔλXi) represents the channel flux at the excitation wavelength interval ΔλXi. The third term on the right hand side of Eq. (26) represents the contribution of each excitation wavelength interval ΔλXi (with central wavelength λXi) to the selected emission wavelength interval.
In Eq. (26) µa and µs are respectively the absorption and scattering coefficient for the selected emission wavelength interval, µa(ΔλXi) is the absorption coefficient at excitation wavelength interval ΔλXi. QE(ΔλXi) expresses the quantum efficiency the excitation wavelength interval. The wavelength ratio is used to convert the number of photons to powerflux.
Since, the emission of the fluorescence is isotropic, the value for the phase function for fluorescence has been taken as 1/2 . Only a part of the fluorescent light generated by each excitation wavelength interval ΔλXi will be emitted within the selected wavelength interval ΔλMj (schematically shown in Fig. 5 ). A weight factor w(ΔλMj) has been introduced to describe the fraction of light generated by an excitation at ΔλXi and emitted within the selected emission wavelength region. The coefficient w(ΔλMj) is normalized so that the integration over the entire emission wavelength region is 1 with Eq. (27):
The same steps to solve the radiative transfer equation at non-emission wavelengths have been used to solve Eq. (26). After discretization and integration over a slab 0-1 with thickness Δd, Eq. (26) can be rewritten as Eq. (28) and Eq. (29) (see appendix B):
The matrix px is defined by Eq. (31):
T(ΔλXi) and R(ΔλXi) are the reflection and transmission matrices calculated at the excitation wavelength region ΔλXi, calculated with the method for non-fluorescent slabs. The matrix Γ is still defined by Eq. (25). A more extensive description of the derivation of the equations can be found in appendix B.
Similar steps can be followed to calculate R01, T10, Rc01(ΔλXi, ΔλMj) and Tc10(ΔλXi, ΔλMj), the results will be the same as Eqs. (32-34) since the slab will behave identical for light entering the slab in the opposite direction.
The initial matrices R10 and T01, describing the behavior of incident light at the emission wavelength, are identical to the matrices in defined Eqs. (23) and (24). Indeed, fluorescence is independent on the absorption and scattering characteristics at the emission wavelengths. The Rc and Tc matrices describing the contribution due to fluorescence are equal for transmission and reflection. This meets the expectations since the emission due to fluorescence is assumed to be isotropic in the single-scatter layer. If there is no fluorescence (QE = 0), the matrices with superindex c become zero and the equations become equal to the equations valid for the non-fluorescent case. The formulas for the non-fluorescent case can thus be considered as special cases of the general formulas taking fluorescence into account.
To validate the extended adding-doubling method for fluorescent applications, the predictions obtained with this method have been compared to simulation results obtained by Monte Carlo ray tracing with TracePro 7.1 (Lambda Research).
The model in the Monte Carlo simulation consists of a slab with thickness of 1 mm and refractive index of 1.5. In the adding-doubling method, the geometry of the material is assumed to be a plane-parallel slab of infinite size with finite thickness. In the Monte Carlo simulation, the size of the sample has been chosen large enough such that the rays escaping through the edges will be negligible. The incident light is chosen perpendicular onto the surface of the slab.
Only absorption by arbitrarily chosen fluorescent particles has been assumed. Gaussian absorption and emission spectra with respectively central wavelength 450 and 600 nm and σ = 12 and 20 nm have been used.
The scattering properties of the fluorescent particle are determined by the scattering coefficient and the asymmetry factor. The spectral characteristics of these parameters have been chosen arbitrarily. The Henyey-Greenstein phase function was used for both calculations and simulations. The spectrum of the incident light has been chosen to be constant between 380 and 780 nm. For the calculations, a wavelength step of 10 nm, corresponding to 41 wavelengths, has been selected .
The angular distribution of the scattered light calculated with the adding-doubling method was validated by calculating the scattered intensities I(θs) for three wavelengths: one excitation wavelength (450 nm), one emission wavelength (600 nm), and a wavelength not participating in the fluorescence (750 nm).
For the adding-doubling calculations, 32 internal channels were used, 16 of them give rise to light escaping from the sample. The light within the other 16 channels is internally reflected. In Fig. 6 , the scattered intensities for incident angle θi = 0° obtained with the Monte Carlo simulations and with the adding- doubling method are shown for the three selected wavelengths.
In Fig. 6 a good agreement between the results obtained from the Monte Carlo simulations and those calculated with the extended adding-doubling method can be observed. The average deviation between the two methods is 2.7%. However, the total time required for the Monte Carlo simulations, tracing 500 000 rays per wavelength, was over 3 hours. Using this number of rays per wavelength, the relative standard deviation of the scattered intensity distribution due to statistical noise was approximately 1%.
The adding-doubling calculations, using 32 channels, took approximately 25 seconds. The extended adding-doubling code was implemented in Matlab 7.10 and run on the same PC as the TracePro simulations.
In Fig. 7 the total hemispherical reflection and transmission is shown for each wavelength. The deviation between the Monte Carlo simulations and the adding-doubling calculations is found to be within 2%.
By extending the theory of the adding-doubling method, a fast and accurate method to obtain the spectral and angular transmission and reflection characteristics of a fluorescent slab has been developed. With this method, accurate reflection and transmittance profiles for infinite plane parallel slabs with finite thickness, optically flat surfaces and an azimuthally symmetric incident beam can be obtained with the same accuracy, but over 400 times faster than with Monte Carlo simulations.
This extended adding-doubling method can be used to get a quick yet accurate estimate of the transmission and reflection characteristics of a simplified sample geometry for fluorescent applications. Due to the fast algorithm, the extended adding-doubling for fluorescence is particularly suitable for optimization purposes.
Rxy and Txy are the reflection and transmission matrices at the selected emission wavelength, representing scattering and absorption. The dependence on wavelength has been omitted from the equations for readability reasons.
The matrices Rcxy(ΔλXi, ΔλMj) and Tcxy(ΔλXi, ΔλMj) represent the conversion of light from an excitation wavelength interval ΔλXi towards a particular emission wavelength interval ΔλMj for each incident and scattered conical segment. Similar as with the adding-doubling method at non-emission wavelengths, the same equations can be written for a slab 1-2, given by Eq. (37) and Eq. (38):
I1-(ΔλXi) can be written in function of I2-(ΔλXi) and I0+(ΔλXi) given by Eq. (40):
Rxy(ΔλXi) and Txy(ΔλXi) are the reflection and transmission matrices at excitation wavelength interval ΔλXi. In a similar way, I1+(ΔλXi) can be written in function of I2-(ΔλXi) and I0+(ΔλXi), given by Eq. (41):
After further manipulation we obtain Eq. (43):
The calculation of the initial matrices starts from the radiative transfer equations for light within a particular emission wavelength interval ΔλMj with central wavelength λMj given by Eq. (48):
In Eq. (48) I(ν) is the channel flux at the selected emission wavelength interval, I(ν’, ΔλXi) represents the channel flux at the excitation wavelength interval ΔλXi (with central wavelength λXi). In Eq. (48) µa, µs and P(ν,ν’) are respectively the absorption and scattering coefficient and phase function for the selected emission wavelength interval, µa(ΔλXi) is the absorption coefficient at excitation wavelength interval ΔλXi. QE(ΔλXi) expresses the quantum efficiency the excitation wavelength interval. The wavelength ratio is used to convert the number of photons to powerflux. A weight factor w(ΔλMj) has been to describe the fraction of light generated by an excitation at ΔλXi and emitted within the selected emission wavelength region.
I+ and I- represent a channel flux at the selected emission wavelength interval in vector notation with the angle between the propagation direction and the positive Z-axis respectively smaller and larger than 90°. I+(ΔλXi) and I-(ΔλXi) represent the fluxes at excitation wavelength interval ΔλXi.
<I±1/2> is defined by Eq. (56):
The “diamond initialization” (DI) is used as the thin layer approximation . The DI assumes a linear decrease of the channel flux within the thin layer. Using the DI approximation, Eq. (55) can be rewritten as Eq. (57):
To further solve this equation we use the relationships at excitation wavelength given by Eqs. (2) and (3). For the initial homogeneous layer, the reflection and transmission characteristics are independent on the irradiation direction, thus R10=R01 and T10=T01. Using this identities in Eq. (63), Eq. (64) is obtained:
With Γ defined by Eq. (65):
The authors would like to thank 'Strategic Initiative Materials' in Flanders (SIM) and the Institute for Innovation through Science and Technology in Flanders (IWT) for their financial support through the Sol-Cap project within the SIBO program. Wouter Saeys is funded as a postdoctoral fellow of the Research Foundation – Flanders (FWO).
References and links
2. A. Liebert, H. Wabnitz, N. Zołek, and R. Macdonald, “Monte Carlo algorithm for efficient simulation of time-resolved fluorescence in layered turbid media,” Opt. Express 16(17), 13188–13202 (2008). [CrossRef] [PubMed]
4. I. Seo, J. Jung, B. J. Oh, and K. Whang, “Improvement of luminance and luminous efficacy of mercury-free, flat fluorescent lamp by optimizing phosphor profile,” IEEE Trans. Plasma Sci. 38(5), 1097–1100 (2010). [CrossRef]
5. J. P. You, N. T. Tran, Y. Lin, Y. He, and F. G. Shi, “Phosphor-concentration-dependent characteristics of white LEDs in different current regulation modes,” J. Electron. Mater. 38(6), 761–766 (2009). [CrossRef]
6. J. H. Park and J. H. Ko, “Optimization of the emitting structure of flat fluorescent lamps for LCD backlight applications,” J. Opt. Soc. Korea 11(3), 118–123 (2007). [CrossRef]
7. P. Chung, H. Chung, and P. H. Holloway, “Phosphor coatings to enhance Si photovoltaic cell performance,” J. Vac. Sci. Technol. A 25(1), 61–66 (2007). [CrossRef]
8. W. G. J. H. M. van Sark, “Enhancement of solar cell performance by employing planar spectral converters,” Appl. Phys. Lett. 87(15), 151117 (2005). [CrossRef]
9. E. Klampaftis and B. S. Richards, “Improvement in multi-crystalline silicon solar cell efficiency via addition of luminescent material to EVA encapsulation layer,” Prog. Photovolt. Res. Appl. 19(3), 345–351 (2011). [CrossRef]
10. D. Bera, S. Maslov, L. Qian, J. S. Yoo, and P. H. Holloway, “Optimization of the yellow phosphor concentration and layer thickness for down-conversion of blue to white light,” J. Disp. Technol. 6(12), 645–651 (2010). [CrossRef]
11. R. G. Young and E. G. F. Arnott, “The effect of phosphor coating weight on the lumen output of luorescent lamps,” J. Electrochem. Soc. 112(10), 982–984 (1965). [CrossRef]
13. Y. Shuai, N. T. Tran, and F. G. Shi, “Nonmonotonic phosphor size dependence of luminous efficacy for typical white LED emitters,” IEEE Photon. Technol. Lett. 23(9), 552–554 (2011). [CrossRef]
14. C. C. Chang, R.-L. Chern, C. C. Chang, C.-C. Chu, J. Y. Chi, J.-C. Su, I.-M. Chan, and J.-F. T. Wang, “Monte Carlo simulation of optical properties of phosphor-screened ultraviolet light in a white light-emitting device,” Jpn. J. Appl. Phys. 44(8), 6056–6061 (2005). [CrossRef]
15. D.-Y. Kang, E. Wu, and D.-M. Wang, “Modeling white light-emitting diodes with phosphor layers,” Appl. Phys. Lett. 89(23), 231102 (2006). [CrossRef]
19. G. G. Stokes, “On the intensity of the light reflected from or transmitted through a pile of plates,” Proc. R. Soc. Lond. 11(0), 545–556 (1860). [CrossRef]
20. Z. Zhang, P. Yang, G. Kattawar, H.-L. Huang, T. Greenwald, J. Li, B. A. Baum, D. K. Zhou, and Y. Hu, “A fast infrared radiative transfer model based on the adding–doubling method for hyperspectral remote-sensing applications,” J. Quant. Spectrosc. Radiat. Transf. 105(2), 243–263 (2007). [CrossRef]
22. P. J. Flatau and G. L. Stephens, “On the fundamental solution of the radiative transfer equation,” J. Geophys. Res. 93(D9), 11037–11050 (1988). [CrossRef]
23. H. G. Völz, Industrial color testing, Fundamentals and Techniques (Wiley-VCH, 2001)
25. W.-F. Cheong, S. A. Prahl, and A. J. Welch, “A review on the optical properties of biological tissues,” IEEE J. Quantum Electron. 26(12), 2166–2185 (1990). [CrossRef]
26. A. Rosema, W. Verhoef, J. Schroote, and J. F. H. Snel, “Simulating fluorescence light-canopy interaction in support of laser-induced fluorescence measurements,” Remote Sens. Environ. 37(2), 117–130 (1991). [CrossRef]
27. C. F. Bohren and D. Huffman, Absorption and Scattering of Light by Small Particles (John Wiley, 1983).
28. W. J. Wiscombe, “On initialization, error and flux conservation in the doubling method,” J. Quant. Spectrosc. Radiat. Transf. 16(8), 637–658 (1976). [CrossRef]
29. L. G. Henyey and J. L. Greenstein, “Diffuse radiation in the galaxy,” Astrophys. J. 93, 70–83 (1941). [CrossRef]
30. J. F. Beek, P. Blokland, P. Posthumus, M. Aalders, J. W. Pickering, H. J. C. M. Sterenborg, and M. J. C. Gemert, “In vitro double-integrating-sphere optical properties of tissues between 630 and 1064 nm,” Phys. Med. Biol. 42(11), 2255–2261 (1997). [CrossRef] [PubMed]
31. W. Saeys, M. A. Velazco-Roa, S. N. Thennadil, H. Ramon, and B. M. Nicolaï, “Optical properties of apple skin and flesh in the wavelength range from 350 to 2200 nm,” Appl. Opt. 47(7), 908–919 (2008). [CrossRef] [PubMed]
32. D. Yudovsky and L. Pilon, “Modeling the local excitation fluence rate and fluorescence emission in absorbing and strongly scattering multilayered media,” Appl. Opt. 49(31), 6072–6084 (2010). [CrossRef]