A Python toolbox for unbiased statistical analysis of fluorescence intermittency of multi-level emitters
AA Python toolbox for unbiased statistical analysisof fluorescence intermittency of multi-levelemitters
Isabelle M. Palstra † , ‡ and A. Femius Koenderink ∗ , ‡ † Institute of Physics, University of Amsterdam, Science Park 904, 1098 XH Amsterdam,The Netherlands ‡ Center for Nanophotonics, AMOLF, Science Park 104, 1098 XG, Amsterdam, TheNetherlands
E-mail: [email protected]
Abstract
We report on a Python-toolbox for unbiased statistical analysis of fluorescence inter-mittency properties of single emitters. Intermittency, i.e., step-wise temporal variationsin the instantaneous emission intensity and fluorescence decay rate properties are com-mon to organic fluorophores, II-VI quantum dots and perovskite quantum dots alike.Unbiased statistical analysis of intermittency switching time distributions, involvedlevels and lifetimes is important to avoid interpretation artefacts. This work providesan implementation of Bayesian changepoint analysis and level clustering applicable totime-tagged single-photon detection data of single emitters that can be applied to realexperimental data and as tool to verify the ramifications of hypothesized mechanisticintermittency models. We provide a detailed Monte Carlo analysis to illustrate thesestatistics tools, and to benchmark the extent to which conclusions can be drawn on the a r X i v : . [ phy s i c s . c o m p - ph ] F e b hotophysics of highly complex systems, such as perovskite quantum dots that switchbetween a plethora of states instead of just two. Introduction
Since the seminal first observation of single molecule emitters in fluorescence microscopythree decades ago, single quantum emitter photophysics has taken center stage in a largebody of research. On one hand, single quantum emitters as single photon sources areheld to be an essential part of quantum communication networks and are deemed essentialfor building optically addressed and cavity-QED based quantum computing nodes. Thishas particularly spurred research in III-V semiconductor quantum dots, color centers indiamond, silicon carbide and 2D materials, and organic molecules at low temperature. Onthe other hand, classical applications of ensembles of emitters for displays, lighting, lasersand as microscopy-tags drive the continuous development of new types of emitters, suchas II-VI self-assembled quantum dots 20 years ago, and inorganic perovskite quantumdots just recently.
For all these systems, understanding the photophysics on the singleemitter level is instrumental, whether the intended use is at the single or ensemble level.A common challenge for almost all types of emitters is that they exhibit intermittency,also known as blinking.
Under constant pumping emitters switch, seemingly at random,between brighter and dimmer states, often corresponding with higher and low quantum yield(QY), and different fluorescence decay rates. Frequently the switching behavior also showspeculiar, power-law distributed, random distributions of durations of events. Determiningthe mechanism through which emitters blink, i.e., the origin of the involved states, thepower-law distribution of residence times, and the cause of switching, have been the topic ofa large number of studies particularly for II-VI quantum dots as recently reviewed by Efrosand Nesbitt. Recent studies on inorganic perovskite quantum dots uncover intermittencybehavior that does not fit common models for intermittency in their II-VI counterparts. ime stampsCPAGrouping Figure 1: A schematic overview of the working of the toolbox. (A) An illustration of the twomethods available to obtain TCSPC data, either by photolumiescence TCSPC measurementsof a single emitter (top panel) or through simulation of a single emitter (bottom panel). Thelatter is provided in the toolbox. Both will result in a stream of timestamps that can then befurther analyzed by the toolkit. (B) Starting out with this stream of (photon) events, usingCPA (C), the changepoints are found, and with these, the instantaneous intensities. (D)Subsequently, these events are grouped in order to find the most likely underlying intensitylevels between the behavior. Following this, a number of analyses can be done, such as(E) ascertaining whether the time between switching events is powerlaw distributed, (F)the visible separation of states in FDID diagrams, and (G) the presence of memory in theswitching behavior. Here we show an example of a simulated, 4 level emitter, with simulationparameters chosen for clarity of illustration. 3n order to quantify intermittent behavior, the simplest and most commonly employedmethod is to subdivide a measurement stream of individual photon-arrival times into shortbins of a few ms, to calculate the intensity (in counts/second) of each bin. Every bin canthen be assigned to a state (on, off or grey) according to its brightness so that on/off timesas well as intensity levels can be defined and analyzed. For pulsed laser excitation, alsoquasi-instantaneous fluorescence decay rates can be obtained.
However, it is well knownthat this method of binning time streams and histogramming binned intensities causes detri-mental artefacts.
Retrieved parameters of the quantum dot behavior often exhibit adependency on the choice of the bin width, which affect estimates of switching time distri-butions and power laws,and also the objective assignment of intensities to intrinsic levels.Narrower bin widths in principle allow better resolution, but run into shot noise limits, whileconversely choosing larger bins suppresses noise, but will render the analysis blind to fastevents. To overcome these issues, Watkins and Yang proposed changepoint analysis (CPA)as a Bayesian statistics approach for the unbiased determination of switching times that isoptimal, i.e., gives the best performance given the constraints of shot noise in the data.With CPA, the data is partitioned in time segments between jumps. Through this, one hasaccess to the best estimate of intensity, decay rates, and clustering of measured intensitiesin distinct levels (states). While such Bayesian statistics methods are widespread in, e.g.,the context of stock market changes and twitter feed analysis, and despite the reported su-perior performance, in the domain of single photon counting data only very few groups haveadopted it. While II-VI quantum dots for which CPA was orginally developed,are generally understood to switch between just two or three states, the problem of accurateanalysis of intermittency is gaining in prominence with the advent of novel emitters, suchas perovskite quantum dots which appear to switch not between just two, but instead amultitude of states.
In this paper, we provide, benchmark, and document a Python toolbox for changepointanalysis, state clustering, and analysis of fluorescence-intensity-decay rate correlations (pro-4ided as supporting information, and posted as a github repository ). This toolbox is bothapplicable to real quantum dot data, and valuable as a testbed for both testing models andanalysis techniques on synthetic, i.e., numerically generated data. Indeed, the code toolboxwe supply includes code to generate numerically random photon time arrival data streamsfor ‘synthetic’ quantum dots that jump between an arbitrary set of intensity levels and decayrates, with jump time statistics and photon budgets that can be set by the user. As results weprovide benchmarks on the performance of CPA for detecting change points in function of thenumber of intensity levels and total photon budget, and we explore the limits to the numberof distinct states that the clustering analysis can reliably separate. Moreover, the toolboxallows to us test the accuracy of jump time statistics, such as power law statistics, for suchmultilevel dots. Finally, the test suite also allows to benchmark the accuracy of fluorescencedecay model fitting with maximum likelihood estimation, and we discuss the constructionof fluorescence-intensity-decay rate correlations from CPA-partitioned data. This paper isstructured as follows. In the Methods section we summarize the Bayesian statistics toolswe implemented to analyze all aspects of our data. Next, we benchmark the performanceof changepoint analysis to pinpoint intensity jumps, and of level clustering to identify thenumber of levels between which a dot switches on the basis of Monte Carlo simulations. Nextwe present considerations on the dependence of on-off time distributions, decay rate fit andso-called ‘fluorescence decay-rate versus intensty diagrams’(FDIDs) on count rates. Methods
In this section we present all the methods implemented in our Python toolbox, as well asthe methods for benchmarking them. Benchmark results are presented in the results section.We refer to the supporting information for a manual to the code and the code itself.5 hangepoint analysis
First, we summarize changepoint analysis (or CPA), a Bayesian statistics method for the un-biased determination of jumps or ‘changepoints’ in time traces of discrete events.
Bayesian statistics is a paradigm that reverses the usual standpoint of probability theory.Usual probability theory views a data set as a random draw from a probability distribution,given a hypothesis on the parameters of the underlying physical process. In this framework,one can calculate the likelihood of drawing the specific measured data set. Bayesian statis-tics, on the other hand, compares the likelihood of distinct hypotheses, given a measureddata set and assumptions on the underlying measurement noise.We consider time-tagged single photon counting data consisting of an ordered list ofmeasured photon arrival times s k , collected over a measurement time T . For a single emitterwith no memory that emits at a count rate of N photons in a time T , the waiting times- i.e., the times between photon arrivals - are exponentially distributed with waiting time τ w = T /N . In order to determine whether there is a changepoint in some segment q , CPAcompares the likelihood of two distinct hypotheses, (1) there is a jump in emission intensity(i.e. the average waiting time τ w jumping from some value to another) against (2) there isthe same intensity throughout the measurement interval. When testing for a jump at photondetection event k at time s k in this trajectory q with time duration T q containing N q photonevents, this leads to a log-likelihood ratio, or ‘Bayes factor’ L k = 2 k ln kV k + 2( N q − k ) ln N q − k − V k − N q ln N q , where V k = s k /T q . Derivation of this log-likelihood ratio involves several steps. First,it incorporates the assumption that in between jumps, the waiting time between photonsis exponentially distributed, on basis of which one can assess the likelihood of measuringthe given data set for a given hypothesis on the exponential waiting time τ w . Second, ituses maximally non-informative priors for L k to compare the hypothesis of presence versus6bsence of a changepoint without further restrictive assumptions on the involved intensitylevels.It should be noted that there are other ways to arrive at the same log-likelihood ratiotest. One alternative starting point is a binary time series in which there is an underlyinguniform and small probability distribution of photon detection per bin (e.g., imagining thetime axis binned in by the timing card resolution (of order 0.1 ns for typical hardware)). Such a uniform distribution would emerge as a direct consequence of exponential waitingtime distributions. In this case one should start from a binomial distribution and ultimatelyarrives at the same formula after application of Stirling’s formula. Another starting point isCPA applied to binned data with wider bins with multiple counts, i.e. to series of Poissondistributed intensities instead of discrete events.
However, the binning would introducean undesirable time scale through the chosen bin width. Of these three methods, workingwith photon arrival times is the most data efficient approach and introduces no artificialpartitioning whatsoever. We refer to Ref. for a derivation of the log-likelihood ratio inall these three scenarios, which includes a precise description of the use of maximally non-informative priors.Following Watkins and Yang and Ensign, the most likely location of a changepoint,if any, is at the k that maximizes the Bayes factor L k . The hypothesis that this most likelychangepoint is indeed a real event is accepted if L k exceeds a critical threshold value for L k or‘skepticism’. This value is chosen to balance false positives against missed events. A full dataset is partitioned recursively, i.e., by recursively checking if data sets between two acceptedchangepoints themselves contain further changepoints. This results in a division of the dataset into segments, each of which starts and ends at an accepted changepoint, and with thelevel of skepticism as stop criterion for the recursion. The resulting segmentation providesthe most likely description of data as consisting of segments within which the intensity isconstant, given the value chosen for the degree of ‘skepticism’, and given the amount ofdata collected. Since the algorithm works with the list of individual photon arrival times7his segmentation entails no arbitrary partitioning. Following Ensign, we use a value ofskepticism of 8. Clustering
Changepoint analysis splits the data in segments separated by jumps (a list of Q jumpsdelineate Q − m r actually underlie the N − I . . . I Q − . Toanswer this question Watkins and Yang proposed a clustering approach. The recent workof Li and Yang provides a detailed explanation of the reasoning involved, though quotingresults for Gaussian instead of Poissonian distributed data. The idea is that with the Q − I q one can use expectationmaximimization to calculate, for a hypothesized and fixed number of levels n G what the mostlikely underlying intensity levels I m are (with m ∈ . . . n G ), and how probable it is thateach segment is ascribed to a given level (probability p mq ). Subsequently Bayesian inferenceis used to establish what the most likely number of levels (i.e. states) m r and associatedintensities I m , with m ∈ { . . . m r } is that describes the data.Following Ref., the expectation minimization in our toolbox is implemented as aninterative algorithm started by a first guess of the segmentation. This guess is obtained by ahierarchical clustering of Q − m = 1 , , . . . Q − m = 1 , , . . . Q − p mq for segment q to belong to the m th level, as well as an estimate of the intensities of these levels I m . In each iteration the8ntensities of all levels are estimated from the level assignment from p m,q . Following this, theprobability distribution p mq is updated to redistribute the segments over the levels. In thiscalculation it is important to understand the type of noise statistics the data obeys. In thecase of single-photon measurements, and for the purpose of this discussion, the intensities arePoisson distributed. The iteration is repeated until p mq converges (practically also cappedby a maximum number of iterations). The final outcome is a most likely assignment ofthe measured segments into n G levels. Next for each value of nG one asseses the ‘Bayesianinformation criterion’ (BIC). This criterion is a measure for how good the description ofthe segmented intensity trace is with n G intensity levels given the assumption of Poissoncounting statistics for each fixed intensity level. Beyond a mere ‘goodness of fit’ metricthat would simply improve with improved number of parameters available to describe thedata, this metric is penalized for the number of parameters to avoid overfitting. For Poissondistributed data the criterion is derived in asBIC = 2 L EM − (2 n G −
1) ln Q − Q ln N where Q again is the number of change points detected, n G is the number of available levels.The term L EM is the log-likelihood function optimized in the expectation maximizationstep, i.e L EM = (cid:80) q (cid:80) n G m =1 p mq ln[ p m P ( I q ; I m )] with P ( x ; λ ) the Poisson probability functionat mean λ , p m the probability of drawing level m . The second term in the BIC s the termpenalizing the BIC for overfitting. The accepted best description of an emitter in n G -levelsis taken to be at the value of n G where the BIC peaks. Intensity cross/autocorrelation and maximum likelihood lifetimefitting
Many single photon counting experiments are set up with pulsed laser excitation for flu-orescence decay rate measurements, and with multiple detectors to collect intensity auto-9orrelations (e.g., to verify antibunching in g (2) ( τ ) for time intervals τ comparable to thefluorescence decay rate, and shorter than the commonly longer detector dead time). In atypical absolute-time tagging set up, this results in multiple data streams S A , S B , S R of timestamps corresponding to the detection events on each detector, and the concomitant laserpulses that created them, respectively. Our Python toolbox contains an implementation ofthe correlation algorithm of Wahl et al. that operates on timestamp series, and returns,for any combination of channels S , S (1 , ∈ { A,B,R } ), the cross correlation C ( τ )∆ τ , i.e.,the number of events in the time series S and the time series S that coincided when shiftedover τ , within a precision ∆ τ .Cross-correlating detected photons and laser arrival times, taking ∆ τ to be the binningprecision of the counting electronics and the range of τ equal to the laser pulse repetitionrate, returns a histogram of the delay times between photon detection events and laserpulses. To obtain g (2) ( τ ) to investigate antibunching, streams of photon events from twodetectors in a Hanbury-Brown Twiss set up are cross-correlated. ∆ τ is taken to be thebinning precision of the counting electronics and the sampled range of τ as an interval istaken symmetrically around τ = 0 and several times the laser pulse interval. Finally auto-or cross correlating detector streams over τ -ranges from nanoseconds to seconds, coarseningboth τ and ∆ τ to obtain equidistant sampling on a logarithmic time axis, results in longtime intensity autocorrelations of use in intermittency analysis. Our toolbox also providesthis logarithmic time-step coarsening version of the correlation algorithm of Wahl et al. Of particular interest for intermittent single emitters is the analysis of fluorescence decayrates in short segments of data as identified by changepoint analysis, that may be so short asto contain only 20 to 1000 photons. For each of the photon detection events in a single CPAsegment, cross-correlation with the laser pulse train yields a histogram of the N q photonsin segment q . In each of the bins (with width ∆ τ ) the photon counts are expected to bePoisson distributed. Therefore the optimum fit procedure to extract decay rates employsthe Maximum Likelihood Estimate procedure for Poisson distributed data, as described by10ajzer et al. In brief, for a decay trace sampled at time points τ i relative to the laserexcitation, with counts per bin D ( τ i ), the merit function reads M = − (cid:88) all data bins i { D ( τ i ) log [ F A ( τ i )] − F A ( τ i ) } . (1)Assuming a chosen fit function F A ( t ) the parameter set A that minimizes this merit functionprovides the parameter values that most likely correspond to the data. The estimated errorsin these parameters then follows from the diagonal elements of the inverse of the Hessian of M relative to the parameters A . Importantly, the fact that the Poisson distribution is tiedto absolute numbers of counts, implies that this approach requires that the data is neither scaled nor background subtracted. Instead, the background should be part of the fit functioneither as a free parameter or a known constant. Furthermore, it should be noted that timebins with zero counts are as informative to the fit as non-empty ones, and should not be leftout. Generating synthetic quantum dot data
To benchmark the CPA and clustering method and to test its limits, our toolbox providesan example routine to generate artificial data mimicking quantum dot intermittency. Toobtain mimicked quantum dot data, we first choose a number m of intensity levels I ,m between which we assume the dot to switch. Next we generate switching times for eachof the states. In this work, we choose all switching times from a power law distribution. For benchmark purposes we will present results with power law exponent α = 1 .
5, thoughany exponent can be set in the code. On the assumption that intensity levels appear in arandom and uncorrelated order, this segments the time axis in a list of switching events T ,j , j = 1 , , . . . m , where for each segment we randomly assign one of the nominal intensities I ,m . Next, to mimic a pulsed excitation experiment, we imagine each of these segmentsto be subdivided in intervals of length τ L equivalent to a laser repetition rate ( τ L = 100 ns11n the examples in this work). We assign each of these intervals to be populated with onephoton at probability p m = I ,m τ L . This ensures that the number of photons N q in everysegment is drawn from a Poisson distribution at mean T ,q I ,m . By removal of all emptybins, the binary list is translated into a list S = ( t , t , t , . . . ) of photon arrival time stampsat resolution τ L to which one can directly apply CPA to attempt a retrieval of switchingtimes, and apply clustering to retrieve the number of states.To also enable fluorescence decay rate analysis, we further refine the photon arrival timelist. Recalling that we have generated switching events T ,q between intensity states I ,m ,we now also assume fluorescence decay rates γ ,m . As each segment q , was already chosen tocorrespond to some level m q , we now impose decay rate γ ,m q on the photon arrival times. Todo so, for each of the photon events k = 1 . . . N q already generated at resolution τ L we nowrandomly draw a delay time ∆ k relative to its exciting laser from an exponential distributioncharacterized by rate γ ,m q . To mimic the behavior of typical TPSPC counting equipment,the delay time is discretized at a finite time resolution ∆ τ (in this work chosen as 165 ps tomatch the hardware in the provided example experimental data measured in our lab (Becker& Hickl DPC-230), though of course in the toolbox the value can be set to match that ofany TCSPC card vendor). Testing of fluorescence decay trace fitting can operate directly onthe generated list of delay times, or alternatively one can synthesize a TCSPC experimentby re-assigning S to represent laser-pulse arrival times S R , and defining photon arrival timesas the events in S A , S B each shifted by its delay time, i.e. S X = ( t + ∆ , t + ∆ , . . . ),withX ∈ A,B. Cross correlation of S R and S X returns the delay time list. We note that althoughour work does not focus on antibunching, our quantum dot simulation routine providesdata distributed over two detector channels, where emission events antibunch, while anuncorrelated background noise level of the detectors can also be set.12 ractical implementation We have implemented the toolbox ingredients in Python 3.8. As timestamp data can besubstantial in size, we use the ‘parquet’ binary format to store timestamps as 64-bit inte-gers. Processing and plotting the data is dependent on Pythons’ standard libraries numpy,matplotlib, while we use Numba, a just-in-time compiler, to accelerate the time-stamp cor-relation algorithms. An example script to generate syntehtic data, and to run the entireworkflow on simulated data is provided. We refer to the supporting information for a guideto the practical implementation and use of our toolbox. The toolbox comes also with exampleexperimental data on single CsPbBr quantum dots from a recent experiment. Results
The remainder of this work is devoted to presenting benchmarks of the provided methods.Benchmarks for emitters with ‘binary’ switching, i.e., two well-separated intensity levels asis typical for II-VI quantum dots, have already been presented in literature.
However,emitters under current study, such as perovskite quantum dots appear to have a multitude,or perhaps even a continuum, of intensity levels. Our tests hence focus on determining theperformance of CPA and level clustering for many-level single photon emitters.
Precision of identifying individual changepoints
Figure 2(A) and (B) show examples of CPA analysis applied to a simulated quantum dotwith a single jump in its behavior, from an intensity level 10 to 10 cts/s resp. from 4 . × to 2 . × , with ∼
900 and ∼
150 photons left and right of the changepoint, respectively.Purely for visualization purposes, the data is plotted in a binned format, as the analysisitself does not make use of any binning. Alongside the binned intensity trace, we also showthe log-likelihood ratio L k . In both cases, the log-likelihood ratio clearly peaks at or closeto the point where there is a changepoint in the data. Since the Bayes factor is actually a13igure 2: Demonstration of changepoint detection applied to a synthesized data set witha single changepoint, with equal photon counts before and after the changepoint. In (A)and (B) the contrast between intensities is a factor of 10 and 2, respectively, while thetotal photon budget is approximately 2000 and 300. The bottom panels show the log-likelihood ratio test, which clearly peaks at the changepoint in both cases. The robustnessof the method is demonstrated in (C) and (D), where we show the likelihood of detecting achangepoint in such a series for different intensity contrasts, and the variance of the foundtimes, respectively, as a function of the total photon numbers. To gather accurate statistics,10 photon traces were generated for each data point. The data is plotted in a binned fashion(ms bins) for visualization purposes only. 14 ogarithmic measure for the comparison of hypotheses, the algorithm indeed identifies thechangepoint with high probability and to within just a few photon events, even where thejump is far smaller than the shot noise in the binned representation in the plot, at a relativeintensity contrast of just a factor of 2. Generally, the probability with which the algorithmidentifies or misses the changepoint is dependent on the total number of photons recordedboth before and after the changepoint, and on the contrast in intensities, consistent with thefindings of Watkins, and Ensign. To identify the limits of CPA we consider the feasibility of identifying changepointsof contrast I /I as function of the total number of photons in the time record. The resultsare shown in Figure 2C for the likelihood of detecting a changepoint, and Figure 2D for theerror in identifying the precise event k at which the changepoint that is identified occurred.Here, we only consider the case where there are roughly an equal number of photon eventsbefore and after the changepoint. This data is obtained by simulating 10 switching eventsof the type as shown in Fig 2(A,B) for each contrast and mean photon count shown. Athigh intensity contrast, exceeding a factor 5, a total photon count as low as 300 is enoughfor near-unity detection. Moreover, for sufficiently high photon count left and right of thechangepoint, even very small changes in intensity have a high likelihood of being accuratelydetected, even if in binned data representations the jump is not visible within the shot noise.Figure 2(D) provides a metric for the accuracy to within which changepoints are pinpointed.Changepoint analysis returns the most likely photon event k in which the jump occurred,which in our analysis can be compared to the actual photon event index k at which weprogrammed the Monte Carlo simulation to show a jump.Figure 2d reports the mean error( (cid:112) var( k − k ) as a metric of accuracy. At jump contrasts above a factor 2, changepoints areidentified to within an accuracy of almost one photon event even with just 10 photon countsin the total event record. At very small contrasts, the error in determining the location of achangepoint is generally on the level of one or two photon events, only worsening when thereare fewer than 200 counts. This observation highlights the fact that if the photon record has15ust a few counts in total, the error in estimating the count rate before and after the jumpbecomes comparable to the magnitude of the jump. Intermittency and on-off time histograms
Figure 3: (A) Typical time trace of a simulated quantum dot. The intensity duty cycleswitches between 0 . × and 2 × counts/s. It shows an on/off input duty cycle generatedwith a power-law distribution (orange), the duty cycle with Poissonian noise (purple), and theretrieved duty cycle (green). Overall, the original intensities and lifetimes (B) are retrievedwell. (C) A histogram of the number of switching events as a function of their duration. Eachdata point represents 10.000, 10 s power-law distributed time traces. The input shows theinitial power-law distribution, the lighter colors show the number of retrieved changepoints,after applying Poisson noise and CPA, at different contrasts, with I = 10 . counts/s. We cansee that even at low contrast, events with long times between switching are retrieved, buteach contrast has a characteristic duration below which changepoints can’t be accuratelyretrieved. This puts a fundamental limit on the information that can be extracted froma given data set. (D) An occupancy diagram illustrating the behavior of the clusteringalgorithm for a system with m =4 intensity levels. The color scale indicates the amount oftime spent in each state m i after the assignment of states for a given number of availablestates n G . We see that when n G > m , effectively all segments are distributed across only n G ≤ m intensity levels.As next step in our Monte Carlo benchmarking we turn to time series with many, insteadof single, jumps. Fig. 3A shows a representative example for a simulated intermittent quan-tum dot with two states, assuming a contrast ratio between states of 2 × and 5 × s − . We generally observe that the recursive CPA algorithm accurately identifies switchingevents, barring a number of missed events of very short duration. From the CPA analysis weretrieve the time duration between switching events. Figure 3C shows a histogram of timedurations, plotted as a probability density function obtained from a whole series of Monte16arlo simulated time traces of varying contrast between states (see legend). Noteably, ifwe simulate quantum dots that have switching times that are power-law distributed, theretrieved distribution indeed follows the assumed power-law particularly for long times. Atshorter times, the histogram remains significantly below the power law, particular at low in-tensity ratios between the two assumed states. This indicates that CPA misses fast switchingevents, and is consistent with the observation from Figure 2C that a minimum photon countis required to observe switching events of a given contrast. As a rule of thumb, usual II-VIcolloidal quantum dots have a contrast between dark and bright states of around 5, meaningthat of order 200 photons are required to detect a change point with near-certainty. At theassumed count rates (2 × s − for the bright state) this means one expects CPA to failfor switching times below 10 ms, where the on-off time histogram indeed shows a distinctroll off. This result suggests one should interpret on-off time histograms from changepointdetection with care: one can generally rely on the long-time tail, but should determine theshortest time scale below which the histogram is meaningless on basis of the intensity levelspresent in the data. Performance of level clustering
Next we consider the performance of the grouping algorithm applied to the segmentation ofsimulated time traces. For reference, Figure 4A shows the Bayesian Information Criterionversus n G for the example of simulated dots with m = 2 , , n G approaches the actual number of levels with which the data wassimulated, and gently decreases once n G exceeds the actual number of levels in the data m . The fast rise indicates that within the assumption of Poisson distributed intensities,the data can not at all be described by fewer than m levels. The slow decrease is dueto the penalization of the BIC by the number of fit parameters. Since the BIC criterionactually relates to the logarithm of the probability with which n G states are the appropriatedescription of the data, even an apparently gentle maximum in BIC actually coincides with17 n G % , & / $ m m P ( m r | m , N ) m r › N fi = 5 · % › N fi = 5 · & › N fi = 5 · ' › N fi = 5 · ( Figure 4: (A) Examples of the BIC for simulated quantum dots, with m = 2 , ,
4. Themost likely number of levels is indicated by a small peak in L k . In the inset, we showthat the distributions indeed peak at their respective m . It should be noted that the BICshows a very sharp rise, and then from m onwards appears almost flat. There is in fact ashallow downward slope. Here one should keep in mind that the BIC is a logarithmic metric.On a linear scale the maximum is significant. (B-E) Likelihood P of finding m r levels fora simulated system, given m initial intensity levels and (cid:104) N (cid:105) the mean number of photoncounts. We see that the photon budget plays a defining role in the total number of statesthhat one can reliably resolve. At low photon budgets, the number of levels is systematicallyunderestimated, whereas at high photon budgets, P ( m i = m ) remains high even at high m . We see that m r is never overestimated.an accurate, unique determination of m .To gauge the accuracy of the retrieval of the number of states for multi-state quan-tum dots, we simulated quantum dot data with power-law distributed ( α = 1 .
5) switchingevents, assuming switching between from ˜ m = 2 to 10 equally likely levels, where we as-sumed intensity levels to be assigned to segments randomly, and where we assumed levelsfor simplicity evenly spaced from dark to bright. Lastly, all segments were reassigned anintensity according to Poisson statistics. In other words, we added shot noise.For many random realizations with different m and (cid:104) N (cid:105) we determine the most likelynumber of states m r (BIC( m r ) = max(BIC)) according to the clustering algorithm, andconstruct histograms of outcomes. The total photon budget is set by the product of as-sumed record length and the mean count rates of the different levels. The outcome of thesecalculations are shown in Figure 4(B-E), where each panel corresponds to a different photonbudget. A plot with only diagonal entries signifies that the number of levels retrieved by theclustering algorithm always correspond to the number of levels assumed, so m r = m . Athigh photon budgets (Figures 4D, E), the retrieval of the number of states is indeed robust,18ven for simulated dots that switch between as many as 10 intensity levels. At low totalphoton budgets (Figures 4B, C), we see that m r is often underestimated. This signifies thatthere is high uncertainty due to shot noise in the assigned intensity levels, so that levels cannot be discriminated within the photon budget.As a secondary metric, additional to the BIC, is in how the clustering algorithm assignsoccupancy to the levels. The clustering algorithm assigns to each data segment the mostlikely intensity level that it was drawn from. Occupancy is a metric for how often each of the nG levels available to the algorithm is actually visited by the measured intensity sequence.We find that if one allows the clustering algorithm to use more levels than originally usedto synthesize the quantum dot data ( n G > m ), the additional levels take essentially nooccupancy. We show this in Figure 3D for an example of a dot assumed to have fourintensity levels with a total photon count of 5 × . As soon as additional states ( m ≥ m . Thus we confirmagain that the grouping algorithm does not over-estimate the number of states. Accuracy of decay rates versus photon budget
Figure 5(A) shows examples of fitted simulated data for slow and fast decays, as examplesof the Monte Carlo simulations we have performed to benchmark the accuracy of decay ratefitting in function of photon budget, and decay rate (panel (B)). We find that the error in γ very roughly follows roughly a power law with an exponent of 0 . − .
1, with slower decayrates showing higher errors. Consistent with Ref. we find by Monte Carlo simulation thatone requires approx. 200 (50) counts to obtain an error below 10% (20%) in decay rate if onefits mono-exponential decay with free parameters. A problem intrinsic to the use of CPA isthat fast switching events may be missed, leading to an averaging of short time intervals withothers. This leads to decay traces that are in fact not attributable to a single exponentialdecay. 19igure 5: (A) Two examples of decay traces ( γ = 0 . , − ) with low photon count ( N =300) fitted to a single-exponential decay. (B) The standard error of the fitted decay ratew.r.t. the input decay rate as a function of the total photon count, for different decay rates.Each data point is the average of 10 simulated decay traces.20 DID diagrams γ Q V − ρ Q R U P I F W V V $ T c ∞ % T c & T c ' T c Figure 6: Four FDID diagrams of a simulated bimodal quantum dot with (
I, γ ) = (2 × s − , 0.1 ns − ) and (0 . × s − , 0.4 ns − ). In each diagram, we apply a different cutoff time T c to the simulated powerlaw at ∞ , 10 s, 1 s, and 0.1 s. We find that a shorter cutoff causesan increasingly strong smearing effect between the two states.Correlative FDID diagrams that plot correlations between intensity levels and fluores-cence lifetime form a powerful visualization of quantum dot photophysics. Essentialto the construction of FDID diagrams is that for each detected photon also the delay timeto the laser pulse that generated it is known so that decay rates can be fitted even to shortsegments of a time trace segmented by CPA.Finally we discuss the construction of FDID diagrams derived from CPA-analysis, againillustrated by examining simulated data for a quantum dot that switches between two statesof distinct intensity and lifetime. FDID diagrams are conventionally constructed from time-binned data, where it is interpreted as a simple histogram in which each time-bin contributesa single histogram count to one single intensity-decay rate bin. It is not trivial to extendthis notion to CPA-segmented data since CPA segments intrinsically have very different timedurations instead of having equal width as in conventional time-binning. We propose twomodifications to the construction of a FDID as a histogram. First, instead of representingFDID entries as a single binary entry in just one histogram bin (one time segment con-tributes one count to a single pixel in a FDID), we propose to incorporate the uncertaintyin intensity and decay rate that is associated with each time segment. To this end, each21egment contributes to the FDID diagram according to a 2D Gaussian function centered atthe CPA-segment decay rate and intensity (total counts C j in segment j divided by segmentduration T j ), where the width is given by the fit error in the decay rate and the shot noiseerror in the segment (cid:112) C j /T j . If one would apply this logic to regular time-binned data,giving each Gaussian contributor the same integrated weight, one obtains a diagram similarto a regular FDID histogram except that the results is smooth and with less dependencyon a chosen histogram bin width. Instead the feature size in FDID represents the actualuncertainties in intensity and rate.As a second modification, we propose to reconsider the weights of the Gaussians — i.e.,the integrated contribution to each entry in the FDID. For time-binned data one assigns eachsegment equal weight so that equal lengths of time contribute equal weight. Since CPA resultsin segments of unequal length, several choices for constructing FDID diagrams are possible,which to our knowledge have not been discussed in CPA literature. Giving equal weight toeach CPA fragment will lead to FDID diagrams from those obtained from binned data, sinceeffectively long time segments are then underrepresented compared to short segments. Thisleads to under representation of states with steeper powerlaw distributions in their switchingtimes. The direct equivalent to regular FDID-weighting for CPA-segmented data is thata segment of duration T i has a weight proportional to T i (henceforth ‘duration-weightedFDIDs’). Alternatively one could argue that since time-averaged intensity and fluorescencedecay traces are rather set by the contribution in emitted photons, one could instead use thetotal number of photons C i in each segment as weight (henceforth ‘count-weighted FDID’).It has been established in a multitude of studies of II-VI quantum dots that the dis-tribution of on/off times follows a power-law (exponents α on , off ) truncated by an exponen-tial with specific time τ l , giving the distribution t − α e − t/τ l . We analyze Monte Carlosimulated FDIDs to establish if there are conditions of the truncation time under whicha two-state quantum dot would appear not as a bimodal distribution. Figure 6 showsduration-weighted FDIDs for simulated quantum dot data. For panel (A) we consider a22wo-state quantum dot with power-law distributed switching times. In panels (B-D) weshow a quantum dot simulated with similar parameters, but with the maximum durationof the segments T c = 10 , , . t − α e − t/τ l ), unless truncation time τ l are as short as 20 ms, so that there are no segments with over approximately 10 counts.This limit of our benchmarking space zooms in on the regime where CPA intrinsically fails(Fig. 2(C)). In this regime, the originally assumed bimodal quantum dot behavior no longerresults in a bi-modal FDID. Instead a significant broadening is evident. We can concludethat for most realistic quantum dot systems, CPA-generated FDIDs will not suffer from thisartificial broadening artefact. Conclusion
In this work, we have provided a Python toolbox for change-point analysis, and for deter-mining the most likely intensity level assignment for intermittent multilevel emitters. Weinvestigated the limits of changepoint analysis and clustering as fundamentally set by pho-ton budget, and for the case of many state emitters. We have shown that for long switchingtimes, the typical power-law behavior of many quantum emitters can be accurately retrieved.We also show that in the case of many-state emitters, the number of intensity levels can beretrieved with high fidelity, provided the photon count is high enough. At low photon counts,the number of states is systematically underestimated. This shows in which way the photonbudget puts a fundamental limit on the amount of information that can be retrieved froma given TCSPC data set. We show that the photon budget also poses a limitation on theaccuracy at which the slope of a single-exponential decay can be retrieved. Additionally, weinvestigate the commonly used intensity-decay time diagrams. We show that with CPA, atwo-state simulated quantum dot is well-represented in an FDID, but when a cutoff is intro-23uced to match that commonly found in literature, the states become increasingly poorlydefined. While the Bayesian inference methods in this toolbox were reported earlier, thistoolbox and its applicability to many-state emitters will, in our view, be of large practical usefor many workers analyzing the complex photophysics of, e.g., inorganic perovskite quantumdots, Also, the toolbox can be used for theory development, following a workflow in whichhypotheses are cast in synthetic photon counting data, which in turn can be subjected tothe CPA analysis suite, to assess how hypothesized mechanisms express in observables, andin how far they are testable within a given photon budget.
Supporting Information Available
Detailed manual to the code and description of the example data set.
Acknowledgement
This work is part of the research program of the Netherlands Organization for Scientific Re-search (NWO). We are grateful to Ilan Shlesinger and Erik Garnett for their critical scrutinyof paper and code. Lastly we would like to express our gratitude to Tom Gregorkiewiczwho passed away in 2019, and whose encouragement and guidance in the early stages of theproject were invaluable.
References (1) Orrit, M.; Bernard, J. Single pentacene molecules detected by fluorescence excitationin a p-terphenyl crystal.
Phys. Rev. Lett. , , 2716–2719.(2) Lounis, B.; Orrit, M. Single-photon sources. Rep. Prog. Phys. , , 1129–1179.(3) Kimble, H. J. The quantum internet. Nature , , 1023–1030.244) Lodahl, P.; Mahmoodian, S.; Stobbe, S. Interfacing single photons and single quantumdots with photonic nanostructures. Rev. Mod. Phys. , , 347–400.(5) Somaschi, N. et al. Near-optimal single-photon sources in the solid state. Nat. Photonics , , 340–345.(6) Doherty, M. W.; Manson, N. B.; Delaney, P.; Jelezko, F.; Wrachtrup, J.; Hollen-berg, C. L. The nitrogen-vacancy colour centre in diamond. Phys. Rep. , ,1 – 45.(7) Castelletto, S.; Boretti, A. Silicon carbide color centers for quantum applications. J.Phys: Photonics , , 022001.(8) Toth, M.; Aharonovich, I. Single Photon Sources in Atomically Thin Materials. Annu.Rev. Phys. Chem. , , 123–142.(9) Toninelli, C. et al. Single organic molecules for photonic quantum technologies. preprint:arXiv:2011.05059 ,(10) Murray, C. B.; Norris, D. J.; Bawendi, M. G. Synthesis and characterization of nearlymonodisperse CdE (E = sulfur, selenium, tellurium) semiconductor nanocrystallites. J.Am. Chem. Soc. , , 8706–8715.(11) Talapin, D. V.; Lee, J.-S.; Kovalenko, M. V.; Shevchenko, E. V. Prospects of ColloidalNanocrystals for Electronic and Optoelectronic Applications. Chem. Rev. , ,389–458.(12) Shirasaki, Y.; Supran, G.; Bawendi, M.; Bulovic, V. Emergence of colloidal quantum-dot light-emitting technologies. Nat. Photonics , , 13–23.(13) Protesescu, L.; Yakunin, S.; Bodnarchuk, M. I.; Krieg, F.; Caputo, R.; Hendon, C. H.;Yang, R. X.; Walsh, A.; Kovalenko, M. V. Nanocrystals of Cesium Lead Halide Per-25vskites (CsPbX 3 , X = Cl, Br, and I): Novel Optoelectronic Materials Showing BrightEmission with Wide Color Gamut. Nano Lett. , , 3692–3696.(14) Swarnkar, A.; Chulliyil, R.; Ravi, V.; Irfanullah, M.; Chowdhury, A.; Nag, A. ColloidalCsPbBr 3 Perovskite Nanocrystals: Luminescence beyond Traditional Quantum Dots. Angew. Chem. Int. Ed. ,(15) Park, Y.-S.; Guo, S.; Makarov, N. S.; Klimov, V. I. Room Temperature Single-PhotonEmission from Individual Perovskite Quantum Dots.
ACS Nano , , 10386–10393.(16) Li, B.; Huang, H.; Zhang, G.; Yang, C.; Guo, W.; Chen, R.; Qin, C.; Gao, Y.; Biju, V. P.;Rogach, A. L.; Xiao, L.; Jia, S. Excitons and Biexciton Dynamics in Single CsPbBr3Perovskite Quantum Dots. J. Phys. Chem. Lett. , , 6934–6940.(17) Gibson, N. A.; Koscher, B. A.; Alivisatos, A. P.; Leone, S. R. Excitation IntensityDependence of Photoluminescence Blinking in CsPbBr3 Perovskite Nanocrystals. J.Phys. Chem. C , , 12106–12113.(18) Seth, S.; Mondal, N.; Patra, S.; Samanta, A. Fluorescence Blinking and Photoactivationof All-Inorganic Perovskite Nanocrystals CsPbBr3 and CsPbBr2I. J. Phys. Chem. Lett. , , 266–271.(19) Yuan, G.; Ritchie, C.; Ritter, M.; Murphy, S.; G´omez, D. E.; Mulvaney, P. The Degra-dation and Blinking of Single CsPbI3 Perovskite Quantum Dots. J. Phys. Chem. C , , 13407–13415.(20) Hou, L.; Zhao, C.; Yuan, X.; Zhao, J.; Krieg, F.; Tamarat, P.; Kovalenko, M. V.;Guo, C.; Lounis, B. Memories in the photoluminescence intermittency of single cesiumlead bromide nanocrystals. Nanoscale , , 6795–6802.(21) Frantsuzov, P.; Kuno, M.; Jank´o, B.; Marcus, R. Universal emission intermittency inquantum dots, nanorods, and nanowires. Nat. Phys. , , 519.2622) Krauss, T. D.; Peterson, J. J. Bright Future for Fluorescence Blinking in SemiconductorNanocrystals. J. Phys. Chem. Lett. , , 1377–1382.(23) Efros, A.; Nesbitt, D. Origin and control of blinking in quantum dots. Nat. Nanotechn. , , 661–671.(24) Galland, C.; Ghosh, Y.; Steinbr¨uck, A.; Sykora, M.; Hollingsworth, J. A.; Klimov, V. I.;Htoon, H. Two types of luminescence blinking revealed by spectroelectrochemistry ofsingle quantum dots. Nature , , 203–207.(25) Rabouw, F. T.; Lunnemann, P.; van Dijk-Moes, R. J. A.; Frimmer, M.; Pietra, F.; Koen-derink, A. F.; Vanmaekelbergh, D. Reduced Auger Recombination in Single CdSe/CdSNanorods by One-Dimensional Electron Delocalization. Nano Lett. , , 4884–4892.(26) Watkins, L. P.; Yang, H. Detection of Intensity Change Points in Time-Resolved Single-Molecule Measurements. J. Phys. Chem. B , , 617–628.(27) Hoogenboom, J. P.; den Otter, W. K.; Offerhaus, H. L. Accurate and unbiased estima-tion of power-law exponents from single-emitter blinking data. J. Chem. Phys , , 204713.(28) Ensign, D. L.; Pande, V. S. Bayesian Single-Exponential Kinetics in Single-MoleculeExperiments and Simulations. J. Phys. Chem. B , , 12410–12423.(29) Ensign, D. L.; Pande, V. S. Bayesian Detection of Intensity Changes in Single Moleculeand Molecular Dynamics Trajectories. J. Phys. Chem. B , , 280–292.(30) Crouch, C. H.; Sauter, O.; Wu, X.; Purcell, R.; Querner, C.; Drndic, M.; Pelton, M.Facts and Artifacts in the Blinking Statistics of Semiconductor Nanocrystals. NanoLett. , , 1692–1698. 2731) Houel, J.; Doan, Q. T.; Cajgfinger, T.; Ledoux, G.; Amans, D.; Aubret, A.; Dom-injon, A.; Ferriol, S.; Barbier, R.; Nasilowski, M.; Lhuillier, E.; Dubertret, B.; Du-jardin, C.; Kulzer, F. Autocorrelation Analysis for the Unbiased Determination ofPower-Law Exponents in Single-Quantum-Dot Blinking. ACS Nano , , 886–893.(32) Bae, Y. J.; Gibson, N. A.; Ding, T. X.; Alivisatos, A. P.; Leone, S. R. Understandingthe Bias Introduced in Quantum Dot Blinking Using Change Point Analysis. J. Phys.Chem. C , , 29484–29490.(33) Ensign, D. L. Bayesian statistics and single-molecule trajectories. Ph.D. thesis, StanfordUniversity, 2010.(34) Schmidt, R.; Krasselt, C.; Von Borczyskowski, C. Change point analysis of matrixdependent photoluminescence intermittency of single CdSe/ZnS quantum dots withintermediate intensity levels. , , 9–14.(35) Palstra, I. M.; Koenderink, A. F. A Python toolbox for unbiased statistical analysis offluorescence intermittency of multi-level emitters. Github repository https: // github.com/ AMOLFResonantNanophotonics/ CPA/ . First commit archived at Zenodo http:// doi. org/ 10. 5281/ zenodo. 4557226 ; 2021.(36) Zhang, K.; Chang, H.; Fu, A.; Alivisatos, A. P.; Yang, H. Continuous Distribution ofEmission States from Single CdSe/ZnS Quantum Dots.
Nano Lett. , , 843–847.(37) G’omez, D. E.; van Embden, J.; Mulvaney, P.; Fern´ee, M. J.; Rubinsztein-Dunlop, H.Exciton-Trion Transitions in Single CdSe-CdS Core–Shell Nanocrystals. ACS Nano , , 2281–2287.(38) Cordones, A. A.; Bixby, T. J.; Leone, S. R. Direct Measurement of Off-State TrappingRate Fluctuations in Single Quantum Dot Fluorescence. Nano Lett. , , 3366–3369. 2839) Schmidt, R.; Krasselt, C.; G¨ohler, C.; von Borczyskowski, C. The Fluorescence Inter-mittency for Quantum Dots Is Not Power-Law Distributed: A Luminescence IntensityResolved Approach. ACS Nano , , 3506–3521.(40) Rabouw, F. T.; Antolinez, F. V.; Brechb¨uhler, R.; Norris, D. J. Microsecond BlinkingEvents in the Fluorescence of Colloidal Quantum Dots Revealed by Correlation Analysison Preselected Photons. J. Phys. Chem. Lett. , , 3732 – 3738.(41) Li, H.; Yang, H. Statistical Learning of Discrete States in Time Series. J. Phys. Chem.B , , 689–701.(42) Wahl, M.; Gregor, I.; Patting, M.; Enderlein, J. Fast calculation of fluorescence corre-lation data with asynchronous time-correlated single-photon counting. Optics Express , , 3583.(43) Bajzer, ˇZ.; Therneau, T. M.; Sharp, J. C.; Prendergast, F. G. Maximum likelihoodmethod for the analysis of time-resolved fluorescence decay curves. Eur. Biophys. J. , , 247–262.(44) Frimmer, M. Spontaneous Emission near Resonant Optical Antennas. Ph.D. thesis,Universiteit van Amsterdam, 2012.(45) Clauset, A.; Rohilla Shalizi, C.; Newman, M. E. J. Power-Law Distributions in Empir-ical Data. SIAM Review , , 661–703.(46) Palstra, I.; Wenniger, I.; Patra, B. K.; Garnett, E. C.; Koenderink, A. F. Intermit-tency of CsPbBr perovskite quantum dots analyzed by an unbiased statistical analysis. Preprint on arxiv. 2102.09333. ,(47) K¨ollner, M.; Wolfrum, J. How many photons are necessary for fluorescence-lifetimemeasurements ?
Chem. Phys. Lett. , , 199–204.29 upporting information - manual to the Python code In this supplement, we will discuss the different uses of the Python toolbox for unbiasedstatistical analysis. The toolkit accepts both data taken in a TCSPC measurement and datacreated in Monte-Carlo (MC) simulations. We will first discuss the data processing of photontime stamp files, and subsequently the provided routine to simulate m -level single photonemitters. Processing toolbox - overview
Figure 7: Flow chart of the processing toolbox described in the main text.Figure 7 shows the different parts of the data analysis tools provivded in the toolbox, anddescribes a workflow to perform a full analysis. The diagram shows main functions as dark30lue ellipses, input and output files as green boxes, with blue frames indicating groupings orrelated tasks the code performs.The starting point for the code is a provided set of three lists of timestamps, whichrepresent three channels, i.e., two photon counting channels, and a pulsed laser referencerecord. These timestamps are provided as sorted lists of 64-bit integers through ”parquet”files. To convert integer timestamps to actual times, the user can specify a discretizationunit, which for experiments would generally be the intrinsic binsize of the TCSPC card used.The data is expected to be in ”reverse start-stop” format as usual for TCSPC experiment,meaning that the pulsed laser reference trace contains time stamps for the laser pulsesimmediately after each photon event. Simulated data is simulated with a user specifiedintrinsic binsize and assumed pulsed excitation repetition rate.The processing code has the following functionality. First, the data is segmented into in-tervals using changepoint analysis. This returns the jumptimes at which the emitter at handjumps in intensity, and between which the dot is constant in intensity, within the CPA thresh-old. Next, for each segment, the detector and photon channels are time correlated (block”segmented correlations”) to obtain fluorescence decay trace photon counting histogramsper segment, to which a decay rate can be fitted with the maximum likelyhood method(MLE). These processes lead to a set of intermediate, ASCII-formatted files that contain thesegmented jumptimes, intensities and decay rates. The reminder of the routines are foundon the right hand side of the diagram. The top block ”timestamp related functions” acts noton the segments, but on the timestamp list as a whole. They convert the timestamps listinto time-correlations between detectors, and detector and laser, across the whole dataset,to provide a decay trace, g (2)( τ ) and long-term autocorrelation trace. The remainder of theblocks post-processes the segmented data. This includes on one hand grouping accordingto Bayesian analysis, and on the other hand plotting of the segmented data as time traces,FDID diagrams, and switching time histograms.The philosophy of the code is that it can interchangeably run on experimental and simu-lated data simply by setting a Boolean switch in the code. The Boolean switch redirects theroutine to take information from the preamble XXX.py files, where XXX points to experimen-tal resp. simulated data. The preamble files contain the input and output data file pathsand required metadata for processing. The input time stamps are taken from the folderslabelled
Data XXX , while the output is saved in the folder
Output XXX with Each functioncan be run independently provided that any required intermediate files have already beengenerated.To demonstrate a complete run through the entire workflow, we provide the script
En-tireWorkflowScript.py . This script is just a sequence of independent function calls, with nopassing of variables between functions. We note that to work with experimental example,the user needs only run or adapt
EntireWorkflowScript.py . For exploring a large variety ofsimulated problems, it suffices to run the simulation routine (see below), and work with
EntireWorkflowScript.py , largely without need to adapt preamble simulated.py .Under the hood, the actual dataprocessing algorithms are contained in the four basicroutines:1.
CPA functions.py — segments the timestamp file using changepoint analysis2. grouping functions.py — clusters lists of jumptimes and counts per segment into the31ost likely set of levels, according to Bayesian information criterion.3. correlate jit.py — contains the timestamp-based channel auto and cross-correlation al-gorithm reported by Wahl et al. It provides decay trace histograms when correlatingphoton and laser channel, used for both segments (”segmented correlations”) and ag-gregate data record (”unsegmented correlation”). It provides g (2) ( τ ), when givingdetector channels A and B as task to correlate, and it provides on a logarithmicallycoarsened time axis the normalized intensity autocorrelation to perform the intermit-tency analysis proposed by Houel et al . We note that for typical time records of10 − events per channel, obtaining the g ( τ ) photon-photon correlations arethe most time-intensive. For this reason this part of the code is accelerated with thePython numba jit just in time compiler, and only called with 64-bit integer datatypes.4. fitexp MLE.py — performs maximum likelihood estimation based fitting of fluorescencedecay traces according to Bajzer et al., which is applied to both segmented correla-tions and the aggregate time trace.The remainder of the routines are for file input-output, bookkeeping and plotting, relyingon the core routines listed above for the actual work. Segmenting the data
Module name:
CPA wrapper.py
Input files: timestamps chX bin.parquet
Output files: metadata dataset.txt , segmentbounds.txt In this script, the data — experimental or generated by the dot-simulation routine — issegmented according to the user’s instructions. The script loads the data from three parquetfiles, which are binary data files containing all the photon event timestamps on the threechannels: the two measurement channels from the Hanbury-Brown Twiss setup of photo-detectors, and the reference channel. This is a ubiquitous setup for single-photon detection,but we note that, in the case of a different number of detectors, the code can be adaptedto accommodate for this, as in most scripts (except for g (2)( τ )), the events from the twophotodetectors are merged into a single list.The script initially retrieves a number of properties about the data, such as the numberof counts in all channels, and the length of the time record (timing of the last events), whichtogether with the timing bin width (specified in preamble.py ) are written to a metadata file.This file also includes a retrieve time offset between channels A, B and R estimated from thedata. That is, an event on a detector and the corresponding reference event will be offset by(1) the time it takes for the emitter to decay, and (2) the path difference between detectorsand reference diode. The quantities Tau0A bin, Tau0B bin refer to the latter. In simulateddata there is no such time offset, but in experimental data there are typically optical andelectronic (cable length) path differences which one needs to correct for before splicing twodetector channels into a single fluorescence decay trace.The script next calls cpa functions.py , which takes the combined photon timestamps fromchannel A and B as input, and returns the timestamps at which a changepoint occurs accord-32ng to Bayesian inference. This routine takes the basic formalism to find a single changepointin a dataset that has photon timestamps distributed according to exponential waiting timestatistics, as specified by Watkins, Ensign, and in the main text. This formalism is ap-plied iteratively to the resulting segments, subdividing them until no more changepoints arefound. The routines returns both the timestamps (jump times t , t , t , ...) and the associ-ated photon event indices, which is equivalent to the cumulative number of counts up to thetime of the jump. The single-changepoint detection routine requires a level of ‘skepticism’to accept or reject changepoints, set in the preamble python file. For an explanation on thewhich values are appropriate, we point the reader to and. Lifetimes per segment
Module name: segmented crosscorrs.py
Input files: timestamps chX bin.parquet , segmentbounds.txt Output files: corr seg.txt , segment fit data.csv This script again imports the photon events from the parquet files, as well as the seg-ment bounds (jumptimes) calculated and exported by
CPA wrapper . Looping through thesegments, the events on either detector are correlated with the events in the reference chan-nels using correlate jit , following the method of. The correlations are calculated over atime range equal to the (assumed) repetition rate of the excitation laser. This gives two2 × Q − × q matrix, with Q the number of found CPA jumps (including time series startand stop) and q the ratio of τ L the time between excitation pulses over the timing cardresolution. Thus each each of the 2( Q −
1) (factor 2 for the two detectors) rows represents afluorescence decay histogram for a segment. Since the correlation calculations can be quitelengthy, they are saved to corr seg.txt, for loading as needed. In addition, we speed up thesecalculations significantly by making use of Numba, which is a so-called just-in-time (JIT)compiler. This allows a specified section of Python code to be accelerated by just in timecompilation. For these segmented correlations, this gives a speedup of between one andtwo orders of magnitude. This JIT compiler is applied for all correlations calculated, othersincluding the g ( τ ), the full decay trace, and the coarsened autocorrelation.Following the interpretation-agnostic conversion from timestamps to decay histogramsper segment, an exponential decay (by default single exponential — several other dependen-cies are preprogrammed and could be chosen by the user) is fitted to each segment. This isdone using the maximum likelihood estimator (MLE) relevant to Poisson statistics, explainedby Bajzer et al. through the routine fitexp MLE . Before fitting, the data is prepped. Thisincludes timeshifting the photon channels such that the peak of the decay histograms landsat τ = 0 (each might have a different shift, determined by the routine as Tau0A bin and
Tau0B bin ), reversing the time axis (exponential decay tail towards positive time, undoingreverse start stop) clipping out negative times (detector rise), and includes an option forclipping out user-specified data intervals. This is useful for experimental data, where elec-tronic ringing in DPC/TCSPC chips sometimes causes artefacts that should not be fitted.Note that excluding this range from the fit is not the same as setting the counts for thisrange to zero, instead removing them entirely.After the data prepping per detector channel, the detector channels are spliced, so that33ecay rates are determined from the full photon record (sum of both detectors). Importantlyfor MLE-fitting of Poisson distributed data in which the probability distribution is linkedto absolute count values, the data is neither scaled, nor background substracted. By fittinga single exponential A exp( − iγt ) + b , three values and their respective errors are obtained.The routine exports A , γ , the error in γ , and b . It should be noted that in the subsequentprocessing routines only γ and δγ are used, as the best estimate for the count rate persegment simply follows from the number of counts in a segment N q divided by the segmentduration T q (error follows from the square root of the number of counts, divided by thesegment duration). The fitted values are saved, and can be reloaded as a pandas dataframein segment fit data.csv . The fit routine starts from an initial guess for parameter values thatare estimated from the data itself.For insightful explanations of the maximum likelihood estimator, we refer the readerto. Grouping, determining most likely levels
Module names: grouping wrapper.py
Input files: segmentbounds.txt
Output files: grouping.txt
This routine operates on the result of
CPA wrapper , i.e., the file segmentbounds.txt pro-duced by it. This consists of a list of jump times ( t , t , .. ) and (cumulative) counts persegment, which convert into segment durations T q = t q − t q − , counts per segment N q andper-segment count rates I q = N q /T q . Following the main text, the segments are clusteredinto a first guess for the actual Bayesian grouping analysis using a recursive hierarchicalclustering algorithm that crawls through a data set of Q − Q − for large data sets, which consists of hierarchically grouping small blocks of300 subsequent segments to 30 levels each, concatenating those, and then grouping. Theresult is groupings of the data in a range of n G levels (in our analysis we didn’t go beyond( n G ) = 20).After this initial guess, the segment allocation to intensity levels, and the level-values areoptimized following Section 2.3.2 of. Here, the segments, divided into exactly n G groups bythe initial hierarchical clustering, are iteratively swapped out, until the most likely groupingon the assumption of exactly n G levels is achieved, with the best estimate for the intensitylevels. Thus for each hypothetical distribution in n G ∈ mingroups , ..., maxgroups levels, themost likely assignment of each segment to an intensity level I , ... I n G is determined. Thisgives a “trajectory” of the emitter through the different levels, and is the which-level infor-mation for each segment, for each of the tested n G . Finally, to determine which value of n G most likely describes the data, the Bayesian Information Criterion (BIC) is calculatedas function of n G . This gives the likelihood of a given n G to be the true, underlying set oflevels. For simulated m -level dots indeed, the BIC peaks at n G = m The grouping routine returns its result as plots and ASCII data in grouping.txt . This34ata folder contains three sets of information. The first is a list of the values of the BayesianInformation Criterion, discussed in the main text and shown in Fig. 4a. The second is asquare matrix, where the values b n G ,m describe the fraction of time spent in state m , giventhat the grouping algoithm had n G intensity levels available to it. This yields figures such asFig.4b in the main text. Lastly, the file contains a 2D matrix that describes the “trajectory”of the emitter through the available levels, given that n G levels are available to it. There isno figure of this in the main text, though Fig. 1D. provides a sketch. Postprocessing segmented data, plot routines
Time trace
Module names: plot timetrace.py
Input files: timestamps chX bin.parquet , segmentbounds.txt , segment fit data.csv Output files: none
This script functions purely as a visualization of the behavior of the emitter as a functionof time. For plotting purposes the photon events are binned in equidistant time bins to obtainan intensity trace in time, which is plotted alongside the intensities and decay rates retrievedby CPA (i.e., without any binning). Additionally, a histogram of the intensity levels in thetime trace is plotted. This is an unweighted histogram, where all CPA-determined segmentshave equal weight, independent of their duration. The routine also plots a time trace of thefitted decay rates with its corresponding histogram besides it. For plotting purposes, theroutine attempts to guess reasonable plot axes from the data.
Switching time statistics
Module names: switchingtimehistogram.py
Input files: segmentbounds.txt
Output files: switchinghist.txt
In this section of the code, the time durations T q of all found CPA segments are takeninto a histogram with logarithmic bin, implementing apropriate normalization to convertthe segment durations to a segment duration probability distribution function distribution.Next, the routine fits a power law, as typically observed for many single emitters follow intheir on/off times. Because CPA is more likely to miss switching events if they are veryshort (segment with <
50 to 100 photons), the found histogram will slope off from the initialpower law at short switching times. The routine will therefore fit a power law to the partof the curve at larger switching times than a chosen threshold value. It uses the minimizemodule from scipy.optimize to do a least-squares fit to the data, reporting as output the fittedpower law exponent and a plot to assess the quality of fit.
FDID diagrams
Module names:
FDID wrapper.py
Input files: segmentbounds.txt, segment fit data.csv
FDID.txt
In this script, the segment intensities and segement-resolved decay rates are compiled intofluorescence decay rate intensity diagrams (FDIDs). For this, the I q , γ q values of all segmentsexported into segment fit data.csv from CPA wrapper.py are used to build a 2D histogram. Forthe FDID, instead of simply compiling a histogram of square bins, each segment contributesa smooth 2 D Gaussian, the widths of which are determined by the error on the segmentdecay rate and the estimated uncertainty in the count rate for that segment, derived fromthe fact that counts per segment are Poisson distributed. The weight of each Gaussian isdefined as its integral. As explained in the main text, three different definitions of weightper segment can be used: (1) all CPA segments have the same weight, (2) the weight ofa segment is determined by the number of photon events in it, (3) the segment weight isdetermined by the time duration of the segment. When applying weighting, it is the volumeunder each 2 D Gaussian that is weighted. A single routine
FDID functions constructs theFDIDs as sum of Gaussians, with the desired choice of weight passed as argument fromthe routine
FDID wrapper that handles file input, and plotting. For plotting purposes, thewrapper routine attempts to guess reasonable plot axes from the data. The resulting FDIDsare output as ASCII files and plots. Output: • FDID prio none, 2D array, the FDID where all CPA segments have the same weight • FDID prio cts, 2D array, the FDID where the CPA segments are weighted by theamount of counts • FDID prio dur, 2D array, the FDID where the CPA segments are weighted by theirduration • FDIDextent, tuple, the (x1, x2, y1, y2) extent of the FDID. Used for plotting
Memory effects
Module names:
FDID wrapper.py
Input files: segmentbounds.txt, segment fit data.csv
Output files: none
This script converts the segmentation results from
CPA wrapper into correlative plotsthat highlight if there are memory effects and correlations between segment observables.The routine reads in segment fit data.csv , i.e., sequencies of intensities I q , the fitted decayrates γ q , and segment durations T q . Due to the large range of segment durations, plots arebased on log T q . The generated plots report on:1. Aging:
Time trace of I , γ and log T q versus ‘wall clock time’. In these figures, theaging of an emitter can be observed, e.g. if the dot becomes dimmer over time, orsegment lengths grow or shrink. 36. Concurrent observable correlations
A set of 2D histograms of the log T q versus(1) intensity I q and (2) the fitted decay rate γ q . This can be used to detect if thereare correlations between these quantities. These diagrams are anaologous to FDID orFLID diagrams, which are just 2D histograms of I q versus γ q .3. Memory
A set of 2 × x n against x n − and x n , with x ∈ I, γ, log T q . These are 2D histograms that examinethe same observable along both their axes, but at subsequent steps in the segmenthistory. For instance, plotting γ n against γ n − highlights if subsequent decay rates areuncorrelated, or correlated.4. Segment sequence autocorrelations of the lists I q , γ q and log T q . Note thatautocorrelating a list such as I q is not equivalent to calculating the intensity autocor-relation, i.e., autocorrelating I ( t ) as provided by Houel autocorr.py , since the segmentintensity list I q contains intensity levels of segments with highly unequal duration,whereas the autocorrelation of I ( t ) depends on the switching dynamics in time. Thesegment sequence autocorrelation brings out the number of switching cycles over whichthe correlated observable shows a memory. Unsegmented data analysis
Decay trace and fit, entire dataset
Module names: unsegmented crosscorrs.py
Input files: timestamps chX bin.parquet
Output files: corr unseg.txt
The entire data set is processed by correlate jit to obtain the correlation between detectorchannel and reference channel, leading to two decay traces (for detector A, and B). Asin the case of the segmented decay trace determination routine, this data corrected byzeroeing the time axis, reversing the start-stop reversed time axis and removing the detectorartefacts (risetime, electronic ringing). The two detector channels are then merged to obtaina fluorescence decay histogram. To this, an n -exponential curve can be fitted. The samePoisson-MLE is used as for the unsegmented data. Using the Minimize function from theScipy Optimize module, we minimize our merit function such that a n-exponential decay F A j ,γ j ,B ( τ ) = b + n (cid:88) j =1 A j exp( γ j τ )is fitted to the calculated decay trace. Here, A j and γ j are the proportionality constantsand the decay rates of the decay j , respectively, and b describes constant background noise.When fitting, we optimize for A j , γ j , B . How the fit merit function is calculated from thedata and the fit function is discussed in the main text.37 utocorrelation g (2) ( τ )Module names: g2 tau.py Input files: timestamps chX bin.parquet
Output files: g2.txt
This function reads the photon detection time stamps from parquet file, and then calcu-lates the unnormalized intensity autocorrelation g (2) ( τ ) by correlating the two detector chan-nels following the method outlined in to obtain g (2) ( τ ). For data taken in a Hanbury-BrownTwiss setup, this should reveal antibunching, as is the case for the example experimentaldata set. The simulated data generated by running the provided routine also antibunches.It should be noted that the routine does not correct for relative constant offsets in channeltimings due to, e.g., differing electronic delays (cable lengths), which, if present, expressesas a small shift of the the g (2) ( τ ) along the the time-axis. Long-time intensity correlations
Module name:
Houel autocorr.py
Input files: timestamps chX bin.parquet
Output files:
Houel autocorr.txt
Houel et al. proposed that intensity autocorrelations on the time scales of milliseconds toseconds express an unbiased assessment of intermittency characteristics, relating at least forbinary on/off behaviors to power law exponents of dwell times. While Houel et al. operatedon intrinsically timebinned data (camera frames), we provide a routine that operates directlyon timestamps. The routine Houel autocorr.py reads in timestamps for two detector channelsfrom Parquet files, splices the two channels in a single array of time stamps, and then calls correlate jit to autocorrelate the timestamp array. In view of the large time span over whichthe correlation is requested compared to the time-discretization intrinsic to the data (secondversus 0.1 ns), a standard correlation as used for g (2) ( τ ) is prohibitively slow. Therefore correlate jit provides a version of the timestamp correlation algorithm with logarithmicallyincreasing steps in τ through a data coarsening approach as outlined in. As opposed to theunnormalized g (2) for antibunching, now the output that we refer to as ACF is normalized.To ACF-1 a powerlaw with exponential roll off At − C exp( − Bt )is fitted where we use the minimize function of Scipy’s Optimize module and a least-squaresfit. According to Houel et al. the fit parameters trace back to power law exponents in thecase of binary (on/off) intermittency. Auxiliary routines
The role of the remaining routines in the processing toolbox is • acc functions.py contains shorthand for auxiliary mathematical functions, and 2D plotcolor map 38 EntireWorkflowScript.py example script executing full sequence of processing steps.Booleans in the file set the requested tasks (e.g., simulated or experimental data),but otherwise no data is passed from function to function. The only dependency isthrough the intermediate data files • loaddata functions.py contains file input output commands for reading and writing thevarious files in the appropriate formats and directories • pypreamble measured.py contains the file paths to the example experimental data, andrequired output data, as well as metadata required for processing. These include theTCSPC discretization bin size of the experiment (0.165 ns, Becker and Hickl DPC-230)and the laser repetition rate (10 MHz) in the experiment that generated the data. • pypreamble simulated.py contains the file paths to simulated emitter data, and re-quired output data, as well as metadata required for processing. It should be notedthat parameters required for processing, i.e. the assumed time bin discretization unitand laser repetition rate need to be specified in the file, and need to be consistent withhow the data was generated. They are not automatically passed from the simulationroutine to maximize the direct exchangeability of experimental and simulated data fortesting the toolbox. Binning instead of using CPA
In the sections above, we have discussed the analysis following segmentation of the data usingCPA. However, the same analysis can be applied following segmentation using equidistantbins. If the user is interesting in applying equidistant binning, this can be done by flippingthe Boolean
CPA insteadof binned in the workflow module. The output from this segmen-tation and analysis will be saved in separate folders, to avoid overwriting CPA-sectioneddata. The only script that will not work when equidistant binning is applied, is of course switchingtimehistogram , as its function starts off with building a histogram of the durationof the segments.
Simulate N-level emitters
Module name: simulate qdot.py
Output files: timestamps chX bin.parquet, nominalsegmentation.csv
We provide a routine to generate time-tagged data of an m -level single emitter for MonteCarlo simulations testing the toolbox, or for testing how assumed photophysics would appearin experiments. This routine is provided in a separate folder from, and is in terms of codeindependent from, the processing toolbox. It produces timestamps in parquet file formatthat, together with an appropriately specified preamble simulated.py file can be processed inexactly the same manner as the example experiment data set .In this section, we will describe the method used to simulate time-tagged data of an m -level single emitter with MC simulations. The entire workflow of generating a dot can39e executed by running entire simulate qot workflow.py . This script generates photon timestamps, and also creates plots and a reference datafile to examine the assumed and simulatedintermittency behavior.The simulation of the emitter behavior happens in a number of steps, as follows: Drawing the segmentation information
First we randomly draw the random segmentation sequence according to which the photontime stamps will be subsequently randomly drawn, using the following procedure:1. The user specifies the number of levels nlevels = m , and the associated nominal in-tensity levels, decay rates, and dwell-time power law exponents through the preamblefile.2. A random list is created with indices that label the order in which the intensity levelsare visited. For each entry in the list, the index is drawn independently and all levelshave equal likelihood. If an entry is a direct repeat of its predecessor, it is removedfrom the list. The segment lengths are drawn using the following formula s = ( a β + ( b β − a ) x ) /β (2)where a is the short time cut off Tshortest ns , b is the long time cut off Tlongest ns ,and β is 1 − α , with alpha the specified power law exponent for the given level. x is a random number generated with numpy’s random.rand module, and has units ofdtau ns. Drawing photon events
The next step is to draw photon timestamps according to the randomly drawn segmentationsequence. The segmentation sequence specifies a list of jump times t , t , ..., t Q and a con-comitant sequence of levels visited, which specifies the nominal count rate, and decay. Wedraw the photon events in a two step process:1. Since we envision simulating a single emitter in a pulsed experiment, every laser repe-tition cycle there can be at most one photon resulting from an excitation event. Thatmeans that for a given nominal count rate there is a known probability p that onephoton is generated. We thus draw photon time stamps with a resolution equal to thelaser repetition time τ L by separating a segment of length T q in T q /τ L bins, drawing foreach bin a random uniformly distributed number between 0 and 1, and conclude thata given cycle corresponds to a photon if the corresponding drawn random number ex-ceeds p . In the provided example, τ L = 100 ns, which gives the granularity with whichthe time stamps are defined at this step. The probability of a photon event is givenby p = I q × τ L . It should be noted that drawing rare events according to uniform dis-tribution results in Poisson distributed count rates upon binning. Programmatically,only the list of timestamps labelling laser cycles with a resulting emitted photon are40igure 8: Flow chart of the routine to simulate multi-level emitters.kept, and stored in list timestamps chR bin . This ultimately is the reference channelfor the simulated two-detector TCSPC measurement.2. Next we fine-grain the photon time stamps, to account for the fluorescent decay trace.For each timestamp in the reference channel, a photon arrival time relative to the laserpulse is drawn through numpy’s random.exponential module, with the appropriatedecay rate from g lst for the relevant segment.3. from a memory management viewpoint, step 1 is prohibitive if applied to the entiredesired time trace (typically 10 laser cycles, for only 10 photon events), so we applyit segment by segment. To avoid the slow python append function, the timestampsare stored in a preallocated, overdimensioned container array, which is clipped to sizeafter generating all the time stamps. Since we provide a simulation toolbox to mimicTCSPCS experiments, all times are expressed as integer multiples of an assumed cardbinsize dtau ns . This allows us to express all time quantities in int64 format.4. uncorrelated dark noise for two detector channels is drawn according to step 1, within step 2 a uniform random arrival time distribution.41 istributing events over channels The outcome of the previous step is a set of time stamps timestamps chR bin correspondingto laser pulses with a concomitant detection, and a list of delay times relative to the laserpulses. These are translated into antibunching detector channels as follows:1. We construct photon arrival time stamps in reverse start stop, by taking timestamps chR bin adding the delay times, and subtracting the duration of exactly one laser period. Thisleads to a list of photon time stamps timestamps bin .2. For each event in timestamps bin we randomly assign it to detector channel A or de-tector channel B, with equal and independently drawn probability.3. Finally we splice the generated uncorrelated detector noise events noiseclicks chX with X ∈ A, B into the data stream. The concomitant laser pulse arrival times was previ-ously sent into timestamps chR bin .These operations ensure that the two simulated detector channels will show antibunching in g (2) ( τ ), with a visibility limited by the signal to background count rate. The backgroundnoise photon events are only added now, as adding them earlier would cause the noise toantibunch as well. Input and output files
The simulation module requires input parameters, which are specified through the mod-ule preamble for simulation.py . The input parameters include parameters of the simulatedemitter (number of levels, concomitant nominal count rates, decay rates, and power lawexponents) and of the simulated measurement (laser repetition time, time discretization binsize, background count rate, and total wall clock length of the measurement record).It is important to note that most of the variables are independent of the variables used forprocessing the data, except the variables that pertain to the properties of the “measurement”itself. These are what in a measurement corresponds to the timing resolution of the countercard and the repetition rate of the excitation, called dtau ns and minitimebin ns respectively.When analyzing simulated data, the values of these two variables must be the same as indata generation.As output three lists of timestamps, with timestamps as 64 bit integers in units of thetiming resolution dtau ns , representing two detectors in a Hanbury-Brown Twiss configu-ration, and a laser pulse reference channel for both. These are saved in a compact binaryformat as parquet files. As a service to the user, the simulation writes out the datafile nomi-nalsegmentation.csv with the jump times, and the concomitant sequence of nominal segmentcount rates and decay rates according to which the data was drawn. This can e.g. be usedto check the accuracy of the CPA, grouping, and lifetime fitting routines.42 xample experimental data
We provide select example data sets with the toolbox that were obtained on a CsPbBr quantum dots, as reported in Ref. (dots 2, 4, 6, 8, 16, with reference to supplement of thatpaper). While the full experimental protocol is reported in that paper, here we report theinformation pertinent for interpreting the time tagged data. The data was obtained usingpulsed excitation at 10 MHz repetition rate (LDH-P-C-450B pulsed laser diode, PicoQuant)at 10 MHz, with confocal excitation and collection through an oil objective (Nikon Plan APOVC, NA=1.4). Detection was on two fiber-coupled avalanche photodiodes (APDs) (SPCM-AQRH-14, Excelitas) in a Hanbury-Brown & Twiss configuration. The APDS were coupledto a photon correlator (Becker & Hickl DPC-230). It uses a 0.165 nm temporal binsize, andsimultaneously records absolute photon arrival times from both APDs, and the arrival timesof laserpulses subsequent to photon detection on any detector. We note that this means ourdecay traces are taken in reverse start-stop configuration, as is standard in TCSPC to reducedata rates. A small time interval centered at around 30 ns is subject to an electronic artefactwhich we attribute to a ringing in the DPC-230 TDS timing chips. This interval is includedin the supplied data - the toolbox is set up to exclude it from the analysis. Our measurementsare taken over 120 seconds of acquisition time. The timetagged data was converted fromthe Becker and Hickl dataformat using C Dependencies
The code is posted on github https://github.com/AMOLFResonantNanophotonics/CPA/ . Thiscode has been developed in Python 3.8. It requires to install the libraries numpy [1.19.2], mat-plotlib [3.3.2], pandas [1.1.3], as well as numba [0.15.2] for jit-speedup of the Wahl algorithm,and pyarrow [2.0.0] for dealing with parquet files. Numbers in square brackets [.] indicatethe version installed when developing the code. The code and the sample experimental dataare archived at Zenodo http://doi.org/10.5281/zenodo.4557226http://doi.org/10.5281/zenodo.4557226