Transient microscopy has emerged as a powerful tool for imaging the diffusion of excitons and free charge carriers in optoelectronic materials. In many excitonic materials, extraction of diffusion coefficients can be simplified because of the linear relationship between signal intensity and local excited state population. However, in materials where transport is dominated by free charge carriers, extracting diffusivities accurately from multidimensional data is complicated by the nonlinear dependence of the measured signal on the local charge carrier density. To obtain accurate estimates of charge carrier diffusivity from transient microscopy data, statistically robust fitting algorithms coupled to efficient 3D numerical solvers that faithfully relate local carrier dynamics to raw experimental measurables are sometimes needed. Here, we provide a detailed numerical framework for modeling the spatiotemporal dynamics of free charge carriers in bulk semiconductors with significant solving speed reduction and for simulating the corresponding transient photoluminescence microscopy data. To demonstrate the utility of this approach, we apply a fitting algorithm using a Markov chain Monte Carlo sampler to experimental data on bulk CdS and methylammonium lead bromide (MAPbBr_{3}) crystals. Parameter analyses reveal that transient photoluminescence microscopy can be used to obtain robust estimates of charge carrier diffusivities in optoelectronic materials of interest, but that other experimental approaches should be used for obtaining carrier recombination constants. Additionally, simplifications can be made to the fitting model depending on the experimental conditions and material systems studied. Our open-source simulation code and fitting algorithm are made freely available to the scientific community.

## INTRODUCTION

Understanding electronic energy transport properties such as carrier diffusivity in materials can enhance optimization of optoelectronic device performance, including solar cells,^{1} displays,^{2} and sensors.^{3} A variety of techniques have been used to study the movement of free charge carriers or excitons (bound electron–hole pairs) in emerging semiconductor materials,^{4} such as Hall-effect measurements,^{5,6} field-effect transistors,^{7–9} terahertz spectroscopy,^{10,11} and photoluminescence quenching.^{12–14} However, in many materials, the nanoscale morphology and structural inhomogeneities have an impact on transport behavior.^{15–19} Recent studies look to observe both the temporal and spatial movement of charge carriers in optoelectronic materials, including transient grating techniques,^{20,21} transient absorption microscopy,^{18,22,23} transient scattering microscopy,^{19} and time-resolved photoluminescence microscopy (TPLM).^{24–32}

Time-resolved photoluminescence microscopy (TPLM) is particularly attractive for luminescent materials because it is a background-free technique that usually requires lower fluence than other approaches.^{4} In this technique, a diffraction-limited population is excited in the sample by a pulsed laser, and an avalanche photodiode (APD) is raster scanned across a magnified image plane to collect temporal and spatial photoluminescence (PL) data [Fig. 1(a)].

A complete description of carrier dynamics in semiconductors depends sensitively on a range of material parameters and experimental conditions, including photoexcited carrier density, doping density, exciton binding energy, temperature, bandgap, and electronic and vibrational band structure, among other factors, as described in detail by others.^{33–35} Here, we highlight a couple limiting examples to illustrate the impact of the nature of carrier dynamics on the interpretation of TPLM data.

In one extreme, electrons and holes move together as bound pairs and the local PL signal is proportional to the local density of excitons. As a result, the measured data can be directly analyzed to obtain the diffusivity [Fig. 1(b)]. For example, applying a 1D diffusion model enables the collected spatiotemporal data (*S*_{PL}) to be normalized at each time point and fit using a Gaussian or related convolution function.^{24} The diffusivity can then be extracted using the mean squared displacement: $MSDt=\sigma 2t\u2212\sigma 20=2Dt\alpha $, where $\sigma 2t$ is the variance at time *t*, *D* is the scaling factor, and *α* is the diffusion exponent. Under normal diffusion conditions (*α* = 1), *D* is the exciton diffusivity within the material.

In materials where transport is dominated by free charge carriers, a more complex approach is required [Fig. 1(c)]. In many cases, the measured photoluminescence signal is no longer proportional to the local population density. If the photogenerated carrier density is sufficiently high, the PL signal is expected to have a quadratic dependence on the charge carrier density. Accordingly, the measured PL signal is not directly proportional to the charge carrier density and additional steps are required in data analysis. In particular, bimolecular and higher-order recombination processes can contribute to spatial broadening in a TPLM experiment—even in the absence of carrier diffusion—because of the nonlinear dependence of PL signal decay rate on local carrier density.^{26,30,32} Hence, trying to fit the time-normalized spatial data to a Gaussian function and extracting the mean squared displacement in the same manner as the excitonic materials result in a curve instead of a straight line, even under normal diffusion conditions [Fig. 1(c)].

Previous works have developed detailed models to describe the behavior of charge carriers in TPLM experiments.^{26,30,32} Sridharan *et al.* presented a 3D model for fitting diffusion and recombination coefficients and applied that model to carrier transport in thin film perovskites.^{26} In addition to accounting for potential artifacts in such measurements along with photon recycling effects, they performed a time-dependent fit of the diffusion coefficient, which suggested evidence for non-diffusive transport. Recently, deQuilettes *et al.* presented an analytical framework for modeling lateral excitonic and free carrier transport.^{32} The effects of photon recycling and grain boundaries were considered, and modified mean squared displacement formulas for extracting diffusivities were derived for specific dominant recombination pathways.

Herein, a detailed framework for numerical 3D modeling of the free charge carriers in TPLM experiments is presented. We begin with established models relating diffusivity to unimolecular and bimolecular recombination of charge carriers. We then devote particular attention to conversion between the experiment and model frame and outline the process for determining the initial condition for the solver. We describe the numerical implementation of the model in MATLAB, employing a Strang splitting algorithm,^{36} vectorization, dimension reduction, and sparse matrix techniques^{37} to speed up the solution time while maintaining reasonable accuracy. We then demonstrate the analysis capabilities enabled by faster solving speed by performing parameter analysis using a Markov Chain Monte Carlo (MCMC) sampler as part of the fitting algorithm on TPLM data collected from methylammonium lead bromide (MAPbBr_{3}) single crystals and bulk CdS. Using the results of these analyses, we conclude that TPLM can be used to obtain statistically confident estimates of ambipolar charge carrier diffusivities but should be coupled with other independent techniques to obtain accurate estimates for unimolecular and bimolecular recombination constants. Furthermore, depending on the experimental conditions chosen (charge carrier densities) and relative magnitudes of material parameters, reductions can be made to the model to simplify analysis. With our open-source code for the model and the fitting algorithm (see the supplementary material), we present a toolbox that can be used by other members of the community to apply to their material systems.

## MODELING THE TPLM EXPERIMENT

A schematic of the experimental setup, along with model coordinates is shown in Fig. 2. The incident pulsed laser is focused to a diffraction-limited spot on the sample and the emitted photoluminescence is detected by raster-scanning an avalanche photodiode (APD) across the magnified image plane to provide the temporal and spatial dimension to the data.

To extract charge carrier diffusivities from TPLM, the experimental data can be modeled using the following partial differential equation (PDE):^{26,30,32}

where *n* = *n*(*x*, *y*, *z*, *t*) is the charge carrier density as a function of space and time (cm^{-3}), *D* is the diffusion constant (cm^{2} s^{−1}), *A* is the first-order recombination constant (including trapping effects) (s^{−1}), *B* is the second-order recombination constant (including radiative recombination effects) (cm^{3} s^{−1}), and *C* is the third-order recombination constant (including Auger recombination effects) (cm^{6} s^{−1}). To simplify calculations, it is possible to neglect the third-order term by experimentally choosing an incident power that is sufficiently low such that *Cn*^{3} is negligible.

Additionally, anisotropy may be incorporated by defining different diffusivities along each spatial direction. If knowledge of electron and hole behavior is known *a priori*, Eq. (1) may be replaced by separate PDEs modeling the electron and hole density separately (in this case, drift terms must also be included and the equations solved self-consistently with Poisson’s equation).^{34,35} In the absence of specific information about the nature of charge carriers within a material, we make the simplifying assumption of ambipolar transport, in which case the diffusivity *D* in Eq. (1) represents the isotropic ambipolar diffusivity, $D=2DeDhDe+Dh$, where the subscripts *e* and *h* denote the electron and hole, respectively.^{35}

To go from the experimental frame to the model frame, we must develop equations that relate the measured PL signal *S*_{PL,measured} to the local charge carrier density *n*, within the material.

From the laser power and repetition rate, we can calculate the total number of charge carriers photoexcited in the sample per laser pulse, *N*_{0,total}, as

where *P* is the laser power (W), *rep* is the laser repetition rate (s^{−1}), *E*_{photon} is the energy of the photon (J), and *γ* is the fraction of reflected light at the surface. This total number of charge carriers is spatially distributed in the sample at time *t* = 0 such that

with *V* as the total volume of the sample (cm^{3}). We assume that the absorbed charge carrier density at *t* = 0 can be approximated using some function that describes the distribution of charge carriers in the $x,y$-plane depending on the shape of the laser excitation pulse used (often Gaussian), and Beer’s law in the *z*-dimension,

where *n*_{0} is the maximum local charge carrier density in the sample at *t* = 0 (cm^{−3}) and *α* is the absorption coefficient of the sample at the excitation wavelength (cm^{−1}).

An expression for the PL emitted across the sample surface can be written as

where *C* is a constant relating charge carrier density to emitted PL signal (obtained from fitting the initial time point signal) and $\alpha \u0304$ is the absorption coefficient at the sample emission wavelength (cm^{−1}). The actual signal that is measured from the sample *S*_{PL,measured} at each time point is related to *I*_{PL,emitted} by

where *w* is the APD sensor width (cm) and *τ* is the integration time used (s). Using these equations and fitting to the initial time data measured $SPL,measuredx,t=0$, the initial condition for the model PDE can be obtained.

Along with the initial condition, boundary conditions also need to be specified to solve the PDE. A point of symmetry at the center of the focal plane yields

As for the edges of the sample, a general boundary condition can be used to account for

where *v*_{s} is the surface recombination velocity. A reflective boundary can be obtained by setting the surface recombination velocity to zero. Thick samples where the absorption depth is much thinner than the thickness of the sample can be considered as semi-infinite.

Photon recycling may affect the lateral energy transport, particularly in excitonic materials that have short radiative lifetimes—or when considering very long-range transport;^{38} additional details have been discussed by deQuilettes *et al.*^{32} Due to the length and time scales considered for the materials of interest in this study (3D perovskites, bulk semiconductors), the effect was considered negligible in our analysis.

## NUMERICAL IMPLEMENTATION

### Explicit finite difference

The PDE in Eq. (1) requires an accurate but fast solver to perform detailed parameter analysis and estimation. Numerical finite differencing techniques can be used to solve the model at discrete space and time points by approximating the derivatives as finite differences.

Using an explicit finite differencing scheme, the derivatives in the PDE can be replaced by the following algebraic expressions:

where the superscript *t* denotes the timepoint discretization and the subscripts *i*, *j*, *k* denote the *x*, *y*, *z* spatial discretization, respectively. Dimensional analysis comparing the different length scales in the problem can be used to determine the best spatial step scale Δ*x*, Δ*y*, Δ*z* to use. From a fast-numerical solving standpoint, choosing Δ*x* to match the experimental length scale avoids inaccuracies from interpolation without increasing the PDE solution time. For Δ*z*, the relevant length scale is $1\alpha $, and depending on the material and experimental wavelength used, Δ*z* should be just under $1\alpha $ to adequately capture the Beer’s law decay in the z-dimension without excessively increasing the solution time.

As for the simulation time scale Δ*t*, the relevant physical time scales in the problem include $1A,N0,totalB,2R2D,Z2D,Zvs$ and $1D\alpha \u03042$, the relative magnitudes of which are dependent on the material system being studied and experimental conditions used. However, another important factor that needs to be considered when choosing Δ*t* is the Courant–Friedrichs–Lewy (CFL) condition for an explicit finite differencing scheme.^{39} When using explicit solving methods, the condition $D\Delta t\Delta x2+D\Delta t\Delta y2+Dz\Delta t\Delta z2<0.5$ must be satisfied for the solution to converge. As a result, a small Δ*t* is usually preferred to obtain accurate solutions, but with a trade-off for long solution times.

### Numerical techniques for speeding up PDE solver

To perform parameter analysis using the numerical model, an accurate but fast solver is required. Using an explicit finite difference method with the full three dimensions led to a PDE single-solve time of more than 2 h on a typical laptop computer, which is impractical for any statistically robust fitting algorithm. Fortunately, numerical techniques can be used to speed up the solver without compromising accuracy.

For isotropic materials, making the assumption of radial symmetry and transforming to cylindrical coordinates for $nr,z,t$ can reduce the solving speed. To go from $nr,z,t$ to $SPL,measuredx,t$, the calculated $IPL,emittedr,t$ should be interpolated at the desired *x*, *y* points and the rest of the equations introduced above can be used. Figure 3 demonstrates an example numerical model 3D solution at different time points and its conversion to the simulated *S*_{PL,measured} for comparison with experimental data.

When a PDE can be represented as a sum of differential operators, the Strang splitting technique^{36} may be beneficial. Strang splitting is particularly effective at increasing solution speed when the time scales of the two operators differ greatly. As shown in Fig. 4(a), the PDE introduced in Eq. (1) can be split into two separate problems *f*_{1} and *f*_{2}, each of which can be solved faster individually than in the combined problem. *f*_{1}, which represents the diffusion problem, can be solved using finite differencing schemes, while a direct analytical solution for *f*_{2} (the recombination problem) can be used instead. By combining the faster solvers for *f*_{1} and *f*_{2} in the presented scheme, a solution of error $\u223cO\Delta t$ is obtained so long as the CFL stability condition is met.

Other finite differencing schemes were also investigated for the *f*_{1} solver, including implicit finite difference and Crank–Nicolson schemes. The advantage of these methods relative to the explicit finite difference method is that there are no stability conditions limiting the choice of Δ*t*. However, the accuracy of the solution is still affected by the size of Δ*t*, and the solution method requires solving ** Ax** =

**type problems involving matrix inversions, which tend to be computationally more expensive than matrix multiplications associated with explicit methods. A comparison of the solution speed with the explicit finite differencing scheme is presented in Fig. S1. With a combination of Strang splitting, sparse matrix formalisms,**

*b*^{37}and vectorization using MATLAB, the fastest solution speed without compromising accuracy was found to involve an explicit finite difference-based solver for

*f*

_{1}combined with an analytical solver for

*f*

_{2}, as demonstrated in Figs. 4(b) and S1. Overall, these approaches led to a reduction in PDE solving time from >2 h to only a few seconds.

## APPLICATION TO EXPERIMENTAL DATA

With an accurate and fast solver, parameter analyses can be conducted on TPLM data collected from real materials. Lead-halide perovskites have attracted a significant amount of research efforts as a candidate for next-generation photovoltaics, among other optoelectronic devices.^{40–42} Consequently, obtaining robust estimates of ambipolar diffusivity is of great interest in this material class.^{5,6,13,16,27,30}

Millimeter-sized single crystals of MAPbBr_{3} were synthesized following the synthetic protocol detailed in Saidaminov *et al.*^{43} with slight modifications [Figs. 5(a)–5(c), more details available in the section titled Methods]. Power-dependent photoluminescence was measured to determine the transport regime [Fig. 5(d)]. There is a transition from a linear to quadratic dependence of emitted photoluminescence with input power as free photogenerated charges become the dominant charge carriers within the semiconductor. To avoid sample degradation effects, we chose to excite the sample with incident power just beyond the transition regime.

Figures 5(e) and 5(f) show the obtained TPLM data for single crystal MAPbBr_{3}. Raw data are presented on the left and intensity-normalized data are shown on the right to emphasize broadening of the spatial distribution. TPLM experiments were also performed on commercially available single crystal CdS (supplementary material, Note 3). Relevant experimental parameters and material constants used for numerical solution of the PDE are presented in Table I.

Simulation model parameters . | P (W)
. | rep (s^{−1})
. | γ
. | σ(t = 0) (nm)
. | α (cm^{−1})
. | $\alpha \u0304$ (cm^{−1})
. | Δx (μm)
. | Δz (μm)
. | Δt (μs)
. | Fitted time range (ns) . |
---|---|---|---|---|---|---|---|---|---|---|

MAPbBr_{3} | 4.75 × 10^{−8} | 2.5 × 10^{6} | 0.2 | 572 | 2.66 × 10^{5}^{44} | 3.25 × 10^{3}^{44} | 0.101 | 0.0505 | 1.28 × 10^{−6} | 6.4 |

CdS | 1.52 × 10^{−6} | 1.25 × 10^{6} | 0.2^{45} | 478 | 1 × 10^{5}^{45} | 3 × 10^{4}^{45} | 0.101 | 0.0505 | 1.28 × 10^{−6} | 3.2 |

Simulation model parameters . | P (W)
. | rep (s^{−1})
. | γ
. | σ(t = 0) (nm)
. | α (cm^{−1})
. | $\alpha \u0304$ (cm^{−1})
. | Δx (μm)
. | Δz (μm)
. | Δt (μs)
. | Fitted time range (ns) . |
---|---|---|---|---|---|---|---|---|---|---|

MAPbBr_{3} | 4.75 × 10^{−8} | 2.5 × 10^{6} | 0.2 | 572 | 2.66 × 10^{5}^{44} | 3.25 × 10^{3}^{44} | 0.101 | 0.0505 | 1.28 × 10^{−6} | 6.4 |

CdS | 1.52 × 10^{−6} | 1.25 × 10^{6} | 0.2^{45} | 478 | 1 × 10^{5}^{45} | 3 × 10^{4}^{45} | 0.101 | 0.0505 | 1.28 × 10^{−6} | 3.2 |

### Markov chain Monte Carlo sampler fitting algorithm

To obtain statistically robust estimates of carrier diffusivity from TPLM experimental data, we apply a Bayesian framework for parameter estimation.^{46–48} This approach has advantages over the more common Levenberg–Marquardt (LM) fitting algorithms deployed in commercial fitting software, including robustness to local optima in parameter space. The method is applied here using a Markov chain Monte Carlo (MCMC) sampler as part of the fitting algorithm. We note that fast solving speed of the model itself is essential for stochastic fitting routines—such as the MCMC sampler used here or for related genetic optimization algorithms^{49}—because the PDE must be solved repeatedly using unique parameter combinations [Eq. (1) is solved ∼10^{5} times for the statistical analysis presented here]. MATLAB code implementing the MCMC sampler and the fast numerical solver is available in the supplementary material.

The MCMC sampling method we applied to analyze TPLM data was adapted from the previously reported methodology by Winslow *et al.*^{48} and Ashner *et al.*^{50} Briefly, in the defined parameter space, a defined number of “walkers” are initialized at random starting points, with each walker representing a distinct combination of fitting parameters (i.e., *A*, *B*, and *D*). Each walker takes random “steps” around the parameter space and the steps are accepted with a probability that depends on the log-likelihood function, which is determined as below.

Using the Bayesian inference framework for the probability distributions, and because we are only interested in the relative probabilities, we use

where *θ* is the set of parameters of our model $A,B,D$. The first term is the ratio of the likelihood functions and the second term is the ratio of the “priors,” in which the prior knowledge about the parameter space before obtaining *S*_{PL,measured} can be included. We assume that the probability of observing a given data point *S*_{PL,measured} can be described by a normal distribution around the true value from the model such that $PSPL,measured|\theta \u223c\u220fx,te\u2212SPLx,t;\theta \u2212SPL,measuredx,t\sigma e(x,t)2$. The natural logarithm of the likelihood function can, then, be expressed as $logPSPL,measured|\theta \u223c\u2212\u2211x,tSPLx,t;\theta \u2212SPL,measuredx,t\sigma ex,t2$. The log-likelihood function includes the variance at $\sigma e2x,t$ at each of the data points, which can be used to give weightings to different data points. [Note that the variance $\sigma e2x,t$ in the likelihood function is different from the variance of the spatial distribution obtained experimentally in a transient microscopy experiment]. As for the prior, we assume a flat distribution within an upper and lower limit based both on what is physically reasonable for our material system (i.e., no negative values for *A*, *B*, or *D*) and on the order of magnitude estimates for parameters obtained previously from literature.

A two-part search process is applied to find the best combination of model parameters $A,B,D$ that maximizes the log-likelihood function in the defined parameter space. In the first part, the walkers are initialized randomly throughout the entire parameter space with the goal of searching for the global optimum and narrowing down the parameter space. In the second part, a new set of walkers are randomly dispersed in a ball of parameter space that is near the maximum likelihood estimate (MLE) determined in the first part. The second part provides statistical information about the parameters including their uncertainties and correlations.

A key decision in any fitting routine is deciding what should be optimized, known as the objective function. To apply the MCMC sampling method to our model, two different objective functions with different values of $\sigma ex,t$ were considered. For objective function Nos. 1 and 2, we used $\sigma ex,t=SPL,measuredx,t$ and $\sigma ex,t=1$, respectively.

Objective function No. 1 is given by

Objective function No. 2 is given by

Objective function No. 1 normalizes the error signal by the measured intensity $SPL,measuredx,t$ at each space and time point, while objective function No. 2 considers the absolute error. In effect, objective function No. 1 considers the fractional error such that all time and space points have equal weightings, while objective function No. 2 gives a greater weighting to the earlier time points where the experimentally measured signal is greatest.

Figure 6 compares the obtained MLE fits for the TPLM data from the MAPbBr_{3} crystals using each objective function. The fit comparisons between the spatially normalized data [Figs. 6(a)–6(c)] reveal that objective function No. 1 is better at capturing the spatial broadening that occurs over time. However, objective function No. 2 is better at capturing the overall signal decay [Fig. 6(d)].

The difference between the two objective functions can be better understood by observing the parameter statistical distributions and cross-correlations obtained from the MCMC analysis (Fig. 7). In Fig. 7(a), we see that objective function No. 1 results in flat distributions for *A* and *B*, while providing a uniquely determined value of *D*. The determined MLE with 95% confidence intervals for the uncertainty of *D*_{MLE,No. 1} = 0.51 ± 0.02 cm^{2} s^{−1}, which is consistent with literature reports of ambipolar diffusivities in MAPbBr_{3} measured using light-induced transient grating experiments^{51} and on the same order of magnitude when calculated from mobilities obtained using Hall mobility measurements.^{52} The flat parameter distributions for *A* and *B* indicate that these parameters are over-specified for the model, resulting in undetermined values. In other words, *A* and *B* could be arbitrarily varied within the indicated range without changing the quality of the fit. This makes intuitive sense as parameters *A* and *B* are primarily responsible for determining the rate of PL signal decay over time. Since objective function No. 1 considers the intensity-normalized error at every time point, this objective function is less sensitive to signal decay, as shown graphically in Fig. 6(d).

Using objective function No. 2, we find that *A* and *B* are uniquely determined and anticorrelated. The obtained MLE values with 95% confidence intervals for the uncertainty were *A*_{MLE,No. 2} = (3.42 ± 0.01) × 10^{8} s^{−1} and $BMLE,2=8.00\xb10.02\xd710\u221210cm3s\u22121$. *B* is on the same order of magnitude with literature measurements using transient photoluminescence^{51} and transient absorption techniques;^{53} however, *A* is 3–4 orders of magnitude larger than the values obtained from the same literature measurements. Since *A* includes the effects of trap recombination, it is sample-dependent and dependent on the growth conditions. However, the obtained MLE *D* value approaches zero and is undetermined. Again, these results make intuitive sense given the choice of objective function. Objective function No. 2 emphasizes the overall signal intensity, which is most sensitively determined by parameters *A* and *B*. However, this accuracy in tracking the overall signal decay [Fig. 6(d)] is obtained at the expense of accuracy in predicting the spatial broadening of the PL signal at later times [Fig. 6(c)], which is mostly determined by the diffusivity *D*. The negative correlation between *A* and *B* can also be understood intuitively—since increasing either *A* or *B* results in faster signal decay, a greater value of *A* would result in a smaller value of *B* for the fit and vice versa. A *D* value near zero could also explain the unreasonably large *A* value found by objective function No. 2: Diffusion acts to lower the local carrier density, leading to a transient decay in the total photoluminescence signal; if *D* ≈ 0, then *A* must increase to account for signal decay.

We also applied the MCMC fitting algorithm to TPLM data experimentally obtained for CdS single crystals (supplementary material, Note 3). Similar to the MAPbBr_{3} analyses, unique determination of *D*_{MLE,No. 1} = 1.61 ± 0.03 cm^{2} s^{−1} was obtained when using objective function No. 1. This value is again comparable with literature reported values of ambipolar diffusivities obtained using light-induced transient grating^{54,55} or calculated from Hall mobilities.^{56,57} Interestingly, with objective function No. 2, unlike in the case for MAPbBr_{3}, *B* was undetermined for CdS. Given the experimental conditions used, it is likely that the contribution from bimolecular recombination is negligible relative to the other terms. The obtained estimates for the “determined” parameters were *A*_{MLE,No. 2} = (1.27 ± 0.01) × 10^{9} s^{−1} and *D*_{MLE,No. 2} = (1.20 ± 0.08) × 10^{−2} cm^{2} s^{−1}; the estimate for *A* is also larger than the literature reported value using single-beam Z-scanning,^{58} while *D* is underestimated by objective function No. 2. The analysis reveals a negative correlation between *A* and *D* for objective function No. 2—intuitively, both *A* and *D* contribute to decay of the measured PL signal; hence, a larger *A* would result in a lower *D* estimate and vice versa. A summary of the obtained maximum likelihood estimates of the parameters and a comparison against literature values are provided in Table II.

. | MAPbBr_{3}
. | CdS . | ||||
---|---|---|---|---|---|---|

Model parameters . | A (s^{−1})
. | B (cm^{3} s^{−1})
. | D (cm^{2} s^{−1})
. | A (s^{−1})
. | B (cm^{3} s^{−1})
. | D (cm^{2} s^{−1})
. |

Objective function | ||||||

no. 1 (this work) | ⋯ | ⋯ | 0.51 ± 0.02 | ⋯ | ⋯ | 1.61 ± 0.03 |

Light-induced | ⋯ | ⋯ | 0.45 ± 0.08^{51} | ⋯ | ⋯ | 1.6^{54} |

transient grating | 1.7 ± 0.2^{55} | |||||

Hall mobility | ||||||

measurements^{a} | $\u22480.26$^{52} | ⋯ | ⋯ | 1.75^{56,57} | ||

Objective function | ||||||

no. 2 (this work) | $3.42\xb10.01$ × 10^{8} | $8.00\xb10.02$ × 10^{−10} | ⋯ | $1.27\xb10.01$ × 10^{9} | ⋯ | $1.20\xb10.08$ × 10^{−2} |

Transient PL | $7\xb11$ × 10^{5}^{51}^{,}^{b} | $6\xb11$ × 10^{−10}^{51} | ⋯ | ⋯ | ⋯ | ⋯ |

Transient absorption | $3.7\xb10.2$ × 10^{4}^{53}^{,}^{b} | $4.9\xb10.2$ × 10^{−10}^{53} | ⋯ | ⋯ | ⋯ | ⋯ |

Single-beam | ||||||

Z-scan technique | ⋯ | ⋯ | ⋯ | $2.8\xb10.7$ × 10^{8}^{58}^{,}^{b} | ⋯ | ⋯ |

. | MAPbBr_{3}
. | CdS . | ||||
---|---|---|---|---|---|---|

Model parameters . | A (s^{−1})
. | B (cm^{3} s^{−1})
. | D (cm^{2} s^{−1})
. | A (s^{−1})
. | B (cm^{3} s^{−1})
. | D (cm^{2} s^{−1})
. |

Objective function | ||||||

no. 1 (this work) | ⋯ | ⋯ | 0.51 ± 0.02 | ⋯ | ⋯ | 1.61 ± 0.03 |

Light-induced | ⋯ | ⋯ | 0.45 ± 0.08^{51} | ⋯ | ⋯ | 1.6^{54} |

transient grating | 1.7 ± 0.2^{55} | |||||

Hall mobility | ||||||

measurements^{a} | $\u22480.26$^{52} | ⋯ | ⋯ | 1.75^{56,57} | ||

Objective function | ||||||

no. 2 (this work) | $3.42\xb10.01$ × 10^{8} | $8.00\xb10.02$ × 10^{−10} | ⋯ | $1.27\xb10.01$ × 10^{9} | ⋯ | $1.20\xb10.08$ × 10^{−2} |

Transient PL | $7\xb11$ × 10^{5}^{51}^{,}^{b} | $6\xb11$ × 10^{−10}^{51} | ⋯ | ⋯ | ⋯ | ⋯ |

Transient absorption | $3.7\xb10.2$ × 10^{4}^{53}^{,}^{b} | $4.9\xb10.2$ × 10^{−10}^{53} | ⋯ | ⋯ | ⋯ | ⋯ |

Single-beam | ||||||

Z-scan technique | ⋯ | ⋯ | ⋯ | $2.8\xb10.7$ × 10^{8}^{58}^{,}^{b} | ⋯ | ⋯ |

^{a}

The ambipolar diffusivity for CdS from Hall mobility measurements was obtained by first calculating the electron and hole diffusivities from the respective mobilities obtained through Hall mobility measurements using *D*_{e/h} = *μ*_{e/h}*k*_{B}*T*/*q*, where *μ* is the mobility, *k*_{B} is the Boltzmann constant, *T* is the temperature, and *q* is the charge of the electron. Then, they were inserted into $D=2DeDhDe+Dh$. For MAPbBr_{3}, the reported charge carrier mobility was used directly.

^{b}

The unimolecular recombination constant was obtained from recombination lifetimes *τ* using *A* = 1/*τ*.

## CONSIDERATIONS WHEN USING THE MODEL

It is important to note that compared to a Levenberg–Marquardt (LM) fitting technique, the MCMC sampling analysis identifies undetermined parameters where the LM fitting technique returns a single best fit value even if the parameters are truly undetermined for the specific objective function. To demonstrate, two built-in functions using the LM algorithm from MATLAB were used to determine the “best fit” parameters for MAPbBr_{3} and compared to the MCMC sampling method (see the supplementary material, Note 4). The LM fitting technique returned “best fit” values for all parameters; however, these values depended on what initial guess was provided to the solver, settling on local optima instead of finding the true global optimum parameter values. In contrast, the MCMC sampling method is robust to wide variation in initial guesses. Parameter variances and cross-correlations could also be obtained from the LM method by calculating the variance–covariance matrix; however, some of the results were unphysical (see the supplementary material). Nevertheless, the LM fitting algorithms present significant speed benefits over the MCMC fitting algorithm, and can provide useful estimates of material parameters if sufficient prior knowledge of the material exists to generate accurate initial guesses.

Given the strong dependence of the fitting results on the choice of objective function, alternative values for the estimated experimental error $\sigma ex,t$ were also considered (supplementary material, Note 5). To estimate experimental error, a combination function for $\sigma ex,t$ involving the sum of the measurement shot noise, $SPL,measured(x,t)$, which scales as the square root of the measured signal, and the detector noise signal (≈1 dark count per bin per acquisition period) was used such that $\sigma ex,t=SPL,measured(x,t)+1$. This places a greater weighting on early time data while also accounting for broadening effects. The resulting parameter analyses led to similar results as with objective function No. 2, with both *A* and *B* determined for MAPbBr_{3}, but *D* undetermined. Alternative scalings for $\sigma ex,t$ were also tested, the most interesting of which, $\sigma ex,t=SPL,measuredx,t0.75+1$, resulted in unique estimates for both *A* and *D*, but *B* remained undetermined for MAPbBr_{3}. While *B* has been shown to contribute to the broadening,^{26,30,32} it is likely that with the experimental conditions used for this sample, the contribution from *B* is less significant than that from *D*.

While unique estimates of *D* from TPLM data were possible, the parameter analyses suggest that *A* and *B* are over-specified if using objective function No. 1. Indeed, performing an order of magnitude analysis reveals that, based on the experimental and material conditions used, the *D*∇^{2}*n* term in the PDE dominates over the other two terms. In fact, for many material systems of interest, simplifications to the fitting model can be made by judiciously selecting a desirable *n*_{0} via the laser power used for the experiment.

Figure 8 shows that depending on the value of the ratio between the diffusivity and the bimolecular recombination terms $Dn0dx2Bn02$, it is possible to make simplifications to the model. For example, if the diffusion term dominates (over an order of magnitude), such as in MAPbBr_{3} and CdS single crystals, the model can be simplified to $\u2202n\u2202t=D\u22072n$ and the diffusivity may be obtained from TPLM data without having to consider the effects of recombination on the charge carrier density. This is consistent with the conclusions drawn by deQuilettes *et al.*^{32} On the other hand, if the diffusion term is negligible, a unique determination of the diffusivity from TPLM data will be challenging. For materials where the diffusion and bimolecular recombination terms are of similar magnitude, it is recommended to use the full model approach as detailed in this paper with the code provided in the supplementary material.

In the case of new material systems where relative estimates of the orders of magnitudes of *A*, *B*, and *D* are unknown, it is recommended to first use alternative techniques, such as thickness-dependent transient photoluminescence with a defocused laser beam,^{26,59,60} to obtain estimates of *A* and *B*, before using TPLM and the full model to estimate *D*.

## METHODS

### MAPbBr_{3} single crystal synthesis

Methylammonium lead bromide (MAPbBr_{3}) single crystals were synthesized following the method of Saidaminov *et al.*^{43} with slight modifications. 1M solution of MABr/PbBr_{2} (1:1 molar ratio) was prepared in *N*,*N*-dimethylformamide (DMF). The solution was filtered using a 0.2 *µ*m PTFE filter, and 2 ml of the filtered solution was placed in a 5 ml vial. The solution in the vial was immersed in an oil bath preheated to 50 °C. The temperature was increased at 10 °C/20 min until 100 °C and left for 2 h. The solution was discarded and the crystals were collected. Then, the crystals were dried under vacuum and stored under inert atmosphere.

### Characterization (PL, Abs, XRD) sample preparation

The characterization sample was prepared by crushing and grounding MAPbBr_{3} crystals between two glass substrates. Grounded crystals on a glass substrate were used for downstream characterization steps.

### Photoluminescence and absorption spectra

For the photoluminescence spectra, a 365 nm fiber-coupled LED (Thorlabs) was used to excite the grounded perovskite crystal sample on a glass substrate. PL spectrum was collected using an Avantes fiber-optic spectrometer. The absorption spectrum was collected by a Cary 5000 UV–vis spectrophotometer. Both were collected in air at room temperature.

### X-ray diffraction (XRD)

Powder XRD was performed using a PANalytical X’Pert Pro operating at 45 kV and 40 mA with a copper radiation source. HighScore software was used for background subtraction. To confirm the crystal structure, collected XRD pattern was compared with cubic MAPbBr_{3} XRD pattern in crystal structure database (ICSD 268785).

### CdS single crystal

Bulk CdS single crystal was purchased from MTI corporation (Richmond, CA). The ⟨0001⟩ oriented crystals were 5 × 5 mm^{2}, 0.5 mm thick, and single-side polished. TPLM was performed on the polished surface.

### Power-dependent photoluminescence

Power-dependent photoluminescence was measured by spectrally integrating the steady-state photoluminescence spectrum under varying excitation laser power. A 405 nm diode laser (LDH-D-C-405M, continuous wave mode) was used and its incident intensity was varied using neutral density filters (Thorlabs). The excitation beam was directed to an inverted microscope (Nikon, Ti-U Eclipse) and focused onto the surface of the crystals using an objective lens (Nikon, 40×, 0.6 numerical aperture). During experiments, all samples were mounted inside an evacuated microscope cryostat to protect from degradation due to air and water. The photoluminescence was collected in the epi configuration and filtered by a dichroic mirror and longpass filter. The emission was then directed into a spectrograph (Princeton Instruments, SP-2500), mounted with a cooled charge-coupled detector (Princeton Instruments, Pixis).

### Time-resolved photoluminescence microscopy

Time-resolved photoluminescence microscopy (TPLM) was performed using a homebuilt microscope. The samples were mounted on a piezo stage (attocube, ANC350) under vacuum inside a microscopy cryostat (Montana Instruments) to protect the samples from degradation due to air and water exposure. Free carriers were generated using a 405 nm pulsed laser source (PDL 800-D, <100 ps pulse width) that was focused down to a near diffraction-limited spot onto the surface of the crystals through an objective (Zeiss EC Epiplan-Neofluar 100×/0.85 NA). Epifluorescence from the sample was collected by the same objective lens and filtered by a dichroic mirror (Semrock, Di02-R405) and a longpass colored glass filter (Thorlabs, FGL435M). The emission then passed through a tube lens (Thorlabs, TTL200-S8) and a telescope (Thorlabs, AC254-030-A and AC254-125-A), resulting in a 495× magnified image plane after the telescope. An avalanche photodiode (APD, Micro Photon Devices, timing resolution ∼50 ps, 50 *µ*m sensor width) was raster scanned across the image plane in one dimension. The APD was mounted on stepper-motor stages (Thorlabs) for controlling the x-y position and connected to timing electronics (PicoQuant PicoHarp 300) such that time-correlated single-photon counting measurements could be taken at each position, providing spatiotemporal photoluminescence data.

## SUPPLEMENTARY MATERIAL

A package of executable MATLAB code and data files is included in the supplementary material, which also contains a comparison of solution speed for different finite differencing schemes; alternate presentations of fit quality; data and analysis results for CdS single crystal; a comparison of the MCMC fitting results to built-in MATLAB nonlinear fitting algorithms; and MCMC results using alternative objective functions.

## ACKNOWLEDGMENTS

This work was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, under Award No. DE-SC0019345. N.N.W. was partially supported by a MathWorks Engineering Fellowship. The opinions and views expressed in this publication are from the authors and not necessarily from MathWorks.

## AUTHOR DECLARATIONS

### Conflict of Interest

Yes, N.N.W. was partially supported by a fellowship from MathWorks, Inc., creator of the MATLAB programming and numerical computing platform used in this work. MathWorks did not review the manuscript or comment on this work prior to submission. The opinions and views expressed in this publication are those of the authors.

### Author Contributions

**Narumi Nagaya Wong**: Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Validation (lead); Writing – original draft (lead). **Seung Kyun Ha**: Data curation (supporting). **Kristopher Williams**: Conceptualization (supporting); Methodology (supporting); Software (supporting). **Wenbi Shcherbakov-Wu**: Data curation (supporting); Methodology (supporting); Supervision (supporting). **James W. Swan**: Formal analysis (supporting); Methodology (supporting); Software (supporting); Supervision (supporting). **William A. Tisdale**: Conceptualization (lead); Funding acquisition (lead); Project administration (lead); Supervision (lead); Writing – review & editing (lead).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material.

## REFERENCES

_{2}

_{2}/WSe

_{2}heterostructures

_{3}and MAPbBr

_{3}crystals measured under one- and two-photon excitations

_{3}NH

_{3}PbBr

_{3}and CH

_{3}NH

_{3}PbI

_{3}perovskite films: Influence of exciton binding energy

_{3}NH

_{3}PbI

_{3}perovskite