## Abstract

The behaviors of the light propagation in the phosphor play a vital role in determining the optical performance of the phosphor-converted light-emitting diodes (pc-LEDs). In this paper, we presented a general model based on the radiative transfer equation integrated with fluorescence (FRTE) to describe the overall light propagating properties in the phosphor layer in terms of light absorption, strong forward scattering, and fluorescence. The model was established by accounting for general boundary conditions including the LED Lambertian incidence, the diffuse reflection at the substrate/reflector, and the Fresnel reflection at the phosphor-air interface. The spectral element method (SEM) was extended to numerically solve FRTE. The radiant intensity at any location and direction of blue and yellow light was iteratively calculated, in which case the angular properties could be further evaluated. The model was validated by comparing the light extraction efficiency (LEE) and angular correlated color temperature (CCT) calculated by the presented model with the experimental results. Good agreements were achieved between model predictions and measurements with the corresponding maximum deviation of 4.9% and 3.7% for LEE and CCT, respectively. We also conducted a comparison between our model and the previous Kubelka-Munk (KM) theory. It has been revealed that the KM theory may overestimate the phosphor heating due to lacking of the blue light scattering effect.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Phosphor, as a wavelength converter, has been widely applied in the lighting industry, such as fluorescent lamps, phosphor-converted light-emitted diodes (pc-LEDs), phosphor-converted laser diodes (pc-LDs), and so on. Among them, pc-LEDs have almost dominated the lighting market with the advantages of energy saving, long lifetime, and high luminous efficiency [1]. In pc-LEDs, cerium-doped yttrium aluminum garnet (YAG:Ce) phosphor is usually applied to absorb part of the blue light from the GaN LED and then generates white light by mixing the transmitted blue light with the converted yellow light [2,3]. It has been reported that the phosphor parameters and layer structures play an important role in determining the overall performances of pc-LEDs, e.g., luminous efficiency, color rendering index, correlated color temperature (CCT) and so forth [3,4].

Light propagating within the phosphor layer in pc-LED exhibits simultaneous absorption, scattering, and fluorescence characteristics [5]. It is hard to observe the internal light propagation properties by experiments. Researchers have turned to developing theoretical and numerical methods to model the phosphor, including the Monte-Carlo ray-tracing simulations [2,6–8], the diffusion-approximation (DA) method [9,10], and a simplified one-dimensional model based on an extended Kubelka–Munk (KM) theory for fluorescence [11–17]. Monte-Carlo simulations can model the phosphor with high accuracy, but it often takes plenty of computational time and resources because usually a million rays need to be traced to get accurate results [18]. The DA method basically applies the first-order spherical harmonics approximation and it fails to model such as media with significant absorption characteristic [19]. The commonly used KM theory solves the forward and backward scattering light intensity analytically by applying a two-flux approximation [11]. It has been a very fast and efficient method compared to the Monte-Carlo simulations. However, it also has some limitations: (1) the scattering effect characterized by the scattering phase function is not taken into account, making it inapplicable when the scattering is strongly anisotropic; (2) the Fresnel reflection at the interface caused by the refractive index mismatch is not considered; (3) only two fluxes travelling in the thickness direction are considered and the angular performance of the phosphor layer cannot be evaluated [20,21]. Alternatively, an extended adding-doubling method for fluorescent materials has been proposed by Leyre *et al.* [22] to calculate the reflection and transmission characteristics of a combination of phosphor layers. In their method, the incident flux is divided into a number of channels and the lights propagating in the forward and backward direction are calculated for each channel. But the Fresnel reflection is also not taken into account. Hence, there is need to develop a general and efficient model for strongly absorbing and strongly anisotropic scattering phosphor with accurate boundary conditions for actual applications in pc-LEDs.

In essence, these methods are all based on the radiative transfer equation (RTE). Light propagation in absorbing and scattering media has been extensively analyzed [23–28] and governed by the steady state RTE. For fluorescing media, e.g. phosphor, the emission wavelength differs from the excitation wavelength. Fluorescence can be interpreted as inelastic volume scattering [29]. In theory, the elastic and inelastic scattering can be described by an extension of RTE integrated with fluorescence (FRTE), which has been used to model the biological tissues [29]. Therefore, it is feasible to apply the FRTE to model the phosphor.

It has been known that the exact analytical solutions of the partial differential-integral equations of RTE only exist in very simple cases and many numerical solutions are usually applied, including the Monte-Carlo method [23], the finite-volume method [24], the finite-element method [25] and the spectral element method (SEM) [26,27]. Among these methods, SEM possesses the advantages of both spectral and the finite-element methods and has been demonstrated to be a very efficient and accurate method to solve RTE [27]. Hence, it is a natural idea to extend this method to solve FRTE.

In this work, a general model for a phosphor is presented by solving FRTE using SEM. The light absorption, anisotropic scattering, and fluorescence were considered. To model the actual phosphor layer in pc-LEDs, the Lambertian incident light intensity distribution for LED, the diffuse reflection at the substrate/reflector, and the Fresnel reflection at the interface between the phosphor and the surrounding air were taken into account. Through discretization in both spatial domain and angular domain, the light intensity at an arbitrary location in any direction was calculated. Then, the light extraction efficiency (LEE) and angular CCT distribution were obtained, followed by experimental validations. The model presented was also compared with the previous KM theory.

## 2. Governing equations

For YAG:Ce phosphor, which absorbs the incident blue light and re-emits yellow light, the extended governing equations FRTEs can be expressed as follows [29]:

For blue light:

*I*(

**R**,

**Ω**) denotes the radiative intensity at spatial location

**R**and direction

**Ω**;

*κ*

_{a}and

*κ*

_{s}are the absorption and scattering coefficients, respectively; the scripts

*B*and

*Y*represent the blue and yellow light, respectively. It should be mentioned that the spectral distribution is ignored, i.e., only the peak wavelengths for excitation and emission are considered.

*p*(

**Ω, Ω**

^{ʹ}) is the scattering phase function, representing the probability of the energy scattering from direction

**Ω**

^{ʹ}to differential solid angle

*d*Ω around direction

**Ω**[28]. For blue light in Eq. (1), the expression of the left side denotes the gradient of the intensity in the specified direction

**Ω**. The two terms on the right denote the intensity decrement due to absorption and outscattering, and intensity increment due to inscattering, respectively [28]. For yellow light in Eq. (2), the similar expressions can be interpreted by the same characteristics except for the contribution of the converted part from the absorbed blue light. This part is defined according to the steady-state rate equation [30], and characterized by the conversion efficiency

*η*

_{con}and the absorption coefficient for blue light. It should be noted that

*η*

_{con}accounts for the phosphor quantum efficiency

*η*

_{qm}as well as the Stokes Shift loss, and it can be expressed as

*η*

_{con}=

*η*

_{qm}·

*λ*

_{Y}/

*λ*

_{B}, where

*λ*

_{B}and

*λ*

_{Y}are wavelengths for excitation blue light and emission yellow light, respectively. It is assumed that the converted yellow light is emitted isotropically [29].

## 3. Model definition

Figure 1 illustrates the schematic of a one-dimensional phosphor model and boundary conditions. It is assumed that the phosphor layer has an infinite length in *x*- and *y*- direction and only thickness (*z*-) direction is considered. Hence, the radiant intensity *I*(**R**, **Ω**) for blue and yellow light to be solved can be simplified as *I*(*z*, **Ω**) or *I*(*z*, *θ*, *φ*), in the units of W·sr^{−1}, where *θ* (0 ≤ *θ* ≤ π) and *φ* (0 ≤ *φ* ≤ 2π) denote the polar and azimuth angle of the solid angle **Ω**, respectively.

For the external LED light source (boundary 1), the typical Lambertian intensity distribution is applied, which can be written as:

*P*

_{in}is the total output optical power from LED and

*θ*denotes the polar angle of the solid angle

**Ω**with respect to z-axis. Here,

**n**

*is the unit inward normal vector at the substrate boundary (*

_{w}*z*= 0) and the incident light travels inward, which satisfies

**Ω**·

**n**

*> 0. It needs to be noted that the arbitrary incident intensity distribution with respect to*

_{w}*θ*and

*φ*can also be easily applied.

When light hits the substrate or the reflector, part of light is absorbed and the other is reflected. For simplicity, the diffuse reflection for both blue and yellow light at the substrate is applied (boundary 2) [12]:

*γ*

_{B}and

*γ*

_{Y}are the diffuse reflectivity at substrate surface for blue and yellow light, respectively.

Considering the refractive index mismatch of the phosphor layer and the surrounding air, light will be specularly reflected back into the layer when attempting to travel out through the phosphor-air boundary [29,31]. However, in previous models [12–17], the refractive indexes of the phosphor and the air are considered identical and equal to one. In this context, light will travel through the layer without deflection and the obtained angular intensity distribution may be misleading. In order to model the accurate condition, by assuming the phosphor-air interface to be a smooth and optically flat surface, the partly-reflecting boundary for blue and yellow light is applied at the phosphor-air interface (boundary 3) [31]:

*d*is the thickness of the phosphor plate,

**n**is the unit inward normal vector at the phosphor-air interface and

**Ω**is the specular reflection of

**Ω**

^{ʹ}, which points outward (

**Ω**

^{ʹ}·

**n**< 0) and satisfies

**Ω**

^{ʹ}=

**Ω**-2(

**Ω**·

**n**)

**n**. The reflectivity

*R*(

**Ω**

^{ʹ}·

**n**) in Eq. (5) is given by [31]:

*θ*

^{ʹ}and the refracted angle

*θ*

_{t}satisfy cos

*θ*

^{ʹ}=

**Ω**

^{ʹ}·

**n**and

*n*

_{p}·sin

*θ*

^{ʹ}=

*n*

_{0}·sin

*θ*

_{t}(Snell’s Law), respectively. Here,

*n*

_{p}and

*n*

_{0}are the refractive index of the phosphor layer and the air. The critical angle

*θ*

_{c}for total internal reflection (

*R*= 1) is defined as

*n*

_{p}·sin

*θ*

_{c}=

*n*

_{0}.

As a whole, the boundary conditions for both blue and yellow light at input and output location can be respectively given by:

To solve the governing Eqs. (1) and (2) together with the boundary conditions Eqs. (3)–(8), the SEM based on the discrete-ordinates equation (DOE) is applied. The implementation of SEM for solving RTE in an absorbing and scattering media has been extensively discussed in the previous literature [27]. A short summary of SEM will be given and extended to solve FRTE.

Figure 2 shows the flowchart of solving FRTE using SEM. First, the input parameters appearing in the governing equations and boundary conditions need to be obtained. Then, the solution domain is discretized using spectral element approximation, and the solution nodes *z*_{i} and corresponding basis functions *ξ*_{i} can be generated using Chebyshev points and polynomials, respectively. Next, the total solid angle (4π sr) is divided uniformly in *θ* and *φ* directions and the corresponding weight *ω* for each angle is derived. This angle discretization method is termed the piecewise constant angular (PCA) approximation and the details can be found in [32]. In the following, the DOEs are built for one-dimensional FRTE described by Eqs. (1) and (2):

*m*denotes a specific angular direction

**Ω**

^{m},

*M*is the total angular number,

*μ*

^{m}is the direction cosine along

*z*direction and

*ω*

^{m}is the angular weight corresponding to the incident direction

**Ω**

^{m}. For simplicity, the in-scattering terms, including the elastic and inelastic scattering, are treated implicitly as part of the source term. In this way, the system of partial differential-integral equations are transformed to a partial differential equations [27]. By using global intensity approximation and integrating Eqs. (9) and (10) over the solution domain, the discretized system of linear equations is obtained as

**K**

^{m}·

**I**

^{m}=

**H**

^{m}. The global approximation intensity

*I*

_{g}can be expressed as:where

*I*

_{j}represents the radiative intensity at solution node

*j*and

*N*

_{sol}is the total number of solution nodes. The detailed expression and assembling of the global stiffness matrix

**K**

^{m}and the vector

**H**

^{m}have been demonstrated in [27].

Then, begin to loop each angular direction for blue light and impose boundary conditions in Eq. (7). The radiative intensity for blue light *I*_{B}(*z*, **Ω**^{m}) on each solution node *z* for each direction **Ω**^{m} can be obtained by solving Eq. (9). If the stop criterion (the maximum relative error of the radiative intensity below 10^{−4}) is satisfied, the iteration process terminates. Otherwise, go back to the loop. After solving the blue light intensity, the same steps can be repeated for yellow light *I*_{Y}(*z*, **Ω**^{m}) by solving Eq. (10) with boundary conditions in Eq. (8).

Finally, the optical performance of the phosphor can be extracted using the obtained radiative intensity of blue and yellow light. The radiant flux for both blue and yellow light can be derived as the integration over the total solid angle in each solution node [19]:

When simulating a phosphor material, the optical performance metrics are always important. When the incident angle *θ* at the phosphor-air interface is below the critical angle *θ*_{c}, the light will transmit the phosphor layer with the transmitted angle *θ*_{t} and the intensity is given by:

In addition, the angular yellow to blue intensity ratio distribution can be applied to determine the angular correlated color temperature (CCT) distribution [33]. The output optical power of blue and yellow light can also be derived as:

*M*

^{’}and the corresponding angular weight is

*ϖ*obtained by the linear interpolation of the original weight

*ω*. Similarly, the ratio of the

*P*

_{out,Y}to

*P*

_{out,B}can also determine the average CCT. The light extraction efficiency (LEE) of the phosphor layer can be determined by the ratio of the total output power

*P*

_{out}to the input power

*P*

_{in}:Moreover, based on the energy conservation law, the phosphor heating power

*Q*

_{ph}can be further obtained by [34]:Considering the light-to-heat conversion process, the heat generation density versus invasion depth

*q*(

*z*) within the phosphor layer can be derived as the sum of the heat generated by the absorbed blue and yellow light:

*q*(

*z*) can be coupled with the heat diffusion equation to further solve the temperature distribution of the phosphor layer.

## 4. Experiments

To validate the one-dimensional model, a remote pc-LED module was fabricated. As shown in the inset in Fig. 3, a commercial LED module bonded onto a substrate was assembled with a designed aluminum reflector with a reflectivity of 0.8 and then a planar phosphor plate was placed onto the top of the reflector. We fabricated two group of phosphor samples with varying thickness of 0.6 mm and 1.0 mm. respectively. For each thickness, there are six different phosphor concentrations varying from 0.05 g/cm^{3} to 0.30 g/cm^{3} with an interval of 0.05 g/cm^{3}. The optical power and average CCT of the module without and with phosphor plate were measured by an integrating sphere (ATA-1000, Everfine Inc.). The corresponding output optical power under driving current of 350 mA were *P*_{opt,1} and *P*_{opt,2}, respectively. The LEE could be derived as the ratio of the *P*_{opt,2} to *P*_{opt,1}. It should be noted that each module was measured for ten times at the same condition. The error was obtained by calculating the average value and standard deviation of these data. In addition, the angular CCT distribution was measured experimentally as depicted in Fig. 3. The pc-LED module was fixed at a designed turntable, which can rotate from 0° to 360°. A spectrophotometer (XYC-I, Enci Co. Ltd.) was applied to measure the color coordinates and CCTs. Before test, the position of the spectrophotometer needs to be adjusted horizontal to the center of the LED. And then the spectrophotometer is fixed at the selected point, which was 1 m away from the module. During test, the CCT distribution was measured by rotating the turntable from 0° to 90° with an increment of 5°.

## 5. Results and discussion

#### 5.1 Calculated results

To implement the model, it is necessary to obtain accurate estimates for the input parameters. The absorption and scattering coefficients of both blue and yellow light are calculated by Mie–Lorenz theory [35]. Figures 4(a) and 4(b) plot the obtained absorption and scattering coefficients versus varying phosphor concentration for both blue and yellow light, respectively. The phosphor quantum efficiency *η*_{qm} is assumed to be 0.9 and the individual wavelength of blue and yellow light are measured to be 445 nm and 558 nm, respectively, corresponding to the conversion efficiency *η*_{con} of 0.718. The diffuse reflectivity *γ* at *z* = 0 for both blue and yellow light is assumed equal to that of the used reflector. The refractive index *n*_{p} of the phosphor plate is set to be 1.53 [3]. In this work, we use the Henyey-Greenstein (HG) phase function [36], which has been widely applied for medical [37] and biological [38] applications as well as the phosphor [39]:

**Ω**·

**Ω**

^{ʹ}. The anisotropy parameter

*g*is set to be 0.82 for YAG:Ce phosphor calculated by Mie–Lorenz theory [35].

After the definition of parameters, the one-dimensional solution domain is discretized with *N*_{sol} Chebyshev-Gauss-Lobatto points. By using the PCA method, the total solid angle is discretized uniformly with *N*_{θ} and *N*_{φ} discrete angles in the polar and azimuth direction, respectively. The implementation of the model was conducted in the commercial software MATLAB using a 3.3-GHz Intel(R) Core(TM) i3-3220 processor. Figure 5 shows the effects of those three discrete numbers on computing time and the accuracy of the results. It should be noted that the accuracy was indicated by the change of the radiant flux at *z* = 0 and *z* = *d*, i.e., Φ(0) and Φ(d), for blue and yellow lights. It can be easily understood that the as the discrete number increases, the computing time keeps rising and the radiant flux tends to be stable. From the figure, we can see that the number of the discrete angles poses a larger effect. As long as *N*_{θ} and *N*_{φ} are respectively greater than 20 and 30, the result can be regarded accurate. Moreover, even when the total number of discrete numbers reaches 1200 (*N*_{θ} × *N*_{φ} = 60 × 20), the computing time is only 45 seconds, which is more efficient than the Monto-Carlo ray-tracing simulations. In the following, we choose those three numbers as *N*_{sol} = 40, *N*_{θ} = 40, and *N*_{φ} = 40. The calculated results will be presented and discussed.

Figures 6(a) and 6(b) show the calculated angular intensity distribution at varying normalized invasion depth *z*/*d* of 0, 0.5 and 1.0 for blue and yellow light, respectively. Considering the symmetry in the azimuth direction of the optical system, only the angular intensities in the polar direction are plotted. It should be noted that *θ* varies from 0° to 180° and the left part of these two figures is just symmetric with the right part. The angle range from 0° to 360° is just labeled to match the usual cases.

For forward-scattering light (*θ* < 90°), the yellow intensity distribution at the incident surface (*z*/*d* = 0) is isotropic, representing the diffuse reflection boundary. However, for blue light, it is caused by the combination of Lambertian incidence and diffuse reflection. With the invasion depth increasing, the blue light intensity decreases in each direction in the similar shape. However, for the yellow light, intensity rises with depth at small angle range, but rises first then drops at high angle range (e.g., *θ* < 75°). For backward-scattering light (*θ* > 90°), the developing trends with depth are just the opposite with that of forward-scattering blue and yellow light. It can be seen that there is a sudden intensity change for both lights. This behavior is related to the partially reflective boundary and can be explained by the plotted angular specular reflectivity at *z*/*d* = 1 shown in Fig. 6(c). When the incident angle *θ* is below the critical angle (*θ*_{c} *=* 41° for *n*_{p} = 1.53), the reflectivity rises very fast as *θ* grows. Otherwise, the reflectivity equals to 1. In this case, the sudden change in backward-scattering intensity at z/d = 1 can be interpreted by the product of the forward-scattering intensity and the angular reflectivity, with a corresponding critical angle of 139°. Figure 6(d) shows the output angular intensity of blue and yellow light obtained by Eq. (13). It can be seen that the blue and yellow light have different intensity distribution. Moreover, the relative intensity distribution is used to determine the CCT distribution to compare with the experimental results in the following.

With the obtained angular intensities for blue and yellow lights in Figs. 6(a) and 6(b), the radiant flux Φ(*z*) and the normalized phosphor heat generation density *q*(*z*)/*P*_{in} can be further calculated using Eqs. (12) and (17), respectively. Figure 7 shows the calculated results. It can be seen that the radiant flux decreases and increases with invasion depth for blue and yellow lights, respectively. It is attributed to the constant absorption of blue light and conversion of yellow light as lights travel in the layer. As for the light-to-heat characteristic, we can see that *q*(*z*)/*P*_{in} shows a monotonic decrease with invasion depth. Because the absorption of blue light accounts for the main part of the heat generation (i.e., *κB a* >> *κY a*), *q*(*z*) has a similar trend with Φ(*z*). This finding is consistent with our previous work [40]. It is worth noting that the temperature field of the phosphor can be obtained by inputting *q*(*z*) into the heat diffusion equation. This coupled radiative and conductive heat transfer problems have been successfully resolved [41,42].

#### 5.2 Comparisons with experiments

Figure 8 shows the comparisons between the experimental results and the model predictions. It should be noted that the fabricated phosphor plate has the diameter of 30 mm and thickness of 0.6 mm and 1.0 mm. In this context, the one-dimensional assumption is feasible. Figure 8(a) plots the comparisons of light extraction efficiency under phosphor plate thickness of 0.6 mm and 1.0 mm versus phosphor concentration varying from 0.05 to 0.30 g/cm^{3} with an interval of 0.05 g/cm^{3}. For both cases, LEE decreases with concentration due to more absorption of light and conversion to heat [17]. It can be seen that LEE also decreases with the thickness at each concentration. For comparison, the LEEs of experiment and model show good agreement in terms of developing trend and absolute value with the maximum deviation of 4.9%.

In addition, we also compare the normalized phosphor heating power *Q*_{ph}/*P*_{in} calculated by our model and KM theory. As shown in Fig. 8(b), the obtained *Q*_{ph}/*P*_{in} by KM are obviously higher than that of our model. It is mainly because that the Lambertian incident light pattern of LED and blue light scattering effect inside the phosphor layer cannot be considered in the KM theory. In this case, the light energy in other directions may not contribute to the light extraction and be dissipated into heat. For further comparison, the angular performance is also compared between experiment and our model. Figures 8(c) and 8(d) plot the comparisons of the angular CCT distribution versus view angle under phosphor concentration of 0.10 g/cm^{3} and 0.15 g/cm^{3}, with corresponding average CCTs of 6643 K and 4551 K, respectively. As a whole, the calculated results agree well with the measured results for both cases, despite the deviation in the low view angle range (0° to 15°) and high view angle range (70° to 90°) with the maximum deviation of 3.7%. The main sources of the deviations may lie in that the input parameters, including the absorption, scattering coefficient, and the phase function, don’t match perfect with the real case. Moreover, the smooth surface assumption may also cause deviation. As a whole, our model is accurate enough to model the light propagation within the phosphor.

## 6. Model applicability and limitations

The presented one-dimensional model can be applied to semi-infinite and uniform phosphor layer with smooth surfaces exposed to LED Lambertian light or other incident light patterns. In addition, this model is feasible for other types of media with characteristics of absorption, scattering and fluorescence. However, there are also some limitations. First, the blue and yellow lights are represented by the corresponding peak wavelengths. To model the actual spectrum, the wavelength-dependent intensity need to be further considered. It may be feasible to discretize the full spectrum into a number of wavelength intervals. And then solve the radiative intensity for each individual wavelength. In this case, the re-absorption effects can be easily included. However, this method may greatly decrease the computational efficiency. This task need to be resolved in the future work. Second, for most commercial pc-LEDs fabricated by direct-coating process, the phosphor is usually in hemispherical shape [3]. In this case, 2-D and 3-D model should be developed in the future work. Third, the surface of the phosphor layer is not actually smooth due to mechanical processing, and the surface topography can be further included to obtain accurate angular properties of phosphor [43]. Alternatively, if a BRDF/BSDF data set is provided, the surface boundary can be imposed accurately by defining the angular reflected and transmitted intensity.

## 7. Summary and conclusion

This work presented a general phosphor model based on the extended radiative transfer equation integrated with fluorescence. The overall light propagating properties were described by taking into account of the light absorption, anisotropic scattering, and fluorescence, simultaneously. We also established accurate boundary conditions by considering the LED Lambertian incident light intensity distribution, the diffuse reflection at the substrate/reflector, and the Fresnel reflection at the phosphor-air interface. The radiant intensity at any location and direction was solved iteratively using spectral element method. For validation, the calculated light extraction efficiency and angular correlated color temperature were verified by experiments and good agreements were reached between the measured and predicted results. The model presented was also compared with the KM theory and it was found that the KM theory may overestimate the phosphor heating due to lacking of the blue light scattering effect. Furthermore, this model could be easily coupled with the heat diffusion equation through the obtained heat generation density to solve the temperature distribution within the phosphor layer. Therefore, the presented efficient and accurate model can be very useful in the optimization processes for better optical and thermal performance of phosphor and other fluorescent media.

## Funding

National Natural Science Foundation of China (51625601, 51576078, 51606074); Creative Research Groups Funding of Hubei Province (2018CFA001); National Key Research and Development Program of China (2016YFB0100901, 2016YFB0400804); Ministry of Science and Technology of the People’s Republic of China (2017YFE0100600).

## References and links

**1. **E. F. Schubert, J. K. Kim, H. Luo, and J. Q. Xi, “Solid state lighting–a benevolent technology,” Rep. Prog. Phys. **69**(12), 3069–3099 (2006). [CrossRef]

**2. **C. H. Hung and C. H. Tien, “Phosphor-converted LED modeling by bidirectional photometric data,” Opt. Express **18**(S3Suppl 3), A261–A271 (2010). [CrossRef] [PubMed]

**3. **X. B. Luo, R. Hu, S. Liu, and K. Wang, “Heat and fluid flow in high-power LED packaging and applications,” Pror. Energy Combust. Sci. **56**, 1–32 (2016). [CrossRef]

**4. **Y. Peng, R. X. Li, H. Zhen Chen, H. Li, and M. X. Chen, “Facile preparation of patterned phosphor-in-glass with excellent luminous properties through screen-printing for high-power white light-emitting diodes,” J. Alloys Compd. **693**, 279–284 (2017). [CrossRef]

**5. **A. Lenef, J. F. Kelso, and A. Piquette, “Light extraction from luminescent light sources and application to monolithic ceramic phosphors,” Opt. Lett. **39**(10), 3058–3061 (2014). [CrossRef] [PubMed]

**6. **Z. Liu, K. Wang, X. Luo, and S. Liu, “Precise optical modeling of blue light-emitting diodes by Monte Carlo ray-tracing,” Opt. Express **18**(9), 9398–9412 (2010). [CrossRef] [PubMed]

**7. **S. W. Jeon, J. H. Noh, K. H. Kim, W. H. Kim, C. Yun, S. B. Song, and J. P. Kim, “Improvement of phosphor modeling based on the absorption of Stokes shifted light by a phosphor,” Opt. Express **22**(S5Suppl 5), A1237–A1242 (2014). [CrossRef] [PubMed]

**8. **A. Correia, P. Hanselaer, and Y. Meuret, “An efficient optothermal simulation framework for optimization of high-luminance white light sources,” IEEE Photonics J. **8**(4), 1–15 (2016). [CrossRef]

**9. **A. Lenef, J. Kelso, M. Tchoul, O. Mehl, J. Sorg, and Y. Zheng, “Laser-activated remote phosphor conversion with ceramic phosphors,” Proc. SPIE **9190**, 919000C (2014).

**10. **A. Lenef, J. Kelso, Y. Zheng, and M. Tchoul, “Radiance limits of ceramic phosphors under high excitation fluxes,” Proc. SPIE **8841**, 884107 (2013).

**11. **T. Shakespeare and J. Shakespeare, “A fluorescent extension to the Kubelka–Munk model,” Color Res. Appl. **28**(1), 4–14 (2003). [CrossRef]

**12. **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]

**13. **R. Hu and X. B. Luo, “A model for calculating the bidirectional scattering properties of phosphor layer in white light-emitting diodes,” J. Lightwave Technol. **30**(21), 3376–3380 (2012). [CrossRef]

**14. **R. Hu, B. Cao, Y. Zou, Y. M. Zhu, S. Liu, and X. B. Luo, “Modeling the light extraction efficiency of bi-layer phosphor in white LEDs,” IEEE Photonics Technol. Lett. **25**(12), 1141–1144 (2013). [CrossRef]

**15. **R. Hu, Y. M. Wang, Y. Zou, X. Chen, S. Liu, and X. B. Luo, “Study on phosphor sedimentation effect in white LED packages by modeling multi-layer phosphors with the modified Kubelka-Munk theory,” J. Appl. Phys. **113**(6), 063108 (2013). [CrossRef]

**16. **R. Hu, H. Zheng, J. Y. Hu, and X. B. Luo, “Comprehensive study on the transmitted and reflected light through the phosphor layer in light-emitting diode packages,” J. Disp. Technol. **9**(6), 447–452 (2013). [CrossRef]

**17. **Y. P. Ma, W. Lan, B. Xie, R. Hu, and X. B. Luo, “An optical-thermal model for laser-excited remote phosphor with thermal quenching,” Int. J. Heat Mass Transf. **116**, 694–702 (2018). [CrossRef]

**18. **J. Chen and X. Intes, “Comparison of Monte Carlo methods for fluorescence molecular tomography-computational efficiency,” Med. Phys. **38**(10), 5788–5798 (2011). [CrossRef] [PubMed]

**19. **T. J. Tarvainen, V. Kolehmainen, S. R. Arridge, and J. P. Kaipio, “Image reconstruction in diffuse optical tomography using the coupled radiative transport–diffusion model,” J. Quant. Spectrosc. **112**(16), 2600–2608 (2011). [CrossRef]

**20. **W. E. Vargas and G. A. Niklasson, “Applicability conditions of the Kubelka-Munk theory,” Appl. Opt. **36**(22), 5580–5586 (1997). [CrossRef] [PubMed]

**21. **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]

**22. **S. Leyre, G. Durinck, B. Van Giel, W. Saeys, J. Hofkens, G. Deconinck, and P. Hanselaer, “Extended adding-doubling method for fluorescent applications,” Opt. Express **20**(16), 17856–17872 (2012). [CrossRef] [PubMed]

**23. **J. R. Howell, “Application of Monte Carlo to heat transfer problems,” Adv. Heat Transf. **5**, 1–54 (1968).

**24. **J. C. Chai, “One-dimensional transient radiation heat transfer modeling using a finite-volume method,” Numer. Heat Transf. B **44**(2), 187–208 (2003). [CrossRef]

**25. **L. H. Liu, “Finite element simulation of radiative heat transfer in absorbing and scattering media,” J. Thermophys. Heat Transfer **18**(4), 555–557 (2004). [CrossRef]

**26. **J. M. Zhao and L. H. Liu, “Discontinuous spectral element method for solving radiative heat transfer in multidimensional semitransparent media,” J. Quant. Spectrosc. **107**(1), 1–16 (2007). [CrossRef]

**27. **J. M. Zhao and L. H. Liu, “Least-squares spectral element method for radiative heat transfer in semitransparent media,” Numer. Heat Transf. B **50**(5), 473–489 (2006). [CrossRef]

**28. **W. A. Fiveland, “Three-dimensional radiative heat-transfer solutions by the discrete-ordinates method,” J. Thermophys. Heat Transf. **2**(4), 309–316 (1988). [CrossRef]

**29. **A. D. Klose, “Radiative transfer of luminescence light in biological tissue,” in *Light Scattering Reviews 4*, (Springer, 2009).

**30. **E. M. Sevick-Muraca and C. L. Burch, “Origin of phosphorescence signals reemitted from tissues,” Opt. Lett. **19**(23), 1928–1930 (1994). [CrossRef] [PubMed]

**31. **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]

**32. **W. A. Fiveland and J. P. Jessee, “Comparison of discrete ordinates formulations for radiative heat transfer in multidimensional geometries,” J. Thermophys. Heat Transf. **9**(1), 47–54 (1995). [CrossRef]

**33. **S. McCamy, “Correlated color temperature as an explicit function of chromaticity coordinates,” Color Res. Appl. **17**(2), 142–144 (1992). [CrossRef]

**34. **Y. P. Ma, R. Hu, X. J. Yu, W. C. Shu, and X. B. Luo, “A modified bidirectional thermal resistance model for junction and phosphor temperature estimation in phosphor-converted light-emitting diodes,” Int. J. Heat Mass Transf. **106**, 1–6 (2017). [CrossRef]

**35. **Z. Liu, S. Liu, K. Wang, and X. Luo, “Measurement and numerical studies of optical properties of YAG:Ce phosphor for white light-emitting diode packaging,” Appl. Opt. **49**(2), 247–257 (2010). [CrossRef] [PubMed]

**36. **L. G. Henyey and J. L. Greenstein, “Diffuse radiation in the galaxy,” Astrophys. J. **93**, 70–83 (1941). [CrossRef]

**37. **J. F. Beek, P. Blokland, P. Posthumus, M. Aalders, J. W. Pickering, H. J. C. M. Sterenborg, and M. J. van 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]

**38. **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]

**39. **S. Leyre, G. Durinck, J. Hofkens, G. Deconinck, and P. Hanselaer, “Experimental determination of the absorption and scattering properties of YAG:Ce phosphor,” in Light, Energy and the Environment, OSA Technical Digest (online) (Optical Society of America, 2014), paper DTu4C.4.

**40. **X. B. Luo and R. Hu, “Calculation of the phosphor heat generation in phosphor-converted light-emitting diodes,” Int. J. Heat Mass Transf. **75**, 213–217 (2014). [CrossRef]

**41. **J. M. Zhao and L. H. Liu, “Spectral element approach for coupled radiative and conductive heat transfer in semitransparent medium,” J. Heat Transfer **129**(10), 1417–1424 (2007). [CrossRef]

**42. **J. Sun, H. L. Yi, and H. P. Tan, “Local RBF meshless scheme for coupled radiative and conductive heat transfer,” Numer. Heat Transf. A **69**(12), 1390–1404 (2016). [CrossRef]

**43. **D. H. Lee, J. Y. Joo, and S. K. Lee, “Modeling of reflection-type laser-driven white lighting considering phosphor particles and surface topography,” Opt. Express **23**(15), 18872–18887 (2015). [CrossRef] [PubMed]