Markov Chain Monte Carlo Sampling for Target Analysis of Transient

Subscriber access provided by Uppsala universitetsbibliotek

Article

Markov Chain Monte Carlo Sampling for Target Analysis of Transient Absorption Spectra Matthew Nickol Ashner, and William A. Tisdale J. Phys. Chem. A, Just Accepted Manuscript • Publication Date (Web): 22 Mar 2019 Downloaded from http://pubs.acs.org on March 22, 2019

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Markov Chain Monte Carlo Sampling for Target Analysis of Transient Absorption Spectra Matthew N. Ashner and William A. Tisdale* Department of Chemical Engineering, Massachusetts Institute of Technology, 77 Massachusetts Ave., Cambridge, MA, 02139, USA *[email protected]

ABSTRACT Global and target analysis techniques are ubiquitous tools for interpreting transient absorption (TA) spectra. However, characterizing uncertainty in the kinetic parameters and component spectra derived from these fitting procedures can be challenging. Furthermore, fitting TA spectra of inorganic nanomaterials where the component spectra of different excited states are nearly or completely overlapped is particularly problematic. Here, we present a target analysis model for extracting excited state spectra and dynamics from TA data using a Markov chain Monte Carlo (MCMC) sampler to visualize and understand uncertainty in the model fits. We demonstrate the utility of this approach by extracting the overlapping component spectra and dynamics of singleand biexciton states in CsPbBr 3 nanocrystals. Significantly, refinement of the component spectra is accomplished by fitting the entire fluence-dependent series of ensemble TA data using the 1 ACS Paragon Plus Environment

The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Poisson statistics of photon absorption, providing multiple checks for internal consistency. The MCMC method itself is highly general and can be applied to any dataset or model framework to accurately characterize uncertainty in the fit and aid model selection when choosing between different models.

INTRODUCTION Transient absorption (TA) spectroscopy is a broadly used technique for characterizing excitedstate dynamics in a large variety of chemical systems. A small sampling of applications include understanding charge carrier transport in photosynthetic complexes,1,2 exciton multiplication processes in small organic molecules,3,4 and optical gain in semiconductor nanomaterials.5–7 Analyzing TA data often involves performing some kind of fitting routine that can vary in complexity depending on the characteristics of the spectra and information desired. In the simplest case, with a single or few well-separated spectral features, the data at selected wavelength points can be fit to one or more exponential decay models. However, TA spectra tend to be broad and have complex line shapes, so the spectra of different components often partially overlap. A common technique to disentangle the overlapping TA spectrum and simultaneously recover the spectra and dynamics of each component is to perform global or target analysis.8 Because the component spectra often have some distinguishing features in organic systems, the model results can be qualitatively verified against the raw data. Inorganic and solid state systems provide a challenge because their excited state characteristics are often encoded in subtle changes of the spectral line shape and complex dynamics. For example, the system considered in this work is the energetics and dynamics of biexciton states in CsPbBr 3 perovskite nanocrystals.9,10 Several studies have used TA spectroscopy to attempt to measure the 2 ACS Paragon Plus Environment

Page 2 of 29

Page 3 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


biexciton lifetimes and binding energies in these materials.7,11–13 The data analysis generally involves invoking a physical argument to propose a feature assignment to obtain the binding energies and fitting the intensity of the single bleach feature to some number of exponential decays to obtain the biexciton lifetimes. In contrast, a global or target analysis technique like that used for organic systems would allow interpretation of the data without having to propose a feature assignment a priori. However, given that the spectral features are nearly entirely overlapping with few to no independent parts, qualitative evaluation of the fit quality is much more difficult, complicating model selection and making it more difficult to evaluate if data have been overfit. Overcoming this challenge requires a method to characterize the uncertainty of the fit in a way that is straightforward to interpret. There are two common methods for performing a global fit, each with their own procedures for characterizing uncertainty. The first uses a linear least-squares algorithm to find the component spectra nested within a nonlinear least-squares optimization to fit a proposed kinetic model.8 Although many nonlinear optimizers will assign goodness-of-fit metrics, these are calculated from the Hessian of the objective function (matrix of second-order partial derivatives), or effectively its curvature. Due to the nested linear least-squares problem, the Hessian cannot be analytically evaluated. Numerical estimations of the Hessian depend on the step size at convergence, and may not accurately capture the uncertainty in the fit parameters if the objective function is nonmonotonic on the relevant scale. Furthermore, the solution yielded by a nonlinear optimizer can vary with the initial guess, especially if multiple disparate solutions give similar amounts of error. The second method employs a singular value decomposition (SVD) of the data yielding a set of component spectra and amplitudes that can be fit to a kinetic model if desired.14–16 SVD is powerful because it requires no assumption of a physical model, and the number of significant components 3 ACS Paragon Plus Environment


can be determined from the decomposition. It is also robust to noise in the data, and there are statistical methods for characterizing uncertainties in the spectra and amplitude vectors produced in the decomposition.8,17 However, SVD is a purely mathematical construction that is not bound by physical constraints. Consequently, SVD can sometimes yield non-physical results. For example, applying SVD to the data in this work produced amplitude profiles for major components that changed sign and did not decay to zero amplitude at long times. Furthermore, a multiple parameter fit to the amplitude vectors suffers from many of the same problems noted above, such as the lack of a reliable way to determine the uncertainty in extracted fit parameters. To overcome these challenges, we require a procedure to efficiently examine the parameter space of a proposed model and independently evaluate the uncertainty in those parameters. Here, we introduce an alternative TA data fitting protocol that can be used to guide model selection and intuitively display the uncertainty of the fits. The method is based on target analysis, where a kinetic model is proposed and fit using the nested optimization method described above. However, instead of using a nonlinear least-squares optimization algorithm as the outer layer, we utilize a Markov chain Monte Carlo (MCMC) algorithm to sample and visualize uncertainty in the parameter space. We demonstrate how the MCMC sampler results can be used to guide decisionmaking when fitting data to complex models, and show that this approach is successful in recovering the overlapping single- and bi-exciton TA spectra in colloidal CsPbBr 3 nanocrystal dispersions.

EXPERIMENTAL METHODS Sample preparation. Colloidal CsPbBr 3 nanocrystals were synthesized by previously reported methods.9 Briefly, lead oleate (5 mL, 2.5 mmol), cesium oleate (4 mL, 1.6 mmol), and 3-(N,N4 ACS Paragon Plus Environment

Page 4 of 29

Page 5 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


dimethyloctadecylammonio)-propanesulfonate (0.215 g, 0.5 mmol) were mixed with 50 mL 1octadecene and heated to 120°C under vacuum, where the atmosphere was changed to argon. The temperature was further elevated to 130°C, and trioctylphosphine bromide was injected (5 mL, 5 mmol of halides). The reaction was cooled down immediately by an ice bath. To purify the nanocrystals, 128 mL of ethyl acetate was added to the 64 mL crude solution, and the nanocrystals were precipitated by centrifugation at 29500g for 10 minutes. The nanocrystals were dispersed in 20 mL of toluene and precipitated with 40 mL of ethyl acetate and centrifugation at 29500g for 1 min. During the second and third purification step solvents were reduced by factor 2 for each step. They were finally dispersed in Toluene at a concentration of up to 100 mg/mL and centrifuged once more at 29500g for 1 min. Transient absorption spectroscopy. For the femtosecond transient absorption experiments, the 1040 nm output from a high repetition-rate Yb amplifier (Spectra-Physics Spirit) operating at 200 kHz was split into two beams. One arm was used to generate a 400 nm pump pulse by frequency doubling the signal from a commercial non-collinear optical parametric amplifier (Spectra-Physics Spirit-NOPA) in a BBO crystal. The other arm was used to generate a broadband probe pulse (480650 nm) by focusing the 1040 nm fundamental into a 4 mm YAG window with a 50 mm focal length lens. The pump pulse was compressed to ~50 fs FWHM using a fused silica prism compressor, sent through a mechanical delay stage (Newport) to control the pump-probe delay, and focused into a mechanical chopper (ThorLabs) to modulate the pump at 5 kHz. The pump beam fluence was controlled with a variable neutral density wheel. The pump and probe beams were focused by a 200 mm focal length concave mirror into the sample, which was a colloidal dispersion held and stirred in a 2 mm quartz cuvette (Spectrocell). The probe beam was coupled into a spectrometer with a high-speed data acquisition system (Ultrafast Systems) to collect spectra 5 ACS Paragon Plus Environment


at 10 kHz. The acquisition and chopping were synchronized with the amplifier repetition rate to ensure a consistent number of laser pulses per spectrometer acquisition. Each measurement was the average of 6 scans with 0.5 s acquisition per time point to ensure that dynamics did not change under laser exposure. The raw data were corrected for chirp in the probe pulse and photoluminescence background. The pump-probe cross correlation was measured using the optical Kerr effect in a 1 mm fused silica window. A wavelength-dependent time zero position was extracted from the cross-correlation and used to correct the data for timing variations due to the probe chirp. After that, the average of 10 spectra from negative time delays was subtracted from the data to remove contributions from spontaneous photoluminescence to yield the processed data.

RESULTS AND DISCUSSION Fluence-Dependent Transient Absorption Spectroscopy. A series of femtosecond TA measurements were performed on the CsPbBr 3 nanocrystal dispersions at varying pump laser powers. The result is a series of two-dimensional traces such as the one shown in Figure 1a, showing the pump-induced change in the absorption spectrum as a function of wavelength and pump-probe time delay. At low fluences, the trace will reflect the energetics and dynamics of the single exciton state. As the fluence is increased, the biexciton state also begins to contribute, with the relative populations of excitons and biexcitons contributing to the ensemble TA spectrum governed by the Poisson statistics of photon absorption. The effect of the biexciton state is captured as changes in the TA line shape at early times and multi-exponential dynamics of the overall signal amplitude, because the biexciton state decays faster than the exciton state. Figure 1b demonstrates the subtlety of the changes in the spectrum due to the biexciton, with the only noticeable changes being the shift from a small induced absorption to a small bleach on the low energy side of the 6 ACS Paragon Plus Environment

Page 6 of 29

Page 7 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


peak and a slight broadening on the high energy side. The changes in the dynamics are more dramatic, with a clear fast component emerging as the pump fluence is increased (Figure 1c). From these data, it is relatively straightforward to estimate dynamics for the biexciton state. A kinetic model can be fit to the data in Figure 1c,18 or an appropriately weighted version of the lowest fluence trace can be subtracted from trace at a higher fluence to display the biexciton dynamics directly.13 Extracting the spectrum and energetics of the biexciton state is significantly more challenging. It is clear from the slight fluence-dependent changes in the line shape that the spectra of the exciton and biexciton states are almost completely overlapped. To overcome this challenge, we require a global model that takes advantage of the full range of information in the dataset, linking the dynamics with the spectral changes. The goal will be to obtain a model for the TA spectral line shapes that captures the features noted above and a more accurate estimation of the dynamics of the biexciton state derived not just from the bleach amplitude dynamics, but also the dynamics of the changes in the spectral line shape.

Figure 1. Sample femtosecond transient absorption data. (a) Two-dimensional color plot displaying the full transient absorption trace at 200 µW of applied pump power. (b) Transient absorption spectra at a pump-probe time delay of 2 ps (the break between the linear and logarithmic scales in panel a) and 30 – 200 µW of applied pump power (2.8 – 18.7 µJ/cm2). (c) Transient absorption intensity at the bleach peak (502 nm, 2.47 eV, dotted white line in panel a) over the same range of pump fluences normalized to the same amplitude at long pump-probe time delays.

7 ACS Paragon Plus Environment


Page 8 of 29

Target Analysis with Markov Chain Monte Carlo Sampling. To fit the data, a model was constructed using the framework of target analysis to fit dynamics and species associated spectra. The goal of this method is to decompose the transient absorption data into the sum of some number of component spectra, si ( λ ) , each weighted by a timedependent intensity, ci ( t , J ) , which depends on the initial power of the pump beam J , ∆A(λ , t , J ) = ∑ ci (t , J )si (λ ) .

(1)

i

Eq. (1) assumes that the component spectrum si ( λ ) of each species in the kinetic model is invariant in time, and that the time dependence of the ensemble transient absorption spectrum can be completely explained by the population dynamics of the different excited state species contributing to the TA spectrum. The first step is choosing a kinetic model. For the dataset shown above, we used a 2-component model to fit the data at time delays >2 ps. The two components are the populations of biexcitons and single excitons in the ensemble. We ignore the effects of hot carrier thermalization since this process is completed faster than 2 ps at the fluences used.19 The kinetic model can be written as

dBX BX = − dt τ BX dX BX X = dt τ BX τ X

(2)

where the rate parameters are expressed as time constants for more intuitive interpretation later.

BX and X are the biexciton and exciton population, respectively. In an ensemble of nanocrystals, the populations of single excitons and biexcitons excited by the laser pulse is determined using the Poisson distribution with an expected value depending on the absorption cross-section multiplied by the photon flux, 8 ACS Paragon Plus Environment

Page 9 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(

)

= P N k= N

k

N e

− N

k!

.

(3)

We can use this relationship to initialize the exciton and biexciton populations and capture the fluence-dependence in the model with a single fit parameter (the absorption cross section). To maintain quantitative accuracy in the intensity, we initialize the biexciton population with the probability of 2 or more excitations per NC, where higher excitations will have small contributions that decay quickly in the fluence ranges used. We express the average excitation as a function of the average pump power, J, and a fit parameter, α, that lumps the absorption cross-section with the spot size, pulse repetition rate, and pump photon energy. The initial conditions are thus expressed as

1− P ( N = 0) − P ( N = 1) = 1 − e −α J (α J + 1) BX (0) = P ( N ≥ 2) = X (0)= P ( N= 1= ) α Je−α J

(4)

The complete kinetic model has 3 unknown parameters: the biexciton and exciton decay time constants (τ BX ,τ X ) and the absorption cross-section parameter (α ) .



Figure 2. Schematic demonstrating the target analysis fitting algorithm. In general terms, the algorithm can be considered to have two independent parts indicated by the large dashed boxes: an optimization algorithm that controls the parameter proposals, and a model evaluation that the optimization operates on. This schematic depicts the details of the target analysis model within the “model evaluation” section. To fit the full time-dependent TA spectral data, we must find the best-fit kinetic model parameters (τ X , τ BX , α) as well as the underlying spectrum for each of the components (X and BX). This dataset has 275 spectral points to give a total of 2x275+3=553 unknown values, far too many to fit efficiently with a nonlinear least-squares algorithm. The typical framework of target analysis overcomes this by taking advantage of the fact that the problem is linear with respect to each point in the component spectra, so the problem can be split into separate linear and nonlinear optimizations. The algorithm, illustrated in Figure 2, and can be considered to have two independent, modular parts illustrated by the dotted boxes. The target analysis framework serves as the model evaluation section, and the details of the framework are depicted in the schematic. 10 ACS Paragon Plus Environment

Page 10 of 29

Page 11 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Each iteration begins with the proposal of a set of kinetic parameters from the optimization algorithm, followed by the analytical solution and evaluation of the kinetic model at each time point and pump power, yielding the predicted component dynamics, c i . The next step is to find the component spectra, s i . To do this, we set up the problem by recasting Equation 1 in matrix form,

ΔA = C ⋅ S .

(5)

ΔA contains the TA spectra at each wavelength, time delay, and pump power, with the data at each pump power concatenated along the time delay dimension. C contains the component weights at each time point and pump power, concatenated as in ΔA, and S contains the component spectra. Inserting the component dynamics evaluated in the previous step as C, and the experimental data as ΔA sets up a linear system with the unknown spectra S. The best fit spectra for a given set of kinetic parameters can be evaluated efficiently using a matrix decomposition method such as the QR algorithm. The recovered spectra, S, are inserted into Equation 3 with the dynamics, C to calculate a fitted set of TA data, and the objective function for the fit is calculated from the sum of squared errors and returned to the optimization algorithm. Although this model makes fewer assumptions than simply fitting dynamics to certain features in the data, there are some assumptions that should be considered. First, it assumes that the exciton and biexciton states have a particular spectrum that does not change with time delay or fluence. Restricting the fit to time delays after hot carrier thermalization is complete is important for this assumption, because the hot carrier spectrum is continuously changing in time until thermalization is complete (~ 2ps for these CsPbBr 3 nanocrystals).19 Second, limiting the model to the biexciton and exciton states assumes that there is no contribution from N=3 or higher states. Limiting the fit to lower fluences where N is less than ~1 limits this contribution, and the excitation quantities



Page 12 of 29

are accounted by using the probability of N≥2 as the initial biexciton population. This is justified by the fast decay of N≥3 states. With the model in place, we must choose an optimization algorithm to propose parameters based on the model output (blue box in Figure 2). The typical choice is a nonlinear least-squares optimizer such as Levenberg-Marquardt or a trust-region algorithm. However, a key drawback of such algorithms applied to the target analysis model is that it can be difficult to know if one has “overfit” the data by including too many unnecessary fitting parameters or too many components in the kinetic model. Although the intensities of the spectra at each wavelength point are independent from each other, the model technically has 553 unknown parameters. It is important to understand how much uncertainty exists in the extracted spectra and model parameters in order to draw any trustworthy conclusions from the results of the model. Because traditional goodnessof-fit metrics are difficult to calculate in this system and may be inaccurate, we require a procedure to efficiently examine and visualize the parameter space and a different method for thinking about uncertainty in the fits. We do this by recasting the problem in the framework of Bayesian inference. We are interested in the conditional probability that a certain set of parameters, θ , describes the system given the observed data, x , referred to as the posterior probability. This can be calculated from Bayes’ Theorem:20

P (θ | x) =

P ( x | θ ) P(θ ) . P( x)

(6)

For this application we are only concerned with the relative probabilities of different sets of parameters given a single set of data, and not the absolute probability. So, the denominator can be ignored, and the equation can be rewritten as a proportionality statement, 12 ACS Paragon Plus Environment

Page 13 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


P (θ | x) ~ P ( x | θ ) P (θ ) ,

(7)

where the first term is the likelihood function, and the second term is the “prior”, or the probability of a parameter set being correct given any prior knowledge about the system before observing the data. In this application, we use an “uninformed prior,” where P(θ ) is a flat probability distribution within reasonable hard bounds and zero outside of those bounds, reflecting limited prior knowledge about the model parameters. In principle, observations from other measurements could be incorporated here as prior knowledge. To calculate the likelihood function, we assume that the probability of observing a given data point is a normal distribution around a true value given by the model. Taking the natural logarithm shows that the log-likelihood function is proportional to the sum of squared error with each term weighted by its variance,

P( x | θ ) ~ ∏ e

 ∆A ( t , λ , J ,θ ) −∆A ( t , λ , J )  −  σ ( t ,λ , J )  

2

t ,λ , J

.

(8)

ln [ P( x | θ ) ] ~ − ∑ ( ∆A(t , λ , J , θ ) − ∆A(t , λ , J ) ) σ (t , λ , J ) 2

−2

t ,λ , J

Since TA spectra are acquired by averaging multiple scans, the sample variance from the set of scans provides a good estimate for the variance at each data point, σ 2 (t , λ , J ) . This inclusion gives another key advantage over typical goodness-of-fit parameters as it explicitly propagates the known uncertainties in the measurement into the uncertainties of the estimated parameters. It also more heavily weights data points with higher certainty over those with lower certainty. Finally, since this method is only concerned with relative probabilities as the fit parameters are varied, the normalization can be neglected. With the derivation of the posterior probability complete, we next require a method to obtain it in a way that yields useful information about the parameter estimation. As with any probability distribution, a Markov chain Monte Carlo (MCMC) method is a suitable tool for this task. Briefly, 13 ACS Paragon Plus Environment


MCMC works by initializing a series of “walkers” in parameter space. Each walker represents a combination of possible model parameters that leads to a unique TA data set. To evaluate the parameter distribution, the walkers are randomly varied and the goodness of fit is evaluated for each walker. Random moves are accepted with some probability depending on the goodness of fit (thus adding it to the Markov chain), and the process is repeated until a suitably large Markov Chain is accumulated. Once the Markov chain has converged, the walkers can be sampled to compute a histogram representing the posterior distribution, as described in the following section. In the context of the larger algorithm in Figure 2, this amounts to using the MCMC sampling algorithm as the optimization algorithm on the left-hand side, proposing sets of parameters and keeping a record of each iteration. Additional information on the MCMC initialization and acceptance procedures can be found in the Supporting Information. Some consideration must be given to choosing an MCMC algorithm because different parameters in the model (and thus dimensions in parameter space) will have different scales which may not be known a priori. The algorithm used in this work is the stretch move with MetropolisHastings rule proposed by Goodman and Weare.21 This algorithm has a useful property known as affine invariance, which states that the performance of the algorithm does not depend on the relative scales or offsets of the different dimensions in parameter space. The key advantage of this algorithm is that it will run efficiently without any manual tuning of the step sizes in each dimension, unlike a simple random walk. Any affine invariant MCMC algorithm will be similarly effective. Executable MATLAB code for performing the MCMC target analysis algorithm described here is available in the online Supporting Information.


Page 14 of 29

Page 15 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3. Corner plot displaying projections of the posterior probability distribution sampled by the MCMC method applied to the model in Equation 2. The diagonal panels show histograms for each of the individual parameters. The off-diagonal panels show the cross-correlations between each pair of parameters produced by applying a kernel smoothing estimation algorithm to the MCMC samples. The contour lines are 10% intervals.

Table 1. Kinetic parameters fit from the target analysis model with 90% credible intervals extracted from the MCMC sampler. Parameter

τ BX τX α

Value 56.7156.77 56.65 ps 4, 9024,906 4,898 ps −1 0.0058360.005839 0.005833 µ W

Visualization of Model Outputs. To display the fit results, histograms are constructed from the samples in the Markov chain to give a representation of the posterior probability distribution, and 15 ACS Paragon Plus Environment


thus the uncertainty in the fitted parameters. These can be displayed in the form of a corner plot (Figure 3). The diagonal plots are histograms of each individual parameters, or the probability distribution projected onto each individual dimension. The results can also be tabulated as shown in Table 1 where the range given is the 90% credible region, or the range containing 90% of the MCMC samples. The off-diagonal plots display cross-correlations between each pair of parameters, or 2-dimensional projections of the probability distribution. The diagonal plots give the desired information about the breadth and shape of the uncertainty in each parameter. For example, if a particular parameter has an asymmetric uncertainty or bimodal distribution, it will appear in the corresponding diagonal plot. The cross-correlations provide the first of several important checks on the validity of the chosen model. If pairs of parameters show strong and unexpected cross-correlations, that can be an indication that the model has too many or redundant parameters and the data are being overfit. The off-diagonal panels of Figure 3 depict almost no cross-correlation between any of the parameters, indicating that the model parameters are independent from each other and there is no redundancy in the chosen model. Another way to visualize the fit quality is to take a random selection of states in the Markov chain, evaluate the model to obtain the fitted spectra for each one, and plot them on top of the data with some transparency. Figures 4a-c show the fits from 100 random points in the Markov chain (solid lines), plotted with the raw spectra at three pump-probe time delays (open circles). In this figure, varying certainty in the fit would be reflected as thickening and blurring of the solid line. The plotted fits appear solid and narrow on the relevant scale, and show a good, but not perfect, reproduction of the line shapes. Figure 4d shows a close-up view of a single data point with its error bars well outside of the field of view that shows this blurring. The extremely small scale of the variations between the 100 sample fits indicates a high confidence that the depicted fit is the 16 ACS Paragon Plus Environment

Page 16 of 29

Page 17 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


best possible fit of the proposed model to the data. Since a key goal of this data fitting procedure is to obtain the component spectra and their uncertainties, the same technique can be applied to those as shown in Figure 4e-f. Once again, the lines are narrow and sharp indicating very little uncertainty in the fit. The key use of these visualizations is that they easily differentiate random from non-random discrepancies between the model and the data. The thin, unblurred appearance of the fit lines and component spectra indicate that any discrepancies are not random, whereas blurred lines that encompass the data points would reflect random error. This allows us to confidently make some conclusions about our data. For example, the shift between the raw data and fits in Figure 2c indicate that there is a peak shift that the model cannot capture as formulated. There is also a jagged region on the blue end of the raw data that appears to be “noise,” but it is reproduced well in the component spectra in Figure 4e-f. Thus, this feature is “real” in the sense that it is accurately reproducing the data. However, since we expect this to be unphysical, we can attribute it to an experimental artifact.



Figure 4. Results of the global fit using the model in Equation 2 and the MCMC sampler. (a)-(c) Comparison of the fitted and raw data. The colored circles are raw TA spectra at a pump-probe time delay of 2, 200, and 2000 ps respectively, and the colors indicate differing pump powers from low (red) to high (blue). The black lines are the fitted spectra extracted from 100 random states of the Markov chain overlaid with some transparency. (d) Close-up view of a data point from panel a demonstrating the blurring of the fit line. The data point is not on the same scale, and its error bar is 1.17x10-5 dOD, well outside of the field of view. (e) Component spectra of the biexciton (green) and exciton (blue) states extracted from 100 random states of the Markov chain overlaid with some transparency. (f) Component spectra normalized to have the same bleach amplitude for visual comparison.

Application of MCMC to Model Selection. In addition to providing accurate characterizations and useful visualizations of the fit uncertainty, the results can provide a wealth of information to aid model selection. For example the bleach peak position at 2 ps and 200 ps (Fig. 4a,b) is wellrepresented by the model fit, but the model and data diverge at a much longer time delay of 2000 ps, as shown in Fig. 4c. This indicates a dynamic change in the line shape on a ns time scale that is not captured in the proposed model. In contrast to plotting the results from a nonlinear least18 ACS Paragon Plus Environment

Page 18 of 29

Page 19 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


squares fit, Figures 4a-c straightforwardly differentiate whether a deficient fit is due to statistical uncertainty or the choice of an inadequate model because they contain a visualization of the model uncertainty. In this case, the thin black lines indicate high certainty in the best fit of the proposed model to the data, so a different model must be proposed to accurately model this shift. To demonstrate the ability of the MCMC algorithm to aid model selection, we propose an alternate model in an attempt to capture this long-time spectral shift. Since nanocrystals can become charged – resulting in trion formation upon single photon excitation – we might consider that there is a second single exciton component with a different spectrum (i.e. the trion, X*) that occurs in some fraction of the population in the system, φ . To test this hypothesis, we run the simulation with a new kinetic model proposed,

dBX BX = − dt τ BX dX BX X = (1 − φ ) dt τ BX τ X dX * BX X * =φ dt τ BX τ X * BX (0) =P ( N ≥ 2 ) =− 1 e −α J (α J + 1) X (0) = (1 − φ ) P ( N = 1) = (1 − φ )α Je −α J X * (0)= φ P ( N= 1= ) φα Je−α J


.

(9)


Figure 5. Example of an over-specified model. Corner plot displaying projections of the posterior probability distribution sampled by the MCMC method applied to the model proposed by Equation 7. The diagonal panels show histograms for each of the individual parameters. The off-diagonal panels show the cross-correlations between each pair of variables produced by applying a kernel smoothing estimation algorithm to the MCMC samples. The contour lines are 10% intervals. The flat posterior distribution of the parameter ϕ indicates an over-specified model.

The corner plot from the fit to Eq. 7 is shown in Figure 5, and the corresponding spectra are shown in Figure 6. Immediately, we recognize problems with the newly proposed model. The histogram in the bottom-right panel of Figure 5 shows that the yield parameter φ can vary freely between 0 and 1, and thus it carries no useful information. Noting that the amplitudes of the sampled exciton component spectra (Figure 6d) vary wildly, while the biexciton spectrum (Figure 6d, inset) does not leads to a hypothesis as to why this occurs: the intensities of the component spectra extracted from the linear least-squares compensate for the different absolute populations 20 ACS Paragon Plus Environment

Page 20 of 29

Page 21 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of each component as the yield varies. It is straightforward to investigate this using the visualizations of the MCMC sampler results. Figures 6d-e provide circumstantial evidence for this by showing that the sampled component spectra show high certainty in their shape when normalized. Figure 7 displays a corner plot comparing the intensities of the component spectra with the split fraction and clearly showing a strong cross-correlation, directly verifying this hypothesis and confirming that the yield fraction parameter should not be used. Furthermore, comparing fits from this model shown in Figure 6a-c to the fits from the 2-component model shown in Figure 4a-c reveals that adding the third component makes only a marginal improvement in the fit. It subtly improves the fit to the line shape at early times, but still does not capture the more obvious transient shift on the ns time scale. Thus, the increased complexity of the 3-component model is not justified by a corresponding increase in fit quality, so this model should not be selected.

Figure 6. Example of an over-specified model. Results of the global fit using the model in Equation 7 and the MCMC sampler. (a)-(c) Comparison of the fitted and raw data. The colored 21 ACS Paragon Plus Environment


circles are raw TA spectra at a pump-probe time delay of 2, 200, and 2000 ps respectively, and the colors indicate differing pump powers from low (red) to high (blue). The black lines are the fitted spectra extracted from 100 random states of the Markov chain overlaid with some transparency. (d) Component spectra of the two exciton states, X (blue) and X* (red), extracted from 100 random states of the Markov chain overlaid with some transparency. (inset) Component spectra for the biexciton state. (e) Component spectra normalized to have the same bleach amplitude for visual comparison. (f) Normalized component spectra zoomed-in for better comparison of the X and BX line shapes.

Figure 7. Corner plot comparing the trion fraction parameter ϕ to the minimum intensity of the X and X* component spectra. The off-diagonal panels show little random spread and a very strong correlation between the yield parameter and both intensities as well as a strong correlation between the 2 intensities.

Self-Consistency of Target Analysis Model. While the MCMC sampling gives a strong indicator of how much to trust the global model, the target analysis model itself as applied to semiconductor nanocrystals includes other checks for self-consistency and agreement with other known material properties in the interpretation. Notably, it is highly confidence-inspiring that the early-time line shapes of the fluence-dependent spectra can be reproduced using a model where the relative biexciton and exciton populations are determined by Poisson statistics with a single fitting parameter (the absorption cross section). As mentioned above, the fitting parameter α can 22 ACS Paragon Plus Environment

Page 22 of 29

Page 23 of 29 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


be converted to an estimate of the true absorption cross section. Using the measured laser beam diameter from our experiment, we retrieve a value of 3.1x10-14 cm2, in excellent agreement with values reported in the literature (

Markov Chain Monte Carlo Sampling for Target Analysis of Transient

Recommend Documents