# Radiative analysis of luminescence in photoreactive systems: Application to photosensitizers for solar fuel production

Caroline Supplis, Jérémi Dauchet, Victor Gattepaille, Fabrice Gros, Thomas Vourc'h, Jean-François Cornet

<figure><img src="image_1.png" /></figure>

## To cite this version:

Caroline Supplis, Jérémi Dauchet, Victor Gattepaille, Fabrice Gros, Thomas Vourc'h, et al., Radiative analysis of luminescence in photoreactive systems: Application to photosensitizers for solar fuel production. PLoS ONE, 2021, 16 (7), pp.e0255002. 10.1371/journal.pone.0255002. hal-03379559

HAL Id: hal-03379559

https://hal.science/hal-03379559

Submitted on 25 Jun 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Radiative analysis of luminescence in photoreactive systems: Application to photosensitizers for solar fuel production

Caroline Supplis, Jérémi Dauchet, Victor Gattepaille, Fabrice Gros, Thomas Vourc'h, Jean-François Cornet

► To cite this version:

Caroline Supplis, Jérémi Dauchet, Victor Gattepaille, Fabrice Gros, Thomas Vourc'h, et al.. Radia- tive analysis of luminescence in photoreactive systems: Application to photosensitizers for solar fuel production. PLoS ONE, 2021, 16 (7), pp.e0255002. 10.1371/journal.pone.0255002 . hal-03379559

HAL Id: hal-03379559

https://hal.science/hal-03379559v1

Submitted on 25 Jun 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

# Radiative analysis of luminescence in photoreactive systems: Application to photosensitizers for solar fuel production

Caroline Supplis<sup>©*</sup>, Jérémi Dauchet<sup>©*</sup>, Victor Gattepaille<sup>‡</sup>, Fabrice Gros<sup>‡</sup>, Thomas Vourc'h<sup>‡</sup>, Jean-François Cornet<sup>‡</sup>

Université Clermont Auvergne, Clermont Auvergne INP, CNRS, Institut Pascal, F-63000 Clermont-Ferrand, France

© These authors contributed equally to this work.

‡ These authors also contributed equally to this work.

* caroline.supplis@uca.fr (CS); jeremi.dauchet@sigma-clermont.fr (JD)

<figure><img src="image_3.png" /></figure>

## Abstract

Most chemical reactions promoted by light and using a photosensitizer (a dye) are subject to the phenomenon of luminescence. Redistribution of light in all directions (isotropic luminescence emission) and in a new spectral range (luminescence emission spectrum) makes experimental and theoretical studies much more complex compared to a situation with a purely absorbing reaction volume. This has a significant impact on the engineering of photo-reactors for industrial applications. Future developments associated with photoreactive system optimization are therefore extremely challenging, and require an in-depth description and quantitative analysis of luminescence. In this study, a radiative model describing the effect of luminescence radiation on the calculation of absorptance is presented and analyzed with the multiple inelastic-scattering approach, using Monte Carlo simulations. The formalism of successive orders of scattering expansion is used as a sophisticated analysis tool which provides, when combined with relevant physical approximations, convenient analytical approximate solutions. Its application to four photosensitizers that are representative of renewable hydrogen production *via* artificial photosynthesis indicates that luminescence has a significant impact on absorptance and on overall quantum yield estimation, with the contribution of multiple scattering and important spectral effects due to inelastic scattering. We show that luminescence cannot be totally neglected in that case, since photon absorption lies at the root of the chemical reaction. We propose two coupled simple and appropriate analytical approximations enabling the estimation of absorptance with a relative error below 6% in every tested situation: the zero-order scattering approximation and the gray single-scattering approximation. Finally, this theoretical approach is used to determine and discuss the overall quantum yield of a bio-inspired photoreactive system with Eosin Y as a photosensitizer, implemented in an experimental setup comprising a photoreactor dedicated to hydrogen production.

<figure><img src="image_4.png" /></figure>

OPEN ACCESS

**Citation:** Supplis C, Dauchet J, Gattepaille V, Gros F, Vourc'h T, Cornet J-F (2021) Radiative analysis of luminescence in photoreactive systems: Application to photosensitizers for solar fuel production. PLoS ONE 16(7): e0255002. https://doi.org/10.1371/journal.pone.0255002

**Editor:** Linda Ramball Jones, College of Charleston, UNITED STATES

**Received:** March 18, 2021

**Accepted:** July 7, 2021

**Published:** July 22, 2021

**Copyright:** © 2021 Supplis et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the manuscript and its Supporting information files.

**Funding:** This work was sponsored by a public grant overseen by the French National Research Agency as part of the Investissements d'Avenir through the IMobS3 Laboratory of Excellence (ANR-10-LABX-0016) and the IDEX-ISITE initiative CAP 20-25 (ANR-16-IDEX-0001).

**Competing interests:** The authors have declared that no competing interests exist.
RESEARCH ARTICLE

# Radiative analysis of luminescence in photoreactive systems: Application to photosensitizers for solar fuel production

Caroline Supplis<sup>©*</sup>, Jérémi Dauchet<sup>©*</sup>, Victor Gattepaille<sup>‡</sup>, Fabrice Gros<sup>‡</sup>, Thomas Vourc'h<sup>‡</sup>, Jean-François Cornet<sup>‡</sup>

Université Clermont Auvergne, Clermont Auvergne INP, CNRS, Institut Pascal, F-63000 Clermont-Ferrand, France

© These authors contributed equally to this work.

‡ These authors also contributed equally to this work.

* caroline.supplis@uca.fr (CS); jeremi.dauchet@sigma-clermont.fr (JD)

<figure><img src="image_2.png" /></figure>

## Abstract

Most chemical reactions promoted by light and using a photosensitizer (a dye) are subject to the phenomenon of luminescence. Redistribution of light in all directions (isotropic luminescence emission) and in a new spectral range (luminescence emission spectrum) makes experimental and theoretical studies much more complex compared to a situation with a purely absorbing reaction volume. This has a significant impact on the engineering of photo-reactors for industrial applications. Future developments associated with photoreactive system optimization are therefore extremely challenging, and require an in-depth description and quantitative analysis of luminescence. In this study, a radiative model describing the effect of luminescence radiation on the calculation of absorptance is presented and analyzed with the multiple inelastic-scattering approach, using Monte Carlo simulations. The formalism of successive orders of scattering expansion is used as a sophisticated analysis tool which provides, when combined with relevant physical approximations, convenient analytical approximate solutions. Its application to four photosensitizers that are representative of renewable hydrogen production *via* artificial photosynthesis indicates that luminescence has a significant impact on absorptance and on overall quantum yield estimation, with the contribution of multiple scattering and important spectral effects due to inelastic scattering. We show that luminescence cannot be totally neglected in that case, since photon absorption lies at the root of the chemical reaction. We propose two coupled simple and appropriate analytical approximations enabling the estimation of absorptance with a relative error below 6% in every tested situation: the zero-order scattering approximation and the gray single-scattering approximation. Finally, this theoretical approach is used to determine and discuss the overall quantum yield of a bio-inspired photoreactive system with Eosin Y as a photosensitizer, implemented in an experimental setup comprising a photoreactor dedicated to hydrogen production.

<figure><img src="image_3.png" /></figure>

OPEN ACCESS

**Citation:** Supplis C, Dauchet J, Gattepaille V, Gros F, Vourc'h T, Cornet J-F (2021) Radiative analysis of luminescence in photoreactive systems: Application to photosensitizers for solar fuel production. PLoS ONE 16(7): e0255002. https://doi.org/10.1371/journal.pone.0255002

**Editor:** Linda Ramball Jones, College of Charleston, UNITED STATES

**Received:** March 18, 2021

**Accepted:** July 7, 2021

**Published:** July 22, 2021

**Copyright:** © 2021 Supplis et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the manuscript and its Supporting information files.

**Funding:** This work was sponsored by a public grant overseen by the French National Research Agency as part of the Investissements d'Avenir through the IMobS3 Laboratory of Excellence (ANR-10-LABX-0016) and the IDEX-ISITE initiative CAP 20-25 (ANR-16-IDEX-0001).

**Competing interests:** The authors have declared that no competing interests exist.
# 1 Introduction

Molecular homogeneous photoreactive systems (i.e. chemical molecules or clusters of molecules which can interact with light and then react) have many applications and significant industrial interest. We can cite first of all applications in industrial photochemistry using artificial light sources [1]. More recently, new solar applications for green photochemistry [2, 3], the production of renewable hydrogen from water [4–6] or carbon dioxyde photoreduction [7–9] involving artificial photosynthesis (solar fuels) are becoming increasingly important. This last application, which motivated the present work, is a considerable challenge in the 21st century. Having efficient and inexpensive catalysts for the photolysis of water, associated with solar processes designed with high thermodynamic efficiency, is the guarantee of being able to easily store solar energy in chemical vectors.

In the quest for cheap catalyst development, molecular catalysts synthesized from relatively abundant elements on Earth (in contrast with particle catalysts based on doped-semiconductors, which are still more effective today [6, 10–12] but out of the scope of this article), are an attractive prospect for the future [13–15]. However, they must always be associated with one or more photosensitizers, which will initiate the reaction by elementary photon absorption leading to an excited state as a singlet or triplet. The use of these photosensitizers goes far beyond the framework of the production of solar energy vectors, since they are at the root of all the applications mentioned above.

One of the characteristic phenomena encountered in the use of photosensitizers is that electrons in the first excited state can spontaneously relax to the ground state by fluorescence (for a singlet) or phosphorescence (for a triplet) emission without a reactive event. These photoluminescence phenomena are both responsible for light emission in the volume, at the expense of photoreaction. We argue that photoluminescence requires a sound description and quantitative analysis when formulating the thermokinetic coupling law and performing the radiative balance in a photoreactor, for two reasons. First, photons emitted by luminescence at one location in the medium can potentially be absorbed at another location and participate in the reaction. Therefore, photoluminescence cannot be treated simply as a loss for the reaction; it must be addressed in the description of photon transport within the photoreactor. Thus, the considerations of luminescence in transmittance and absorptance calculations or measurements may be of crucial importance as regards the above-mentioned photochemical and artificial photosynthesis applications. Second, photoluminescence is responsible for light emission in all directions (isotropic emission) and with a spectral distribution (emission spectrum) that is different from that of the incident radiation. The radiative balance in the photoreactor then becomes extremely challenging without appropriate modeling and analysis tools.

In more detail, from the engineering point of view, and considering that photoreactive processes are controlled at different scales by radiative transfer (or photon transport) [16], an understanding of the photon absorption process (leading to a reactive event) is essential to the radiative analysis, and thus to the optimization of any photoreactive process at the application scale. It is therefore necessary to gain further insights into this physical phenomenon, which can be quantified properly by determining the local $A(\mathbf{r})$ and spatial (mean) $\langle A \rangle$ volumetric rates of photons absorbed (MVRPA), where we have introduced for convenience the bracket notation $\langle \cdot \rangle = \frac{1}{V} \iiint_V \cdot dV$. This quantification also serves as a basis for the formulation of coupling laws between the photon absorption rate and the local chemical reaction rate $r_i(\mathbf{r})$. The physical or engineering treatment of any photoreactive system relies on a local thermokinetic coupling involving an overall quantum yield $\varphi$ with the general form [16, 17]:

$$r_i(\mathbf{r}) = \varphi A(\mathbf{r}) \qquad (1)$$
# 1 Introduction

Molecular homogeneous photoreactive systems (i.e. chemical molecules or clusters of molecules which can interact with light and then react) have many applications and significant industrial interest. We can cite first of all applications in industrial photochemistry using artificial light sources [1]. More recently, new solar applications for green photochemistry [2, 3], the production of renewable hydrogen from water [4–6] or carbon dioxyde photoreduction [7–9] involving artificial photosynthesis (solar fuels) are becoming increasingly important. This last application, which motivated the present work, is a considerable challenge in the 21st century. Having efficient and inexpensive catalysts for the photolysis of water, associated with solar processes designed with high thermodynamic efficiency, is the guarantee of being able to easily store solar energy in chemical vectors.

In the quest for cheap catalyst development, molecular catalysts synthesized from relatively abundant elements on Earth (in contrast with particle catalysts based on doped-semiconductors, which are still more effective today [6, 10–12] but out of the scope of this article), are an attractive prospect for the future [13–15]. However, they must always be associated with one or more photosensitizers, which will initiate the reaction by elementary photon absorption leading to an excited state as a singlet or triplet. The use of these photosensitizers goes far beyond the framework of the production of solar energy vectors, since they are at the root of all the applications mentioned above.

One of the characteristic phenomena encountered in the use of photosensitizers is that electrons in the first excited state can spontaneously relax to the ground state by fluorescence (for a singlet) or phosphorescence (for a triplet) emission without a reactive event. These photoluminescence phenomena are both responsible for light emission in the volume, at the expense of photoreaction. We argue that photoluminescence requires a sound description and quantitative analysis when formulating the thermokinetic coupling law and performing the radiative balance in a photoreactor, for two reasons. First, photons emitted by luminescence at one location in the medium can potentially be absorbed at another location and participate in the reaction. Therefore, photoluminescence cannot be treated simply as a loss for the reaction; it must be addressed in the description of photon transport within the photoreactor. Thus, the considerations of luminescence in transmittance and absorptance calculations or measurements may be of crucial importance as regards the above-mentioned photochemical and artificial photosynthesis applications. Second, photoluminescence is responsible for light emission in all directions (isotropic emission) and with a spectral distribution (emission spectrum) that is different from that of the incident radiation. The radiative balance in the photoreactor then becomes extremely challenging without appropriate modeling and analysis tools.

In more detail, from the engineering point of view, and considering that photoreactive processes are controlled at different scales by radiative transfer (or photon transport) [16], an understanding of the photon absorption process (leading to a reactive event) is essential to the radiative analysis, and thus to the optimization of any photoreactive process at the application scale. It is therefore necessary to gain further insights into this physical phenomenon, which can be quantified properly by determining the local $A(\mathbf{r})$ and spatial (mean) $\langle A \rangle$ volumetric rates of photons absorbed (MVRPA), where we have introduced for convenience the bracket notation $\langle \cdot \rangle = \frac{1}{V} \iiint_V \cdot dV$. This quantification also serves as a basis for the formulation of coupling laws between the photon absorption rate and the local chemical reaction rate $r_i(\mathbf{r})$. The physical or engineering treatment of any photoreactive system relies on a local thermokinetic coupling involving an overall quantum yield $\varphi$ with the general form [16, 17]:

$$r_i(\mathbf{r}) = \varphi A(\mathbf{r}) \qquad (1)$$
Local coupling is required because the local volumetric rate of photons absorbed (LVRPA) field $A(r)$ is highly heterogeneous when most of the incident radiation is absorbed in the vessel, which is the aim of photoreactive processes. Nevertheless, if the quantum yield $\phi$ does not depend on the radiation field, which is often the case for molecular photoreactive systems as studied in this article, averaging these two rates at the spatial scale of the enclosure (the scale of physical observable and measured quantities) is straightforward and leads to a linear thermokinetic coupling law [18, 19]:

$$ \langle r_i \rangle = \phi \langle A \rangle \qquad (2) $$

In these situations, the kinetic and thermodynamic formulations of a knowledge model for any photoreactor geometry are fully provided by a thorough evaluation of the mean spatial volumetric rate of radiant energy absorbed $\langle A \rangle$ [16, 19, 20]. This description of luminescence has two important consequences in the thermokinetic coupling law (Eq 2):

* Luminescence cannot be considered as a loss in the thermokinetic coupling. Thus the luminescence quantum yield is not included in the overall quantum yield $\phi$,

* Luminescence is included in the radiative model. It enables us to account for the multiple absorption and emission of luminescent radiation when estimating $\langle A \rangle$.

Even though the present study is of great interest for all applications using homogeneous (molecular) photoreactive systems, here we focused on four standard photosensitizers. They were chosen because they are widely used for the production of solar fuels, in association with a suitable catalyst and a sacrificial electron donor. The luminescent molecules studied in this paper are either fluorescent molecules such as Eosin Y [21–24], derived from triazatriangulennium (TATA$^+$) [25] and Rhodamine B [26, 27] or phosphorescent molecules such as tris(2,2’-bipyridine)ruthenium (Ru[Bpy]$^{2+}$) [28, 29]. They were treated using a unified approach, requiring only basic information (experimental data or DFT calculations, for example) such as spectral extinction coefficient, luminescence emission spectrum and luminescence quantum yield.

Regarding the application, we considered a homogeneous reaction medium with luminescent species implemented in a photoreactor modeled as a slab whose radiative configuration and boundary conditions are detailed in section 2.1. In particular, heterogeneous photoreactive systems with particles leading to elastic scattering are not discussed in this article. First we will focus on situations where the luminescent photosensitizer is the only species that interacts with radiation—there is no other absorbing species—and results will subsequently be extended to configurations where the catalyst also absorbs radiation. Extension is quite straightforward and the conclusions are not modified, since absorption by catalyst is not desirable (it does not lead to reaction) and is therefore low when working with efficient molecular photoreactive systems.

The corresponding steady-state Radiative Transfer Equation (RTE) will be presented in section 2.2 and is rigorously solved using the Monte Carlo method. Efforts are then essentially focused on the analysis of luminescence effects on MVRPA, and as a result on the construction of relevant analytical approximate solutions. The collision term that describes luminescence in the radiative transfer equation can either be interpreted as an emission source in the volume [30, 31], or as an inelastic and isotropic scattering phenomenon [32–34]. The second approach enables an intuitive understanding of the photon transport physics and was chosen here. Moreover, we took the decision to describe the MVRPA solution as successive orders of scattering for several reasons: i) it provides intuitive physical pictures when describing the absorption of luminescence radiation, ii) it is a common method used in radiative analysis to develop
Local coupling is required because the local volumetric rate of photons absorbed (LVRPA) field $A(r)$ is highly heterogeneous when most of the incident radiation is absorbed in the vessel, which is the aim of photoreactive processes. Nevertheless, if the quantum yield $\phi$ does not depend on the radiation field, which is often the case for molecular photoreactive systems as studied in this article, averaging these two rates at the spatial scale of the enclosure (the scale of physical observable and measured quantities) is straightforward and leads to a linear thermokinetic coupling law [18, 19]:

$$ \langle r_i \rangle = \phi \langle A \rangle \qquad (2) $$

In these situations, the kinetic and thermodynamic formulations of a knowledge model for any photoreactor geometry are fully provided by a thorough evaluation of the mean spatial volumetric rate of radiant energy absorbed $\langle A \rangle$ [16, 19, 20]. This description of luminescence has two important consequences in the thermokinetic coupling law (Eq 2):

* Luminescence cannot be considered as a loss in the thermokinetic coupling. Thus the luminescence quantum yield is not included in the overall quantum yield $\phi$,

* Luminescence is included in the radiative model. It enables us to account for the multiple absorption and emission of luminescent radiation when estimating $\langle A \rangle$.

Even though the present study is of great interest for all applications using homogeneous (molecular) photoreactive systems, here we focused on four standard photosensitizers. They were chosen because they are widely used for the production of solar fuels, in association with a suitable catalyst and a sacrificial electron donor. The luminescent molecules studied in this paper are either fluorescent molecules such as Eosin Y [21–24], derived from triazatriangulennium (TATA$^+$) [25] and Rhodamine B [26, 27] or phosphorescent molecules such as tris(2,2’-bipyridine)ruthenium (Ru[Bpy]$^{2+}$) [28, 29]. They were treated using a unified approach, requiring only basic information (experimental data or DFT calculations, for example) such as spectral extinction coefficient, luminescence emission spectrum and luminescence quantum yield.

Regarding the application, we considered a homogeneous reaction medium with luminescent species implemented in a photoreactor modeled as a slab whose radiative configuration and boundary conditions are detailed in section 2.1. In particular, heterogeneous photoreactive systems with particles leading to elastic scattering are not discussed in this article. First we will focus on situations where the luminescent photosensitizer is the only species that interacts with radiation—there is no other absorbing species—and results will subsequently be extended to configurations where the catalyst also absorbs radiation. Extension is quite straightforward and the conclusions are not modified, since absorption by catalyst is not desirable (it does not lead to reaction) and is therefore low when working with efficient molecular photoreactive systems.

The corresponding steady-state Radiative Transfer Equation (RTE) will be presented in section 2.2 and is rigorously solved using the Monte Carlo method. Efforts are then essentially focused on the analysis of luminescence effects on MVRPA, and as a result on the construction of relevant analytical approximate solutions. The collision term that describes luminescence in the radiative transfer equation can either be interpreted as an emission source in the volume [30, 31], or as an inelastic and isotropic scattering phenomenon [32–34]. The second approach enables an intuitive understanding of the photon transport physics and was chosen here. Moreover, we took the decision to describe the MVRPA solution as successive orders of scattering for several reasons: i) it provides intuitive physical pictures when describing the absorption of luminescence radiation, ii) it is a common method used in radiative analysis to develop
analytical approximate solutions [35], and iii) it is convenient to study the luminescence effect on $\langle A \rangle$, as will be presented in section 3.

The overall effect of luminescence on absorptance is quantified, indicating that luminescence can hardly be neglected in applications concerning solar fuel production. We therefore analyzed the respective influence of two characteristics of luminescence: multiple scattering and inelasticity. For this purpose, the weight of the successive scattering orders was first analyzed, and then the effect of inelasticity in scattering was quantified by comparing the reference results with those obtained thanks to equivalent gray or elastic scattering models. From these developments important conclusions describing the impact of luminescence on MVRPA quantification are given, and rules of thumb for reasoning from well-defined optical thicknesses are proposed for any application using photosensitizers. Various useful analytical approximations for MVRPA assessment, also valid for any photosensitized reactive process in the absence of particles, are fully developed and compared to rigorous calculations in section 4 in order to discuss their validity in detail.

Finally, in section 5 the implementation of the proposed method will make it possible to revisit results already published regarding the use of a bio-inspired catalyst for the production of H₂ in photoreactors, where the fluorescence of Eosin Y was neglected [18]. This will enable us to accurately estimate the overall quantum yields $\varphi$ and thus to refine the performances achievable by this kind of photoreactive system.

# 2 Radiative model

## 2.1 Radiative configuration

Photon transport was modeled in a one dimensional slab of thickness L with incident-normal collimated radiation on the left-hand side (see Fig 1). The incident monochromatic photon flux density inside the medium, at $x = 0$, is $q_{0,\lambda} = q_0 p^i(\lambda)$, where $q_0 = \int_0^{+\infty} d\lambda \, q_{0,\lambda}$ is the total photon flux density and $p^i(\lambda)$ is the incident spectral probability density function, hereafter called *incident spectrum*: $\int_0^{+\infty} d\lambda \, p^i(\lambda) = 1$. The interfaces at both sides of the slab are assumed to be transparent in order to analyze the roles of luminescence and separate them from the impact of reflection-refraction at the boundaries.

The slab is filled with a homogeneous photocatalytic system solution including luminescent photosensitizer at concentration C (there is no particle leading to elastic scattering in such systems). In the present configuration the photosensitizer is the only molecule that interacts with

<figure><img src="image_5.png" /><figcaption>**Fig 1. Radiative configuration.** (a) Incident-normal collimated radiation at $x = 0$ with spectrum $p^i(\lambda)$ and surface flux density $q_0$. (b) Slab of thickness L containing a luminescent photosensitizer solution at concentration C. The phase space ($\mathbf{r}$, $\mathbf{u}$), with $\mathbf{u}$ the propagation direction, can be reduced to ($x$, $\mu$) in such a one-dimensional slab. In this case, the isotropic distribution $p(\mathbf{u}) = \frac{1}{4\pi}$ becomes $p(\mu) = \frac{1}{2}$ as $\int_{4\pi} d\mathbf{u} \frac{1}{4\pi} = \int_0^{2\pi} \frac{d\varphi}{2\pi} \int_0^{\pi} d\theta \frac{1}{2} \sin\theta = \int_{-1}^1 d\mu \, p(\mu)$ with $\mu = \cos\theta$ [35]. (c) Examples of molar extinction cross-section $E_\lambda$ and luminescence spectra $p^L(\lambda)$ of the photosensitizer.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0255002.g001
analytical approximate solutions [35], and iii) it is convenient to study the luminescence effect on $\langle A \rangle$, as will be presented in section 3.

The overall effect of luminescence on absorptance is quantified, indicating that luminescence can hardly be neglected in applications concerning solar fuel production. We therefore analyzed the respective influence of two characteristics of luminescence: multiple scattering and inelasticity. For this purpose, the weight of the successive scattering orders was first analyzed, and then the effect of inelasticity in scattering was quantified by comparing the reference results with those obtained thanks to equivalent gray or elastic scattering models. From these developments important conclusions describing the impact of luminescence on MVRPA quantification are given, and rules of thumb for reasoning from well-defined optical thicknesses are proposed for any application using photosensitizers. Various useful analytical approximations for MVRPA assessment, also valid for any photosensitized reactive process in the absence of particles, are fully developed and compared to rigorous calculations in section 4 in order to discuss their validity in detail.

Finally, in section 5 the implementation of the proposed method will make it possible to revisit results already published regarding the use of a bio-inspired catalyst for the production of H₂ in photoreactors, where the fluorescence of Eosin Y was neglected [18]. This will enable us to accurately estimate the overall quantum yields $\varphi$ and thus to refine the performances achievable by this kind of photoreactive system.

# 2 Radiative model

## 2.1 Radiative configuration

Photon transport was modeled in a one dimensional slab of thickness L with incident-normal collimated radiation on the left-hand side (see Fig 1). The incident monochromatic photon flux density inside the medium, at $x = 0$, is $q_{0,\lambda} = q_0 p^i(\lambda)$, where $q_0 = \int_0^{+\infty} d\lambda \, q_{0,\lambda}$ is the total photon flux density and $p^i(\lambda)$ is the incident spectral probability density function, hereafter called incident spectrum: $\int_0^{+\infty} d\lambda \, p^i(\lambda) = 1$. The interfaces at both sides of the slab are assumed to be transparent in order to analyze the roles of luminescence and separate them from the impact of reflection-refraction at the boundaries.

The slab is filled with a homogeneous photocatalytic system solution including luminescent photosensitizer at concentration C (there is no particle leading to elastic scattering in such systems). In the present configuration the photosensitizer is the only molecule that interacts with

<figure><img src="image_4.png" /><figcaption>**Fig 1. Radiative configuration.** (a) Incident-normal collimated radiation at $x = 0$ with spectrum $p^i(\lambda)$ and surface flux density $q_0$. (b) Slab of thickness L containing a luminescent photosensitizer solution at concentration C. The phase space ($\mathbf{r}$, $\mathbf{u}$), with $\mathbf{u}$ the propagation direction, can be reduced to ($x$, $\mu$) in such a one-dimensional slab. In this case, the isotropic distribution $p(\mathbf{u}) = \frac{1}{4\pi}$ becomes $p(\mu) = \frac{1}{2}$ as $\int_{4\pi} d\mathbf{u} \frac{1}{4\pi} = \int_0^{2\pi} \frac{d\varphi}{2\pi} \int_0^{\pi} d\theta \frac{1}{2} \sin\theta = \int_{-1}^1 d\mu \, p(\mu)$ with $\mu = \cos\theta$ [35]. (c) Examples of molar extinction cross-section $E_\lambda$ and luminescence spectra $p^L(\lambda)$ of the photosensitizer.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0255002.g001
radiation in the spectral range of interest (UV-visible). Some catalysts also absorb radiation and we will show in section 5 that the extension to this situation is straightforward.

Luminescence can either be interpreted as an emission source in the volume or as an inelastic scattering phenomenon, as will be discussed in Section 2.2. In this work, the luminescence of photosensitizers is approached from the multiple inelastic-scattering point of view, because of the intuitive physical pictures it provides when describing the absorption of luminescence radiation. There is therefore an important distinction between 1) extinction cross-section $E_\lambda$, which characterizes all interactions with radiation, 2) scattering cross-section $E_{S,\lambda}$, which characterizes interactions leading to luminescence emission only and 3) absorption cross-section $E_{A,\lambda}$, which characterizes interactions leading to other phenomena, including photocatalytic reactions (and other relaxation phenomena), with $E_\lambda = E_{S,\lambda} + E_{A,\lambda}$. The luminescence quantum yield $\Phi_\lambda$ is the proportion of interactions that leads to luminescence emission, and therefore $\Phi_\lambda = E_{S,\lambda}/E_\lambda$ [36–38]. From the photon transport point of view, $\Phi_\lambda$ is the single-scattering albedo: if a photon interacts with a photosensitizer molecule, it is scattered with probability $\Phi_\lambda$ or absorbed with probability $(1 - \Phi_\lambda)$. In the following, the radiative properties of photosensitizers are provided by the molar extinction cross-section $E_\lambda$ and luminescence quantum yield $\Phi_\lambda$, since absorption and scattering cross-sections are easily deduced from these quantities: $E_{A,\lambda} = (1 - \Phi_\lambda)E_\lambda$ and $E_{S,\lambda} = \Phi_\lambda E_\lambda$. In the same way, the extinction coefficient of the solution is $k_\lambda = CE_\lambda$, the absorption coefficient is $(1 - \Phi_\lambda)k_\lambda$ and the scattering coefficient is $\Phi_\lambda k_\lambda$. Scattering events redistribute propagation directions according to isotropic phase function [38] and wavelengths according to the luminescence spectral probability density function $p^L(\lambda|\lambda')$, hereafter called *luminescence spectrum*: $\int_0^{+\infty} d\lambda p^L(\lambda|\lambda') = 1$. The luminescence spectrum is the third radiative property of the photosensitizer, together with $E_\lambda$ and $\Phi_\lambda$, that may vary depending on the photocatalytic system. In the following, Kasha's rule is assumed, and therefore the luminescence quantum yield $\Phi$ and spectrum $p^L$ do not depend on the excitation radiation wavelength: $\Phi_\lambda \equiv \Phi$ and $p^L(\lambda|\lambda') \equiv p^L(\lambda)$ [39]. In other words, the spectral distribution of scattered photons is independent of the spectral distribution of photons incident on the molecule.

For the purpose of radiative analysis, we define the following dimensionless optical thicknesses:

* spectral extinction optical thickness

$$ \tau_\lambda = CLE_\lambda \qquad (3) $$

* mean absorption optical thickness for incident radiation

$$ \tau_A^i = (1 - \Phi)CL\bar{E}^i \qquad (4) $$

with the mean molar extinction cross-section over incident spectrum (gray approximation)

$$ \bar{E}^i = \int_0^{+\infty} d\lambda p^i(\lambda) E_\lambda \qquad (5) $$

* mean scattering optical thickness for luminescence radiation

$$ \tau_S^L = \Phi CL\bar{E}^L \qquad (6) $$

with the mean molar extinction cross-section over luminescence spectrum (gray
radiation in the spectral range of interest (UV-visible). Some catalysts also absorb radiation and we will show in section 5 that the extension to this situation is straightforward.

Luminescence can either be interpreted as an emission source in the volume or as an inelastic scattering phenomenon, as will be discussed in Section 2.2. In this work, the luminescence of photosensitizers is approached from the multiple inelastic-scattering point of view, because of the intuitive physical pictures it provides when describing the absorption of luminescence radiation. There is therefore an important distinction between 1) extinction cross-section $E_\lambda$, which characterizes all interactions with radiation, 2) scattering cross-section $E_{S,\lambda}$, which characterizes interactions leading to luminescence emission only and 3) absorption cross-section $E_{A,\lambda}$, which characterizes interactions leading to other phenomena, including photocatalytic reactions (and other relaxation phenomena), with $E_\lambda = E_{S,\lambda} + E_{A,\lambda}$. The luminescence quantum yield $\Phi_\lambda$ is the proportion of interactions that leads to luminescence emission, and therefore $\Phi_\lambda = E_{S,\lambda}/E_\lambda$ [36–38]. From the photon transport point of view, $\Phi_\lambda$ is the single-scattering albedo: if a photon interacts with a photosensitizer molecule, it is scattered with probability $\Phi_\lambda$ or absorbed with probability $(1 - \Phi_\lambda)$. In the following, the radiative properties of photosensitizers are provided by the molar extinction cross-section $E_\lambda$ and luminescence quantum yield $\Phi_\lambda$, since absorption and scattering cross-sections are easily deduced from these quantities: $E_{A,\lambda} = (1 - \Phi_\lambda)E_\lambda$ and $E_{S,\lambda} = \Phi_\lambda E_\lambda$. In the same way, the extinction coefficient of the solution is $k_\lambda = CE_\lambda$, the absorption coefficient is $(1 - \Phi_\lambda)k_\lambda$ and the scattering coefficient is $\Phi_\lambda k_\lambda$. Scattering events redistribute propagation directions according to isotropic phase function [38] and wavelengths according to the luminescence spectral probability density function $p^L(\lambda|\lambda')$, hereafter called *luminescence spectrum*: $\int_0^{+\infty} d\lambda p^L(\lambda|\lambda') = 1$. The luminescence spectrum is the third radiative property of the photosensitizer, together with $E_\lambda$ and $\Phi_\lambda$, that may vary depending on the photocatalytic system. In the following, Kasha's rule is assumed, and therefore the luminescence quantum yield $\Phi$ and spectrum $p^L$ do not depend on the excitation radiation wavelength: $\Phi_\lambda \equiv \Phi$ and $p^L(\lambda|\lambda') \equiv p^L(\lambda)$ [39]. In other words, the spectral distribution of scattered photons is independent of the spectral distribution of photons incident on the molecule.

For the purpose of radiative analysis, we define the following dimensionless optical thicknesses:

* spectral extinction optical thickness

$$ \tau_{\lambda} = CLE_{\lambda} \qquad (3) $$

* mean absorption optical thickness for incident radiation

$$ \tau_{A}^{i} = (1 - \Phi)CL\bar{E}^{i} \qquad (4) $$

with the mean molar extinction cross-section over incident spectrum (gray approximation)

$$ \bar{E}^{i} = \int_{0}^{+\infty} d\lambda p^{i}(\lambda) E_{\lambda} \qquad (5) $$

* mean scattering optical thickness for luminescence radiation

$$ \tau_{S}^{L} = \Phi CL\bar{E}^{L} \qquad (6) $$

with the mean molar extinction cross-section over luminescence spectrum (gray
approximation)

$$ \bar{E}^L = \int_0^{+\infty} d\lambda p^L(\lambda) E_\lambda \qquad (7) $$

The extinction optical thickness for luminescence radiation $\tau^L = CL\bar{E}^L$ will be also used to formulate approximations in Section 4.

In Section 3 $\tau_A^i$ and $\tau_S^L$ are used to characterize incident and luminescence photon transport, respectively. We choose to focus on the absorption of incident radiation because absorbing and converting incident radiation is the aim of photoreactive systems. On the other hand, we choose to focus on luminescence radiation scattering because the aim of this study is to analyze the impact of multiple inelastic scattering in such systems. However, note that it is straightforward to compute the scattering optical thickness for incident radiation $\tau_S^i = \Phi/(1-\Phi)\tau_A^i$ and the absorption optical thickness for luminescence radiation $\tau_A^L = (1-\Phi)/\Phi\tau_S^L$ from the data provided in Table 2, since the previous definitions of cross-sections lead to $\tau = \tau_S + \tau_A$ and $\Phi = \tau_S/\tau$.

## 2.2 Radiative transfer equation

The Radiative Transfer Equation (RTE) for a one-dimensional slab filled with a luminescent medium is

$$ \mu \frac{dI_{\lambda}(x, \mu)}{dx} = -k_{\lambda} I_{\lambda}(x, \mu) + \frac{1}{2} \int_{0}^{+\infty} d\lambda' p^{L}(\lambda) \int_{-1}^{1} d\mu' \Phi k_{\lambda'} I_{\lambda'}(x, \mu') \qquad (8) $$

where $I_\lambda$ is the intensity, $\mu$ is the dot product between the propagation direction and $e_x$ (see Fig 1). The extinction coefficient $k_\lambda$, the single-scattering albedo $\Phi$ (the luminescent quantum yield) and the luminescence spectrum $p^L(\lambda)$ were defined in the previous paragraphs. The collision term in Eq 8 describes luminescence. When it is interpreted as an emission source in the volume, it can be read as follows:

* $\int_0^{+\infty} d\lambda' \int_{-1}^1 d\mu' k_{\lambda'} I_{\lambda'}(x, \mu')$ is the volumetric rate at which photons with all wavelength $\lambda'$, propagating in all directions $\mu'$, are absorbed at location $x$,

* multiplying the above term by the luminescence quantum yield $\Phi$ leads to the rate of luminescence emission,

* this emission is spectrally distributed according to $p^L(\lambda)$, with isotropic distribution $1/2$ of emission directions $\mu$ (see Fig 1) [35].

When it is interpreted as an inelastic scattering phenomenon, it can be read as follows:

* $\int_{-1}^1 d\mu' \Phi k_{\lambda'} I_{\lambda'}(x, \mu')$ is the volumetric rate at which photons with wavelength $\lambda'$, propagating in all directions $\mu'$, are scattered at location $x$,

* multiplying the above term by the isotropic phase function (see Fig 1) $p(\mu|\mu') = \frac{1}{2}$ leads to the rate at which photons are scattered in direction $\mu$ (per unit $d\mu$),

* $d\lambda'p^L(\lambda)$ is the probability that a photon with wavelength in $d\lambda'$ around $\lambda'$ is scattered with wavelength in $d\lambda$ around $\lambda$,

* Integration over $[0, +\infty[$ sums the contributions of all wavelengths $\lambda'$.

The analysis in the article will be mainly based on this second interpretation.
approximation)

$$ \bar{E}^L = \int_{0}^{+\infty} d\lambda p^L(\lambda) E_{\lambda} \qquad (7) $$

The extinction optical thickness for luminescence radiation $\tau^L = CL\bar{E}^L$ will be also used to formulate approximations in Section 4.

In Section 3 $\tau_A^i$ and $\tau_S^L$ are used to characterize incident and luminescence photon transport, respectively. We choose to focus on the absorption of incident radiation because absorbing and converting incident radiation is the aim of photoreactive systems. On the other hand, we choose to focus on luminescence radiation scattering because the aim of this study is to analyze the impact of multiple inelastic scattering in such systems. However, note that it is straightforward to compute the scattering optical thickness for incident radiation $\tau_S^i = \Phi/(1-\Phi)\tau_A^i$ and the absorption optical thickness for luminescence radiation $\tau_A^L = (1-\Phi)/\Phi\tau_S^L$ from the data provided in Table 2, since the previous definitions of cross-sections lead to $\tau = \tau_S + \tau_A$ and $\Phi = \tau_S/\tau$.

## 2.2 Radiative transfer equation

The Radiative Transfer Equation (RTE) for a one-dimensional slab filled with a luminescent medium is

$$ \mu \frac{dI_{\lambda}(x, \mu)}{dx} = -k_{\lambda} I_{\lambda}(x, \mu) + \frac{1}{2} \int_{0}^{+\infty} d\lambda' p^{L}(\lambda) \int_{-1}^{1} d\mu' \Phi k_{\lambda'} I_{\lambda'}(x, \mu') \qquad (8) $$

where $I_\lambda$ is the intensity, $\mu$ is the dot product between the propagation direction and $e_x$ (see Fig 1). The extinction coefficient $k_\lambda$, the single-scattering albedo $\Phi$ (the luminescent quantum yield) and the luminescence spectrum $p^L(\lambda)$ were defined in the previous paragraphs. The collision term in Eq 8 describes luminescence. When it is interpreted as an emission source in the volume, it can be read as follows:

* $\int_{0}^{+\infty} d\lambda' \int_{-1}^{1} d\mu' k_{\lambda'} I_{\lambda'}(x, \mu')$ is the volumetric rate at which photons with all wavelength $\lambda'$, propagating in all directions $\mu'$, are absorbed at location $x$,

* multiplying the above term by the luminescence quantum yield $\Phi$ leads to the rate of luminescence emission,

* this emission is spectrally distributed according to $p^L(\lambda)$, with isotropic distribution $1/2$ of emission directions $\mu$ (see Fig 1) [35].

When it is interpreted as an inelastic scattering phenomenon, it can be read as follows:

* $\int_{-1}^{1} d\mu' \Phi k_{\lambda'} I_{\lambda'}(x, \mu')$ is the volumetric rate at which photons with wavelength $\lambda'$, propagating in all directions $\mu'$, are scattered at location $x$,

* multiplying the above term by the isotropic phase function (see Fig 1) $p(\mu|\mu') = \frac{1}{2}$ leads to the rate at which photons are scattered in direction $\mu$ (per unit $d\mu$),

* $d\lambda' p^L(\lambda)$ is the probability that a photon with wavelength in $d\lambda'$ around $\lambda'$ is scattered with wavelength in $d\lambda$ around $\lambda$,

* Integration over $[0, +\infty[$ sums the contributions of all wavelengths $\lambda'$.

The analysis in the article will be mainly based on this second interpretation.
For our configuration with normal collimated incident radiation, boundary conditions for $\mu > 0$ are:

$$I_{\lambda}(0, \mu = 1) = q_{0,\lambda} = q_{0} p^{i}(\lambda) \qquad (9)$$

$$I_{\lambda}(0, \mu \neq 1) = 0 \qquad (10)$$

## 2.3 Successive orders of scattering expansion

The analysis in this paper is based on expansion into successive orders of scattering (numerical solutions use the Monte Carlo method, see Section 2.5).

**Intensity** obeys the linear transport equation in Eq 8 and can therefore be expanded into successive orders of scattering [40–46]. The intensity $I$ of the entire photon population is formulated as the sum of the intensities $I^{(j)}$ corresponding to photons that have undergone $j$ scattering events:

$$I_{\lambda}(x, \mu) = \sum_{j=0}^{+\infty} I_{\lambda}^{(j)}(x, \mu) \qquad (11)$$

**Ballistic intensity** $I_{\lambda}^{(0)}(x, \mu)$ is due to photons that come directly from the incident-normal collimated source at the boundary $x = 0$, and is only attenuated within the medium:

$$\mu \frac{dI_{\lambda}^{(0)}(x, \mu)}{dx} = -k_{\lambda} I_{\lambda}^{(0)}(x, \mu) \qquad (12)$$

with the boundary condition

$$I_{\lambda}^{(0)}(0, \mu = 1) = q_{0,\lambda} \qquad (13)$$

$$I_{\lambda}^{(0)}(0, \mu \neq 1) = 0 \qquad (14)$$

**Order $j$ intensity** $I_{\lambda}^{(j)}(x, \mu)$ is due to photons that have undergone $j$ scatterings before reaching phase-space location $(x, \mu)$: $I^{(1)}$ accounts for photons that have undergone one scattering event, $I^{(2)}$ accounts for photons that have undergone two scattering events, and so on. Each intensity $I^{(j>0)}$ obeys a radiative transfer equation of its own, in which the source term corresponds to lower-order photons $I^{(j-1)}$ that are scattered locally and move from population $(j-1)$ to population $(j)$:

$$\mu \frac{dI_{\lambda}^{(j)}(x, \mu)}{dx} = -k_{\lambda} I_{\lambda}^{(j)}(x, \mu) + \frac{1}{2} \int_{0}^{+\infty} d\lambda' p^{I}(\lambda) \int_{-1}^{1} d\mu' \Phi k_{\lambda'} I_{\lambda'}^{(j-1)}(x, \mu') \qquad (15)$$

with the boundary condition

$$I_{\lambda}^{(j)}(0, \mu) = 0 \qquad (16)$$

since no photon having undergone $j > 0$ scattering events is incident on the slab.

**Luminescence intensity** $I_{\lambda}^{(S)}$ sums the contributions $j > 0$ of all scattered photons:

$$I_{\lambda}^{(S)} = \sum_{j=1}^{+\infty} I_{\lambda}^{(j)} \qquad (17)$$
For our configuration with normal collimated incident radiation, boundary conditions for $\mu > 0$ are:

$$I_{\lambda}(0, \mu = 1) = q_{0,\lambda} = q_{0} p^{i}(\lambda) \qquad (9)$$

$$I_{\lambda}(0, \mu \neq 1) = 0 \qquad (10)$$

## 2.3 Successive orders of scattering expansion

The analysis in this paper is based on expansion into successive orders of scattering (numerical solutions use the Monte Carlo method, see Section 2.5).

**Intensity** obeys the linear transport equation in Eq 8 and can therefore be expanded into successive orders of scattering [40–46]. The intensity $I$ of the entire photon population is formulated as the sum of the intensities $I^{(j)}$ corresponding to photons that have undergone $j$ scattering events:

$$I_{\lambda}(x, \mu) = \sum_{j=0}^{+\infty} I_{\lambda}^{(j)}(x, \mu) \qquad (11)$$

**Ballistic intensity** $I_{\lambda}^{(0)}(x, \mu)$ is due to photons that come directly from the incident-normal collimated source at the boundary $x = 0$, and is only attenuated within the medium:

$$\mu \frac{dI_{\lambda}^{(0)}(x, \mu)}{dx} = -k_{\lambda} I_{\lambda}^{(0)}(x, \mu) \qquad (12)$$

with the boundary condition

$$I_{\lambda}^{(0)}(0, \mu = 1) = q_{0,\lambda} \qquad (13)$$

$$I_{\lambda}^{(0)}(0, \mu \neq 1) = 0 \qquad (14)$$

**Order $j$ intensity** $I_{\lambda}^{(j)}(x, \mu)$ is due to photons that have undergone $j$ scatterings before reaching phase-space location $(x, \mu)$: $I^{(1)}$ accounts for photons that have undergone one scattering event, $I^{(2)}$ accounts for photons that have undergone two scattering events, and so on. Each intensity $I^{(j>0)}$ obeys a radiative transfer equation of its own, in which the source term corresponds to lower-order photons $I^{(j-1)}$ that are scattered locally and move from population $(j-1)$ to population $(j)$:

$$\mu \frac{dI_{\lambda}^{(j)}(x, \mu)}{dx} = -k_{\lambda} I_{\lambda}^{(j)}(x, \mu) + \frac{1}{2} \int_{0}^{+\infty} d\lambda' p^{I}(\lambda) \int_{-1}^{1} d\mu' \Phi k_{\lambda'} I_{\lambda'}^{(j-1)}(x, \mu') \qquad (15)$$

with the boundary condition

$$I_{\lambda}^{(j)}(0, \mu) = 0 \qquad (16)$$

since no photon having undergone $j > 0$ scattering events is incident on the slab.

**Luminescence intensity** $I_{\lambda}^{(S)}$ sums the contributions $j > 0$ of all scattered photons:

$$I_{\lambda}^{(S)} = \sum_{j=1}^{+\infty} I_{\lambda}^{(j)} \qquad (17)$$
leading to (see Eq 11)

$$I_{\lambda}(x, \mu) = I_{\lambda}^{(0)}(x, \mu) + I_{\lambda}^{(S)}(x, \mu) \qquad (18)$$

$I_{\lambda}^{(S)}$ obeys the following RTE, that is the sum of the RTEs in Eq 15 for all $j > 0$:

$$\mu \frac{dI_{\lambda}^{(S)}(x, \mu)}{dx} = -k_{\lambda} I_{\lambda}^{(S)}(x, \mu) + \frac{1}{2} \int_{0}^{+\infty} d\lambda' p^{L}(\lambda) \int_{-1}^{1} d\mu' \Phi k_{\lambda'} I_{\lambda'}^{(S)}(x, \mu') + \frac{1}{2} \int_{0}^{+\infty} d\lambda' p^{L}(\lambda) \int_{-1}^{1} d\mu' \Phi k_{\lambda'} I_{\lambda'}^{(0)}(x, \mu') \qquad (19)$$

with the boundary condition

$$I_{\lambda}^{(S)}(0, \mu) = 0 \qquad (20)$$

This is the same RTE as that for intensity $I_{\lambda}$ except that the source is not at the boundary because luminescence emission stimulated by ballistic photons $I^{(0)}$ within the volume is the source of scattered photons, as formulated in the last term of Eq 19.

A meaningful property is due to the fact that the expansion has a zero-order closure: $I_{\lambda}^{(0)}$ has a closed-form expression, independent of the $I_{\lambda}^{(j>0)}$, and each higher order $j > 0$ is a function of $I_{\lambda}^{(j-1)}$ only. Therefore, truncating the expansion at order $q$ means that the contribution of photons that have undergone more than $q$ scattering events is neglected. As a result, the $q$-th order expansion in successive orders of scattering systematically underestimates intensity.

## 2.4 Radiative quantities of interest

The Mean Volumetric Rate of Photons Absorbed $\langle A \rangle$ (MVRPA) is the key radiative quantity in the study of photoreactive systems [16]. It is obtained by integrating the intensity $I$ over directions, locations and wavelengths:

$$\langle A \rangle^{(\bullet)} = \frac{1}{L} \int_{0}^{+\infty} d\lambda \int_{0}^{L} dx \int_{-1}^{1} d\mu (1 - \Phi) k_{\lambda} I_{\lambda}^{(\bullet)}(x, \mu) \qquad (21)$$

where superscript $(\bullet)$ is dropped in the formulation of the total rate of absorption $\langle A \rangle$; otherwise, take $\bullet \equiv j$ for the rate $\langle A \rangle^{(j)}$ at which photons that have undergone $j$ scattering events are absorbed, and $\bullet \equiv S$ for the rate $\langle A \rangle^{(S)} = \sum_{j=1}^{+\infty} \langle A \rangle^{(j)}$ at which scattered (i.e. luminescence) photons are absorbed:

$$\langle A \rangle = \sum_{j=0}^{+\infty} \langle A \rangle^{(j)} = \langle A \rangle^{(0)} + \langle A \rangle^{(S)} \qquad (22)$$

Absorptance $P_{A}$ is the key quantity in our radiative analysis. It can be interpreted as the proportion of photons that are absorbed in the medium, leading in our 1D cartesian geometry to:

$$\langle A \rangle = \frac{q_0}{L} P_{A} \qquad (23)$$

The maximum photon absorption rate $\langle A \rangle_{max}$ is then obtained with $P_{A,max}$, which is lower than 1 if any scattering phenomenon exists in the medium. $P_{A,max} = 1$ would be only achieved
leading to (see Eq 11)

$$I_{\lambda}(x, \mu) = I_{\lambda}^{(0)}(x, \mu) + I_{\lambda}^{(S)}(x, \mu) \quad (18)$$

$I_{\lambda}^{(S)}$ obeys the following RTE, that is the sum of the RTEs in Eq 15 for all $j > 0$:

$$\mu \frac{dI_{\lambda}^{(S)}(x, \mu)}{dx} = -k_{\lambda} I_{\lambda}^{(S)}(x, \mu) + \frac{1}{2} \int_{0}^{+\infty} d\lambda' p^{L}(\lambda) \int_{-1}^{1} d\mu' \Phi k_{\lambda'} I_{\lambda'}^{(S)}(x, \mu') + \frac{1}{2} \int_{0}^{+\infty} d\lambda' p^{L}(\lambda) \int_{-1}^{1} d\mu' \Phi k_{\lambda'} I_{\lambda'}^{(0)}(x, \mu') \quad (19)$$

with the boundary condition

$$I_{\lambda}^{(S)}(0, \mu) = 0 \quad (20)$$

This is the same RTE as that for intensity $I_{\lambda}$ except that the source is not at the boundary because luminescence emission stimulated by ballistic photons $I^{(0)}$ within the volume is the source of scattered photons, as formulated in the last term of Eq 19.

**A meaningful property** is due to the fact that the expansion has a zero-order closure: $I_{\lambda}^{(0)}$ has a closed-form expression, independent of the $I_{\lambda}^{(j>0)}$, and each higher order $j > 0$ is a function of $I_{\lambda}^{(j-1)}$ only. Therefore, truncating the expansion at order $q$ means that the contribution of photons that have undergone more than $q$ scattering events is neglected. As a result, the $q$-th order expansion in successive orders of scattering systematically underestimates intensity.

## 2.4 Radiative quantities of interest

**The Mean Volumetric Rate of Photons Absorbed** $\langle A \rangle$ (MVRPA) is the key radiative quantity in the study of photoreactive systems [16]. It is obtained by integrating the intensity $I$ over directions, locations and wavelengths:

$$\langle A \rangle^{(\bullet)} = \frac{1}{L} \int_{0}^{+\infty} d\lambda \int_{0}^{L} dx \int_{-1}^{1} d\mu (1 - \Phi) k_{\lambda} I_{\lambda}^{(\bullet)}(x, \mu) \quad (21)$$

where superscript $(\bullet)$ is dropped in the formulation of the total rate of absorption $\langle A \rangle$; otherwise, take $\bullet \equiv j$ for the rate $\langle A \rangle^{(j)}$ at which photons that have undergone $j$ scattering events are absorbed, and $\bullet \equiv S$ for the rate $\langle A \rangle^{(S)} = \sum_{j=1}^{+\infty} \langle A \rangle^{(j)}$ at which scattered (i.e. luminescence) photons are absorbed:

$$\langle A \rangle = \sum_{j=0}^{+\infty} \langle A \rangle^{(j)} = \langle A \rangle^{(0)} + \langle A \rangle^{(S)} \quad (22)$$

**Absorptance** $P_{A}$ is the key quantity in our radiative analysis. It can be interpreted as the proportion of photons that are absorbed in the medium, leading in our 1D cartesian geometry to:

$$\langle A \rangle = \frac{q_0}{L} P_A \quad (23)$$

The maximum photon absorption rate $\langle A \rangle_{max}$ is then obtained with $P_{A,max}$, which is lower than 1 if any scattering phenomenon exists in the medium. $P_{A,max} = 1$ would be only achieved
if ballistic and scattered radiation could be completely absorbed. We chose to base the radiative analysis in Section 3 mainly on absorptance rather than MVRPA, because it enables us to clearly distinguish the impact of:

* photon transport, which is contained in $P_A$ only and depends on optical thicknesses,

* dimension $L$ of the system, which depends on the experimental setup or the engineering specifications when designing units for industrial-scale production (note that for a given dimension $L$, optical thickness can be controlled by adjusting the photosensitizer concentration $C$),

* incident photon flux $q_0$, which depends on the source and fluctuates in solar conditions.

Substituting Eq 21 in Eq 23 leads to:

$$P_A^{(\bullet)} = \frac{1}{q_0} \int_0^{+\infty} d\lambda \int_0^L dx \int_{-1}^1 d\mu (1-\Phi) k_\lambda I_\lambda^{(\bullet)}(x, \mu) \quad (24)$$

Overall, expansion into successive orders of scattering is:

$$P_A = \sum_{j=0}^{+\infty} P_A^{(j)} = P_A^{(0)} + P_A^{(S)} \quad (25)$$

The distribution $P^{(j)}$ of scattering orders within the expansion indicates the respective weight of each scattering order:

$$P^{(j)} = \frac{P_A^{(j)}}{P_A} \quad (26)$$

The cumulative $P^{(j \le q)} = \sum_{j=0}^q P^{(j)}$ is the proportion of absorption that is described by the $q$-th order expansion in successive orders of scattering. Alternatively, the complementary (tail) cumulative $P^{(j > q)} = \sum_{j=q+1}^{+\infty} P^{(j)}$ is the relative error on both absorptance $P_A$ and MVRPA $\langle A \rangle$ when the expansion is truncated at the $q$-th order. Of course, $\sum_{j=0}^{+\infty} P^{(j)} = 1$ and $P^{(j > q)} = 1 - P^{(j \le q)}$ according to Eq 25. In addition, the contribution of luminescence radiation is $P^{(S)} = P^{(j > 0)} = P_A^{(S)}/P_A$, and we have $P^{(0)} + P^{(S)} = 1$. Note that $P^{(j)}$ can also be interpreted as the proportion of photons absorbed after $j$ scattering events with respect to all the absorbed photons, or as the probability that a photon is absorbed after $j$ scattering events.

## 2.5 Monte Carlo algorithm

An analytical solution for the above multiple inelastic-scattering radiative model is not available. Our study will therefore use a Monte Carlo numerical solution, which has three main advantages: 1) the algorithms are intuitive and easily modified to evaluate physical hypotheses, as will be presented in Section 3.3; 2) it provides a reference solution including statistical error estimation and 3) even if the present work is conducted in a one-dimensional cartesian configuration, it is straightforward to extend implementation to complex geometry in subsequent studies [47].

The Monte Carlo algorithm presented below was used to estimate the absorptance $P_A$ and each scattering order $P_A^{(j)}$ for a given concentration $C$, slab thickness $L$ and radiative properties $\Phi$ and $E_\lambda$. It consists in sampling $N$ independent realizations $w_i$, $i = 1, \dots, N$ with the following sampling procedure:
if ballistic and scattered radiation could be completely absorbed. We chose to base the radiative analysis in Section 3 mainly on absorptance rather than MVRPA, because it enables us to clearly distinguish the impact of:

* photon transport, which is contained in $P_A$ only and depends on optical thicknesses,

* dimension $L$ of the system, which depends on the experimental setup or the engineering specifications when designing units for industrial-scale production (note that for a given dimension $L$, optical thickness can be controlled by adjusting the photosensitizer concentration $C$),

* incident photon flux $q_0$, which depends on the source and fluctuates in solar conditions.

Substituting Eq 21 in Eq 23 leads to:

$$P_A^{(\bullet)} = \frac{1}{q_0} \int_0^{+\infty} d\lambda \int_0^L dx \int_{-1}^1 d\mu (1-\Phi) k_\lambda I_\lambda^{(\bullet)}(x, \mu) \quad (24)$$

Overall, expansion into successive orders of scattering is:

$$P_A = \sum_{j=0}^{+\infty} P_A^{(j)} = P_A^{(0)} + P_A^{(S)} \quad (25)$$

The distribution $P^{(j)}$ of scattering orders within the expansion indicates the respective weight of each scattering order:

$$P^{(j)} = \frac{P_A^{(j)}}{P_A} \quad (26)$$

The cumulative $P^{(j \le q)} = \sum_{j=0}^q P^{(j)}$ is the proportion of absorption that is described by the $q$-th order expansion in successive orders of scattering. Alternatively, the complementary (tail) cumulative $P^{(j > q)} = \sum_{j=q+1}^{+\infty} P^{(j)}$ is the relative error on both absorptance $P_A$ and MVRPA $\langle A \rangle$ when the expansion is truncated at the $q$-th order. Of course, $\sum_{j=0}^{+\infty} P^{(j)} = 1$ and $P^{(j > q)} = 1 - P^{(j \le q)}$ according to Eq 25. In addition, the contribution of luminescence radiation is $P^{(S)} = P^{(j > 0)} = P_A^{(S)}/P_A$, and we have $P^{(0)} + P^{(S)} = 1$. Note that $P^{(j)}$ can also be interpreted as the proportion of photons absorbed after $j$ scattering events with respect to all the absorbed photons, or as the probability that a photon is absorbed after $j$ scattering events.

## 2.5 Monte Carlo algorithm

An analytical solution for the above multiple inelastic-scattering radiative model is not available. Our study will therefore use a Monte Carlo numerical solution, which has three main advantages: 1) the algorithms are intuitive and easily modified to evaluate physical hypotheses, as will be presented in Section 3.3; 2) it provides a reference solution including statistical error estimation and 3) even if the present work is conducted in a one-dimensional cartesian configuration, it is straightforward to extend implementation to complex geometry in subsequent studies [47].

The Monte Carlo algorithm presented below was used to estimate the absorptance $P_A$ and each scattering order $P_A^{(j)}$ for a given concentration $C$, slab thickness $L$ and radiative properties $\Phi$ and $E_\lambda$. It consists in sampling $N$ independent realizations $w_i$, $i = 1, \dots, N$ with the following sampling procedure:
* **Step 1:** Initialization of Monte Carlo weights $w_i = 0$ and $w_i^{(q)} = 0$ for $q \in \mathbb{N}$, scattering counter $j = 0$, emission location $\mathbf{r} = \mathbf{0}$ and emission direction $\mathbf{u} = \mathbf{e}_x$ (see Fig 1).

* **Step 2:** Wavelength $\lambda$ is sampled according to the incident spectrum $p^i(\lambda)$, and the extinction coefficient $k_\lambda = CE_\lambda$ is interpolated from spectral database.

* **Step 3:** A first extinction length $l$ is sampled over $[0, +\infty[$ according to the Bouguer extinction probability density function $p_L(l) = k_\lambda e^{-k_\lambda l}$ and the location is updated: $\mathbf{r} = \mathbf{r} + l \mathbf{u}$.

* **while** $0 \le \mathbf{r} \cdot \mathbf{e}_x \le L$ (location is inside the slab) **do**

**Step 4:** Bernoulli test: uniform sampling of a realization $\xi$ over $[0, 1]$

* **if** $\xi < \Phi$ (scattering event) **then**

**Step 4a.1:** Scattering counter is incremented: $j = j + 1$.

**Step 4a.2:** Scattering direction $\mathbf{u}$ is sampled according to the isotropic phase function.

**Step 4a.3:** Wavelength $\lambda$ is sampled according to the luminescence spectrum $p^L(\lambda)$ and the extinction coefficient $k_\lambda = CE_\lambda$ is interpolated from the spectral database.

**Step 4a.4:** Extinction length $l$ is sampled over $[0, +\infty[$ according to the Bouguer extinction probability density function $p_L(l) = k_\lambda e^{-k_\lambda l}$ and the location is updated: $\mathbf{r} = \mathbf{r} + l \mathbf{u}$.

**end**

* **Else**

**Step 4b.1:** Path sampling is terminated due to absorption:

* weights are computed:

$$w_i^{(j)} = 1$$

$$w_i = 1$$

* go to the end of Monte Carlo realization.

**end**

**end**

Absorptance is estimated as:

$$P_A \approx \frac{1}{N} \sum_{i=1}^{N} w_i \qquad (27)$$

with standard error:

$$\sigma(P_A) = \frac{1}{\sqrt{N-1}} \sqrt{\frac{1}{N} \sum_{i=1}^{N} w_i^2 - \left(\frac{1}{N} \sum_{i=1}^{N} w_i\right)^2} \qquad (28)$$

In Tables 2 and 3, the number $N$ of Monte Carlo realizations has been adjusted to obtain a relative uncertainty $\Delta P_A$ which is always smaller than 1%, with $\Delta P_A$ defined as the 95% confidence interval: $\Delta P_A = t_{0.05, N-1} \sigma(P_A)$, where $t_{0.05, N-1}$ is the Student quantile ($t_{0.05, N-1} \approx 2$ for a typical number $N$ of Monte Carlo realizations). Scattering orders $P_A^{(j)}$ are obtained in the same way, but replacing $w_i$ by $w_i^{(j)}$ in Eqs 27 and 28.

## 3 A potpourri of photosensitizers

Radiative analysis is performed by computing absorptance for incident solar radiation and four commonly-used photosensitizers for solar fuel production: Ru[Bpy]$_3^{2+}$, TATA$^+$, Eosin Y and Rhodamine B. Radiative properties (molar extinction cross-section $E_\lambda$, luminescence emission spectra $p^L(\lambda)$ and quantum yield $\Phi$) obtained from the literature are presented in
* **Step 1:** Initialization of Monte Carlo weights $w_i = 0$ and $w_i^{(q)} = 0$ for $q \in \mathbb{N}$, scattering counter $j = 0$, emission location $\mathbf{r} = \mathbf{0}$ and emission direction $\mathbf{u} = \mathbf{e}_x$ (see Fig 1).

* **Step 2:** Wavelength $\lambda$ is sampled according to the incident spectrum $p^i(\lambda)$, and the extinction coefficient $k_\lambda = CE_\lambda$ is interpolated from spectral database.

* **Step 3:** A first extinction length $l$ is sampled over $[0, +\infty[$ according to the Bouguer extinction probability density function $p_L(l) = k_\lambda e^{-k_\lambda l}$ and the location is updated: $\mathbf{r} = \mathbf{r} + l \mathbf{u}$.

* **while** $0 \le \mathbf{r} \cdot \mathbf{e}_x \le L$ (location is inside the slab) **do**

**Step 4:** Bernoulli test: uniform sampling of a realization $\xi$ over $[0, 1]$

* **if** $\xi < \Phi$ (scattering event) **then**

**Step 4a.1:** Scattering counter is incremented: $j = j + 1$.

**Step 4a.2:** Scattering direction $\mathbf{u}$ is sampled according to the isotropic phase function.

**Step 4a.3:** Wavelength $\lambda$ is sampled according to the luminescence spectrum $p^L(\lambda)$ and the extinction coefficient $k_\lambda = CE_\lambda$ is interpolated from the spectral database.

**Step 4a.4:** Extinction length $l$ is sampled over $[0, +\infty[$ according to the Bouguer extinction probability density function $p_L(l) = k_\lambda e^{-k_\lambda l}$ and the location is updated: $\mathbf{r} = \mathbf{r} + l \mathbf{u}$.

**end**

* **Else**

**Step 4b.1:** Path sampling is terminated due to absorption:

* weights are computed:

$$w_i^{(j)} = 1$$

$$w_i = 1$$

* go to the end of Monte Carlo realization.

**end**

**end**

Absorptance is estimated as:

$$P_A \approx \frac{1}{N} \sum_{i=1}^{N} w_i \qquad (27)$$

with standard error:

$$\sigma(P_A) = \frac{1}{\sqrt{N-1}} \sqrt{\frac{1}{N} \sum_{i=1}^{N} w_i^2 - \left(\frac{1}{N} \sum_{i=1}^{N} w_i\right)^2} \qquad (28)$$

In Tables 2 and 3, the number $N$ of Monte Carlo realizations has been adjusted to obtain a relative uncertainty $\Delta P_A / P_A$ which is always smaller than 1%, with $\Delta P_A$ defined as the 95% confidence interval: $\Delta P_A = t_{0.05, N-1} \sigma(P_A)$, where $t_{0.05, N-1}$ is the Student quantile ($t_{0.05, N-1} \approx 2$ for a typical number $N$ of Monte Carlo realizations). Scattering orders $P_A^{(j)}$ are obtained in the same way, but replacing $w_i$ by $w_i^{(j)}$ in Eqs 27 and 28.

### 3 A potpourri of photosensitizers

Radiative analysis is performed by computing absorptance for incident solar radiation and four commonly-used photosensitizers for solar fuel production: Ru[Bpy]$_3^{2+}$, TATA$^+$, Eosin Y and Rhodamine B. Radiative properties (molar extinction cross-section $E_\lambda$, luminescence emission spectra $p^L(\lambda)$ and quantum yield $\Phi$) obtained from the literature are presented in
<figure><img src="image_6.png" /><figcaption>**Fig 2. Photosensitizer radiative properties.** In all plots, plain lines stand for UV-visible molar extinction cross-section spectrum $E_\lambda$ and dashed lines for luminescence emission spectrum $p^l(\lambda)$. Luminescent quantum yield $\Phi$ is provided in a box. (a) Ru[Bpy]₃Cl₂ (40 µM) in aqueous buffer (pH 6.8) solution [48]. (b) TATA(Cl) (0.5 mM) UV-visible molar extinction cross-section, TATA(Cl) (10 µM) fluorescence emission spectrum and quantum yield under argon in acetate buffer (1 M) at pH 4.5 [25] (c) Eosin Y (24 µM) in a mixture of triethylamine (10% vol.) in water at pH 10.5 [18, 21] (d) Rhodamine B in ethanol [26, 49]. Extinction cross-sections are obtained assuming that absorbance $A_\lambda$ measured by the authors cited above using spectrophotometry experiments can be interpreted as $A_\lambda = 1 - \exp(-E_\lambda CL)$. Overlapping between molar extinction cross-section and luminescence emission spectrum is represented for each photosensitizer (except Ru[Bpy]₃Cl₂ for which no overlapping is observed) by 2 parallel lines linked by a double arrow.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0255002.g002

**Fig 2.** Since these properties depend on the solvent, working conditions were selected to cover a wide range of situations. Incident solar AM1.5 spectrum $p^l(\lambda)$ is used, in the spectral range [280, 650 nm] where the studied photosensitizers are active (see Fig 3). It represents 24.4% of the complete AM1.5 spectrum (in moles of photon distribution, see Fig 3). TATA⁺ is studied in the spectral range [335, 650 nm] only, because molar extinction cross-section is not available in the literature for [280, 335 nm], to the best of our knowledge. In this particular situation it represents a proportion of 24.1% of the total AM1.5 spectrum. Corresponding mean molar extinction cross-sections are provided in Table 1; they are required to evaluate absorption and scattering optical thicknesses with Eqs 4–7. Four incident absorption optical thicknesses were selected: $\tau_A^i = 0.1; 0.3; 1.4; 4.6$. $\tau_A^i = 4.6$ is a photochemistry standard corresponding to 99%
<figure><img src="image_5.png" /><figcaption>**Fig 2. Photosensitizer radiative properties.** In all plots, plain lines stand for UV-visible molar extinction cross-section spectrum $E_\lambda$ and dashed lines for luminescence emission spectrum $p^l(\lambda)$. Luminescent quantum yield $\Phi$ is provided in a box. (a) Ru[Bpy]₃Cl₂ (40 µM) in aqueous buffer (pH 6.8) solution [48]. (b) TATA(Cl) (0.5 mM) UV-visible molar extinction cross-section, TATA(Cl) (10 µM) fluorescence emission spectrum and quantum yield under argon in acetate buffer (1 M) at pH 4.5 [25] (c) Eosin Y (24 µM) in a mixture of triethylamine (10% vol.) in water at pH 10.5 [18, 21] (d) Rhodamine B in ethanol [26, 49]. Extinction cross-sections are obtained assuming that absorbance $A_\lambda$ measured by the authors cited above using spectrophotometry experiments can be interpreted as $A_\lambda = 1 - \exp(-E_\lambda CL)$. Overlapping between molar extinction cross-section and luminescence emission spectrum is represented for each photosensitizer (except Ru[Bpy]₃Cl₂ for which no overlapping is observed) by 2 parallel lines linked by a double arrow.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0255002.g002

**Fig 2.** Since these properties depend on the solvent, working conditions were selected to cover a wide range of situations. Incident solar AM1.5 spectrum $p^l(\lambda)$ is used, in the spectral range [280, 650 nm] where the studied photosensitizers are active (see Fig 3). It represents 24.4% of the complete AM1.5 spectrum (in moles of photon distribution, see Fig 3). TATA⁺ is studied in the spectral range [335, 650 nm] only, because molar extinction cross-section is not available in the literature for [280, 335 nm], to the best of our knowledge. In this particular situation it represents a proportion of 24.1% of the total AM1.5 spectrum. Corresponding mean molar extinction cross-sections are provided in Table 1; they are required to evaluate absorption and scattering optical thicknesses with Eqs 4–7. Four incident absorption optical thicknesses were selected: $\tau_A^i = 0.1; 0.3; 1.4; 4.6$. $\tau_A^i = 4.6$ is a photochemistry standard corresponding to 99%
<figure><img src="image_7.png" /><figcaption>**Fig 3.** Solar AM1.5 incident spectrum in moles of photons. Probability density function $p^i(\lambda)$ between $\lambda_{min}^i = 280$ nm and $\lambda_{max}^i = 65$ nm corresponding to an incident photon flux $q_0 = 1743 \ \mu\text{mol}_{\text{hv}}\text{m}^{-2}\text{s}^{-1}$ (for TATA$^+$, we work between $\lambda_{min}^i = 335$ nm and $\lambda_{max}^i = 650$ nm and therefore $q_0 = 1721 \ \mu\text{mol}_{\text{hv}}\text{m}^{-2}\text{s}^{-1}$). Note that $p^i(\lambda) = 0$ for $\lambda \notin [\lambda_{min}^i, \lambda_{max}^i]$ and spectral integration over incident radiation therefore gives $\int_0^{+\infty} d\lambda \, p^i(\lambda) = \int_{\lambda_{min}^i}^{\lambda_{max}^i} d\lambda \, p^i(\lambda)$.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0255002.g003

of incident photons absorbed in gray and purely absorbing media [1] and other values are chosen to investigate thinner optical thicknesses. In practice, the product CL leading to targeted $\tau_A^i$ values is determined with Eq 4 where $\bar{E}^i$, $\bar{E}^L$ and $\Phi$ are known for each photosensitizer. Then the same CL value is used to compute the scattering optical thickness $\tau_S^L$ with Eq 6. Results gathered in Tables 2 and 3 include:

<table><caption>**Table 1.** The table summarizes mean molar extinction cross-sections for each photosensitizer. $\bar{E}^i$ is weighted by the incident AM1.5 source in the spectral range [280 nm, 650 nm] for Ru[Bpy]$_3^{2+}$, Eosin Y and Rhodamine B, and in the spectral range [335 nm, 650 nm] for TATA$^+$ due to extinction coefficient data availability (see Eq 5 and Fig 3). $\bar{E}^L$ is weighted by the luminescence emission spectra (see Eq 7). Note that $p^i(\lambda) = 0$ for $\lambda \notin [\lambda_{min}^L, \lambda_{max}^L]$, with $\lambda_{min}^L = 500$ nm and $\lambda_{max}^L = 800$ nm for Ru[Bpy]$_3^{2+}$, $\lambda_{min}^L = 550$ nm and $\lambda_{max}^L = 775$ nm for TATA$^+$, $\lambda_{min}^L = 475$ nm and $\lambda_{max}^L = 675$ nm for Eosin Y and $\lambda_{min}^L = 500$ nm and $\lambda_{max}^L = 800$ nm for Rhodamine B. Thus spectral integration over luminescence spectrum gives $\int_0^{+\infty} d\lambda \, p^L(\lambda) = \int_{\lambda_{min}^L}^{\lambda_{max}^L} d\lambda \, p^L(\lambda)$. $\bar{E}^L = 0$ for Ru[Bpy]$_3^{2+}$ because no overlapping is observed between extinction cross section and luminescence spectra (see Fig 2).</caption><thead><tr><th>Photosensitizers</th><th>Ru[Bpy]<sub>3</sub><sup>2+</sup></th><th>TATA<sup>+</sup></th><th>Eosin Y</th><th>Rhodamine B</th></tr></thead><tbody><tr><td>E&#773;<sup>i</sup> (m<sup>2</sup>.mol<sup>-1</sup>)</td><td>564</td><td>7957</td><td>3454</td><td>18337</td></tr><tr><td>E&#773;<sup>L</sup> (m<sup>2</sup>.mol<sup>-1</sup>)</td><td>0</td><td>2833</td><td>2728</td><td>24943</td></tr></tbody></table>

https://doi.org/10.1371/journal.pone.0255002.t001
<figure><img src="image_6.png" /><figcaption>**Fig 3.** Solar AM1.5 incident spectrum in moles of photons. Probability density function $p^i(\lambda)$ between $\lambda_{min}^i = 280$ nm and $\lambda_{max}^i = 65$ nm corresponding to an incident photon flux $q_0 = 1743 \mu mol_{hv} \cdot m^{-2} \cdot s^{-1}$ (for TATA+, we work between $\lambda_{min}^i = 335$ nm and $\lambda_{max}^i = 650$ nm and therefore $q_0 = 1721 \mu mol_{hv} \cdot m^{-2} \cdot s^{-1}$). Note that $p^i(\lambda) = 0$ for $\lambda \notin [\lambda_{min}^i, \lambda_{max}^i]$ and spectral integration over incident radiation therefore gives
$$ \int_{0}^{+\infty} d\lambda \, p^i(\lambda) = \int_{\lambda_{min}^i}^{\lambda_{max}^i} d\lambda \, p^i(\lambda). $$</figcaption></figure>

https://doi.org/10.1371/journal.pone.0255002.g003

of incident photons absorbed in gray and purely absorbing media [1] and other values are chosen to investigate thinner optical thicknesses. In practice, the product CL leading to targeted $\tau_A^i$ values is determined with Eq 4 where $\bar{E}^i$, $\bar{E}^L$ and $\Phi$ are known for each photosensitizer. Then the same CL value is used to compute the scattering optical thickness $\tau_S^L$ with Eq 6. Results gathered in Tables 2 and 3 include:

<table><caption>**Table 1.** The table summarizes mean molar extinction cross-sections for each photosensitizer. $\bar{E}^i$ is weighted by the incident AM1.5 source in the spectral range [280 nm, 650 nm] for Ru[Bpy]$_3^{2+}$, Eosin Y and Rhodamine B, and in the spectral range [335 nm, 650 nm] for TATA+ due to extinction coefficient data availability (see Eq 5 and Fig 3). $\bar{E}^L$ is weighted by the luminescence emission spectra (see Eq 7). Note that $p^i(\lambda) = 0$ for $\lambda \notin [\lambda_{min}^L, \lambda_{max}^L]$, with $\lambda_{min}^L = 500$ nm and $\lambda_{max}^L = 800$ nm for Ru[Bpy]$_3^{2+}$, $\lambda_{min}^L = 550$ nm and $\lambda_{max}^L = 775$ nm for TATA+, $\lambda_{min}^L = 475$ nm and $\lambda_{max}^L = 675$ nm for Eosin Y and $\lambda_{min}^L = 500$ nm and $\lambda_{max}^L = 800$ nm for Rhodamine B. Thus spectral integration over luminescence spectrum gives $\int_{0}^{+\infty} d\lambda \, p^L(\lambda) = \int_{\lambda_{min}^L}^{\lambda_{max}^L} d\lambda \, p^L(\lambda)$. $\bar{E}^L = 0$ for Ru[Bpy]$_3^{2+}$ because no overlapping is observed between extinction cross section and luminescence spectra (see Fig 2).</caption><thead><tr><th>Photosensitizers</th><th>Ru[Bpy]<sub>3</sub><sup>2+</sup></th><th>TATA<sup>+</sup></th><th>Eosin Y</th><th>Rhodamine B</th></tr></thead><tbody><tr><td>E&#773;<sup>i</sup> (m<sup>2</sup>.mol<sup>-1</sup>)</td><td>564</td><td>7957</td><td>3454</td><td>18337</td></tr><tr><td>E&#773;<sup>L</sup> (m<sup>2</sup>.mol<sup>-1</sup>)</td><td>0</td><td>2833</td><td>2728</td><td>24943</td></tr></tbody></table>

https://doi.org/10.1371/journal.pone.0255002.t001
<table><caption>Table 2. Results for Ru[Bpy]$_{3}^{2+}$, TATA$^{+}$, Eosin Y and Rhodamine B photosensitizers (radiative properties are in Fig 2) for the incident solar spectrum in Fig 3, at selected optical thicknesses $\tau_{A}^{l}$ and $\tau_{S}^{l}$ (see Eqs 4–7 and Table 1): Reference absorptance $P_{A}$, reference MVRPA $\langle A \rangle$ (slab thickness $L = 1$ cm; incident flux density $q_{0}$ = 1743 $\mu mol_{h,v} \cdot m^{-2} \cdot s^{-1}$ and 1721 $\mu mol_{h,v} \cdot m^{-2} \cdot s^{-1}$ for TATA$^{+}$), relative difference $\Delta(X)$ between these values and those obtained with a model denoted X (see Eq 29). $X \equiv (\Phi = 0)$ when luminescence is neglected, $X \equiv (j = 0)$ for zero-order scattering expansion ($P_{A}^{X} = P_{A}^{(0)}$) as presented in Section 4.1, $X \equiv (j \le 1)$ for first-order scattering expansion ($P_{A}^{X} = P_{A}^{(0)} + P_{A}^{(1)}$) as presented in Section 4.3, $X \equiv (elastic)$ and $X \equiv (gray)$ for elastic and gray models presented in Section 3.3, $X \equiv (\overline{P1})$ and $X \equiv (\overline{P1})$ for elastic and gray P1 approximations presented in Section 4.2 and finally, $X \equiv (\overline{SS})$ for the gray single-scattering approximation presented in Section 4.3. Results for Ru[Bpy]$_{3}^{2+}$ do not include all X models, since zero-order expansion already provides a reference solution. The relative uncertainty of the Monte Carlo absorptance estimation is always lower than 1%; the precision indicated accounts for case-by-case absolute uncertainty.</caption><tbody><tr><th colspan="5">Photosensitizer</th></tr><tr><td colspan="5">Ru[Bpy]$_{3}^{2+}$ (Φ = 0.062)</td></tr><tr><td>τ<sub>A</sub><sup>l</sup> | τ<sub>S</sub><sup>l</sup></td><td>0.1 | 0</td><td>0.3 | 0</td><td>1.4 | 0</td><td>4.6 | 0</td></tr><tr><td>P<sub>A</sub></td><td>0.0831</td><td>0.1805</td><td>0.3033</td><td>0.3362</td></tr><tr><td>⟨A⟩ [mol<sub>h,v</sub>·m<sup>-3</sup>·s<sup>-1</sup>]</td><td>1.45 × 10<sup>-2</sup></td><td>3.15 × 10<sup>-2</sup></td><td>5.28 × 10<sup>-2</sup></td><td>5.86 × 10<sup>-2</sup></td></tr><tr><td>Δ(Φ = 0)</td><td>6.6%</td><td>6.6%</td><td>6.6%</td><td>6.6%</td></tr><tr><td>Δ(j = 0)</td><td>0%</td><td>0%</td><td>0%</td><td>0%</td></tr><tr><th colspan="5">Photosensitizer</th></tr><tr><td colspan="5">TATA<sup>+</sup> (Φ = 0.15)</td></tr><tr><td>τ<sub>A</sub><sup>l</sup> | τ<sub>S</sub><sup>l</sup></td><td>0.1 | 0.006</td><td>0.3 | 0.019</td><td>1.4 | 0.088</td><td>4.6 | 0.289</td></tr><tr><td>P<sub>A</sub></td><td>0.0909</td><td>0.2265</td><td>0.5164</td><td>0.6832</td></tr><tr><td>⟨A⟩ [mol<sub>h,v</sub>·m<sup>-3</sup>·s<sup>-1</sup>]</td><td>1.56 × 10<sup>-2</sup></td><td>3.89 × 10<sup>-2</sup></td><td>8.89 × 10<sup>-2</sup></td><td>1.17 × 10<sup>-1</sup></td></tr><tr><td>Δ(Φ = 0)</td><td>16.6%</td><td>15.4%</td><td>13.1%</td><td>11.4%</td></tr><tr><td>Δ(j = 0)</td><td>-0.9%</td><td>-1.7%</td><td>-3.8%</td><td>-5.2%</td></tr><tr><td>Δ(j ≤ 1)</td><td>0%</td><td>0%</td><td>0%</td><td>-0.2%</td></tr><tr><td>Δ(elastic)</td><td>0.2%</td><td>0.5%</td><td>0.4%</td><td>0.4%</td></tr><tr><td>Δ(ΓP1)</td><td>0.5%</td><td>0.5%</td><td>0.2%</td><td>0%</td></tr><tr><td>Δ(gray)</td><td>0.3%</td><td>1.1%</td><td>3.4%</td><td>5.8%</td></tr><tr><td>Δ(ΓP1)</td><td>-0.5%</td><td>-0.3%</td><td>2.3%</td><td>5.5%</td></tr><tr><td>Δ(SS)</td><td>0.2%</td><td>0.9%</td><td>2.9%</td><td>4.6%</td></tr><tr><th colspan="5">Photosensitizer</th></tr><tr><td colspan="5">Eosin Y (Φ = 0.35)</td></tr><tr><td>τ<sub>A</sub><sup>l</sup> | τ<sub>S</sub><sup>l</sup></td><td>0.1 | 0.043</td><td>0.3 | 0.128</td><td>1.4 | 0.595</td><td>4.6 | 1.956</td></tr><tr><td>P<sub>A</sub></td><td>0.0818</td><td>0.1751</td><td>0.3397</td><td>0.4934</td></tr><tr><td>⟨A⟩ [mol<sub>h,v</sub>·m<sup>-3</sup>·s<sup>-1</sup>]</td><td>1.43 × 10<sup>-2</sup></td><td>3.05 × 10<sup>-2</sup></td><td>5.92 × 10<sup>-2</sup></td><td>8.59 × 10<sup>-1</sup></td></tr><tr><td>Δ(Φ = 0)</td><td>46.5%</td><td>40.2%</td><td>29.5%</td><td>21.5%</td></tr><tr><td>Δ(j = 0)</td><td>-5.3%</td><td>-9.3%</td><td>-15.9%</td><td>-21.1%</td></tr><tr><td>Δ(j ≤ 1)</td><td>-0.6%</td><td>-1.1%</td><td>-2.6%</td><td>-4.7%</td></tr><tr><td>Δ(elastic)</td><td>0.5%</td><td>1.1%</td><td>1.7%</td><td>1.5%</td></tr><tr><td>Δ(ΓP1)</td><td>-1.5%</td><td>-1.2%</td><td>-0.2%</td><td>1.1%</td></tr><tr><td>Δ(gray)</td><td>1.7%</td><td>4.7%</td><td>11.4%</td><td>11.4%</td></tr><tr><td>Δ(ΓP1)</td><td>-1.2%</td><td>0.8%</td><td>9.9%</td><td>10.7%</td></tr><tr><td>Δ(SS)</td><td>0.7%</td><td>2.0%</td><td>2.5%</td><td>-3.0%</td></tr><tr><th colspan="5">Photosensitizer</th></tr><tr><td colspan="5">Rhodamine B (Φ = 0.7)</td></tr><tr><td>τ<sub>A</sub><sup>l</sup> | τ<sub>S</sub><sup>l</sup></td><td>0.1 | 0.317</td><td>0.3 | 0.952</td><td>1.4 | 4.444</td><td>4.6 | 14.600</td></tr><tr><td>P<sub>A</sub></td><td>0.0752</td><td>0.1379</td><td>0.2288</td><td>0.2935</td></tr><tr><td>⟨A⟩ [mol<sub>h,v</sub>·m<sup>-3</sup>·s<sup>-1</sup>]</td><td>1.31 × 10<sup>-2</sup></td><td>2.40 × 10<sup>-2</sup></td><td>3.99 × 10<sup>-2</sup></td><td>5.12 × 10<sup>-2</sup></td></tr><tr><td>Δ(Φ = 0)</td><td>170%</td><td>142%</td><td>116.5%</td><td>105.5%</td></tr><tr><td>Δ(j = 0)</td><td>-19.1%</td><td>-27.4%</td><td>-35.0%</td><td>-38.3%</td></tr><tr><td>Δ(j ≤ 1)</td><td>-3.7%</td><td>-7.8%</td><td>-12.7%</td><td>-15.2%</td></tr><tr><td>Δ(elastic)</td><td>7.7%</td><td>16.8%</td><td>24.2%</td><td>26.4%</td></tr><tr><td>Δ(ΓP1)</td><td>3.3%</td><td>13.4%</td><td>22.2%</td><td>26.3%</td></tr><tr><td>Δ(gray)</td><td>14.6%</td><td>37.7%</td><td>68.2%</td><td>70.3%</td></tr><tr><td>Δ(ΓP1)</td><td>5.5%</td><td>30.6%</td><td>65.9%</td><td>69.6%</td></tr><tr><td>Δ(SS)</td><td>4.6%</td><td>6.2%</td><td>2.6%</td><td>-1.1%</td></tr></tbody></table>

https://doi.org/10.1371/journal.pone.0255002.t002
<table><caption>Table 3. Successive orders of scattering for Ru[Bpy]₃⁺, TATA⁺, Eosin Y and Rhodamine B photosynthetizers (radiative properties are presented in Fig 2) at selected optical thicknesses τᵢₐ and τᵢₛ (see Eqs 4–7 and Fig 2): distribution P(j) = Pⱼₐ / Pₐ of scattering orders (or the probability that an absorbed photon is absorbed after j scattering events) and cumulative tail distribution P(j-q) = Σⱼ=q+1 P(j) are indicated with a gray background (see Eq 26). P(j>q) is the relative error on both absorbance Pₐ and MVRPA ⟨A⟩ when the expansion in successive scattering orders is truncated at the q-th order. The weight of luminescence radiation is P(j>0) = P(S).</caption><thead><tr><th></th><th colspan="4">Ru[Bpy]<sub>3</sub><sup>2+</sup> (Φ = 0.062)</th></tr></thead><tbody><tr><td><strong>Photosensitizer</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>τ<sub>A</sub><sup>i</sup> | τ<sub>S</sub><sup>L</sup></td><td>0.1 | 0</td><td>0.3 | 0</td><td>1.4 | 0</td><td>4.6 | 0</td></tr><tr><td>P<sup>(0)</sup></td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>P<sup>(j&gt;0)</sup> = P<sup>(S)</sup></td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td><strong>Photosensitizer</strong></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td colspan="4">TATA<sup>+</sup> (Φ = 0.15)</td></tr><tr><td>τ<sub>A</sub><sup>i</sup> | τ<sub>S</sub><sup>L</sup></td><td>0.1 | 0.006</td><td>0.3 | 0.019</td><td>1.4 | 0.088</td><td>4.6 | 0.289</td></tr><tr><td>P<sup>(0)</sup></td><td>0.991</td><td>0.983</td><td>0.962</td><td>0.948</td></tr><tr><td>P<sup>(1)</sup></td><td>0.009</td><td>0.019</td><td>0.037</td><td>0.049</td></tr><tr><td>P<sup>(j&gt;0)</sup> = P<sup>(S)</sup></td><td>0.009</td><td>0.017</td><td>0.038</td><td>0.052</td></tr><tr><td>P<sup>(j&gt;1)</sup></td><td>0</td><td>0</td><td>0</td><td>0.002</td></tr><tr><td><strong>Photosensitizer</strong></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td colspan="4">Eosin Y (Φ = 0.35)</td></tr><tr><td>τ<sub>A</sub><sup>i</sup> | τ<sub>S</sub><sup>L</sup></td><td>0.1 | 0.043</td><td>0.3 | 0.128</td><td>1.4 | 0.595</td><td>4.6 | 1.956</td></tr><tr><td>P<sup>(0)</sup></td><td>0.947</td><td>0.907</td><td>0.841</td><td>0.789</td></tr><tr><td>P<sup>(1)</sup></td><td>0.048</td><td>0.082</td><td>0.133</td><td>0.164</td></tr><tr><td>P<sup>(2)</sup></td><td>0.002</td><td>0.007</td><td>0.022</td><td>0.036</td></tr><tr><td>P<sup>(j&gt;0)</sup> = P<sup>(S)</sup></td><td>0.053</td><td>0.093</td><td>0.159</td><td>0.211</td></tr><tr><td>P<sup>(j&gt;1)</sup></td><td>0.006</td><td>0.011</td><td>0.026</td><td>0.047</td></tr><tr><td>P<sup>(j&gt;2)</sup></td><td>0.003</td><td>0.004</td><td>0.004</td><td>0.011</td></tr><tr><td><strong>Photosensitizer</strong></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td colspan="4">Rhodamine B (Φ = 0.7)</td></tr><tr><td>τ<sub>A</sub><sup>i</sup> | τ<sub>S</sub><sup>L</sup></td><td>0.1 | 0.317</td><td>0.3 | 0.952</td><td>1.4 | 4.444</td><td>4.6 | 14.600</td></tr><tr><td>P<sup>(0)</sup></td><td>0.809</td><td>0.726</td><td>0.650</td><td>0.617</td></tr><tr><td>P<sup>(1)</sup></td><td>0.154</td><td>0.196</td><td>0.223</td><td>0.232</td></tr><tr><td>P<sup>(2)</sup></td><td>0.030</td><td>0.054</td><td>0.080</td><td>0.091</td></tr><tr><td>P<sup>(3)</sup></td><td>0.006</td><td>0.016</td><td>0.029</td><td>0.036</td></tr><tr><td>P<sup>(j&gt;0)</sup> = P<sup>(S)</sup></td><td>0.191</td><td>0.274</td><td>0.350</td><td>0.383</td></tr><tr><td>P<sup>(j&gt;1)</sup></td><td>0.037</td><td>0.078</td><td>0.127</td><td>0.152</td></tr><tr><td>P<sup>(j&gt;2)</sup></td><td>0.007</td><td>0.024</td><td>0.047</td><td>0.061</td></tr><tr><td>P<sup>(j&gt;3)</sup></td><td>0</td><td>0.008</td><td>0.018</td><td>0.025</td></tr></tbody></table>

https://doi.org/10.1371/journal.pone.0255002.t003

* absorptance Pₐ provided by the reference model in Section 2, solved with the Monte Carlo algorithm presented in Section 2.5,

* MVRPA ⟨A⟩ obtained from Eq 23, for a slab of thickness L = 1 cm and incident photon flux density q₀ = 1743 µmol hv.m⁻².s⁻¹, which matches with common laboratory configurations (q₀ = 1721 µmol hv.m⁻².s⁻¹ will be used for TATA⁺, due to its slightly narrower spectral range). ⟨A⟩ for any other configuration can easily be obtained by the reader from Eq 23,

* relative difference

$$ \Delta(X) = \frac{P_{A}^{X} - P_{A}}{P_{A}} \qquad (29) $$

between reference results Pₐ and Pₐₓ values obtained with a model denoted X. Note that the relative difference for ⟨A⟩ is equal to that for Pₐ (see Eq 23),
<table><caption>Table 3. Successive orders of scattering for Ru[Bpy]₃⁺, TATA⁺, Eosin Y and Rhodamine B photosynthetizers (radiative properties are presented in Fig 2) at selected optical thicknesses τᵢₐ and τᵢₛ (see Eqs 4–7 and Fig 2): distribution P(j) = Pⱼⱼₐ / Pₐ of scattering orders (or the probability that an absorbed photon is absorbed after j scattering events) and cumulative tail distribution P(j-q) = Σⱼ=q+1 P(j) are indicated with a gray background (see Eq 26). P(j>q) is the relative error on both absorbance Pₐ and MVRPA ⟨A⟩ when the expansion in successive scattering orders is truncated at the q-th order. The weight of luminescence radiation is P(j>0) = P(S).</caption><thead><tr><th></th><th colspan="4">Ru[Bpy]<sub>3</sub><sup>2+</sup> (Φ = 0.062)</th></tr></thead><tbody><tr><td><strong>Photosensitizer</strong></td><td></td><td></td><td></td><td></td></tr><tr><td>τ<sub>A</sub><sup>i</sup> | τ<sub>S</sub><sup>L</sup></td><td>0.1 | 0</td><td>0.3 | 0</td><td>1.4 | 0</td><td>4.6 | 0</td></tr><tr><td>P<sup>(0)</sup></td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>P<sup>(j&gt;0)</sup> = P<sup>(S)</sup></td><td>0</td><td>0</td><td>0</td><td>0</td></tr><tr><td><strong>Photosensitizer</strong></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td colspan="4">TATA<sup>+</sup> (Φ = 0.15)</td></tr><tr><td>τ<sub>A</sub><sup>i</sup> | τ<sub>S</sub><sup>L</sup></td><td>0.1 | 0.006</td><td>0.3 | 0.019</td><td>1.4 | 0.088</td><td>4.6 | 0.289</td></tr><tr><td>P<sup>(0)</sup></td><td>0.991</td><td>0.983</td><td>0.962</td><td>0.948</td></tr><tr><td>P<sup>(1)</sup></td><td>0.009</td><td>0.019</td><td>0.037</td><td>0.049</td></tr><tr><td>P<sup>(j&gt;0)</sup> = P<sup>(S)</sup></td><td>0.009</td><td>0.017</td><td>0.038</td><td>0.052</td></tr><tr><td>P<sup>(j&gt;1)</sup></td><td>0</td><td>0</td><td>0</td><td>0.002</td></tr><tr><td><strong>Photosensitizer</strong></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td colspan="4">Eosin Y (Φ = 0.35)</td></tr><tr><td>τ<sub>A</sub><sup>i</sup> | τ<sub>S</sub><sup>L</sup></td><td>0.1 | 0.043</td><td>0.3 | 0.128</td><td>1.4 | 0.595</td><td>4.6 | 1.956</td></tr><tr><td>P<sup>(0)</sup></td><td>0.947</td><td>0.907</td><td>0.841</td><td>0.789</td></tr><tr><td>P<sup>(1)</sup></td><td>0.048</td><td>0.082</td><td>0.133</td><td>0.164</td></tr><tr><td>P<sup>(2)</sup></td><td>0.002</td><td>0.007</td><td>0.022</td><td>0.036</td></tr><tr><td>P<sup>(j&gt;0)</sup> = P<sup>(S)</sup></td><td>0.053</td><td>0.093</td><td>0.159</td><td>0.211</td></tr><tr><td>P<sup>(j&gt;1)</sup></td><td>0.006</td><td>0.011</td><td>0.026</td><td>0.047</td></tr><tr><td>P<sup>(j&gt;2)</sup></td><td>0.003</td><td>0.004</td><td>0.004</td><td>0.011</td></tr><tr><td><strong>Photosensitizer</strong></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td colspan="4">Rhodamine B (Φ = 0.7)</td></tr><tr><td>τ<sub>A</sub><sup>i</sup> | τ<sub>S</sub><sup>L</sup></td><td>0.1 | 0.317</td><td>0.3 | 0.952</td><td>1.4 | 4.444</td><td>4.6 | 14.600</td></tr><tr><td>P<sup>(0)</sup></td><td>0.809</td><td>0.726</td><td>0.650</td><td>0.617</td></tr><tr><td>P<sup>(1)</sup></td><td>0.154</td><td>0.196</td><td>0.223</td><td>0.232</td></tr><tr><td>P<sup>(2)</sup></td><td>0.030</td><td>0.054</td><td>0.080</td><td>0.091</td></tr><tr><td>P<sup>(3)</sup></td><td>0.006</td><td>0.016</td><td>0.029</td><td>0.036</td></tr><tr><td>P<sup>(j&gt;0)</sup> = P<sup>(S)</sup></td><td>0.191</td><td>0.274</td><td>0.350</td><td>0.383</td></tr><tr><td>P<sup>(j&gt;1)</sup></td><td>0.037</td><td>0.078</td><td>0.127</td><td>0.152</td></tr><tr><td>P<sup>(j&gt;2)</sup></td><td>0.007</td><td>0.024</td><td>0.047</td><td>0.061</td></tr><tr><td>P<sup>(j&gt;3)</sup></td><td>0</td><td>0.008</td><td>0.018</td><td>0.025</td></tr></tbody></table>

https://doi.org/10.1371/journal.pone.0255002.t003

* absorptance Pₐ provided by the reference model in Section 2, solved with the Monte Carlo algorithm presented in Section 2.5,

* MVRPA ⟨A⟩ obtained from Eq 23, for a slab of thickness L = 1 cm and incident photon flux density q₀ = 1743 µmol<sub>hv</sub>.m⁻².s⁻¹, which matches with common laboratory configurations (q₀ = 1721 µmol<sub>hv</sub>.m⁻².s⁻¹ will be used for TATA⁺, due to its slightly narrower spectral range). ⟨A⟩ for any other configuration can easily be obtained by the reader from Eq 23,

* relative difference

$$ \Delta(X) = \frac{P_{A}^{X} - P_{A}}{P_{A}} \qquad (29) $$

between reference results Pₐ and Pₐₓ values obtained with a model denoted X. Note that the relative difference for ⟨A⟩ is equal to that for Pₐ (see Eq 23),
* probability $P(j)$ that an absorbed photon is absorbed after $j$ scattering events (in Table 3), calculated according to Eq 26 where $P_A$ and $P_A^{(j)}$ are reference solutions obtained with the Monte Carlo algorithm presented in Section 2.5.

The conclusions of this analysis provide guidelines to construct relevant radiative transfer approximations in Section 4.

## 3.1 Overall luminescence effect

Here we compare absorptance values $P_A^{(\Phi=0)}$ obtained when luminescence is neglected, i.e. when $\Phi$ is set to zero in our model, with $P_A$ values obtained from the reference model. The results in Table 2 present the relative difference $\Delta(\Phi = 0)$ computed with Eq 29, where $X \equiv (\Phi = 0)$.

Neglecting luminescence always leads to an overestimation of absorptance: $\Delta(\Phi = 0) > 0$. Indeed when luminescence is not taken into account, each interaction with photosensitizers necessarily leads to absorption, whereas luminescence interaction has a probability $\Phi$ of leading to scattered photons that might exit the medium. We record $\Delta(\Phi = 0)$ values between 10% and 20% for TATA$^+$, between 20% and 50% for Eosin Y and larger than 100% for Rhodamine B. The only photosensitizer for which $\Delta(\Phi = 0)$ is constant and lower than 10% is $Ru[Bpy]_{3}^{2+}$, but we will see in the next paragraphs and in Section 4 that strictly accounting for luminescence is very simple in this particular case. As the absorption optical thickness $\tau_A^i$ increases, absorptance given by both models increases (see $P_A$ values in Table 2) and the error $\Delta(\Phi = 0)$ decreases. Nevertheless, $\Delta(\Phi = 0)$ does not approach 0 when photosensitizer concentration goes to infinity. Indeed, reference $P_A$ approaches a limit value lower than that for $P_A^{(\Phi=0)}$, because luminescence radiation cannot be completely absorbed, for two reasons. The first reason is caused by inelastic scattering: part of the luminescence radiation belongs to a spectral range where the extinction cross section is equal to zero (see Fig 2 and Section 2.3) and therefore photons exit the medium without interacting with photosensitizer. The second is driven by multiple scattering: some photons are backscattered in the vicinity of $x = 0$ and exit the medium through the front face of the reactor (even when concentration goes to infinity, or equivalently, for a semi-infinite medium).

Overall, these results indicate significant luminescence effects that should not be neglected when studying photosensitized systems.

### 3.1.1 Remark on spectral effects that are not due to luminescence.
In Table 2, for absorption optical thickness $\tau_A^i = 4.6$ it could be expected that $P_A^{(\Phi=0)} = 0.99$ [1] but this is not the case: $P_A^{(\Phi=0)} = 0.358$ for $Ru[Bpy]_{3}^{2+}$, $P_A^{(\Phi=0)} = 0.761$ for TATA$^+$, $P_A^{(\Phi=0)} = 0.599$ for Eosin Y and $P_A^{(\Phi=0)} = 0.603$ for Rhodamine B. Indeed, $\tau_A^i = 4.6$ leads to 99% of incident photons being absorbed only for gray media (or monochromatic incident radiation). Here we choose to compute absorptance in the spectral range [280 nm, 650 nm] (or [335 nm, 650 nm] for TATA$^+$), but the molar extinction coefficient is equal to zero over part of this spectral range (see Fig 2). Therefore, even when $\tau_A^i$ approaches infinity, absorptance is significantly lower than 0.99. Nevertheless, when $\tau_A^i = 4.6$ most of the incident radiation belonging to a spectral range for which $E_\lambda \neq 0$ interacts with the medium. If for each photosensitizer absorptance $P_A^{(\Phi=0)}$ is calculated in the spectral range where $E_\lambda \neq 0$, we obtain: 0.96 for $Ru[Bpy]_{3}^{2+}$, 0.86 for TATA$^+$, 0.82 for Eosin Y and 0.86 for Rhodamine B. The difference with 0.99 is due to second-order spectral effects that are not related to luminescence. Therefore, those effects will not be discussed further in this article that focuses on luminescence.
* probability $P(j)$ that an absorbed photon is absorbed after $j$ scattering events (in Table 3), calculated according to Eq 26 where $P_A$ and $P_A^{(j)}$ are reference solutions obtained with the Monte Carlo algorithm presented in Section 2.5.

The conclusions of this analysis provide guidelines to construct relevant radiative transfer approximations in Section 4.

## 3.1 Overall luminescence effect

Here we compare absorptance values $P_A^{(\Phi=0)}$ obtained when luminescence is neglected, i.e. when $\Phi$ is set to zero in our model, with $P_A$ values obtained from the reference model. The results in Table 2 present the relative difference $\Delta(\Phi = 0)$ computed with Eq 29, where $X \equiv (\Phi = 0)$.

Neglecting luminescence always leads to an overestimation of absorptance: $\Delta(\Phi = 0) > 0$. Indeed when luminescence is not taken into account, each interaction with photosensitizers necessarily leads to absorption, whereas luminescence interaction has a probability $\Phi$ of leading to scattered photons that might exit the medium. We record $\Delta(\Phi = 0)$ values between 10% and 20% for TATA$^+$, between 20% and 50% for Eosin Y and larger than 100% for Rhodamine B. The only photosensitizer for which $\Delta(\Phi = 0)$ is constant and lower than 10% is $Ru[Bpy]_{3}^{2+}$, but we will see in the next paragraphs and in Section 4 that strictly accounting for luminescence is very simple in this particular case. As the absorption optical thickness $\tau_A^i$ increases, absorptance given by both models increases (see $P_A$ values in Table 2) and the error $\Delta(\Phi = 0)$ decreases. Nevertheless, $\Delta(\Phi = 0)$ does not approach 0 when photosensitizer concentration goes to infinity. Indeed, reference $P_A$ approaches a limit value lower than that for $P_A^{(\Phi=0)}$, because luminescence radiation cannot be completely absorbed, for two reasons. The first reason is caused by inelastic scattering: part of the luminescence radiation belongs to a spectral range where the extinction cross section is equal to zero (see Fig 2 and Section 2.3) and therefore photons exit the medium without interacting with photosensitizer. The second is driven by multiple scattering: some photons are backscattered in the vicinity of $x = 0$ and exit the medium through the front face of the reactor (even when concentration goes to infinity, or equivalently, for a semi-infinite medium).

Overall, these results indicate significant luminescence effects that should not be neglected when studying photosensitized systems.

### 3.1.1 Remark on spectral effects that are not due to luminescence.
In Table 2, for absorption optical thickness $\tau_A^i = 4.6$ it could be expected that $P_A^{(\Phi=0)} = 0.99$ [1] but this is not the case: $P_A^{(\Phi=0)} = 0.358$ for $Ru[Bpy]_{3}^{2+}$, $P_A^{(\Phi=0)} = 0.761$ for TATA$^+$, $P_A^{(\Phi=0)} = 0.599$ for Eosin Y and $P_A^{(\Phi=0)} = 0.603$ for Rhodamine B. Indeed, $\tau_A^i = 4.6$ leads to 99% of incident photons being absorbed only for gray media (or monochromatic incident radiation). Here we choose to compute absorptance in the spectral range [280 nm, 650 nm] (or [335 nm, 650 nm] for TATA$^+$), but the molar extinction coefficient is equal to zero over part of this spectral range (see Fig 2). Therefore, even when $\tau_A^i$ approaches infinity, absorptance is significantly lower than 0.99. Nevertheless, when $\tau_A^i = 4.6$ most of the incident radiation belonging to a spectral range for which $E_\lambda \neq 0$ interacts with the medium. If for each photosensitizer absorptance $P_A^{(\Phi=0)}$ is calculated in the spectral range where $E_\lambda \neq 0$, we obtain: 0.96 for $Ru[Bpy]_{3}^{2+}$, 0.86 for TATA$^+$, 0.82 for Eosin Y and 0.86 for Rhodamine B. The difference with 0.99 is due to second-order spectral effects that are not related to luminescence. Therefore, those effects will not be discussed further in this article that focuses on luminescence.
## 3.2 Successive orders of scattering analysis

Here we analyze the effect of multiple luminescence scattering. In the previous paragraphs we discussed the fact that due to luminescence, interaction has a probability $\Phi$ of leading to scattered photons that might exit the medium. Therefore $P_A$ is decreased compared to situations without luminescence, leading to a loss with respect to the photoreaction. However, those scattered photons emitted by luminescence do not necessarily exit the medium. They can interact (again) with photosensitizers and be absorbed—in this case they participate in the reaction—or scattered, and so on.

Table 3 presents the distribution $P^{(j)}$ of scattering orders and its tail cumulative $P^{(j>q)}$ (see Eq 26): $P^{(0)}$ is the contribution of ballistic photons, $P^{(1)}$ is the contribution of photons absorbed after one scattering event (first order), $P^{(2)}$ is the contribution of photons absorbed after two scattering events (second order), etc., with decreasing contributions $P^{(j)} > P^{(j+1)}$ and normalization $\sum_{j=0}^{+\infty} P^{(j)} = 1$. The key parameter in analyzing this distribution is the scattering optical thickness $\tau_S^L$. The absorption of luminescence radiation increases with optical thickness, leading to a higher luminescence contribution $P^{(j>0)} = P^{(S)}$ and lower ballistic contribution $P^{(0)}$ (the distribution is normalized). Furthermore, as regards luminescence contributions, the higher the $\tau_S^L$, the larger the contribution of high scattering orders.

The following case-by-case analysis is organized according to $\tau_S^L$ value, which is strongly influenced by the overlap between extinction cross-section and luminescence emission spectra via $\bar{E}^L$ (see Eqs 6 and 7 and Fig 2).

### 3.2.1 Case 1: $\tau_S^L = 0$ because extinction and luminescence spectra do not overlap

This case is illustrated here with Ru[Bpy]$_3^{2+}$ photosensitizer (see Fig 2). The photosensitizer does not interact with scattered photons and there is no multiple scattering. Luminescence is therefore a pure loss for the reaction, since only incident solar radiation participates in absorption. In terms of successive orders of scattering, it translates into $P^{(0)} = 1$ and $P^{(j>0)} = 0$ in Table 3, or equivalently, $P_A^{(S)} = 0$ and therefore $P_A = P_A^{(0)}$ in Eq 25: zero-order expansion is sufficient to rigorously describe luminescence. Indeed, $\Delta(j=0) = 0$ for Ru[Bpy]$_3^{2+}$. Thus strictly accounting for luminescence is very simple in this particular case, as it will be discussed in Section 4.1.

### 3.2.2 Case 2: $\tau_S^L \neq 0$ because extinction and luminescence spectra overlap

This is the case for most fluorescent molecules, as illustrated here with TATA$^+$, Eosin Y and Rhodamine B. In this situation the photosensitizer interacts with scattered luminescence radiation that can be absorbed after one, two, ... $j$ scattering events:

* Looking at the TATA$^+$ case, ballistic photons $P^{(0)}$ are responsible for at least 95% of absorption. The proportion of photons $P^{(j>1)}$ absorbed after more than one scattering event is lower than 0.002 in the studied cases, indicating that ballistic photons and photons absorbed after one scattering event $P^{(1)}$ represent about 99% of photon absorption.

* In the case of Eosin Y, the photon proportion absorbed after more than one scattering event $P^{(j>1)}$ is lower than 0.05 in the studied cases, indicating that ballistic photons and photons absorbed after one scattering event contribute to 95% of absorption. For $\tau_S^L = 0.043$, we are in the same situation as for TATA$^+$: ballistic photons alone are responsible for 95% of absorption, and adding the first scattering order $P^{(1)}$ leads to 99%. For the higher optical thicknesses, ballistic photons represent 80 to 90% of absorption. For $\tau_S^L = 0.128$, again, ballistic and first-order $P^{(0)} + P^{(1)}$ are responsible for 99% of absorption. For $\tau_S^L = 0.595$ and 1.956, ballistic and first-order contribute to more than 95% of absorption. For these two last optical thicknesses, accounting for photons absorbed after 0, 1 and 2 scattering events (second-order expansion) is required to describe 99% of absorption.
## 3.2 Successive orders of scattering analysis

Here we analyze the effect of multiple luminescence scattering. In the previous paragraphs we discussed the fact that due to luminescence, interaction has a probability $\Phi$ of leading to scattered photons that might exit the medium. Therefore $P_A$ is decreased compared to situations without luminescence, leading to a loss with respect to the photoreaction. However, those scattered photons emitted by luminescence do not necessarily exit the medium. They can interact (again) with photosensitizers and be absorbed—in this case they participate in the reaction—or scattered, and so on.

Table 3 presents the distribution $P(j)$ of scattering orders and its tail cumulative $P^{(j>q)}$ (see Eq 26): $P^{(0)}$ is the contribution of ballistic photons, $P^{(1)}$ is the contribution of photons absorbed after one scattering event (first order), $P^{(2)}$ is the contribution of photons absorbed after two scattering events (second order), etc., with decreasing contributions $P^{(j)} > P^{(j+1)}$ and normalization $\sum_{j=0}^{+\infty} P^{(j)} = 1$. The key parameter in analyzing this distribution is the scattering optical thickness $\tau_S^L$. The absorption of luminescence radiation increases with optical thickness, leading to a higher luminescence contribution $P^{(j>0)} = P^{(S)}$ and lower ballistic contribution $P^{(0)}$ (the distribution is normalized). Furthermore, as regards luminescence contributions, the higher the $\tau_S^L$, the larger the contribution of high scattering orders.

The following case-by-case analysis is organized according to $\tau_S^L$ value, which is strongly influenced by the overlap between extinction cross-section and luminescence emission spectra via $\bar{E}^L$ (see Eqs 6 and 7 and Fig 2).

### 3.2.1 Case 1: $\tau_S^L = 0$ because extinction and luminescence spectra do not overlap

This case is illustrated here with Ru[Bpy]$_3^{2+}$ photosensitizer (see Fig 2). The photosensitizer does not interact with scattered photons and there is no multiple scattering. Luminescence is therefore a pure loss for the reaction, since only incident solar radiation participates in absorption. In terms of successive orders of scattering, it translates into $P^{(0)} = 1$ and $P^{(j>0)} = 0$ in Table 3, or equivalently, $P_A^{(S)} = 0$ and therefore $P_A = P_A^{(0)}$ in Eq 25: zero-order expansion is sufficient to rigorously describe luminescence. Indeed, $\Delta(j = 0) = 0$ for Ru[Bpy]$_3^{2+}$. Thus strictly accounting for luminescence is very simple in this particular case, as it will be discussed in Section 4.1.

### 3.2.2 Case 2: $\tau_S^L \neq 0$ because extinction and luminescence spectra overlap

This is the case for most fluorescent molecules, as illustrated here with TATA$^+$, Eosin Y and Rhodamine B. In this situation the photosensitizer interacts with scattered luminescence radiation that can be absorbed after one, two, ... $j$ scattering events:

* Looking at the TATA$^+$ case, ballistic photons $P^{(0)}$ are responsible for at least 95% of absorption. The proportion of photons $P^{(j>1)}$ absorbed after more than one scattering event is lower than 0.002 in the studied cases, indicating that ballistic photons and photons absorbed after one scattering event $P^{(1)}$ represent about 99% of photon absorption.

* In the case of Eosin Y, the photon proportion absorbed after more than one scattering event $P^{(j>1)}$ is lower than 0.05 in the studied cases, indicating that ballistic photons and photons absorbed after one scattering event contribute to 95% of absorption. For $\tau_S^L = 0.043$, we are in the same situation as for TATA$^+$: ballistic photons alone are responsible for 95% of absorption, and adding the first scattering order $P^{(1)}$ leads to 99%. For the higher optical thicknesses, ballistic photons represent 80 to 90% of absorption. For $\tau_S^L = 0.128$, again, ballistic and first-order $P^{(0)} + P^{(1)}$ are responsible for 99% of absorption. For $\tau_S^L = 0.595$ and 1.956, ballistic and first-order contribute to more than 95% of absorption. For these two last optical thicknesses, accounting for photons absorbed after 0, 1 and 2 scattering events (second-order expansion) is required to describe 99% of absorption.
* Regarding Rhodamine B, scattering optical thicknesses are significantly higher than for other photosensitizers, due to higher values for luminescence quantum yield and overlapping. Consequently, the contribution of photons absorbed after more than one scattering event $P^{(j>1)}$ becomes significant. For $\tau_S^L = 0.317$ and 0.952, first-order expansion is sufficient to describe 95% of absorption but the second order is required to describe 99% of absorption. For $\tau_S^L = 4.444$ and 14.6, second-order expansion describes 95% of absorption but the third and fourth orders (photons absorbed after three and four scattering events) are required to describe 99% of absorption.

In summary, these results indicate that in the case investigated here 95% of photons absorbed are:

* ballistic photons only when $\tau_S^L < 0.05$,

* ballistic photons and photons absorbed after one scattering event when $0.05 < \tau_S^L < 1$,

* ballistic photons and photons absorbed after one and two scattering events when $\tau_S^L > 1$,

### 3.3 Inelastic scattering sensitivity

Inelastic scattering is another conceptual and numerical difficulty when describing photon transport in photosensitized and photocatalytic systems, in addition to multiple scattering. Indeed, wavelengths are redistributed according to the luminescence spectrum at each scattering event. Here we analyze the sensitivity of absorptance and MVRPA to inelastic scattering effects in order to determine if this difficulty must be tackled or can be partially bypassed. For this purpose, we constructed two equivalent models for luminescence radiation that we compared with the reference solution (see $\Delta(\text{gray})$ and $\Delta(\text{elastic})$ in Table 2). In the first one, gray approximation is used in the luminescence spectral range and in the second, the inelastic collision term in the RTE for luminescence intensity is simply replaced by an elastic collision term. In both cases, the spectral distribution for ballistic photons (zero-th order) is preserved; only luminescence intensity is affected by the approximations.

A gray model for luminescence radiation is constructed based on the distinction between ballistic and scattered intensity (see Eq 18). Gray approximation is used in the RTE for scattered (i.e. luminescence) intensity $I_\lambda^{(S)}(x, \mu)$, which turns the inelastic collision term into an elastic collision term.

The RTE in Eq 19 is integrated over wavelengths $\lambda$ and then the following gray approximation is used:

$$ \int_{0}^{+\infty} d\lambda k_{\lambda} I_{\lambda}^{(S)}(x, \mu; k_{\lambda}) \approx \bar{k} I^{(S)}(x, \mu; \bar{k}) \quad (30) $$

where $I^{(S)} = \int_{0}^{+\infty} d\lambda I_{\lambda}^{(S)}$ and $\bar{k} = \int_{0}^{+\infty} d\lambda p^{L}(\lambda) k_{\lambda} = \bar{E}^{L}C$ (see Eq 7). Using the normalization $\int_{0}^{+\infty} d\lambda p^{L}(\lambda) = 1$ and the notation $\bar{I}^{(S)}$ for this gray approximation leads to

$$ \mu \frac{d\bar{I}^{(S)}(x, \mu)}{dx} = -\bar{k} \bar{I}^{(S)}(x, \mu) + \frac{\Phi \bar{k}}{2} \int_{-1}^{1} d\mu' \bar{I}^{(S)}(x, \mu') + \int_{0}^{+\infty} d\lambda' \frac{\Phi k_{\lambda'}}{2} \int_{-1}^{1} d\mu' I_{\lambda'}^{(0)}(x, \mu') \quad (31) $$

with the boundary condition

$$ \bar{I}^{(S)}(0, \mu) = 0 \quad (32) $$

Note that $I_\lambda^{(0)}$ in Eq 31 is the solution of $\text{Ru}[\text{Bpy}]_3^{2+}$ Eq 12: ballistic intensity is not affected by the present approximation. Then the contribution of scattered photons to absorptance is
* Regarding Rhodamine B, scattering optical thicknesses are significantly higher than for other photosensitizers, due to higher values for luminescence quantum yield and overlapping. Consequently, the contribution of photons absorbed after more than one scattering event $P^{(j>1)}$ becomes significant. For $\tau_S^L = 0.317$ and 0.952, first-order expansion is sufficient to describe 95% of absorption but the second order is required to describe 99% of absorption. For $\tau_S^L = 4.444$ and 14.6, second-order expansion describes 95% of absorption but the third and fourth orders (photons absorbed after three and four scattering events) are required to describe 99% of absorption.

In summary, these results indicate that in the case investigated here 95% of photons absorbed are:

* ballistic photons only when $\tau_S^L < 0.05$,

* ballistic photons and photons absorbed after one scattering event when $0.05 < \tau_S^L < 1$,

* ballistic photons and photons absorbed after one and two scattering events when $\tau_S^L > 1$,

### 3.3 Inelastic scattering sensitivity

Inelastic scattering is another conceptual and numerical difficulty when describing photon transport in photosensitized and photocatalytic systems, in addition to multiple scattering. Indeed, wavelengths are redistributed according to the luminescence spectrum at each scattering event. Here we analyze the sensitivity of absorptance and MVRPA to inelastic scattering effects in order to determine if this difficulty must be tackled or can be partially bypassed. For this purpose, we constructed two equivalent models for luminescence radiation that we compared with the reference solution (see $\Delta(\text{gray})$ and $\Delta(\text{elastic})$ in Table 2). In the first one, gray approximation is used in the luminescence spectral range and in the second, the inelastic collision term in the RTE for luminescence intensity is simply replaced by an elastic collision term. In both cases, the spectral distribution for ballistic photons (zero-th order) is preserved; only luminescence intensity is affected by the approximations.

A gray model for luminescence radiation is constructed based on the distinction between ballistic and scattered intensity (see Eq 18). Gray approximation is used in the RTE for scattered (i.e. luminescence) intensity $I_\lambda^{(S)}(x, \mu)$, which turns the inelastic collision term into an elastic collision term.

The RTE in Eq 19 is integrated over wavelengths $\lambda$ and then the following gray approximation is used:

$$ \int_{0}^{+\infty} d\lambda k_{\lambda} I_{\lambda}^{(S)}(x, \mu; k_{\lambda}) \approx \bar{k} I^{(S)}(x, \mu; \bar{k}) \quad (30) $$

where $I^{(S)} = \int_{0}^{+\infty} d\lambda I_{\lambda}^{(S)}$ and $\bar{k} = \int_{0}^{+\infty} d\lambda p^{L}(\lambda) k_{\lambda} = \bar{E}^{L}C$ (see Eq 7). Using the normalization $\int_{0}^{+\infty} d\lambda p^{L}(\lambda) = 1$ and the notation $\bar{I}^{(S)}$ for this gray approximation leads to

$$ \mu \frac{d\bar{I}^{(S)}(x, \mu)}{dx} = -\bar{k} \bar{I}^{(S)}(x, \mu) + \frac{\Phi \bar{k}}{2} \int_{-1}^{1} d\mu' \bar{I}^{(S)}(x, \mu') + \int_{0}^{+\infty} d\lambda' \frac{\Phi k_{\lambda'}}{2} \int_{-1}^{1} d\mu' I_{\lambda'}^{(0)}(x, \mu') \quad (31) $$

with the boundary condition

$$ \bar{I}^{(S)}(0, \mu) = 0 \quad (32) $$

Note that $I_\lambda^{(0)}$ in Eq 31 is the solution of $\text{Ru}[\text{Bpy}]_3^{2+}$ Eq 12: ballistic intensity is not affected by the present approximation. Then the contribution of scattered photons to absorptance is
obtained by applying the gray approximation Eq 30 to the definition in Eq 24:

$$ \bar{P}_{A}^{(S)} = \frac{1}{q_0} \int_{0}^{L} dx \int_{-1}^{1} d\mu (1 - \Phi) \bar{k} \bar{I}^{(S)}(x, \mu) \quad (33) $$

Finally, our gray model for absorptance is obtained by summing the contributions of ballistic and scattered photons:

$$ \mathcal{P}_{A}^{(\text{gray})} = \mathcal{P}_{A}^{(0)} + \bar{\mathcal{P}}_{A}^{(S)} \quad (34) $$

where $\mathcal{P}_{A}^{(0)}$ is not modified by the approximation (see its definition in Eq 24).

The $\mathcal{P}_{A}^{(\text{gray})}$ values presented in Table 2 were obtained with the Monte Carlo algorithm presented in Section 2.5, except that **Step 4a.3** was replaced by:

* **Step 4a.3:** $k_{\lambda} = \bar{k}$

Whatever the wavelength, the extinction coefficient for scattered photons is set to $\bar{k} = \bar{E}^L C$ (see Fig 2).

The relative difference $\Delta(\text{gray})$ between $\mathcal{P}_{A}^{(\text{gray})}$ and the reference $\mathcal{P}_{A}$ is presented in Table 2. Results do not include Ru[Bpy]$^{2+}$, since luminescence radiation does not participate in absorption in this case (see Section 2.3). For other photosensitizers, the fact that $\Delta(\text{gray}) > 0$ (always) indicates that gray approximation for luminescence radiation overestimates absorptance (and indeed, $f(k_{\lambda}) = \int_{0}^{L} dx \int_{-1}^{1} d\mu (1 - \Phi)k_{\lambda}I_{\lambda}^{(S)}(x, \mu)$ in Eq 24 is concave).

$\Delta(\text{gray})$ rapidly increases with the scattering optical thickness $\tau_{S}^{L}$, because the average number of wavelength redistributions due to inelastic scattering events increases with $\tau_{S}^{L}$. Indeed, the analysis in Section 2.3 insists on the significant contribution of scattering orders $j = 1, 2, 3$ in some tested cases, involving $j$ wavelength redistribution events (see Table 3). For TATA$^{+}$, $\Delta(\text{gray}) \lesssim 5$ since at least 95% of absorption is due to ballistic radiation that is not affected by inelastic scattering. For Eosin Y, the first and second scattering orders contribute up to 20% and $\Delta(\text{gray})$ reaches 11%. For Rhodamine B, we record a significant contribution up to the third scattering order, and $\Delta(\text{gray})$ reaches 70%.

$\Delta(\text{gray})$ always increases with $\tau_{S}^{L}$ for a given set of radiative properties, but the magnitude of that increase depends on the extinction cross section and luminescence spectra. For example, if the extinction cross section is gray within the luminescence spectral range (i.e. $E_{\lambda}$ does not vary with $\lambda$ in the spectral range $[\lambda_{min}^{L}, \lambda_{max}^{L}]$ where $p^{L}(\lambda) \neq 0$), then the gray model and reference model are exactly identical, and therefore $\Delta(\text{gray}) = 0 \forall \tau_{S}^{L}$. In general, we expect that the higher the spectral variation of $E_{\lambda}$ on $[\lambda_{min}^{L}, \lambda_{max}^{L}]$ is, the faster $\Delta(\text{gray})$ increases with $\tau_{S}^{L}$. We cannot illustrate this behavior here because the typical photosensitizers we selected have comparable spectral variations over the luminescence spectral range (see Fig 2).

To conclude, for the typical cases tested in this paper, wavelength redistribution due to luminescence can hardly be put aside once scattering optical thickness exceeds 0.1. However, the results in Table 3 indicate that the main contribution to luminescence radiation is due to first-order terms in scattering expansion. Therefore we hereafter study a model that preserves luminescence spectral redistribution during the first scattering event only, and neglects the subsequent ones.

**The elastic-scattering model for luminescence radiation** is simply obtained by replacing the inelastic collision term in the RTE for luminescence intensity (see Eq 19) by an elastic
obtained by applying the gray approximation Eq 30 to the definition in Eq 24:

$$ \bar{P}_{A}^{(S)} = \frac{1}{q_0} \int_{0}^{L} dx \int_{-1}^{1} d\mu (1 - \Phi) \bar{k} \bar{I}^{(S)}(x, \mu) \quad (33) $$

Finally, our gray model for absorptance is obtained by summing the contributions of ballistic and scattered photons:

$$ \mathcal{P}_{A}^{(\text{gray})} = \mathcal{P}_{A}^{(0)} + \bar{P}_{A}^{(S)} \quad (34) $$

where $\mathcal{P}_{A}^{(0)}$ is not modified by the approximation (see its definition in Eq 24).

The $\mathcal{P}_{A}^{(\text{gray})}$ values presented in Table 2 were obtained with the Monte Carlo algorithm presented in Section 2.5, except that **Step 4a.3** was replaced by:

* **Step 4a.3:** $k_{\lambda} = \bar{k}$

Whatever the wavelength, the extinction coefficient for scattered photons is set to $\bar{k} = \bar{E}^L C$ (see Fig 2).

The relative difference $\Delta(\text{gray})$ between $\mathcal{P}_{A}^{(\text{gray})}$ and the reference $\mathcal{P}_{A}$ is presented in Table 2. Results do not include Ru[Bpy]$^{2+}$, since luminescence radiation does not participate in absorption in this case (see Section 2.3). For other photosensitizers, the fact that $\Delta(\text{gray}) > 0$ (always) indicates that gray approximation for luminescence radiation overestimates absorptance (and indeed, $f(k_{\lambda}) = \int_{0}^{L} dx \int_{-1}^{1} d\mu (1 - \Phi)k_{\lambda}I_{\lambda}^{(S)}(x, \mu)$ in Eq 24 is concave).

$\Delta(\text{gray})$ rapidly increases with the scattering optical thickness $\tau_{S}^{L}$, because the average number of wavelength redistributions due to inelastic scattering events increases with $\tau_{S}^{L}$. Indeed, the analysis in Section 2.3 insists on the significant contribution of scattering orders $j = 1, 2, 3$ in some tested cases, involving $j$ wavelength redistribution events (see Table 3). For TATA$^{+}$, $\Delta(\text{gray}) \lesssim 5$ since at least 95% of absorption is due to ballistic radiation that is not affected by inelastic scattering. For Eosin Y, the first and second scattering orders contribute up to 20% and $\Delta(\text{gray})$ reaches 11%. For Rhodamine B, we record a significant contribution up to the third scattering order, and $\Delta(\text{gray})$ reaches 70%.

$\Delta(\text{gray})$ always increases with $\tau_{S}^{L}$ for a given set of radiative properties, but the magnitude of that increase depends on the extinction cross section and luminescence spectra. For example, if the extinction cross section is gray within the luminescence spectral range (i.e. $E_{\lambda}$ does not vary with $\lambda$ in the spectral range $[\lambda_{min}^{L}, \lambda_{max}^{L}]$ where $p^{L}(\lambda) \neq 0$), then the gray model and reference model are exactly identical, and therefore $\Delta(\text{gray}) = 0 \forall \tau_{S}^{L}$. In general, we expect that the higher the spectral variation of $E_{\lambda}$ on $[\lambda_{min}^{L}, \lambda_{max}^{L}]$ is, the faster $\Delta(\text{gray})$ increases with $\tau_{S}^{L}$. We cannot illustrate this behavior here because the typical photosensitizers we selected have comparable spectral variations over the luminescence spectral range (see Fig 2).

To conclude, for the typical cases tested in this paper, wavelength redistribution due to luminescence can hardly be put aside once scattering optical thickness exceeds 0.1. However, the results in Table 3 indicate that the main contribution to luminescence radiation is due to first-order terms in scattering expansion. Therefore we hereafter study a model that preserves luminescence spectral redistribution during the first scattering event only, and neglects the subsequent ones.

The elastic-scattering model for luminescence radiation is simply obtained by replacing the inelastic collision term in the RTE for luminescence intensity (see Eq 19) by an elastic
collision term:

$$ \mu \frac{d\hat{I}_{\lambda}^{(S)}(x, \mu)}{dx} = -k_{\lambda} \hat{I}_{\lambda}^{(S)}(x, \mu) + \frac{\Phi k_{\lambda}}{2} \int_{-1}^{1} d\mu' \hat{I}_{\lambda}^{(S)}(x, \mu') \\ + p^{L}(\lambda) \int_{0}^{+\infty} d\lambda' \frac{\Phi k_{\lambda'}}{2} \int_{-1}^{1} d\mu' \hat{I}_{\lambda'}^{(L)}(x, \mu') \quad (35) $$

Luminescence absorptance $\hat{P}_{A}^{(S)}$ is defined in Eq 24 with $I_{\lambda}^{(\bullet)} = \hat{I}^{(S)}$ and

$$ \mathcal{P}_{A}^{(\text{elastic})} = \mathcal{P}_{A}^{(0)} + \hat{P}_{A}^{(S)} \quad (36) $$

where $\mathcal{P}_{A}^{(0)}$ is not affected by the approximation.

The $\mathcal{P}_{A}^{(\text{elastic})}$ values presented in Table 2 were obtained using the Monte Carlo algorithm presented in Section 2.5, except that **Step 4a.3** was replaced by:

**Step 4a.3:**

if $j = 1$ then

The wavelength $\lambda$ is sampled according to the luminescence spectrum $p^{L}(\lambda)$ and the extinction coefficient $k_{\lambda} = CE_{\lambda}$ is interpolated.

end

The wavelength $\lambda$ is sampled during the first scattering event only (i.e. for $j = 1$).

Two spectral integrations are required to solve this model—one over the incident spectrum and one over the luminescence spectrum—whereas only spectral integration over the incident spectrum was considered in the previous gray model. This additional spectral dimension is required to take into account the luminescence spectrum when ballistic photons are scattered and move from population (0) to population (1). The comparison of $\Delta$ (elastic) and $\Delta$ (gray) in Table 2 indicates that it is a significant improvement in the description of absorption, and yet it corresponds to a considerable conceptual and practical simplification compared to the reference model, which involves an infinite number of spectral integrations (one at each scattering order). This simplified description of inelastic scattering provides results with a relative difference $\Delta$ (elastic) below 2% for TATA+ and Eosin Y. This is due to the fact that scattering orders $j > 1$ have little impact with these photosensitizers (see Table 3). For Rhodamine, however, these scattering orders contribute up to 15% of absorbance, in addition to the strong spectral variation in the extinction cross section in the luminescence spectral range [$\lambda_{min}^{L}, \lambda_{max}^{L}$] (see Fig 2d). In this case, inelastic scattering should be described in all its complexity.

## 4 Analytical approximate solutions

Based on the above radiative analysis, we propose in the present section to develop five approximations and discuss their behavior. We aim at providing meaningful analytical solutions enabling a simple evaluation and analysis of luminescence radiation in most cases, without having to use advanced numerical methods. This work is based on zero-order scattering, P1 and single-scattering approximations. But first let us clarify a common practice in photochemistry, which consists in neglecting luminescence.

Neglecting luminescence is equivalent to taking $\Phi = 0$ in the reference RTE Eq 8, with the same boundary conditions as in Eq 9:

$$ \mu \frac{dI_{\lambda}^{(\Phi=0)}(x, \mu)}{dx} = -k_{\lambda} I_{\lambda}^{(\Phi=0)}(x, \mu) \quad (37) $$
collision term:

$$ \mu \frac{d\hat{I}_{\lambda}^{(S)}(x, \mu)}{dx} = -k_{\lambda} \hat{I}_{\lambda}^{(S)}(x, \mu) + \frac{\Phi k_{\lambda}}{2} \int_{-1}^{1} d\mu' \hat{I}_{\lambda}^{(S)}(x, \mu') \\ + p^{L}(\lambda) \int_{0}^{+\infty} d\lambda' \frac{\Phi k_{\lambda'}}{2} \int_{-1}^{1} d\mu' \hat{I}_{\lambda'}^{(L)}(x, \mu') \quad (35) $$

Luminescence absorptance $\hat{P}_{A}^{(S)}$ is defined in Eq 24 with $I_{\lambda}^{(\bullet)} = \hat{I}^{(S)}$ and

$$ \mathcal{P}_{A}^{(\text{elastic})} = \mathcal{P}_{A}^{(0)} + \hat{P}_{A}^{(S)} \quad (36) $$

where $\mathcal{P}_{A}^{(0)}$ is not affected by the approximation.

The $\mathcal{P}_{A}^{(\text{elastic})}$ values presented in Table 2 were obtained using the Monte Carlo algorithm presented in Section 2.5, except that **Step 4a.3** was replaced by:

**Step 4a.3:**

if $j = 1$ then

The wavelength $\lambda$ is sampled according to the luminescence spectrum $p^{L}(\lambda)$ and the extinction coefficient $k_{\lambda} = CE_{\lambda}$ is interpolated.

end

The wavelength $\lambda$ is sampled during the first scattering event only (i.e. for $j = 1$).

Two spectral integrations are required to solve this model—one over the incident spectrum and one over the luminescence spectrum—whereas only spectral integration over the incident spectrum was considered in the previous gray model. This additional spectral dimension is required to take into account the luminescence spectrum when ballistic photons are scattered and move from population (0) to population (1). The comparison of $\Delta$ (elastic) and $\Delta$ (gray) in Table 2 indicates that it is a significant improvement in the description of absorption, and yet it corresponds to a considerable conceptual and practical simplification compared to the reference model, which involves an infinite number of spectral integrations (one at each scattering order). This simplified description of inelastic scattering provides results with a relative difference $\Delta$ (elastic) below 2% for TATA+ and Eosin Y. This is due to the fact that scattering orders $j > 1$ have little impact with these photosensitizers (see Table 3). For Rhodamine, however, these scattering orders contribute up to 15% of absorbance, in addition to the strong spectral variation in the extinction cross section in the luminescence spectral range [$\lambda_{min}^{L}, \lambda_{max}^{L}$] (see Fig 2d). In this case, inelastic scattering should be described in all its complexity.

## 4 Analytical approximate solutions

Based on the above radiative analysis, we propose in the present section to develop five approximations and discuss their behavior. We aim at providing meaningful analytical solutions enabling a simple evaluation and analysis of luminescence radiation in most cases, without having to use advanced numerical methods. This work is based on zero-order scattering, P1 and single-scattering approximations. But first let us clarify a common practice in photochemistry, which consists in neglecting luminescence.

Neglecting luminescence is equivalent to taking $\Phi = 0$ in the reference RTE Eq 8, with the same boundary conditions as in Eq 9:

$$ \mu \frac{dI_{\lambda}^{(\Phi=0)}(x, \mu)}{dx} = -k_{\lambda} I_{\lambda}^{(\Phi=0)}(x, \mu) \quad (37) $$
leading to a simple Bouguer exponential attenuation along straight paths:

$$I_{\lambda}^{(\Phi=0)}(x, \mu = 1) = q_{0,\lambda} e^{-k_{\lambda} x} \qquad (38)$$

$$I_{\lambda}^{(\Phi=0)}(x, \mu \neq 1) = 0 \qquad (39)$$

Substituting this expression into Eq 24 with $\Phi = 0$ leads to

$$\mathcal{P}_{A}^{(\Phi=0)} = \int_{0}^{+\infty} d\lambda \, p^{i}(\lambda) (1 - e^{-k_{\lambda} L}) \qquad (40)$$

With this model, the extinction coefficient $k_{\lambda}$ is equal to the absorption coefficient, and radiation is collimated along the $e_x$ axis (no deviation due to scattering here).

As discussed in Section 3.1, we observe serious $\Delta(\Phi = 0)$ errors when neglecting luminescence (see Table 2). A straightforward and significant improvement is obtained by using instead the 0th order scattering approximation.

## 4.1 Zero-order scattering approximation

In Section 3.2 we observed that the main contribution to absorption is due to ballistic radiation. The expansion in successive orders of scattering is therefore truncated here at order zero: $\mathcal{P}_A \simeq \mathcal{P}_A^{(0)}$. In comparison with the reference solution in Eq 25, the absorption of luminescence radiation $\mathcal{P}_A^{(S)}$ is not taken into account, leading to quite a simple analytic expression: the solution of the RTE Eq 12 for ballistic intensity $I_{\lambda}^{(0)}$ is

$$I_{\lambda}^{(0)}(x, \mu = 1) = q_{0,\lambda} e^{-k_{\lambda} x} \qquad (41)$$

$$I_{\lambda}^{(0)}(x, \mu \neq 1) = 0 \qquad (42)$$

and substituting this expression into Eq 24 leads to

$$\mathcal{P}_{A}^{(0)} = (1 - \Phi) \int_{0}^{+\infty} d\lambda \, p^{i}(\lambda) (1 - e^{-k_{\lambda} L}) = (1 - \Phi) \mathcal{P}_{A}^{(\Phi=0)} \qquad (43)$$

Let us emphasize that $\mathcal{P}_A^{(0)}$ is very easily obtained by multiplying $\mathcal{P}_A^{(\Phi=0)}$, which is routinely calculated or measured in photochemistry practice, by $(1 - \Phi)$.

The solution for $I_{\lambda}^{(0)}$ is the same as for $I_{\lambda}^{(\Phi=0)}$, but a crucial difference is that the luminescence quantum yield $\Phi$ is here taken into account in the expression of $\mathcal{P}_A$ and $\langle A \rangle$. When luminescence is neglected, every photon that interacts with the photosensitizer is counted as being absorbed, leading to an overestimation of $\mathcal{P}_A: \Delta(\Phi = 0) > 0$ in Table 2. Here, in contrast, luminescence scattering is modeled but only ballistic photon absorption is accounted for, leading to the prefactor $(1 - \Phi)$ in Eq 43. Since the absorption of luminescence radiation is ignored, the zero-order scattering approximation underestimates absorbance: $\Delta(j = 0) < 0$ in Table 2. The impact of truncating the expansion in successive orders of scattering has already been discussed in Sections 2.3 and 3.2, and indeed, $\Delta(j = 0)$ in Table 2 is by definition equal to $-P$ ($j > 0$) in Table 3. As seen in previous discussions, the error associated with zero-order scattering approximation is governed by the scattering optical thickness $\tau_S^L$: the higher $\tau_S^L$ is, the higher the error $\Delta(j = 0)$.

Of course, when $\tau_S^L = 0$, the zero-order scattering approximation gives exact results: $I_{\lambda} = I_{\lambda}^{(0)}$, $\langle A \rangle = \langle A \rangle^{(0)}$ and $\mathcal{P}_A = \mathcal{P}_A^{(0)}$. This situation is encountered in two cases. First, for photosensitizers that do not emit luminescence radiation (i.e. when $\Phi = 0$). In this case
leading to a simple Bouguer exponential attenuation along straight paths:

$$I_{\lambda}^{(\Phi=0)}(x, \mu = 1) = q_{0,\lambda} e^{-k_{\lambda} x} \qquad (38)$$

$$I_{\lambda}^{(\Phi=0)}(x, \mu \neq 1) = 0 \qquad (39)$$

Substituting this expression into Eq 24 with $\Phi = 0$ leads to

$$\mathcal{P}_{A}^{(\Phi=0)} = \int_{0}^{+\infty} d\lambda \, p^{i}(\lambda) (1 - e^{-k_{\lambda}L}) \qquad (40)$$

With this model, the extinction coefficient $k_{\lambda}$ is equal to the absorption coefficient, and radiation is collimated along the $e_x$ axis (no deviation due to scattering here).

As discussed in Section 3.1, we observe serious $\Delta(\Phi = 0)$ errors when neglecting luminescence (see Table 2). A straightforward and significant improvement is obtained by using instead the 0th order scattering approximation.

## 4.1 Zero-order scattering approximation

In Section 3.2 we observed that the main contribution to absorption is due to ballistic radiation. The expansion in successive orders of scattering is therefore truncated here at order zero: $\mathcal{P}_A \simeq \mathcal{P}_A^{(0)}$. In comparison with the reference solution in Eq 25, the absorption of luminescence radiation $\mathcal{P}_A^{(S)}$ is not taken into account, leading to quite a simple analytic expression: the solution of the RTE Eq 12 for ballistic intensity $I_{\lambda}^{(0)}$ is

$$I_{\lambda}^{(0)}(x, \mu = 1) = q_{0,\lambda} e^{-k_{\lambda}x} \qquad (41)$$

$$I_{\lambda}^{(0)}(x, \mu \neq 1) = 0 \qquad (42)$$

and substituting this expression into Eq 24 leads to

$$\mathcal{P}_{A}^{(0)} = (1 - \Phi) \int_{0}^{+\infty} d\lambda \, p^{i}(\lambda) (1 - e^{-k_{\lambda}L}) = (1 - \Phi) \mathcal{P}_{A}^{(\Phi=0)} \qquad (43)$$

Let us emphasize that $\mathcal{P}_A^{(0)}$ is very easily obtained by multiplying $\mathcal{P}_A^{(\Phi=0)}$, which is routinely calculated or measured in photochemistry practice, by $(1 - \Phi)$.

The solution for $I_{\lambda}^{(0)}$ is the same as for $I_{\lambda}^{(\Phi=0)}$, but a crucial difference is that the luminescence quantum yield $\Phi$ is here taken into account in the expression of $\mathcal{P}_A$ and $\langle A \rangle$. When luminescence is neglected, every photon that interacts with the photosensitizer is counted as being absorbed, leading to an overestimation of $\mathcal{P}_A: \Delta(\Phi = 0) > 0$ in Table 2. Here, in contrast, luminescence scattering is modeled but only ballistic photon absorption is accounted for, leading to the prefactor $(1 - \Phi)$ in Eq 43. Since the absorption of luminescence radiation is ignored, the zero-order scattering approximation underestimates absorbance: $\Delta(j = 0) < 0$ in Table 2. The impact of truncating the expansion in successive orders of scattering has already been discussed in Sections 2.3 and 3.2, and indeed, $\Delta(j = 0)$ in Table 2 is by definition equal to $-P$ ($j > 0$) in Table 3. As seen in previous discussions, the error associated with zero-order scattering approximation is governed by the scattering optical thickness $\tau_S^L$: the higher $\tau_S^L$ is, the higher the error $\Delta(j = 0)$.

Of course, when $\tau_S^L = 0$, the zero-order scattering approximation gives exact results: $I_{\lambda} = I_{\lambda}^{(0)}$, $\langle A \rangle = \langle A \rangle^{(0)}$ and $\mathcal{P}_A = \mathcal{P}_A^{(0)}$. This situation is encountered in two cases. First, for photosensitizers that do not emit luminescence radiation (i.e. when $\Phi = 0$). In this case
$I_{\lambda} = I_{\lambda}^{(\Phi=0)} = I_{\lambda}^{(0)}$, $\langle A \rangle = \langle A \rangle^{(\Phi=0)} = \langle A \rangle^{(0)}$ and $P_A = P_A^{(\Phi=0)} = P_A^{(0)}$. Second, and more importantly, for photosensitizers with extinction cross-section and luminescence emission spectra that do not overlap, as illustrated with $Ru[Bpy]_{3}^{2+}$ in Section 3.1 (see Case 1). In this case, the error obtained with a model neglecting luminescence is $\Delta(\Phi = 0) = (P_A^{(\Phi=0)} - P_A^{(0)})/P_A^{(0)} = \Phi/(1 - \Phi)$. In Table 2, indeed, $\Delta(\Phi = 0) = 6.6\%$ for $Ru[Bpy]_{3}^{2+}$ with luminescence quantum yield $\Phi = 0.062$. Note that a photosensitizer with a higher $\Phi$ could lead to a significantly greater error: e.g. $\Delta(\Phi = 0) = 25\%$ for $\Phi = 0.2$. In these situations where $\tau_S^L = 0$, the very simple solution provided in Eq 43 should be used.

In Table 2, zero-order scattering approximation is also tested for three photosensitizers with $\tau_S^L \neq 0$: TATA+, Eosin Y and Rhodamine B. As discussed in Section 3.1, the luminescence effect is significant in these cases and should not be neglected ($\Delta(\Phi = 0) = 10\%$ to $170\%$):

* For TATA+, the zero-order scattering approximation provides accurate results with error $\Delta(j = 0) \le 5\%$ in every tested situation. Indeed, Table 3 indicates that at least 95% of absorption is due to ballistic photons. The relative error is 2 to 20 times lower than the error $\Delta(\Phi = 0)$ obtained when ignoring luminescence.

* For Eosin Y and Rhodamine B, even if the zero-order scattering approximation is not sufficient to describe photon transport with accuracy below 5%, we nevertheless record a significant improvement in relative error compared to $\Delta(\Phi = 0)$.

Overall, we record the following accuracy levels for the zero-order scattering approximation as a function of scattering optical thickness $\tau_S^L$:

* $\Delta(j = 0) < 5\%$ for $\tau_S^L < 0.05$,

* $\Delta(j = 0) < 10\%$ for $0.05 < \tau_S^L < 0.3$,

* $\Delta(j = 0) < 30\%$ for $0.3 < \tau_S^L < 1$.

In conclusion, zero-order scattering is a relevant approximation that provides accurate results for optical thicknesses $\tau_S^L < 0.05$ and for TATA+ in general; it includes half the cases investigated in this paper. In addition, this approximation is exceedingly simple because it only requires multiplying $P_A^{(\Phi=0)}$, which is routinely calculated or measured in photochemistry practice, by $(1 - \Phi)$.

## 4.2 P1 approximations for luminescence radiation

For optical thicknesses $\tau_S^L > 0.05$, the absorption of luminescence radiation must be described to improve accuracy compared to zero-order approximation. Here, the luminescence radiation $P_A^{(S)}$ in Eq 25 is described by diffusion equations derived thanks to P1 approximation [35, 50]. But a simple analytical solution is only accessible for elastic scattering processes, and therefore P1 approximation is applied to the gray and elastic models presented in Section 3.3.

### 4.2.1 P1 approximation applied to the gray model for luminescence radiation.
The gray approximation presented in Section 3.3 produces fairly accurate results for TATA+ and Eosin Y (see $\Delta(\text{gray})$ in Table 2) and is also expected to be relevant in the case of photosensitizers with small spectral variations in extinction cross section over the luminescence spectral range. Therefore the P1 approximation is applied to the gray RTE in Eq 31. The trivial reference solution for ballistic photons (order zero) is retained, including incident spectral distribution, but wavelength redistribution due to luminescence is not taken into account: $P_A^{(S)} \approx \bar{P}_A^{(S)}$. Hereafter, a simple analytic expression for P1 approximation $P_A^{(S)} \approx \bar{P}_A^{(S,P1)}$ is obtained and we
$I_{\lambda} = I_{\lambda}^{(\Phi=0)} = I_{\lambda}^{(0)}$, $\langle A \rangle = \langle A \rangle^{(\Phi=0)} = \langle A \rangle^{(0)}$ and $P_A = P_A^{(\Phi=0)} = P_A^{(0)}$. Second, and more importantly, for photosensitizers with extinction cross-section and luminescence emission spectra that do not overlap, as illustrated with $Ru[Bpy]_{3}^{2+}$ in Section 3.1 (see Case 1). In this case, the error obtained with a model neglecting luminescence is $\Delta(\Phi = 0) = (P_A^{(\Phi=0)} - P_A^{(0)})/P_A^{(0)} = \Phi/(1 - \Phi)$. In Table 2, indeed, $\Delta(\Phi = 0) = 6.6\%$ for $Ru[Bpy]_{3}^{2+}$ with luminescence quantum yield $\Phi = 0.062$. Note that a photosensitizer with a higher $\Phi$ could lead to a significantly greater error: e.g. $\Delta(\Phi = 0) = 25\%$ for $\Phi = 0.2$. In these situations where $\tau_S^L = 0$, the very simple solution provided in Eq 43 should be used.

In Table 2, zero-order scattering approximation is also tested for three photosensitizers with $\tau_S^L \neq 0$: TATA+, Eosin Y and Rhodamine B. As discussed in Section 3.1, the luminescence effect is significant in these cases and should not be neglected ($\Delta(\Phi = 0) = 10\%$ to $170\%$):

* For TATA+, the zero-order scattering approximation provides accurate results with error $\Delta(j = 0) \le 5\%$ in every tested situation. Indeed, Table 3 indicates that at least 95% of absorption is due to ballistic photons. The relative error is 2 to 20 times lower than the error $\Delta(\Phi = 0)$ obtained when ignoring luminescence.

* For Eosin Y and Rhodamine B, even if the zero-order scattering approximation is not sufficient to describe photon transport with accuracy below 5%, we nevertheless record a significant improvement in relative error compared to $\Delta(\Phi = 0)$.

Overall, we record the following accuracy levels for the zero-order scattering approximation as a function of scattering optical thickness $\tau_S^L$:

* $\Delta(j = 0) < 5\%$ for $\tau_S^L < 0.05$,

* $\Delta(j = 0) < 10\%$ for $0.05 < \tau_S^L < 0.3$,

* $\Delta(j = 0) < 30\%$ for $0.3 < \tau_S^L < 1$.

In conclusion, zero-order scattering is a relevant approximation that provides accurate results for optical thicknesses $\tau_S^L < 0.05$ and for TATA+ in general; it includes half the cases investigated in this paper. In addition, this approximation is exceedingly simple because it only requires multiplying $P_A^{(\Phi=0)}$, which is routinely calculated or measured in photochemistry practice, by $(1 - \Phi)$.

## 4.2 P1 approximations for luminescence radiation

For optical thicknesses $\tau_S^L > 0.05$, the absorption of luminescence radiation must be described to improve accuracy compared to zero-order approximation. Here, the luminescence radiation $P_A^{(S)}$ in Eq 25 is described by diffusion equations derived thanks to P1 approximation [35, 50]. But a simple analytical solution is only accessible for elastic scattering processes, and therefore P1 approximation is applied to the gray and elastic models presented in Section 3.3.

### 4.2.1 P1 approximation applied to the gray model for luminescence radiation.
The gray approximation presented in Section 3.3 produces fairly accurate results for TATA+ and Eosin Y (see $\Delta(\text{gray})$ in Table 2) and is also expected to be relevant in the case of photosensitizers with small spectral variations in extinction cross section over the luminescence spectral range. Therefore the P1 approximation is applied to the gray RTE in Eq 31. The trivial reference solution for ballistic photons (order zero) is retained, including incident spectral distribution, but wavelength redistribution due to luminescence is not taken into account: $P_A^{(S)} \approx \bar{P}_A^{(S)}$. Hereafter, a simple analytic expression for P1 approximation $P_A^{(S)} \approx \bar{P}_A^{(S,P1)}$ is obtained and we
study:

$$P_A \approx P_A^{(\text{P1})} = P_A^{(0)} + \bar{P}_A^{(\text{S,P1})} \quad (44)$$

where $P_A^{(0)}$ is the reference zero-order absorptance given in Eq 43 and $\bar{P}_A^{(\text{S,P1})}$ is the luminescence absorptance resulting from both gray and P1 approximations (see Eq 53).

First the expression of $I_\lambda^{(0)}$ in Eq 41 is substituted into the gray RTE Eq 31:

$$\mu \frac{d\bar{I}^{(S)}}{dx} = -\bar{k} \bar{I}^{(S)}(x, \mu) + \frac{\Phi \bar{k}}{2} \int_{-1}^{1} d\mu' \bar{I}^{(S)}(x, \mu') + \int_{0}^{+\infty} d\lambda' \frac{\Phi k_{\lambda'}}{2} q_{0,\lambda'} e^{-k_{\lambda'} x} \quad (45)$$

The first moment of this equation is integrated over $\mu \in [-1, 1]$ to obtain the local macroscopic balance:

$$\frac{d\bar{q}^{(S)}}{dx} = -\bar{k}(1-\Phi)\bar{G}^{(S)}(x) + \int_{0}^{+\infty} d\lambda' \Phi k_{\lambda'} q_{0,\lambda'} e^{-k_{\lambda'} x} \quad (46)$$

where $\bar{G}^{(S)}(x) = \int_{-1}^{1} d\mu \bar{I}^{(S)}(x, \mu)$ is the irradiance and $\bar{q}^{(S)}(x) = \int_{-1}^{1} d\mu \bar{I}^{(S)}(x, \mu)$ $\mu$ is the surface flux density at location $x$. Afterwards, Eq 45 is multiplied by $\mu$ and then integrated over $\mu \in [-1, 1]$:

$$\int_{-1}^{1} d\mu \mu^2 \frac{d\bar{I}^{(S)}}{dx} = -\bar{k} \bar{q}^{(S)}(x) \quad (47)$$

P1 approximation assumes linear angular dependence of intensity:

$$\bar{I}^{(S)}(x, \mu) = a(x) + b(x) \mu \quad (48)$$

This functional form of intensity is substituted into the left-hand term of Eq 47, and integration leads to Fick's law of diffusion [35, 50]:

$$\bar{q}^{(\text{S,P1})}(x) = -\frac{1}{3\bar{k}} \frac{d\bar{G}^{(\text{S,P1})}}{dx} \quad (49)$$

Finally, substituting Eq 49 in the local balance Eq 46 gives the diffusion equation:

$$\frac{d^2 \bar{G}^{(\text{S,P1})}}{dx^2} - 3 \bar{k}^2 (1 - \Phi) \bar{G}^{(\text{S,P1})}(x) = - \int_{0}^{+\infty} d\lambda \, 3 \bar{\Phi} \bar{k} \, k_{\lambda} q_{0,\lambda} e^{-k_{\lambda} x} \quad (50)$$

In this work, the Marshak boundary conditions detailed in [50] are used:

$$\bar{G}^{(\text{S,P1})}(0) - \frac{2}{3\bar{k}} \frac{d\bar{G}^{(\text{S,P1})}}{dx}(0) = 0 \quad (51)$$

$$\bar{G}^{(\text{S,P1})}(L) + \frac{2}{3\bar{k}} \frac{d\bar{G}^{(\text{S,P1})}}{dx}(L) = 0 \quad (52)$$

The solution to this equation is derived in S1 Appendix and used in the definition of luminescence absorptance (Eq 33). After symbolic integration, we obtain the expression of total absorptance in Eq 44 where:

$$\bar{P}_A^{(\text{S,P1})} = \int_0^{+\infty} d\lambda p^\prime(\lambda) \Phi(c_{0,\lambda}(1 - e^{-\tau_\lambda}) \\ + c_{1,\lambda}(1 - e^{-m \tau^L}) + c_{2,\lambda}(e^m \tau^L - 1)) \quad (52)$$

where $\tau_\lambda = k_\lambda L = E_\lambda CL$ and $\tau^L = \bar{k}L = \bar{E}^L CL$ are the spectral and mean extinction optical
study:

$$ \mathcal{P}_A \approx \mathcal{P}_A^{(\text{P1})} = \mathcal{P}_A^{(0)} + \bar{\mathcal{P}}_A^{(\text{S,P1})} \quad (44) $$

where $\mathcal{P}_A^{(0)}$ is the reference zero-order absorptance given in Eq 43 and $\bar{\mathcal{P}}_A^{(\text{S,P1})}$ is the luminescence absorptance resulting from both gray and P1 approximations (see Eq 53).

First the expression of $I_{\lambda}^{(0)}$ in Eq 41 is substituted into the gray RTE Eq 31:

$$ \mu \frac{d\bar{\bar{I}}^{(S)}}{dx} = -\bar{k} \bar{\bar{I}}^{(S)}(x, \mu) + \frac{\Phi \bar{k}}{2} \int_{-1}^{1} d\mu' \bar{\bar{I}}^{(S)}(x, \mu') + \int_{0}^{+\infty} d\lambda' \frac{\Phi k_{\lambda'}}{2} q_{0,\lambda'} e^{-k_{\lambda'} x} \quad (45) $$

The first moment of this equation is integrated over $\mu \in [-1, 1]$ to obtain the local macroscopic balance:

$$ \frac{d\bar{q}^{(S)}}{dx} = -\bar{k}(1-\Phi)\bar{G}^{(S)}(x) + \int_{0}^{+\infty} d\lambda' \Phi k_{\lambda'} q_{0,\lambda'} e^{-k_{\lambda'} x} \quad (46) $$

where $\bar{G}^{(S)}(x) = \int_{-1}^{1} d\mu \bar{\bar{I}}^{(S)}(x, \mu)$ is the irradiance and $\bar{q}^{(S)}(x) = \int_{-1}^{1} d\mu \bar{\bar{I}}^{(S)}(x, \mu)$ $\mu$ is the surface flux density at location $x$. Afterwards, Eq 45 is multiplied by $\mu$ and then integrated over $\mu \in [-1, 1]$:

$$ \int_{-1}^{1} d\mu \mu^2 \frac{d\bar{\bar{I}}^{(S)}}{dx} = -\bar{k} \bar{q}^{(S)}(x) \quad (47) $$

P1 approximation assumes linear angular dependence of intensity:

$$ \bar{\bar{I}}^{(S)}(x, \mu) = a(x) + b(x) \mu \quad (48) $$

This functional form of intensity is substituted into the left-hand term of Eq 47, and integration leads to Fick's law of diffusion [35, 50]:

$$ \bar{q}^{(\text{S,P1})}(x) = -\frac{1}{3\bar{k}} \frac{d\bar{G}^{(\text{S,P1})}}{dx} \quad (49) $$

Finally, substituting Eq 49 in the local balance Eq 46 gives the diffusion equation:

$$ \frac{d^2 \bar{G}^{(\text{S,P1})}}{dx^2} - 3 \bar{k}^2 (1 - \Phi) \bar{G}^{(\text{S,P1})}(x) = - \int_{0}^{+\infty} d\lambda \, 3 \bar{\Phi} \bar{k} \, k_{\lambda} q_{0,\lambda} e^{-k_{\lambda} x} \quad (50) $$

In this work, the Marshak boundary conditions detailed in [50] are used:

$$ \bar{G}^{(\text{S,P1})}(0) - \frac{2}{3\bar{k}} \frac{d\bar{G}^{(\text{S,P1})}}{dx}(0) = 0 \quad (51) $$

$$ \bar{G}^{(\text{S,P1})}(L) + \frac{2}{3\bar{k}} \frac{d\bar{G}^{(\text{S,P1})}}{dx}(L) = 0 \quad (52) $$

The solution to this equation is derived in S1 Appendix and used in the definition of luminescence absorptance (Eq 33). After symbolic integration, we obtain the expression of total absorptance in Eq 44 where:

$$ \bar{\mathcal{P}}_{A}^{(\text{S,P1})} = \int_{0}^{+\infty} d\lambda p^{i}(\lambda) \Phi(c_{0,\lambda}(1 - e^{-\tau_{\lambda}}) + c_{1,\lambda}(1 - e^{-m \tau^{L}}) + c_{2,\lambda}(e^{m \tau^{L}} - 1)) \quad (52) $$

where $\tau_{\lambda} = k_{\lambda} L = E_{\lambda} CL$ and $\tau^{L} = \bar{k}L = \bar{E}^{L}CL$ are the spectral and mean extinction optical
thicknesses, respectively (see Section 2.1 and Fig 2), $p^i(\lambda)$ is the incident spectrum (see Section 2.1 and Fig 3), $m = \sqrt{3(1-\Phi)}$, $b = \frac{1+2/3 m}{1-2/3 m}$, $d_\lambda = \tau_\lambda/(m\tau^L)$ and

$$ 
\begin{aligned}
c_{0,\lambda} &= \frac{1}{1-d_{\lambda}^2} \\
c_{1,\lambda} &= -c_{0,\lambda} d_{\lambda} \frac{e^{-(\tau_{\lambda}+m \tau^L)}(1-2/3 d_{\lambda} m) - b(1+2/3 d_{\lambda} m)}{e^{-2 m \tau^L}(1-2/3 m) - b(1+2/3 m)} \\
c_{2,\lambda} &= -c_{1,\lambda} b - c_{0,\lambda} d_{\lambda} \frac{1+2/3 d_{\lambda} m}{1-2/3 m}
\end{aligned}
$$

The results presented in Table 2 indicate that this approximation does not improve accuracy either significantly or systematically compared to zero-order approximation (see $\Delta(\overline{P1})$). For TATA+, which is accurately described by zero-order approximation and the gray model, relative errors are comparable: $|\Delta(\overline{P1})| \simeq |\Delta(j=0)| \simeq |\Delta(\text{gray})|$. For Eosin Y, accuracy is slightly improved, but the relative error exceeds 10% as soon as $\tau_S^L > 0.6$. For Rhodamine B, the relative error is even greater compared to zero-order expansion only.

The key point here is that the $\overline{P1}$ approximation gives results that approach those provided by the gray equivalent model presented in Section 3.3 when scattering optical thickness increases ($\Delta(\overline{P1}) \to \Delta(\text{gray})$ as $\tau_S^L \to \infty$). This is due to the fact that the angular distribution of intensity comes closer to the functional form assumed by the P1 approximation in Eq 48, due to isotropic scattering. As a result, multiple scattering is well described as soon as $\tau_S^L \ge 0.5$ and $\Delta(\overline{P1}) \simeq \Delta(\text{gray})$ in this case. But on the other hand, the higher the $\tau_S^L$, the greater the inelastic scattering effects, which are neglected in the gray model at the root of the $\overline{P1}$ approximation (see Section 3.3). Overall, the $\overline{P1}$ approximation never finds its niche in the present case study. However, we kept this approximation because it would provide accurate results for photosensitizers with small spectral variations over the luminescence spectral range, since the gray model would be accurate in that case (see Section 3.3). We cannot demonstrate this statement with the typical photosensitizers that we selected because extinction cross section has large spectral variations over the luminescence spectral range. However, we can demonstrate it in the following virtual case: we use the same properties as Rhodamine B in Fig 2, except that the extinction cross section is gray and equal to $\overline{E}^L = 24943 \text{ m}^2.\text{mol}^{-1}$ (see Table 1). Comparing the Monte Carlo reference solution and the $\overline{P1}$ solution for $\tau_A^i = 4.6$ we obtain $\Delta(\overline{P1}) = 0.01$, indicating a very good accuracy.

**4.2.2 P1 approximation applied to an elastic-scattering model for luminescence radiation.** In order to improve the description of inelastic scattering effects that prevented the previous $\overline{P1}$ approximation from providing significant results, here we apply the P1 approximation to the elastic-scattering model presented in Section 3.3:

$$ \mathcal{P}_A \approx \mathcal{P}_A^{(\overline{P1})} = \mathcal{P}_A^{(0)} + \hat{\mathcal{P}}_A^{(S,P1)} \quad (53) $$

where $\mathcal{P}_A^{(0)}$ is the reference zero-order absorptance given in Eq 43 and $\hat{\mathcal{P}}_A^{(S,P1)}$ is the luminescence absorptance resulting from both elastic and P1 approximations (see Eq 55).

An analytical expression of $\hat{\mathcal{P}}_A^{(S,P1)}$ is obtained by following exactly the same steps as in Eqs 45 to 51 and S1 Appendix but working with the elastic RTE Eq 35 instead of the gray RTE
thicknesses, respectively (see Section 2.1 and Fig 2), $p^i(\lambda)$ is the incident spectrum (see Section 2.1 and Fig 3), $m = \sqrt{3(1-\Phi)}$, $b = \frac{1+2/3 m}{1-2/3 m}$, $d_\lambda = \tau_\lambda/(m\tau^L)$ and

$$ 
\begin{aligned}
c_{0,\lambda} &= \frac{1}{1-d_{\lambda}^2} \\
c_{1,\lambda} &= -c_{0,\lambda} d_{\lambda} \frac{e^{-(\tau_{\lambda}+m \tau^L)}(1-2/3 d_{\lambda} m) - b(1+2/3 d_{\lambda} m)}{e^{-2 m \tau^L}(1-2/3 m) - b(1+2/3 m)} \\
c_{2,\lambda} &= -c_{1,\lambda} b - c_{0,\lambda} d_{\lambda} \frac{1+2/3 d_{\lambda} m}{1-2/3 m}
\end{aligned}
$$

The results presented in Table 2 indicate that this approximation does not improve accuracy either significantly or systematically compared to zero-order approximation (see $\Delta(\overline{P1})$). For TATA+, which is accurately described by zero-order approximation and the gray model, relative errors are comparable: $|\Delta(\overline{P1})| \simeq |\Delta(j=0)| \simeq |\Delta(\text{gray})|$. For Eosin Y, accuracy is slightly improved, but the relative error exceeds 10% as soon as $\tau_S^L > 0.6$. For Rhodamine B, the relative error is even greater compared to zero-order expansion only.

The key point here is that the $\overline{P1}$ approximation gives results that approach those provided by the gray equivalent model presented in Section 3.3 when scattering optical thickness increases ($\Delta(\overline{P1}) \to \Delta(\text{gray})$ as $\tau_S^L \to \infty$). This is due to the fact that the angular distribution of intensity comes closer to the functional form assumed by the P1 approximation in Eq 48, due to isotropic scattering. As a result, multiple scattering is well described as soon as $\tau_S^L \ge 0.5$ and $\Delta(\overline{P1}) \simeq \Delta(\text{gray})$ in this case. But on the other hand, the higher the $\tau_S^L$, the greater the inelastic scattering effects, which are neglected in the gray model at the root of the $\overline{P1}$ approximation (see Section 3.3). Overall, the $\overline{P1}$ approximation never finds its niche in the present case study. However, we kept this approximation because it would provide accurate results for photosensitizers with small spectral variations over the luminescence spectral range, since the gray model would be accurate in that case (see Section 3.3). We cannot demonstrate this statement with the typical photosensitizers that we selected because extinction cross section has large spectral variations over the luminescence spectral range. However, we can demonstrate it in the following virtual case: we use the same properties as Rhodamine B in Fig 2, except that the extinction cross section is gray and equal to $\overline{E}^L = 24943 \text{ m}^2.\text{mol}^{-1}$ (see Table 1). Comparing the Monte Carlo reference solution and the $\overline{P1}$ solution for $\tau_A^i = 4.6$ we obtain $\Delta(\overline{P1}) = 0.01$, indicating a very good accuracy.

### 4.2.2 P1 approximation applied to an elastic-scattering model for luminescence radiation.
In order to improve the description of inelastic scattering effects that prevented the previous $\overline{P1}$ approximation from providing significant results, here we apply the P1 approximation to the elastic-scattering model presented in Section 3.3:

$$ \mathcal{P}_A \approx \mathcal{P}_A^{(\overline{P1})} = \mathcal{P}_A^{(0)} + \hat{\mathcal{P}}_A^{(S,P1)} \quad (53) $$

where $\mathcal{P}_A^{(0)}$ is the reference zero-order absorptance given in Eq 43 and $\hat{\mathcal{P}}_A^{(S,P1)}$ is the luminescence absorptance resulting from both elastic and P1 approximations (see Eq 55).

An analytical expression of $\hat{\mathcal{P}}_A^{(S,P1)}$ is obtained by following exactly the same steps as in Eqs 45 to 51 and S1 Appendix but working with the elastic RTE Eq 35 instead of the gray RTE
Eq 31, and using the definition of absorptance given in Eq 24 instead of its gray version Eq 33:

$$ \hat{P}_A^{(S,P1)} = \int_0^{+\infty} d\lambda p^i(\lambda) \int_0^{+\infty} d\lambda' p^L(\lambda') \Phi(\hat{c}_{0,\lambda,\lambda'}(1 - e^{-\tau_\lambda}) + \hat{c}_{1,\lambda,\lambda'}(1 - e^{-m\tau_\lambda'}) + \hat{c}_{2,\lambda,\lambda'}(e^{m\tau_\lambda'} - 1)) \quad (54) $$

where $\tau_\lambda = k_\lambda L = E_\lambda CL$ is the extinction optical thickness (see Section 2.1 and Fig 2), $p^i(\lambda)$ and $p^L(\lambda)$ are the incident and luminescence spectrum respectively (see Section 2.1 and Fig 3), $m = \sqrt{3(1 - \Phi)}$, $b = \frac{1+2/3 m}{1-2/3 m}$, $d_{\lambda,\lambda'} = \tau_\lambda/(m\tau_\lambda)$ and

$$ \hat{c}_{0,\lambda,\lambda'} = \frac{1}{1 - d_{\lambda,\lambda'}^2} $$

$$ \hat{c}_{1,\lambda,\lambda'} = -\hat{c}_{0,\lambda,\lambda'} d_{\lambda,\lambda'} \frac{e^{-(\tau_\lambda + m \tau_{\lambda'})} (1 - 2/3 d_{\lambda,\lambda'} m) - b (1 + 2/3 d_{\lambda,\lambda'} m)}{e^{-2 m \tau_{\lambda'}} (1 - 2/3 m) - b (1 + 2/3 m)} $$

$$ \hat{c}_{2,\lambda,\lambda'} = -\hat{c}_{1,\lambda,\lambda'} b - \hat{c}_{0,\lambda,\lambda'} d_{\lambda,\lambda'} \frac{1 + 2/3 d_{\lambda,\lambda'} m}{1 - 2/3 m} $$

The relative error $\Delta(\hat{P1})$ obtained with this model is reported in Table 2. As expected from previous discussions (including Section 3.3), accuracy is significantly improved compared to $\hat{P1}$ approximation, in particular for optical thicknesses $\tau^L > 0.1$ involving the strong impact of inelastic scattering. Double spectral integration in Eq 55 is the price to pay for this improvement. As a result, $\hat{P1}$ approximation is highly accurate for low-to-intermediate optical thicknesses, generating errors which are always below 2% for TATA+ and Eosin Y. Indeed, the two validity conditions for $\hat{P1}$ approximation are satisfied in these situations: first, the elastic-scattering model presented Section 3.3 is accurate ($\Delta(\text{elastic}) < 2\%$) and second, P1 approximation well describes multiple scattering in the framework of this elastic model ($\Delta(\hat{P1}) \approx \Delta(\text{elastic})$, see the discussion on P1 approximation validity in paragraphs below Eq 53). For Rhodamine B, on the other hand, the error increases by up to 26% because, even if the second condition is satisfied ($\Delta(\hat{P1}) \approx \Delta(\text{elastic})$ for $\tau^L > 1$), the elastic-scattering model at the root of the $\hat{P1}$ approximation is not relevant ($\Delta(\text{elastic}) > 10\%$). This is due to high scattering optical thicknesses combined with significant inelastic effects, as discussed in Section 3.3.

## 4.3 Single-scattering approximation

The $\hat{P1}$ approximation developed in the previous section fails to describe absorptance for Rhodamine B photosensitizer, with $\Delta(\hat{P1})$ errors higher than the contribution of scattering orders $j > 1$ (see Table 3). Therefore, the description of luminescence radiation thanks to first-order scattering expansion is investigated here.

### 4.3.1 Single-scattering approximation applied to the reference model.
In the previous section, we were able to accurately account for multiple scattering, but important effects due to inelastic scattering were missing. Here we aim to describe inelastic scattering more accurately, but to do so we have to make a concession in the description of multiple scattering. In Section 3.2 we noticed that most of the luminescence contribution is due to photons absorbed after one scattering event only. Therefore, the expansion in successive orders of scattering is here truncated to the first order, which leads to a simple analytic expression [35, 50]:

$$ P_A \approx P_A^{(j \le 1)} = P_A^{(0)} + P_A^{(1)} \quad (55) $$

where $P_A^{(0)}$ is the trivial zero-order absorptance given in Eq 43 and $P_A^{(1)}$ is given in Eq 58.
Eq 31, and using the definition of absorptance given in Eq 24 instead of its gray version Eq 33:

$$ \hat{P}_A^{(S,P1)} = \int_0^{+\infty} d\lambda p^i(\lambda) \int_0^{+\infty} d\lambda' p^L(\lambda') \Phi(\hat{c}_{0,\lambda,\lambda'}(1 - e^{-\tau_\lambda}) + \hat{c}_{1,\lambda,\lambda'}(1 - e^{-m\tau_\lambda'}) + \hat{c}_{2,\lambda,\lambda'}(e^{m\tau_\lambda'} - 1)) \quad (54) $$

where $\tau_\lambda = k_\lambda L = E_\lambda CL$ is the extinction optical thickness (see Section 2.1 and Fig 2), $p^i(\lambda)$ and $p^L(\lambda)$ are the incident and luminescence spectrum respectively (see Section 2.1 and Fig 3), $m = \sqrt{3(1 - \Phi)}$, $b = \frac{1+2/3 m}{1-2/3 m}$, $d_{\lambda,\lambda'} = \tau_\lambda/(m\tau_\lambda)$ and

$$ \hat{c}_{0,\lambda,\lambda'} = \frac{1}{1 - d_{\lambda,\lambda'}^2} $$

$$ \hat{c}_{1,\lambda,\lambda'} = -\hat{c}_{0,\lambda,\lambda'} d_{\lambda,\lambda'} \frac{e^{-(\tau_\lambda + m \tau_{\lambda'})} (1 - 2/3 d_{\lambda,\lambda'} m) - b (1 + 2/3 d_{\lambda,\lambda'} m)}{e^{-2m \tau_{\lambda'}} (1 - 2/3 m) - b (1 + 2/3 m)} $$

$$ \hat{c}_{2,\lambda,\lambda'} = -\hat{c}_{1,\lambda,\lambda'} b - \hat{c}_{0,\lambda,\lambda'} d_{\lambda,\lambda'} \frac{1 + 2/3 d_{\lambda,\lambda'} m}{1 - 2/3 m} $$

The relative error $\Delta(\hat{P}1)$ obtained with this model is reported in Table 2. As expected from previous discussions (including Section 3.3), accuracy is significantly improved compared to $\hat{P}1$ approximation, in particular for optical thicknesses $\tau^L > 0.1$ involving the strong impact of inelastic scattering. Double spectral integration in Eq 55 is the price to pay for this improvement. As a result, $\hat{P}1$ approximation is highly accurate for low-to-intermediate optical thicknesses, generating errors which are always below 2% for TATA+ and Eosin Y. Indeed, the two validity conditions for $\hat{P}1$ approximation are satisfied in these situations: first, the elastic-scattering model presented Section 3.3 is accurate ($\Delta(\text{elastic}) < 2\%$) and second, P1 approximation well describes multiple scattering in the framework of this elastic model ($\Delta(\hat{P}1) \approx \Delta(\text{elastic})$, see the discussion on P1 approximation validity in paragraphs below Eq 53). For Rhodamine B, on the other hand, the error increases by up to 26% because, even if the second condition is satisfied ($\Delta(\hat{P}1) \approx \Delta(\text{elastic})$ for $\tau^L > 1$), the elastic-scattering model at the root of the $\hat{P}1$ approximation is not relevant ($\Delta(\text{elastic}) > 10\%$). This is due to high scattering optical thicknesses combined with significant inelastic effects, as discussed in Section 3.3.

## 4.3 Single-scattering approximation

The $\hat{P}1$ approximation developed in the previous section fails to describe absorptance for Rhodamine B photosensitizer, with $\Delta(\hat{P}1)$ errors higher than the contribution of scattering orders $j > 1$ (see Table 3). Therefore, the description of luminescence radiation thanks to first-order scattering expansion is investigated here.

### 4.3.1 Single-scattering approximation applied to the reference model.
In the previous section, we were able to accurately account for multiple scattering, but important effects due to inelastic scattering were missing. Here we aim to describe inelastic scattering more accurately, but to do so we have to make a concession in the description of multiple scattering. In Section 3.2 we noticed that most of the luminescence contribution is due to photons absorbed after one scattering event only. Therefore, the expansion in successive orders of scattering is here truncated to the first order, which leads to a simple analytic expression [35, 50]:

$$ \mathcal{P}_A \approx \mathcal{P}_A^{(j \le 1)} = \mathcal{P}_A^{(0)} + \mathcal{P}_A^{(1)} \quad (55) $$

where $\mathcal{P}_A^{(0)}$ is the trivial zero-order absorptance given in Eq 43 and $\mathcal{P}_A^{(1)}$ is given in Eq 58.
The expression of $I_{\lambda}^{(0)}$ given in Eq 41 is substituted in the RTE Eq 15 for first-order intensity $I_{\lambda}^{(1)}$, leading to the ordinary linear differential equation

$$ \mu \frac{dI_{\lambda}^{(1)}(x, \mu)}{dx} = -k_{\lambda} I_{\lambda}^{(1)}(x, \mu) + \int_{0}^{+\infty} d\lambda' p^{L}(\lambda) \frac{\Phi k_{\lambda'}}{2} q_{0,\lambda'} e^{-k_{\lambda'} x} \quad (56) $$

whose solution is

$$ I_{\lambda}^{(1)}(x, \mu) = \begin{cases} \int_{0}^{x} \frac{dx'}{\mu} e^{-\frac{k_{\lambda}(x-x')}{\mu}} p^{L}(\lambda) \int_{0}^{+\infty} d\lambda' \frac{\Phi k_{\lambda'}}{2} q_{0,\lambda'} e^{-k_{\lambda'} x'} & \text{for } \mu > 0 \\ \int_{x}^{L} \frac{dx'}{-\mu} e^{-\frac{k_{\lambda}(x-x')}{\mu}} p^{L}(\lambda) \int_{0}^{+\infty} d\lambda' \frac{\Phi k_{\lambda'}}{2} q_{0,\lambda'} e^{-k_{\lambda'} x'} & \text{for } \mu < 0 \end{cases} \quad (57) $$

Substituting this expression in Eq 24 and performing symbolic integration leads to (see details in S2 Appendix):

$$ \mathcal{P}_A^{(1)} = \int_0^{+\infty} d\lambda \, p^i(\lambda) \int_0^{+\infty} d\lambda' \, p^L(\lambda') \frac{\tau_{\lambda'} \Phi(1-\Phi)}{2} \left\{ \frac{1}{\tau_{\lambda}} [Ei(-\tau_{\lambda'}) (1 + e^{-\tau_{\lambda}}) \right. \\ - Ei(-\tau_{\lambda'} - \tau_{\lambda}) + \frac{1}{2} \ln \left( \frac{\tau_{\lambda} + \tau_{\lambda'}}{\tau_{\lambda'}} \right)^2 + e^{-\tau_{\lambda}} \left( \frac{1}{2} \ln \left( \frac{\tau_{\lambda} - \tau_{\lambda'}}{\tau_{\lambda'}} \right)^2 - Ei(\tau_{\lambda} - \tau_{\lambda'}) \right) \left. \right\} \quad (58) \\ + (e^{-\tau_{\lambda}} - 1) \left( Ei(-\tau_{\lambda'}) + \frac{1}{\tau_{\lambda'}} (e^{-\tau_{\lambda'}} - 1) \right) \right\} $$

where $Ei = \int_{-\infty}^{x} \frac{e^t}{t} dt$ is the exponential integral function (which is available in most scientific computing tools and applications).

Note that applying the single-scattering approximation to the reference model or to the elastic-scattering model presented in Section 3.3 leads to the exact same expression. Indeed, these two models lead to identical zero- and first-order terms in the scattering expansion (higher-order terms are different) because our elastic-scattering model includes spectral redistribution when ballistic photons move from population (0) to population (1). Thus the numerical evaluation of Eq 58 requires double integration, just as for the $\hat{P}1$ approximation in Eq 55.

The relative error $\Delta(j \le 1)$ recorded with this model is reported in Table 2. Note that $\Delta(j \le 1) = -P^{(j>1)}$ in Table 3 and that first-order scattering expansion was already discussed in Section 3.2.

Single-scattering approximation leads to an error below or equal to 5% for TATA⁺ and Eosin Y, for every tested optical thickness. It is slightly less accurate than $\hat{P}1$ approximation, but still a valuable approximation. At the contrary, as expected, accuracy for Rhodamine B is better with single-scattering approximation than with $\hat{P}1$ approximation. Nevertheless, Rhodamine B is definitely challenging: the relative error $\Delta(j \le 1)$ rapidly exceeds 10% when absorption optical thickness increases. Overall, single-scattering approximation is fairly accurate when $\tau^L < 1$, which corresponds to a significant improvement compared to zero-order scattering approximation.

**4.3.2 Single-scattering approximation applied to the gray model for luminescence radiation.** The above approximation provides valuable results but requires double numerical integration. Here, the gray approximation presented in Section 3.3 is applied to the first-order term, in order to eliminate spectral integration over the luminescence spectrum in $\mathcal{P}_A^{(1)}$ (see Eq 58): $\mathcal{P}_A^{(1)} \approx \bar{\mathcal{P}}_A^{(1)}$. In this way we obtain an expression for the single-scattering approximation
The expression of $I_{\lambda}^{(0)}$ given in Eq 41 is substituted in the RTE Eq 15 for first-order intensity $I_{\lambda}^{(1)}$, leading to the ordinary linear differential equation

$$ \mu \frac{dI_{\lambda}^{(1)}(x, \mu)}{dx} = -k_{\lambda} I_{\lambda}^{(1)}(x, \mu) + \int_{0}^{+\infty} d\lambda' p^{L}(\lambda) \frac{\Phi k_{\lambda'}}{2} q_{0,\lambda'} e^{-k_{\lambda'} x} \quad (56) $$

whose solution is

$$ I_{\lambda}^{(1)}(x, \mu) = \begin{cases} \int_{0}^{x} \frac{dx'}{\mu} e^{-\frac{k_{\lambda}(x-x')}{\mu}} p^{L}(\lambda) \int_{0}^{+\infty} d\lambda' \frac{\Phi k_{\lambda'}}{2} q_{0,\lambda'} e^{-k_{\lambda'} x'} & \text{for } \mu > 0 \\ \int_{x}^{L} \frac{dx'}{-\mu} e^{-\frac{k_{\lambda}(x-x')}{\mu}} p^{L}(\lambda) \int_{0}^{+\infty} d\lambda' \frac{\Phi k_{\lambda'}}{2} q_{0,\lambda'} e^{-k_{\lambda'} x'} & \text{for } \mu < 0 \end{cases} \quad (57) $$

Substituting this expression in Eq 24 and performing symbolic integration leads to (see details in S2 Appendix):

$$ \mathcal{P}_A^{(1)} = \int_0^{+\infty} d\lambda \, p^i(\lambda) \int_0^{+\infty} d\lambda' \, p^L(\lambda') \frac{\tau_{\lambda'} \Phi(1-\Phi)}{2} \left\{ \begin{aligned} & \frac{1}{\tau_{\lambda}} [Ei(-\tau_{\lambda'}) (1 + e^{-\tau_{\lambda}}) \\ & - Ei(-\tau_{\lambda'} - \tau_{\lambda}) + \frac{1}{2} \ln \left( \frac{\tau_{\lambda} + \tau_{\lambda'}}{\tau_{\lambda'}} \right)^2 + e^{-\tau_{\lambda}} \left( \frac{1}{2} \ln \left( \frac{\tau_{\lambda} - \tau_{\lambda'}}{\tau_{\lambda'}} \right)^2 - Ei(\tau_{\lambda} - \tau_{\lambda'}) \right) ] \\ & + (e^{-\tau_{\lambda}} - 1) \left( Ei(-\tau_{\lambda'}) + \frac{1}{\tau_{\lambda'}} (e^{-\tau_{\lambda'}} - 1) \right) \end{aligned} \right\} \quad (58) $$

where $Ei = \int_{-\infty}^{x} \frac{e^t}{t} dt$ is the exponential integral function (which is available in most scientific computing tools and applications).

Note that applying the single-scattering approximation to the reference model or to the elastic-scattering model presented in Section 3.3 leads to the exact same expression. Indeed, these two models lead to identical zero- and first-order terms in the scattering expansion (higher-order terms are different) because our elastic-scattering model includes spectral redistribution when ballistic photons move from population (0) to population (1). Thus the numerical evaluation of Eq 58 requires double integration, just as for the $\hat{P}1$ approximation in Eq 55.

The relative error $\Delta(j \le 1)$ recorded with this model is reported in Table 2. Note that $\Delta(j \le 1) = -P^{(j>1)}$ in Table 3 and that first-order scattering expansion was already discussed in Section 3.2.

Single-scattering approximation leads to an error below or equal to 5% for TATA⁺ and Eosin Y, for every tested optical thickness. It is slightly less accurate than $\hat{P}1$ approximation, but still a valuable approximation. At the contrary, as expected, accuracy for Rhodamine B is better with single-scattering approximation than with $\hat{P}1$ approximation. Nevertheless, Rhodamine B is definitely challenging: the relative error $\Delta(j \le 1)$ rapidly exceeds 10% when absorption optical thickness increases. Overall, single-scattering approximation is fairly accurate when $\tau^L < 1$, which corresponds to a significant improvement compared to zero-order scattering approximation.

### 4.3.2 Single-scattering approximation applied to the gray model for luminescence radiation.
The above approximation provides valuable results but requires double numerical integration. Here, the gray approximation presented in Section 3.3 is applied to the first-order term, in order to eliminate spectral integration over the luminescence spectrum in $\mathcal{P}_A^{(1)}$ (see Eq 58): $\mathcal{P}_A^{(1)} \approx \bar{\mathcal{P}}_A^{(1)}$. In this way we obtain an expression for the single-scattering approximation
$\mathcal{P}_A^{(\overline{SS})}$ that is even more straightforward to evaluate:

$$ \mathcal{P}_A \approx \mathcal{P}_A^{(\overline{SS})} = \mathcal{P}_A^{(0)} + \bar{\mathcal{P}}_A^{(1)} \qquad (59) $$

where $\mathcal{P}_A^{(0)}$ is given in Eq 43 and $\bar{\mathcal{P}}_A^{(1)}$ is given in Eq 62.

The same gray equivalent medium representation as in Section 3.3 is applied to Eq 15 with $j = 1$, leading to the following RTE for first-order intensity $\bar{I}_\lambda^{(1)} \approx \int_0^{+\infty} d\lambda I_\lambda^{(1)}$:

$$ \mu \frac{d\bar{I}^{(1)}}{dx} = -\bar{k}\bar{I}^{(1)}(x, \mu) + \int_{0}^{+\infty} d\lambda' \frac{\Phi k_{\lambda'}}{2} \int_{-1}^{1} d\mu' I_{\lambda'}^{(0)}(x, \mu') \qquad (60) $$

with the boundary condition

$$ \bar{I}_{\lambda}^{(1)}(x, \mu) = 0 \qquad (61) $$

Then the same approach as in Eqs 56 to 58 and S2 Appendix leads to

$$ \bar{\mathcal{P}}_A^{(1)} = \int_0^{+\infty} d\lambda p^i(\lambda) \frac{\tau^L \Phi(1-\Phi)}{2} \left\{ \begin{aligned} & \frac{1}{\tau_\lambda} [Ei(-\tau^L)(1+e^{-\tau_\lambda}) - Ei(-\tau^L - \tau_\lambda)] \\ & + \frac{1}{2} \ln \left( \frac{\tau_\lambda + \tau^L}{\tau^L} \right)^2 + e^{-\tau_\lambda} \left( \frac{1}{2} \ln \left( \frac{\tau_\lambda - \tau^L}{\tau^L} \right)^2 - Ei(\tau_\lambda - \tau^L) \right) \end{aligned} \right\} \qquad (62) \\ + (e^{-\tau_\lambda} - 1) \left( Ei(-\tau^L) + \frac{1}{\tau^L} (e^{-\tau^L} - 1) \right) \} $$

where $\tau^L = \bar{k}L = \bar{E}^LCL$ is the mean extinction optical thickness (see Section 2.1), and other definitions are given below Eq 58.

The results in Table 2 present the relative error $\Delta(\overline{SS})$ obtained with this model. Surprisingly, it provides better accuracy than the previous model, even though inelastic spectral effects are here neglected: $|\Delta(\overline{SS})| < |\Delta(j \le 1)|$. The results are remarkably accurate in every tested configuration, with a relative error below 5% (except Rhodamine B for $\tau_S^L = 0.952$ where $\Delta(\overline{SS}) = 6.2\%$).

This remarkable accuracy results in fact from the compensation of two errors:

* On the one hand, truncating scattering expansion always underestimates absorptance, as discussed in Section 2.3. The magnitude of this underestimation is provided by $\Delta(j \le 1)$ in Table 2.

* On the other hand, gray approximation for luminescence radiation always overestimates absorptance (concavity property of the integrand), as discussed in Section 3.3.

Therefore, even if this approximation is of great interest for the study of photosensitized and photocatalytic systems, we cannot be sure that these two antagonistic effects always have similar absolute values and compensate for each other. For example, if a photosensitizer were studied with a gray extinction cross section on the luminescence spectral range (see Section 3.3), then $\Delta(\overline{SS}) = \Delta(j \le 1)$ and accuracy would be significantly degraded at high optical thicknesses. Such situations should rather be addressed using the $\bar{\mathcal{P}}\bar{\mathcal{I}}$ approximation, which would be accurate when $\tau_S^L > 1$ (see Section 4.2).

To conclude, the use of this gray single-scattering approximation requires some precautions, but provides accurate results for every typical case investigated in this paper.
$\mathcal{P}_A^{(\overline{SS})}$ that is even more straightforward to evaluate:

$$ \mathcal{P}_A \approx \mathcal{P}_A^{(\overline{SS})} = \mathcal{P}_A^{(0)} + \bar{\mathcal{P}}_A^{(1)} \qquad (59) $$

where $\mathcal{P}_A^{(0)}$ is given in Eq 43 and $\bar{\mathcal{P}}_A^{(1)}$ is given in Eq 62.

The same gray equivalent medium representation as in Section 3.3 is applied to Eq 15 with $j=1$, leading to the following RTE for first-order intensity $\bar{I}_\lambda^{(1)} \approx \int_0^{+\infty} d\lambda I_\lambda^{(1)}$:

$$ \mu \frac{d\bar{I}^{(1)}}{dx} = -\bar{k}\bar{I}^{(1)}(x, \mu) + \int_{0}^{+\infty} d\lambda' \frac{\Phi k_{\lambda'}}{2} \int_{-1}^{1} d\mu' I_{\lambda'}^{(0)}(x, \mu') \qquad (60) $$

with the boundary condition

$$ \bar{I}_{\lambda}^{(1)}(x, \mu) = 0 \qquad (61) $$

Then the same approach as in Eqs 56 to 58 and S2 Appendix leads to

$$ \bar{\mathcal{P}}_A^{(1)} = \int_0^{+\infty} d\lambda p^i(\lambda) \frac{\tau^L \Phi(1-\Phi)}{2} \left\{ \begin{aligned} & \frac{1}{\tau_\lambda} [Ei(-\tau^L)(1+e^{-\tau_\lambda}) - Ei(-\tau^L - \tau_\lambda)] \\ & + \frac{1}{2} \ln \left( \frac{\tau_\lambda + \tau^L}{\tau^L} \right)^2 + e^{-\tau_\lambda} \left( \frac{1}{2} \ln \left( \frac{\tau_\lambda - \tau^L}{\tau^L} \right)^2 - Ei(\tau_\lambda - \tau^L) \right) \end{aligned} \right\} \qquad (62) \\ + (e^{-\tau_\lambda} - 1) \left( Ei(-\tau^L) + \frac{1}{\tau^L} (e^{-\tau^L} - 1) \right) \right\} $$

where $\tau^L = \bar{k}L = \bar{E}^LCL$ is the mean extinction optical thickness (see Section 2.1), and other definitions are given below Eq 58.

The results in Table 2 present the relative error $\Delta(\overline{SS})$ obtained with this model. Surprisingly, it provides better accuracy than the previous model, even though inelastic spectral effects are here neglected: $|\Delta(\overline{SS})| < |\Delta(j \le 1)|$. The results are remarkably accurate in every tested configuration, with a relative error below 5% (except Rhodamine B for $\tau_S^L = 0.952$ where $\Delta(\overline{SS}) = 6.2\%$).

This remarkable accuracy results in fact from the compensation of two errors:

* On the one hand, truncating scattering expansion always underestimates absorptance, as discussed in Section 2.3. The magnitude of this underestimation is provided by $\Delta(j \le 1)$ in Table 2.

* On the other hand, gray approximation for luminescence radiation always overestimates absorptance (concavity property of the integrand), as discussed in Section 3.3.

Therefore, even if this approximation is of great interest for the study of photosensitized and photocatalytic systems, we cannot be sure that these two antagonistic effects always have similar absolute values and compensate for each other. For example, if a photosensitizer were studied with a gray extinction cross section on the luminescence spectral range (see Section 3.3), then $\Delta(\overline{SS}) = \Delta(j \le 1)$ and accuracy would be significantly degraded at high optical thicknesses. Such situations should rather be addressed using the $\bar{\mathcal{P}}\bar{\mathcal{I}}$ approximation, which would be accurate when $\tau_S^L > 1$ (see Section 4.2).

To conclude, the use of this gray single-scattering approximation requires some precautions, but provides accurate results for every typical case investigated in this paper.
**5 Practical implementation for overall quantum yield estimation**

In a previous publication [18], the overall quantum yield of bio-inspired H₂ production imple-
mented in a benchmark photoreactor was estimated without taking into account lumines-
cence. Hereafter we revisit these results and discuss the effects of luminescence on quantum
yield measurements. In particular, the method established in the previous sections will be
implemented in order to identify a relevant analytical approximate solution based on optical
thickness analysis.

The photocatalytic system is an aqueous solution of an diiron-thiolate (FeFe) complex used as a proton reduction catalyst, fluorescent Eosin Y (EY²⁻) as the photosensitizer and triethylamine (Et₃N) as a sacrificial electron donor. Such a photoreactive system can operate in water with the mechanism (see Fig 4) proposed as follows [21]:

1. When EY²⁻ absorbs photon energy, electrons are excited to a higher electronic level before rapidly relaxing to the first electronic state (the singlet ¹*EY*²⁻),

2. In the most favorable case for hydrogen production, electrons go through a spin conversion by intersystem crossing (ISC) into a triplet state (3*EY*²⁻),

3. Next, an electron is transferred to the diiron catalyst,

4. To close the Eosin Y cycle, EY⁻ produced by the oxydoreduction reaction between 3*EY*²⁻ and catalyst, EY²⁻ is regenerated by an electron transfer from Et₃N,

5. Finally, a succession of oxidation-reduction and protonation with the catalyst, extensively described and analyzed, arises to produce hydrogen [51, 52].

<figure><img src="image_8.png" /><figcaption>**Fig 4. Proposed mechanism.** Hydrogen production with the fluorescent photocatalytic system in homogeneous phase [21, 51, 52].</figcaption></figure>

https://doi.org/10.1371/journal.pone.0255002.g004
**5 Practical implementation for overall quantum yield estimation**

In a previous publication [18], the overall quantum yield of bio-inspired H₂ production imple-
mented in a benchmark photoreactor was estimated without taking into account lumines-
cence. Hereafter we revisit these results and discuss the effects of luminescence on quantum
yield measurements. In particular, the method established in the previous sections will be
implemented in order to identify a relevant analytical approximate solution based on optical
thickness analysis.

The photocatalytic system is an aqueous solution of an diiron-thiolate (FeFe) complex used as a proton reduction catalyst, fluorescent Eosin Y (EY²⁻) as the photosensitizer and triethylamine (Et₃N) as a sacrificial electron donor. Such a photoreactive system can operate in water with the mechanism (see Fig 4) proposed as follows [21]:

1. When EY²⁻ absorbs photon energy, electrons are excited to a higher electronic level before rapidly relaxing to the first electronic state (the singlet ¹*EY*²⁻),

2. In the most favorable case for hydrogen production, electrons go through a spin conversion by intersystem crossing (ISC) into a triplet state (3*EY*²⁻),

3. Next, an electron is transferred to the diiron catalyst,

4. To close the Eosin Y cycle, EY⁻ produced by the oxydoreduction reaction between 3*EY*²⁻ and catalyst, EY²⁻ is regenerated by an electron transfer from Et₃N,

5. Finally, a succession of oxidation-reduction and protonation with the catalyst, extensively described and analyzed, arises to produce hydrogen [51, 52].

<figure><img src="image_7.png" /><figcaption>**Fig 4. Proposed mechanism.** Hydrogen production with the fluorescent photocatalytic system in homogeneous phase [21, 51, 52].</figcaption></figure>

https://doi.org/10.1371/journal.pone.0255002.g004
<figure><img src="image_9.png" /><figcaption>**Fig 5.** **Materials.** (a) Photocatalytic system radiative properties—left: Eosin Y (24 µM) UV-visible molar extinction cross-section *E*<sub>λ</sub>, fluorescence emission spectrum *p*<sup>1</sup>(*λ*) and quantum yield Φ = 0.35 in a mixture of triethylamine (10% vol.) in water at pH 10.5 [18, 21]; right: Catalyst (1.66 mM) molar absorption cross section in methanol. UV-visible spectra were obtained using a Shimadzu UV-visible UV-160 A spectrophotometer and emission spectra from a Shimadzu RF-1501 spectrofluorimeter [18]—(b) Example of experimental pressure time course presenting transient and linear regimes with the slope *dP/dt* in the linear regime—(c) wavelength probability density function incident source i.e. LED panel emission spectrum *p*′(*λ*)—(d) photograph of the photoreactor containing the fluorescent photocatalytic system illuminated by the blue LED panel—(e) Example of the wavelength probability density function spectrum *p*′′(*λ*) of photons measured at the rear of the reactor composed of transmitted LED panel photons and fluorescent photons.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0255002.g005

Mechanisms unfavorable for hydrogen production (called loss mechanisms) can be described as follows:

* $ ^1\text{EY}^{2-} $ can relax by non-radiative emission (see the NR₁ mechanism in Fig 4), thermalization or quenching. These phenomena cause a reduction in singlet state concentration.

* Assuming $ ^1\text{EY}^{2-} $ is converted into $ ^3\text{EY}^{2-} $ by ISC, relaxing to the ground state by phosphorescence does not occur at 25°C [53], but non-radiative emission (see the NR₂ mechanism in Fig 4) still exists, resulting in a reduction in the triplet state concentration.

* Photons absorbed by the catalyst do not lead to hydrogen production. Therefore, the study of this specific catalyst implies taking into account a second absorbing species in the models. This will be presented in Section 5.2.

The singlet state of Eosin Y $ ^1\text{EY}^{2-} $ can drop to the ground state by fluorescence emission (see the F mechanism in Fig 4), and this was observed as shown in Fig 5e. As previously discussed, this mechanism cannot simply be described as a loss mechanism at the spatial scale of a photoreactor, since for Eosin Y an overlap exists between molar extinction cross section and emission spectrum. This means that photons emitted by fluorescence at one location, with a wavelength belonging to the overlapping, can be absorbed by Eosin Y at another location and potentially participate in the hydrogen production mechanism [54].
<figure><img src="image_8.png" /><figcaption>**Fig 5.** **Materials.** (a) Photocatalytic system radiative properties—left: Eosin Y (24 µM) UV-visible molar extinction cross-section $E_λ$, fluorescence emission spectrum $p^λ(λ)$ and quantum yield $Φ = 0.35$ in a mixture of triethylamine (10% vol.) in water at pH 10.5 [18, 21]; right: Catalyst (1.66 mM) molar absorption cross section in methanol. UV-visible spectra were obtained using a Shimadzu UV-visible UV-160 A spectrophotometer and emission spectra from a Shimadzu RF-1501 spectrofluorimeter [18]—(b) Example of experimental pressure time course presenting transient and linear regimes with the slope $dP/dt$ in the linear regime—(c) wavelength probability density function incident source i.e. LED panel emission spectrum $p'(λ)$—(d) photograph of the photoreactor containing the fluorescent photocatalytic system illuminated by the blue LED panel—(e) Example of the wavelength probability density function spectrum $p'(λ)$ of photons measured at the rear of the reactor composed of transmitted LED panel photons and fluorescent photons.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0255002.g005

Mechanisms unfavorable for hydrogen production (called loss mechanisms) can be described as follows:

* $^1$EY$^{2-}$ can relax by non-radiative emission (see the NR$_1$ mechanism in Fig 4), thermalization or quenching. These phenomena cause a reduction in singlet state concentration.

* Assuming $^1$EY$^{2-}$ is converted into $^3$EY$^{2-}$ by ISC, relaxing to the ground state by phosphorescence does not occur at 25°C [53], but non-radiative emission (see the NR$_2$ mechanism in Fig 4) still exists, resulting in a reduction in the triplet state concentration.

* Photons absorbed by the catalyst do not lead to hydrogen production. Therefore, the study of this specific catalyst implies taking into account a second absorbing species in the models. This will be presented in Section 5.2.

The singlet state of Eosin Y $^1$EY$^{2-}$ can drop to the ground state by fluorescence emission (see the F mechanism in Fig 4), and this was observed as shown in Fig 5e. As previously discussed, this mechanism cannot simply be described as a loss mechanism at the spatial scale of a photoreactor, since for Eosin Y an overlap exists between molar extinction cross section and emission spectrum. This means that photons emitted by fluorescence at one location, with a wavelength belonging to the overlapping, can be absorbed by Eosin Y at another location and potentially participate in the hydrogen production mechanism [54].
**5.1 Materials and methods**

**5.1.1 Photocatalytic system preparation.** The catalyst synthesis and the photocatalytic system solution preparation protocols are detailed in [18, 21].

**5.1.2 Radiative properties.** Eosin Y UV-visible extinction cross section $E_\lambda$ and catalyst UV-visible absorption cross section $E_{cat,\lambda}$ spectra were measured with a UV-160A Shimadzu spectrophotometer (see Fig 5a) [18]. The Eosin Y fluorescent emission spectrum was measured using a Shimadzu RF-1501 spectrofluorimeter and is presented in Fig 5a. The Eosin Y overall fluorescence quantum yield $\Phi = 0.35$ was identified by minimising the mean squared errors on transmittance between measurements and Monte Carlo algorithm results [54].

**5.1.3 Photoreactor and incident light source.** The reaction was implemented in a laboratory-scale photoreactor extensively described and characterized elsewhere [18, 19, 55]. The gas tight photoreactor of 190 mL is a square-section torus of thickness L = 2.5 cm with two translucent glass faces at the front and the rear (see Fig 5d). It was filled with a liquid volume $V_L = 160$ mL; consequently the headspace gas volume $V_G$ was equal to 30 mL. The liquid photocatalytic system was composed of fluorescent Eosin Y at concentration C and the catalyst at concentration $C_{cat}$, in a mixture of triethylamine (10% vol.) and water [18]. Eosin Y and the catalyst absorb radiation (see Fig 5a). The photoreactor was perfectly mixed and the temperature was controlled at 25.0 ± 0.1°C by water circulation in the photoreactor vessel jacket [55].

Incident radiation was obtained from a 25-LED panel positioned at a distance of 15 cm from the reactor, providing quasi-collimated blue light with a spectrum $p'(\lambda)$ centered at a wavelength of 457 nm (see Fig 5c). The photon flux density entering the medium at the front face of the photoreactor, accurately controlled and easily modifiable, was measured via actinometry experiments, as described in [19].

In such a configuration, the radiative transfer theory can be approximated as a one-dimensional problem.

**5.1.4 Spectral light flux density measurements at the rear of the photoreactor.** To carry out the radiative analysis, we measured the mean photon flux density exiting the reactor $q_r$ measured at the rear glass window using a LI-COR quantum sensor (LI-190) connected to a LI-189 portable LI-COR Quantum meter/Radiometer/Photometer. This mean photon flux density at the rear of the photoreactor was obtained by measuring the photon flux density at 10 different positions. However, this device counts only the total hemispherical photon flux density and was not able to determine spectral distribution. Thus our experimental setup was complemented with a spectrometer (Ocean Optics USB 2000+) equipped with an optical fiber (QP 400-2-SR) and a cosine corrector (Ocean Optics CC-3) coupled to the fiber. The spectrometer was linked via USB connection to a computer running Oceanview software. Spectral distribution at the rear of the reactor $p'(\lambda)$ was thus obtained (see Fig 5e).

**5.1.5 Hydrogen reaction rate measurement.** Hydrogen was produced in the gas tight reactor under irradiation; this was indicated by an increase in pressure (see Fig 5b) measured using a pressure sensor (Keller PA 33X) located in the headspace of the reactor. To analyze this phenomenon, a complete mass analysis of the hydrogen was performed, taking into account a transient regime followed by a linear regime. During the linear regime, the mean volumetric rate of hydrogen production $\langle r_{H_2} \rangle$ can be calculated from the linear variation in pressure P with time t on the basis of the slope $\frac{dp}{dt}$ [55]. A complete analysis of the gas phase was routinely performed with a micro gas chromatograph to demonstrate that only hydrogen is produced.

**5.1.6 A typical experiment.** For a given photocatalytic system (Eosin Y and catalyst concentrations), the incident photon flux density $q_0$ in front of the photoreactor is fixed. The photon flux exiting at the rear of the photoreactor as well as the rear emission spectrum and the