Markov Chains for Modeling Complex Luminescence, Absorption, and Scattering in Nanophotonic Systems

We develop a method based on Markov chains to model fluorescence, absorption, and scattering in nanophotonic systems. We show that the method reproduces Beer-Lambert's Law and Kirchhoff's Law, but also can be used to analyze deviations from these laws when some of their assumptions are violated. We show how to use the method to analyze a luminescent solar concentrator (LSC) based on semiconductor nanocrystals.


Introduction
Recent progress in the study of nanophotonics, encompassing research topics such as metamaterials, plasmonics, and quantum confined matter has established that the optical interaction between light and matter is broadly tailorable across the electromagnetic spectrum via precise control of sub-wavelength scale optical structures [1,2,3,4,5]. Developments in radiation theory, aided by sophisticated computational methods, have informed experimental studies in nanophotonics, for example, as encapsulated in the inverse design approach [6] (and references therein). Complex features of the optical response of nanoscale structures including scattering, absorption, emission, and nonlinear behavior are increasingly well understood [7,8,9,10,11,12,13,14,15,16,17].
However, a significant outstanding challenge is the behavior of macroscopic systems composed of large assemblies of sub-wavelength nanophotonic elements. In such systems complex interactions between individual elements on the nano-and mesoscale give rise to the macroscopic optical response. The challenge is compounded when these systems are optically pumped and driven far from equilibrium, so that simplifying assumptions commonly employed, such as Kirchhoff's law for equating absorptivity and emissivity, are not a priori guaranteed to be applicable [18,19]. Further, standard analytical or computational strategies, such as finite element methods, may be difficult or intractable to extend beyond nanoscale volumes or without simplifying boundary conditions [20].
This work is focused on demonstrating how a computational strategy based on statistical analysis using Markov chains provides a powerful tool for modeling complex interactions between assemblies of nanophotonic elements that give rise to a macroscopic optical response. These interactions are mathematically described as pathways of radiative and nonradiative energy exchange within the optical system and between itself and the environment, including anisotropic absorption, photoluminescence, scattering, and loss. We show how this approach gives deep insight into the steady state behavior of the distribution of light and energy throughout the system, as well as other important optical properties, such as the angle-dependent emission or scattering from the assembly as a whole, even when the system is driven out of equilibrium.
Markov chains have been used to model lightmatter interactions before, particularly in the context of radiative transfer, for example, see [21,22]. Mathematically, Markov chains also share some similarities with the more commonly used computational approach of Monte Carlo ray tracing. While Monte Carlo ray tracing is popular, it is computationally taxing and it has been shown that the use of Markov chains applied to the same problems is computationally more efficient [23,24]. Monte Carlo methods approximate a distribution by sampling many instances from it, while the Markov Chain method allows us to calculate the distribution directly.
Moreover, we show how modeling nanophotonic assemblies using a Markov chain is also conceptually appealing. The strategy is based on identifying the probabilities that connect the pathways of energy exchange between the nanoscale elements of the system, for example the angle-dependent emission or absorption profile of individual emitters. These descriptions of the optical behavior of individual elements in the ensemble can be obtained using other methods, for example Mie theory [25] or extended Mie theory [26], finite element modeling [27] or finite-difference timedomain [28]. Once the probabilities are identified the steady state behavior of the system is obtained simply by solving one eigenvalue problem.
This work is organized as follows: Section 2 describes the basics of a Markov chain for the kinds of systems we will consider. Section 3 fills in the details of these Markov chains to analyze the optical behavior of three specific systems. The first two systems provide an analysis of Beer-Lambert's law and Kirchhoff's law, respectively. The third models the behavior of a luminescent solar concentrator (LSC) comprised of semiconductor nanocrystals. Section 4 concludes the paper.

Methods
The primary method used to model optical interactions in this work is a probabilistic model known as a Markov chain. First, we give a mathematical introduction to discrete-time Markov chains, Sec. 2.1. Then in Sec. 2.2, we set up the specific Markov chain for the optical systems we analyze.

Markov Chains
Using terms defined below, a discrete time, timehomogeneous Markov chain is an appropriate modeling device for a "memoryless" system in which the transitions at any moment depend only on the current state of the system, and not how the system has been at various points in the past, or how long it has been operating. Under certain simple conditions, the probability distribution of the system converges to a unique stationary distribution, that can be used to find the steady state of the system it represents. Our introduction to Markov chains loosely follows [29,30,31].
A Markov chain is formally defined as a stochastic process that obeys the Markov property. A stochastic process is a collection of random variables X = {X t : t ∈ T } indexed by time. Here, we consider a discrete set of times T . This "time" variable merely counts the number of elapsed transitions between states. We do not assume any correspondence of the number of time steps to a precise measurement of physical time in an experiment. The point of the model will not be to determine the precise time at which any event occurs, but rather the steady state distribution of the photons in the system. The development here assumes that the lifetime of each state is roughly equal, and that transition times are negligible, but the model can be modified to accommodate known divergences from these assumptions.
The random variables X t all take values in a finite set S of possible states. Since the set of states is finite, the probability distribution for a random variable X t is determined by the probabilities of the individual states, P (X t = i). For convenience, we denote i , and the vector of probabilities for all states as If the variables are all independent of each other, then the joint distribution P ( i2 . . . . But a general stochastic process may have arbitrary dependencies among the variables, so that computing this joint distribution would require information about the conditional probabilities of each variable on every combination of the earlier ones. However, if the variables obey the Markov property, then these conditional probabilities are determined by the value of the previous variable, and are independent of the values of all earlier variables. That is, Furthermore, if the process is time-homogeneous, then this conditional probability is independent of t, and we can define This p ji , the transition probability, is the probability of the system transitioning from state i to state j.
By the Law of Total Probability, or d If we think of d (t) as a column vector and let P be the matrix with p ji in the jth row and ith column, called the transition matrix, then this means Since the columns of P sum to 1 (as they represent the probabilities of all values of X t conditional on a particular value of X t−1 ), it is said to be a stochastic matrix.
A Markov chain is said to be reducible if there are states i and j such that, if one of the variables ever takes value i, then no later variable ever takes the value j, and irreducible otherwise. It is said to be periodic if there is some integer n > 1 such that X t = X t only if |t − t | is a multiple of n, and aperiodic otherwise. All our Markov chains are irreducible and aperiodic.
Each vector d (t) represents the probability distribution of the system at a time. If there is a distribution then it is said to be a stationary distribution of the system. For an irreducible, aperiodic Markov chain, there is always a unique stationary distribution [29,30,31], and furthermore, for any initial distribution Thus, for such a Markov chain, the stationary distribution represents the steady state behavior of the system, and we can find the stationary distribution by solving Eq. (7). Note that this is an eigenvalue equation with an eigenvalue of 1.

Model
We first introduce a system with some simplifying assumptions in order to establish a concrete connection between the Markov chain formalism (Sec. 2.1) and fluorescence, scattering, and absorption. In this work, a system comprises a medium and an environment, which interact with each other via photon exchange. The environment sends photons into the medium at a constant rate and angular profile.
The medium could be a solid such as a semiconducting slab (Sec. 3.1), a colloidal suspension of nanoparticles (Sec. 3.1), or nanoparticles embedded in a solid polymer matrix (Sec. 3.2 and 3.3). In this section, the medium is a solid matrix of uniform refractive index n m that hosts uniformly distributed fluorophores. The environment around the medium has refractive index n e .
The emission of the fluorophores is assumed to be azimuthally symmetric. Thus, we represent angle of emission just by the zenith angle θ, and assume emission into any angle is proportional to some function f (θ).
We consider two kinds of fluorophores and label them as isotropic or dipolar. By an isotropic fluorophore, we mean that the fluorophore emits radiation with equal probability in all directions; we also assume that there is no anisotropy in its absorption cross-section σ. Thus, for isotropic fluorophores, f (θ) = 1. In contrast, the dipolar fluorophores display anisotropy in the emission and absorption-cross section, specifically f (θ) = sin 2 θ, see Ref. [5]; and σ(θ) ∝ sin 2 θ, see Ref. [32]. Further, all of the dipolar fluorophores emit radiation with the dipole oriented along the same axis, normal to the top surface of the medium. We will later draw a parallel between this behavior and ensembles of homeotropically aligned nanorod-shaped semiconductor nanocrystals in LSClike devices.
We also assume that the light emitted by either type of fluorophores is monochromatic. At first we assume light coming from the environment is the same wavelength as the fluorophore emission, but in Sec. 3.3 this assumption is relaxed. The host material of the medium is non-absorbing, and the real part of the fluorophores' refractive index is matched to that of the host medium. We assume the fluorophores are far enough apart so that near-field effects and resonant energy transfer can be neglected in the probabilistic description for energy transfer. Finally, we assume that the medium is of infinite extent in the x-and y-directions and has a finite thickness L in the z-direction. Now we are in a position to describe the state space S, random variables X t with distributions d (t) , and transition matrix P, as they relate to this system. The Markov chain tracks the movement of an individual photon through the system. For a system with multiple photons, the stationary distribution of each photon will be the same, and the observed distribu-tion of photons in the steady state will be proportional to this distribution.
The state space S, is constructed by stratifying the medium into M equal thickness parallel slabs. We call these slabs layers. For a medium of thickness L, a layer has thickness ∆z = L/M . Each layer is treated as a state representing a possible location for the photon at a time. In all our simulations reported here, we use M = 1000. We found that increasing M to 10, 000 only made a difference to the results on the order of 10 −8 .
X t represents the "position" of a photon at a time t, where positions are layers. Note that we only track the position of a photon along the z-axis since X t takes values from S which is created by discretizing the system along z. It does not account for the location of the photon in the x and y directions. Our standard assumption is that the systems are large enough and homogeneous enough in these horizontal dimensions that we can treat them as infinite, though see further discussion in Sec. 3.3.
The probability d (t) i = P (X t = i), represents the probability that a photon had its most recent absorption event at layer i at time t, that is, the probability that a photon is "located at" this layer. In describing a system with multiple photons, we will sometimes refer to d A transition probability p ji of the system is identified as the probability that a photon from layer i at time t − 1 moves to layer j at time t. To move from layer i to layer j, a photon must be emitted from layer i, transmitted through all intervening layers, and absorbed at layer j. In this simplified model the probabilities of these three events depend only on the zenith angle at which the photon is emitted; thus, we calculate each of these three probabilities at each angle, multiply them, and integrate over all possible angles of travel between these two layers. (As discussed in Sec. 2.2.3, these are the angles from 0 to π/2.) Denoting the emission probability density as dp (e) (θ), the transmission probability as p (t) (θ, i, j) and the absorption probability as p (a) (θ), the transi-tion probability p ji can be written as We describe dp ( Fig. 1b. Importantly, since the emission probabilities are assumed to be proportional to f (θ), and the transmission probability depends only on the path traversed from one layer to the next, these transition probabilities do not depend on any aspect of the system other than what layer the photon was most recently absorbed by, so the system obeys the Markov condition.
In the physical system we are modeling, photons are sent into the system at a constant rate, move according to the transition probabilities from layer to layer within the medium, and eventually escape. While it does have a unique steady state, such a system is not strictly a Markov process, because of the entry and escape of photons from the set of states. However, for any such system, there is a related fictitious system with the same transition probabilities among the states in the medium, but with an extra dummy state. Photons entering the medium transition from the dummy state to the medium, while photons escaping the medium transition from the medium to this dummy state. That is, in this fictitious system the photons entering the medium are the photons that just escaped. This fictitious system is a Markov process because the photon number is conserved, and therefore has a unique, identifiable steady state distribution.
This fictitious system does not behave like the physical system outside the steady state; in the fictitious system the number of photons entering the system at a time t is the number that escaped at time t−1, while in the physical system it is constant. However, the steady state distributions must be the same. Both systems have the same transition probabilities among the states within the medium. Furthermore, in the steady state of the physical system (assuming for now no non-radiative loss, but see Sec. 2.2.7), the number of escaping photons will equal the number of entering photons, just like in the fictitious system. Thus, our analysis focuses on finding the steady state of this fictitious Markov process, which is also the steady state of the physical system.
We refer to the dummy state as "the environment". The entire set of states is S = {0, 1, . . . , M }, where layer 0 is the environment, while the remaining layers, 1, . . . , M , belong to the medium, see Fig. 1. Transition probabilities from and to the "environment" layer are discussed in Secs. 2.2.5 and 2.2.6 respectively.
The steady state distributions d (s) are eigenvectors, and thus can be scaled by any constant. We will always display steady state distributions scaled so that the population of the "environment" layer is 1. This can be interpreted physically as treating all modeled systems as having the same intensity of light incident on them, and comparing the relative populations of the layers of the different systems under this condition. The populations of the layers are sensitive to how thinly the medium is sliced into these discrete layers, but the trend will be the same. Thus, the normalized steady state distribution is given by

Emission Probability Density
We develop the emission probability density by first considering emission from a layer inside the medium. We will discuss emission from the environment to the medium on a case-by-case basis as different scenarios are considered in Secs. 3.1, 3.2, and 3.3.
For a layer of the medium, the emission probability density dp (e) (θ) is defined as the probability that a photon is emitted at an angle whose zenith is infinitesimally close to θ. If f (θ) is proportional to the probability of emission in solid angle dΩ = sin θdθdφ, then integrating over azimuth angle φ yields dp (e) (θ) = 2πνf (θ) sin θdθ, where ν is a normalization constant to ensure that the total probability is 1. That is, For isotropic fluorophores, f (θ) = 1 and ν = 1/4π. For dipolar fluorophores, f (θ) = sin 2 θ and ν = 3/8π.

Absorption Probability
For a photon traveling at angle θ that starts outside a layer and then reaches it, the probability of being absorbed by that layer can be determined from the well-known Beer-Lambert law [33] (see Ref. [34] for a derivation of Beer-Lambert's law based on stochastics). Recall that ∆z is the thickness of the layer, so the path length through the layer along this angle is ∆z/| cos θ|. Let α(θ) = ρσ(θ), where ρ is the fluorophore concentration and σ(θ) is the absorption cross-section (which has angle dependence in the case of anisotropic fluorophores). Then, the probability of transmitting through the layer is e −α(θ)∆z/| cos θ| while the probability of being absorbed by the layer is given by

Transmission Probability
For a photon that has been emitted by one layer i at an angle θ, to find the probability of it being transmitted to another layer j, we calculate the probability of it not being absorbed along the path from layer i to layer j. The length of the full path the photon travels before reaching layer j depends on where precisely in layer i the photon was emitted. We approximate by assuming that all emissions occur from the center of the layer, and divide the medium into a large number of layers to ensure that this is a good approximation.
The path from layer i to layer j can be direct, or can involve reflections. In all our models, the top surface of layer 1 is an interface between the medium and the environment, and reflection occurs if θ is outside the escape cone; the escape cone is formed by the critical angle θ c which is determined from Snell's law. In Sec. 3.1, the bottom surface of layer M is also an interface between the medium and the environment that reflects at these angles, while in Sec. 3.2 and 3.3 the bottom surface is an ideal mirror that reflects at any angle.
These reflections change the angle of travel from θ to π − θ and vice versa. But in all our systems, dp (e) (θ) = dp (e) (π − θ) and p (a) (θ) = p (a) (π − θ). Thus, in these calculations, we will always treat θ as a number between 0 and π/2, understanding that it might represent the angle from the vertical in either direction, so that reflections do not change it.
Because direct paths and reflected paths are two separate ways for a photon emitted at angle θ from layer i to reach layer j, the transmission probability is the sum of the direct and reflection transmission probabilities Figure 2a depicts the paths associated with the direct and reflection (one reflection) transmission probabilities. Figure 2b shows a probability tree for the transmission probability. The total transmission probability Eq. (14) is what goes into the computing the transition probability from layer i to j using Eq. (9). However, as we observe at the end of this section, much of the reflection term is negligible in cases where the medium is thick and highly absorbing.
Under the assumption that emission occurs at the center of layer i, the length of the direct path from layer i to layer j is ( Fig. 1a. Thus the direct transmission probability from layer i to layer j along an angle θ is (Since θ is now assumed to be between 0 and π/2, we no longer need to take the absolute value of cos θ).
Calculating p (t) R (θ, i, j) involves considering several cases, depending on whether a mirror is present on The transition probability p ji is obtained by integrating the transition probability density over θ from 0 to π/2. (b) A transition diagram for the system. The rectangles represent the states and the arrows represent the transition probabilities.
the bottom layer, and whether θ is within the escape cone for medium-environment interfaces. For angles within the escape cone θ c , there is only one possible reflected path, which occurs only if there is a mirror on the bottom layer, and the photon initially travels downward from layer i and is reflected back up to layer j. But for angles outside the escape cone, the light can initially travel either upwards or downwards, and can reflect any number of times, regardless of whether the top and bottom are mirrors or interfaces with the environment.
Thus, for θ < θ c , the probability of transmission by reflection, p only reflection is a single reflection of a downwards path, But if θ ≥ θ c , a photon can travel from i to j by means of any number of reflections of either an upwards or downwards path. It is convenient to write p (t) R as the sum of several terms: p ↓O R , representing transmission probability along any path that starts by moving downwards and has an odd number of reflections; p ↓E R , representing transmission probability along any path that starts by moving downwards and has an even number of reflections; p ↑O R , representing transmission probability along any path that starts by moving upwards and has an odd number of reflections; and p ↑E R , representing transmission probability along any path that starts by moving upwards and has an even number of reflections.
For each of these terms, there is a shortest path, with either one or two reflections. The longer paths in each term differ in length by 2M ∆z/ cos θ, since they cross the entire thickness of the medium twice to get two more reflections. Thus, the transmission proba-bility along each path is e −2α(θ)M ∆z/ cos θ times the probability along the path with two fewer reflections. The sum of these transmission probabilities is thus a geometric series, whose value can be calculated by multiplying the first term by 1/(1−e −2α(θ)M ∆z/ cos θ ). To find the first term in the series, we calculate the shortest path length starting with an emission in the relevant direction and either one or two reflections. The shortest path for a photon that starts downward and has an odd number of reflections was given in equation (16), and the others are calculated similarly: We can input these expressions into Eq. (14) and see the following. When there is no mirror at the bottom and θ < θ c , then p (t) (θ, i, j) is just given by the term in Eq. (15). When there is a mirror at the bottom and θ < θ c , then p (t) (θ, i, j) is the sum of the terms in equations (15) and (16). When θ ≥ θ c , and there is either a mirror or an environmental interface at the bottom, then p (t) (θ, i, j) is the sum of the terms in Eqs. (15), (17), (18), (19), and (20).
Note that many of these terms depend on e −α(θ)M ∆z/ cos θ . This is the transmission probability of a photon going through the entire thickness of the medium at angle θ. When the medium is thick and highly absorbing at angle θ, this probability is very close to 0. In this case it is a good approximation to treat the denominator of the fraction as 1 (that is, to take just the first term in these infinite series of reflections, rather than the sum of all terms). It is also a good approximation to ignore the even terms entirely (since even two reflections involves passing completely through the thickness of the medium). Additionally, the odd reflection terms can largely be ignored except for the ↑ O term when i and j are both small, or the ↓ O term when i and j are both close to M .

Self-absorption
In Secs. 2.2.2 and 2.2.3, we made the assumption that i = j. This was essential because it allowed us to multiply a single absorption probability by the sum of the transmission probabilities along all paths. But when i = j, the probability of absorption on the direct path is different from the probability of absorption on reflected paths. On the direct path, a photon emitted from the center of the layer only has half of the thickness of the layer in which it can be absorbed, instead of the full thickness of the layer, though it can be emitted in the upward or downward direction on this angle. But on the reflected paths, the photon can be absorbed anywhere along the full thickness of the layer. To find the probability of self-absorption on the direct path, p Thus, the transition matrix element p ii is given by

Transitions from Environment to Medium
In all systems we consider, photons enter the system from a source above layer 1, with some characteristic distribution for that system. Thus, we model transitions from the "environment" layer to the medium by Although the integration goes up to the critical angle, the emission density may be zero for some range of angles less than the critical angle. Specifically, if light from the environment impinges on the medium at an angle θ e then dp (e) is zero in the range θ m < θ ≤ θ c .
Here θ m = sin −1 ([n e /n m ] sin(θ e )). The amount of light, or the emission density dp (e) env , from the environment will differ in the different systems we consider, and will thus be described in each of the relevant sections. p (a) (θ) is the same as in Sec. 2.2.2.
However, p (t) (θ, 0, j) has some more significant differences from p (t) (θ, i, j). Light from the environment is incident precisely at the border of the medium, rather than in the middle of a layer, so the distance traversed is an integer number of layers, rather than a half integer. Furthermore, light from the environment only moves in the downward direction.
Transmission probability from the environment to layer j along the direct path is given by For θ within the escape cone, with no mirror on the bottom, the direct term is the only term, so p (t) (θ, 0, j) = p (t) D (θ, 0, j). For θ within the escape cone, in cases with a mirror on the bottom, Light from outside the system only enters at angles within the escape cone, so we do not consider transmission probabilities at angles outside the escape cone.

Transitions to Environment
Photons that escape the system are said to transition to the "environment" layer. While other transitions can happen at any angle from 0 to π/2, these transitions can only happen at angles within the escape cone. Thus, dp (e) is just as in Sec. 2.2.1. We no longer need to account for absorption, because photons that reach the interface with the environment within the escape cone automatically escape.
When there is no mirror on the bottom layer, there are two direct paths from layer i to the environment -an upward and a downward path. Thus, When there is a mirror on the bottom layer, there is a direct upwards path from layer i to the environment, and a downward path with one reflection. Thus, Finally, we must consider transitions directly from the environment to itself. p 00 = θc 0 dp (e) env (θ)p (t) (θ, 0, 0).
As above, we only consider angles within the escape cone, and there is no need for an absorption term. As in Sec. 2.2.5, dp (e) env (θ) depends on the system. As mentioned in Sec. 2.2.5 dp (e) may be zero for some angles.
p (t) (θ, 0, 0) consists only of a single term, and is either e −α(θ)2M ∆z/ cos θ or e −α(θ)M ∆z/ cos θ , depending on whether there is or is not a mirror at the bottom of the medium.

Variable Quantum Yield
So far, we have assumed that every photon that is absorbed by a layer is re-emitted by that layer. That is, we have assumed there is no non-radiative loss, and in particular that the quantum yield of the fluorophores η Q = 1. However, when there are such losses, these can be modeled as just another sort of transition from each layer of the medium to the dummy "environment" layer. We will consider nonradiative losses due to non-unity quantum yield of the fluorophores.
To transform the transition matrix P assuming unity quantum yield to the matrix P accounting for non-unity quantum yield, we need to slightly modify the transitions. All transition probabilities due to emission of photons from the medium are scaled down by the quantum yield, while the remaining probability is added to the probability of transition to the "environment". Transitions from the environment are unchanged.
The transformation described by Eqs. (28), (29) and (30) preserves the stochastic nature of the transition matrix and the total number of photons in the fictitious system will be conserved. Its steady state will still correspond to the steady state of the physical system, in which the sum of all losses (radiative and non-radiative) is equal to the number of photons entering the system.

Deviations from Beer-Lambert's Law
The well-known Beer-Lambert Law allows one to analyze how a non-fluorescent and non-scattering body absorbs and transmits radiation. A standard experimental geometry usually consists of a collimated light source incident on a sample with known thickness, so the path length of the light as well as the incident and transmitted intensity can be quantified. From this, one can determine the absorption coefficient of a sample, or if the absorption coefficient is known, one can predict the transmitted intensity.
However, when the sample is strongly fluorescent and/or scatters light at the same wavelength as it absorbs or emits, it is not straightforward to analyze the optical response. Photons can be absorbed, reemitted, and/or scattered multiple times before exiting the material, and the overall behavior is not well described by the Beer-Lambert Law. By explicitly accounting for these more complex interactions, our model can predict the transmitted intensity; angledependent emission and scattering profile into the environment; and the light intensity, i.e. photon population, as a function of position in the medium. The latter information can be used to calculate local variations in the photo-induced chemical potential.
For analyzing these deviations from Beer-Lambert's law we consider two systems: one is a slab of the photoluminescent semiconducting material gallium arsenide (GaAs) and the other is a colloidal dispersion of gold nanoparticles with strong optical scattering. In both systems, the medium is sandwiched between two environment layers (no mirror at the bottom), see Fig. 3a and Fig. 4a. Monochromatic light from the environment impinges upon the media at normal incidence i.e. θ e = 0.
Because light from the environment is received only at a single angle, the probability density is given by the Dirac delta function dp and

Gallium Arsenide Slab
We model GaAs as an isotropic emitter. We take the absorption coefficient and the refractive index of GaAs at a wavelength corresponding to its bandgap energy at room temperature, i.e. α = 14247 cm −1 , and n m = 3.63. The GaAs slab is 10 −3 cm (or 10 microns) thick, as such the entire slab absorbs 99.99994% of light impinging upon it at normal incidence, assuming the light is monochromatic. Fig. 3b plots the steady state population as a function of depth for different cases that analyze how its behavior depends on changes in quantum yield η Q with the incident light intensity held constant.
The case of η Q = 0 corresponds to when Beer-Lambert's law should be strictly valid (no emission, only absorption). This case is given by the last trace in Fig. 3b. The constant negative slope corresponds to exponential decay, consistent with usual Beer-Lambert behavior. This log-scaled plot shows the steady state photon population as a function of depth for various fluorescent quantum yields in a 10 µm thick slab of GaAs illuminated at its band gap energy at normal incidence.
Deviations from the Beer-Lambert law occur when emission is "turned on", i.e. η Q = 0. As η Q increases, the steady state population of each layer increases monotonically. At any η Q , population decreases with depth -it is always higher towards the top, where light is incident, and lower towards the bottom, where light is only emitted to the environment. The decrease of population with depth is less pronounced for higher η Q , but even with η Q = 1, the steady state distribution is not uniform. Most significantly, the Markov model enables us to calculate the deviation from strict exponential decay for any η Q > 0. Such deviations can be clearly seen in Fig. 3. Upon close inspection, one observes a subtle tilde-shaped features to plots in Fig. 3 especially for 0.5 ≤ η Q ≤ 0.99. These features are due to total internal reflections. The internally reflected light, reflecting off the top and bottom interfaces has a higher probability of being absorbed in the layers nearest to the interface where the photons were reflected. This causes an increase in the population near the top and bottom of the medium leading the to curvy features observed in Fig. 3.

Colloidal Gold Nanoparticles
We consider a colloid composed of gold nanospheres in water. For this model, the "absorption" and "fluorescence" terms developed above are interpreted to represent incoming and outgoing light, respectively, that is scattered by the nanoparticles. As a simple approximation we assume that the gold spheres scatter light isotropically. A more realistic model would involve Rayleigh-Gans or Lorenz-Mie phase functions; the scattering phase function specifies the angular distribution of the scattered light [35]. Since the angular distribution of scattered light depends on the incident angle, a Markov model for these systems would need to track not just the location of the most recent interaction, but also the incident angle, in its state space. Here we use isotropic scattering as a first rough approximation that can be calculated with the simpler state space. Models where the phase function is employed for light scattering using Markov chains are considered in Refs. [22,23,36].
For comparison with common experimental conditions, the number density of the gold particles is 7.15 × 10 15 m −3 , the extinction cross section is 2.93 × 10 −15 m 2 . The thickness of the medium is 1 cm and the refractive index is taken to be that of water 1.33. For our analysis we consider behavior that results when the extinction cross-section σ e , which is the sum of the absorption cross-section σ and the scattering cross-section σ s , remains fixed as the absorption and scattering cross-section are varied. We keep track of only the variation of the scattering cross-section in terms of the extinction crosssectionω = σ s /σ e , whereω is called the scattering albedo. The scattering albedo is varied from zero to one. For a scattering albedo of zero the particles only absorb light (no scattering), whereas for a scattering albedo of one, the particle only scatters light (no absorption).
We assume that the particles and incident light satisfy the conditions required for elastic and isotropic scattering. In the case of zero scattering albedo, one sees from Fig. 4 that there is an exponential decay of photons (the vertical axis is logarithmic), as expected. As the scattering albedo increases, the photons have a chance of being scattered. It is apparent from Fig. 4 that the attenuation increases as the scattering albedo decreases.

Kirchhoff 's Law and Violations of Detailed Balance
In this section we reproduce Kirchhoff's law of radiation, and analyze deviations from it. Kirchhoff's law of radiation states that the emissivity of a medium at any angle is equal to its absorptivity at that angle, in thermal equilibrium. This is a stricter condition than a steady state, in which it is only required that the total flow of photons from the medium to the environment is equal to the total flow of photons from the environment to the medium.
Kirchhoff's law requires that in equilibrium there is also a detailed balance, meaning that the flow of photons from the medium to the environment at an angle is equal to the flow of photons from the environment to the medium at that same angle. We show that although this detailed balance holds in equilibrium systems, and is very well approximated in a non-equilibrium steady state system with isotropic fluorophores, it is clearly violated in a non-equilibrium steady state system with aligned dipolar fluorophores.
We consider two types of environmental radiation interacting with a medium. The first is an isotropic radiation field characteristic of thermal equilibrium, e.g. the interior radiation field of a blackbody cavity. Under these conditions, light is incident on the medium with a Lambertian angular profile [37], with dp (e) env (θ) = 2 sin 2 θ c cos θ sin θdθ, for θ < θ c . The second radiation source we consider is light received from the solar disk, i.e. sunlight, at normal incidence to the medium. We use dp (e) env (θ) = ν sin θ cos θdθ for θ ≤ 0.267 • and dp (e) env (θ) = 0 for θ > 0.267 • , where ν is a normalization constant equal to n 2 m /(πn 2 e sin 2 (0.267 • )). Since our medium emits at all angles, but absorbs only at angles corresponding to the solar disk, its steady state is not an equilibrium.
We also consider two versions of the medium, one with isotropic fluorophores and the other with dipolar fluorophores. The absorption, transmission, and emission probabilities for these are defined in Sec. 2.2. For the isotropic case, we consider a slab of GaAs with parameters (thickness, refractive index, absorption coefficient and emission/absorption wavelength of light) as described in Sec. 3.1.1. For the dipolar case, we consider cadmium selenide/cadmium sulfide (CdSe/CdS) dot-in-rod crystals that are in a chlo-roform solution, and are aligned with their long axis normal to the surface of the medium. The concentration of the CdSe/CdS dot-in-rod fluorophores is 10 −6 molar concentration or molarity (M), extinction coefficient of the fluorophores is 5 × 10 6 cm −1 M −1 , the wavelength of emission/absorption is taken as 615 nm. The refractive index of the medium is that of chloroform at 615 nm, i.e. 1.445. The thickness of the medium is 2 cm.
In each of the four cases (Lambertian or solar incident radiation; isotropic or dipolar fluorophores) we always have a mirror at the bottom of layer M , and assume the refractive index of the environment is n e = 1.
To probe the angle distribution of photons entering and leaving the system, we partition the range 0 ≤ θ ≤ θ c into a number b max of bins, each of size ∆θ = θ c /b max . We index these bins by a value b ranging from 1 to b max . The bin b consists of angles from (b − 1)∆θ to b∆θ. Note that for the equilibrium case the light entering the medium, upon refraction, spans the range θ ≤ θ c . However, for the non-equilibrium case, the angular range of light entering the medium is less than θ c . In particular, it spans the range θ ≤ sin −1 ([n e /n m ] sin(0.267 • )). However, we still use the same bin range while keeping in mind that the absorption corresponding to bins in the range sin −1 ([n e /n m ] sin(0.267 • )) < θ ≤ θ c is zero.
We analyze transitions in greater detail now by identifying the transition probabilities at angles within each bin. That is, we break up the definitions of the transition probabilities p 0i and p j0 so that instead of integrating over all angles, we integrate just over angles within the bin. We will have p 0i = bmax b=1 p 0i (b) and p j0 = bmax b=1 p j0 (b). For i > 0, p 0i (b) is the probability of a photon in layer i to leave the medium at an angle in bin b: as in Eq. (26). Similarly, for j > 0, p j0 (b) is the probability of a photon from the environment to reach layer j at an angle in bin b: as in Eq. (23). p 00 (b) is the probability of light from the environment entering the medium, being reflected off the mirror on the bottom, and emerging again without being absorbed in the medium. That is: as in Eq. (27).
Once we have these detailed probabilities, we can analyze the total absorption and emission of the medium within each of these bins by summing over all layers. In the steady state, the angle-dependent absorption of the medium at angles in bin b is and the angle-dependent emission of the medium at angles in b is Note that a(b) can be simplified to b∆θ (b−1)∆θ dp (e) env (θ), because the sum of p (t) (θ, 0, j)p (a) (θ) is 1 (since all light must end up somewhere), though there is no comparable simplification of e(b). Because this is the steady state, it is clear that Eqs. (37) and (38) are plotted for the isotropic fluorophores and the dipole fluorophores in Figs. 5b and 5d, respectively, for the case of an isotropic incident radiation field (i.e. the equilibrium state). One sees that the angle-dependent absorption and the angledependent emission are coincident which implies that Kirchhoff's law is upheld. In contrast, the non-equilibrium case shows that the angular distributions of the incoming and outgoing light does not match, i.e. a(b) = e(b). This can be seen in Fig. 6 where the incoming light spans 0 ≤ θ ≤ 0.267 • (not pictured). Light exiting the medium for the isotropic case is shown in Fig. 6b while the exiting light for the aligned dipolar fluorophores it is shown in Fig. 6d. It is commonly assumed that most bulk materials fluoresce in a Lambertian manner, like equilibrium radiation fields. We can see that this happens when the system is composed of isotropic fluorophores, but not when it is composed of aligned dipolar ones, at least with the thickness of medium that we consider. Interestingly, the anisotropy of the emission in the dipolar case is driven more by the lack of reabsorption along certain angles, rather than by the excess emission along other angles.

Light Trapping in Luminescent Solar Concentrators
Luminescent solar concentrators (LSCs) have been proposed as an alternative to conventional lens-based optical concentrators for improving the performance and lowering the cost of high efficiency photovoltaic (PV) solar converters, see Ref. [17] and references therein. An LSC consists of fluorophores such as dye molecules or semiconductor nanocrystals embedded in a dielectric medium. Light is absorbed by the fluorophores and re-emitted. Refractive index contrast between the medium and environment means that some re-emitted light is trapped by total internal reflection into "waveguide modes". PV elements are then placed at the edges of the LSC, where light in the waveguide modes is directed.
The performance of an LSC, and the consequent PV conversion efficiency, depends on the ratio of light occupying waveguide modes compared to the incident sunlight. A useful proxy for this value is the fraction of light each fluorophore emits into the waveguide mode, η trap , which for isotropic emitters is 1 − (n e /n m ) 2 or ∼ 74% for a typical organic polymeric waveguide medium (n m /n e = 1.5) [38].
Anisotropic fluorophores, such as aligned semiconductor nanorods, have been proposed as fluorophores because they emit preferentially into the waveguide modes, and can thus provide a higher η trap . However, a full analysis of LSC performance based on either fluoruphore should also account for the efficiency of absorption of incident light, as well as further reabsorption and re-emission events. Monte Carlo simulations of three-dimensional structures have been used to account for these complexities, as well as scattering and non-unity quantum yield [39,40,41]. Markov models of the sort described in this paper can make this fuller analysis more computationally tractable.
An LSC can be described with the same geometry analyzed in Sec. 3.2, with a mirror on the bottom of the medium. We assume that the fluorophores are suspended in a solid matrix of refractive index 1.49, while the refractive index of the environment is 1. We consider both isotropic fluorophores, as in Fig. 5a, and dipolar as in Fig. 5c. We assume diffuse incident light with a Lambertian profile. We do not track the wavelength of light specifically, but to account for a Stokes shift, we use an absorption coefficient for incident solar radiation of 5 × 10 8 cm −1 M −1 , while we use 5 × 10 5 cm −1 M −1 ( ∼ 1000x lower cross section) for light emitted from the fluorophores, in agreement with common values for semiconductor nanocrystals [17].
We can use the previous analysis to find the steady state population of excited fluorophores in the medium, M i=1 d i /d 0 . The trapping ratio of re-emitted light is η trap = π/2 θc dp (e) (θ) (which is 1 − (n e /n m ) 2 for isotropic emitters). With quantum yield η Q , the net occupation of the waveguide modes, which we will call the "occupation factor", is given by the product of these three factors, Because our model assumes an infinite horizontal extent of the medium, it does not directly track how much light reaches the photovoltaics at the edges. In-stead it identifies how much light occupies the waveguide modes at any point in the medium, just as the analysis based on η trap alone does. To fully track the three-dimensional motion of the photons from initial incidence until PV conversion, the Markov model would require states representing differential volume elements, rather than thin layers, as well as transition probabilities between them, but the process of finding the stationary distribution would be the same. Figure 7 plots the occupation factor with dipolar and isotropic fluorophores and variable quantum yield, as a function of the concentration of nanocrystals in the polymer matrix of the LSC, spanning a range commonly studied experimentally. Values greater than 1 can be achieved because a single photon may be absorbed and re-emitted several times. Note that dipolar fluorophores achieve higher occupation factors when concentration is high, due in part to higher η trap , while they have lower occupation factor at low concentrations due to absorbing less incident light than isotropic fluorophores.

Conclusion
We have demonstrated the use of a Markov model to study the steady state behavior of several optical systems. The method works by finding an eigenvector of the matrix of transition probabilities for photons moving through the system. We have shown that this method yields the expected results in systems that allow for an analytical solution (Beer-Lambert's law or Kirchhoff's law). This model also gives a description of steady states for which analytical solutions are not available. We propose that this method could be particularly useful for describing complex mesoscale interactions between nanophotonic elements, such as in an LSC.
The particular model we developed assumes translational symmetry in two dimensions, and only tracks the movement of photons in the third dimension. This model also assumes that the angle of emission of a photon is independent of the angle of absorption, which makes less accurate predictions when applied to analysis of scattering. However, more detailed Markov models can track these additional variables to give more precise results. Transition probabilities can also be modified to take into account additional features of radiation, such as near-field interactions, wavelength dependence, and deviation from perfect isotropic or sine-squared emission from the fluorophores.