Abstract:
A method of extracting circuit parameters from picosecond-scale photon timing measurements is disclosed. In one embodiment, the method is implemented by a system that comprises a photomultiplier, a data acquirer, and a processing module. The photomultiplier detects photons emitted from current-carrying channels in an integrated circuit, and associates a detection position and a detection time with each detected photon. The data acquirer receives position and time signals from the photomultiplier, and further receives a trigger signal. The data acquirer determines a relative detection time for each photon by combining the time and trigger signals. The data acquirer gradually compiles the photon detection data and makes it available to the processing module. The processing module responsively determines optimal values for a parameterized model of the data. The model is preferably based on non-homogeneous Poisson process statistics, and may employ a maximum likelihood approach to estimating the optimal values.

Description:
BACKGROUND  
         [0001]    [0001]FIG. 1 shows an illustrative cross-section of an integrated circuit (not to scale). An integrated circuit is typically fabricated on the “front” side  102  of a silicon wafer  104 . The integrated circuit is formed in stages, beginning with the careful introduction of impurities into localized regions of the silicon. The regions of impure (“doped”) silicon are termed “wells”. Small, heavily-doped wells are used as connection terminals of individual integrated device structures. The structure of a metal-oxide-semiconductor field effect transistor (MOSFET) includes two closely-spaced wells, with a conductive (often metallic) “gate” layer deposited on the front side  106  of the wafer  104  between the wells. (A thin oxide layer separates the gate from the silicon, hence the description “metal-oxide-semiconductor”). The conductivity of the silicon between the closely-spaced wells (the “channel”) can be changed by applying a voltage to the gate; this provides the controlled switching behavior that allows the device to operate as a transistor.  
           [0002]    Once the transistors and other circuit components have been formed, a series of layers  108  are added to interconnect the circuit components. These layers include conductors that “wire” the components together, and insulation that protects and separates the conductors. The conductive materials in the layers  108  are opaque, but it has been discovered that the integrated circuit can be studied from the “back” side  106  of the wafer  104 . Material may be ground away from the back side  106  to thin the wafer  104 , thereby increasing its transparency.  
           [0003]    Integrated circuits have microscopic dimensions, which make them difficult to test and characterize. Tsang, Kash, and Vallett present a new technique for measuring integrated circuit parameters in “Picosecond imaging circuit analysis,”  IBM J Res. Develop.,  Vol. 44, No. 4, July 2000, which is hereby incorporated by reference. Some applications of what Tsang, et al, call their “PICA” technique are outlined by Lundquist and McManus in “Characterize gate-level transistor performance with PICA,”  Semiconductor International,  July 2001, which is also incorporated herein by reference. Among these applications are the measurement of propagation delays and relative timing parameters in circuits. Such applications may prove invaluable in quick pinpointing of the precise nature and location of faults.  
           [0004]    The PICA technique relies on a combination of physical phenomena that are present in modem digital circuits. First, the MOSFETs employed in digital circuits typically operate in the saturation region of their current-voltage curves when “ON”. In this region, very high electric fields exist in the channel. Charge carriers (i.e., electrons and holes) can quickly gain a significant amount of kinetic energy in such electric fields, and indeed, many “hot” carriers are generated in this manner when current flows through the channel. A variety of scattering and recombination mechanisms may strip the energy from “hot” carriers, and in so doing, may trigger the emission of a photon of light. The light is emitted over a wide range of frequencies, but the infrared band of the spectrum is particularly significant because silicon is relatively transparent there.  
           [0005]    The upshot is that current-carrying transistors in modem digital circuits emit infrared light. FIG. 1 represents the emission of infrared light as wavy lines originating from the transistor channels  110 ,  112 .  
           [0006]    One of the key factors to the success of highly integrated digital circuits is the use of complementary metal-oxide-semiconductor (CMOS) technology. CMOS technology consumes a dramatically smaller amount of current (and hence power) than other transistor technologies. The reduction in current is a consequence of pairing p-type MOSFETs with n-type MOSFETs, so that one of the transistors is “OFF” when the other is “ON”.  
           [0007]    See, for example, FIG. 2, which shows a CMOS inverter. When the input voltage VIN is high (close to the supply voltage VDD), this puts the p-type MOSFET  202  in the “OFF” state, and puts the n-type MOSFET  204  in the “ON” state. The output voltage VOUT is pulled low (close to ground). Conversely, when VIN is low, transistor  202  is “ON” and transistor  204  is “OFF”, causing VOUT to be pulled high. Except briefly in the midst of VIN transitions between high and low voltages, one of the transistors is always “OFF”, so there is no steady-state current consumption. Only during transitions between states does any current flow occur. Because current flow is indicative of a switching transition, and because infrared light emission is indicative of current flow, we can use photon detection to measure the timing parameters of switching in digital circuits.  
           [0008]    However, accurate extraction of timing parameters from PICA measurements is not straightforward. Photon emission is a statistical process, and it is expected that the probability of detected photon emission for current VLSI technology is on the order of  10 - 9  infrared photons per switching event per micron of channel width. To get a statistically useful measurement, many switching events must be generated, and the timing of the separately detected photons must be correlated with a repeatable trigger. The relative timing of the photons can be plotted as a histogram, and the peak of the histogram may be used to determine the timing of the transistor&#39;s switching transition relative to the trigger. Unfortunately, the histograms are “noisy” as a result of the granularity of the events (i.e., either a photon is detected, or it is not detected) and the rarity of events. The noise makes the measurements difficult to interpret. The histograms gradually approach a smooth curve if the measurements are taken over a sufficiently long period of time, but this is not expected to be feasible for small width devices or devices that switch infrequently relative to the clock frequency. Accordingly, an objective method of extracting circuit parameters from picosecond-scale photon timing measurements is desirable.  
         SUMMARY  
         [0009]    Accordingly, an objective method of extracting circuit parameters from picosecond scale photon timing measurements is disclosed herein. In one embodiment, the method is implemented by a system that comprises an imaging photomultiplier, a data acquirer, and a processing module. The imaging photomultiplier detects photons emitted from current-carrying channels in an integrated circuit, and associates a detection position and a detection time with each detected photon. The data acquirer receives position and time signals from the imaging photomultiplier, and further receives a trigger signal. The data acquirer determines a relative detection time for each photon by combining the time and trigger signals. The data acquirer gradually compiles the photon detection data and makes it available to the processing module. The processing module receives data indicative of a photon intensity profile, and responsively determines optimal values for a parameterized model of the data. The model is preferably based on non-homogeneous Poisson process statistics, and may employ a maximum likelihood approach to estimating the optimal values. The intensity function may be modeled using one or more triangular peaks with a constant for background intensity, and timing jitter can be accounted for. Other intensity function models may alternatively be employed.  
       
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS  
       [0010]    A better understanding of the present invention can be obtained when the following detailed description is considered in conjunction with the following drawings, in which:  
         [0011]    [0011]FIG. 1 shows an illustrative integrated circuit device;  
         [0012]    [0012]FIG. 2 shows a schematic of a CMOS inverter;  
         [0013]    [0013]FIG. 3 shows an exemplary picosecond-scale photon detection system;  
         [0014]    [0014]FIG. 4 shows an exemplary functional block diagram of an apparatus for analyzing photon detection data;  
         [0015]    [0015]FIG. 5 shows a hypothetical graph of photon detection data;  
         [0016]    [0016]FIG. 6 shows one parameterized model of a triangular intensity peak;  
         [0017]    [0017]FIG. 7 shows a hypothetical graph of a maximum likelihood estimate of a triangular intensity function with a constant background intensity for the photon detection data;  
         [0018]    [0018]FIG. 8 shows a hypothetical graph of a maximum likelihood estimate of a time-jittered triangular intensity function with a constant background intensity for the photon detection data; and  
         [0019]    [0019]FIG. 9 shows a flowchart of a preferred method of extracting circuit parameters from photon detection data. 
     
    
       [0020]    While the invention is susceptible to various modifications and alternative forms, specific embodiments thereof are shown by way of example in the drawings and will herein be described in detail. It should be understood, however, that the drawings and detailed description thereto are not intended to limit the invention to the particular form disclosed, but on the contrary, the intention is to cover all modifications, equivalents and alternatives falling within the spirit and scope of the present invention as defined by the appended claims.  
       DETAILED DESCRIPTION  
       [0021]    [0021]FIG. 3 shows an illustrative functional block diagram of an apparatus for making picosecond-scale photon timing measurements. The apparatus includes a test pattern generator  302  which generates a sequence of input signal vectors (“test vectors”) that are designed to make the integrated circuit operate in the desired manner. The desired operation depends on the circumstances. For example, when a fault has been detected in the circuit under certain operating conditions, the test pattern generator may be programmed to duplicate those operating conditions. When it is desired to measure a propagation delay between circuit components, the test vectors may be designed to repeatedly cause transitions in those components.  
         [0022]    The test pattern generator  302  provides the test vectors to the device under test  304 , which is supported in a test stand  306 . The test stand  306  accepts input signals from the test pattern generator  302 , and provides the input signals to the device under test  304  along with power and any other needed signals (e.g. clock signals, bias voltages, etc.). The test stand  306  may also accept output signals from device  304  and provide them to another functional block.  
         [0023]    It is expected that device  304  will typically be fully packaged integrated circuit chip which has been modified to make the back side of the silicon wafer viewable. In addition the wafer may have been thinned to increase the transparency of the wafer to infrared light.  
         [0024]    As the device  304  operates, currents flow through various components, and in particular, through metal-oxide-semiconductor field effect transistors (MOSFETs) undergoing switching transitions. The current-carrying components emit photons of infrared light that pass through the wafer and radiate from the back side. A microscope lens  308  captures and redirects the photons into an imaging photomultiplier  310 .  
         [0025]    Photomultiplier  310  performs two functions: amplification and measurement. A single photon impacting on the photomultiplier  310  knocks an electron free from a surface. The photomultiplier amplifies the effect of the freed electron to create a cascade of electrons that form a detectable current pulse. The photomultiplier  310  also measures the position of the incident photon in an image grid by determining where the current pulse was detected. The image grid may be divided into 512 or more rows and 512 or more columns. Photomultiplier  310  preferably provides an X signal indicating the row, and a Y signal indicating the column in which the photon is detected. In addition, the transient current pulse is preferably conditioned to form a T signal that indicates the time when the photon is detected, to within picosecond-scale resolution. These signals may or may not be in digital form.  
         [0026]    The final functional block of the measurement apparatus is the Data Acquirer and Analyzer (DAA)  312 . The DAA receives the X, Y, and T signals from the imaging photomultiplier  310 , and receives a trigger signal from the test pattern generator  302 . The trigger signal may alternatively be provided from the device under test  304  and/or the test stand  306 . The trigger signal is a reference signal, and typically it is correlated with the switching event(s) being measured on device  304 . The DAA  312  measures the T signal relative to the trigger signal. If the signals are digital, the DAA  312  calculates a difference between the T signal and the trigger signal. If the signals are not digital, it determines the photon detection time by measuring an interval between a transition of the trigger signal and a transition of the T signal. This relative time value is hereafter denoted “R”.  
         [0027]    If the X and Y and R values are not in digital form, the DAA  312  preferably converts them to digital form. The DAA  312  then preferably stores the (X,Y,R) triplets for later analysis, although in an alternative embodiment, the DAA  312  constructs a histogram and then discards the triplets. The histogram is a three dimensional array of possible (X,Y,R) triplets, in which a count is associated with each of the values. When a particular (X,Y,R) triplet is received, the count associated with that triplet is incremented. Note that histogram is not necessarily stored in the form of a three dimensional array, but may alternatively be represented by a more compact data structure.  
         [0028]    The data collection performed by the DAA  312  may be divided into sessions. At the beginning of a session, a new data file is created and/or the histogram is initialized with zeros. Data collection then commences and continues until the session is temporarily suspended or is ended. At the end of a session, the data collection is halted and the data file and/or histogram is exported or stored in long term storage. The data can then be analyzed in a variety of ways by DAA  312  or by some external system.  
         [0029]    [0029]FIG. 4 shows an illustrative functional block diagram of a system  400  for analyzing data obtained by the apparatus of FIG. 3. The system  400  may be contained within DAA  312 , or alternatively, it may be an external system coupled to DAA  312 . In yet another embodiment, it is an entirely separate system to which the data from DAA  312  is transported.  
         [0030]    System  400  includes a module  402  which receives the photon measurement data. Module  402  may simply be a port through which the data enters the system, or it may be a more complex module for conditioning and acquiring data signals. The data is stored in a storage module  404 , which maintains a copy of the data for access by the processing module  406 . Storage module may take the form of dynamic random access memory (DRAM), static random access memory (SRAM), Flash ROM (read only memory), optical storage media, magnetic storage media, or any combination thereof. Note that storage module  404  (or a portion thereof) may be removable so as to be portable between systems. Storage module  404  may further maintain a copy of software for analyzing the data as described further below.  
         [0031]    Processing module  406  accepts instructions from a user via input module  408 , and responsively operates on the data from storage module  404  to extract circuit timing parameters. The data and/or parameters may be displayed in some form to the user by display module  410 .  
         [0032]    The photon detection data (or a subset thereof) may be displayed in a variety of forms. Of particular utility to the present application is form shown in FIG. 5. In FIG. 5, the number of photons detected at a given pixel or at a given set of pixel is plotted as a function of relative time. The time axis is divided into intervals, and the data shows the number of photons detected in each interval. Of particular interest is the precise interval between peaks, although other parameters may also be of interest, such as relative amplitude, peak width, and rise and fall times.  
         [0033]    Although a variety of heuristic techniques may be used to estimate these values, a preferred method of extracting these parameters involves Maximum Likelihood Estimation (MLE). MLE is appealing for a number of reasons. The parameters estimates obtained using MLE are asymptotically consistent, i.e., they converge to a determinable value as the sample size increases. The parameters estimates are also asymptotically unbiased, i.e. the determinable value is the true value of the parameters. Further, the variance of the parameter estimates approaches the theoretical minimum for all possible estimation methods for large sample sizes. In addition, the MLE estimates are asymptotically fully efficient, i.e. the mean square error of the estimates approach the theoretical minimum for large sample sizes. Still further, the estimates, when normalized, follow (approximately) a multivariate normal distribution, which allows for construction of approximate confidence bounds. Finally, the estimation equations are generally well-behaved and amenable to solution by numerical means.  
         [0034]    The MLE approach involves first defining a parameterized model, then constructing a likelihood equation for the model, and finally determining the parameters that maximize the likelihood equation for a given data set.  
         [0035]    Model for Photon Arrival Times (Without Jitter)  
         [0036]    Consider an imaging session where photons are detected over some set S of pixels. As relative time evolves from the triggered start time T start =0 to stop time T stop , photons are generated quantum-mechanically by hot carriers in the channel(s) imaged by pixel set S. Let 0&lt;X 1 &lt;X 2  &lt; . . . &lt;X r &lt;T stop  be the true sequence of successive arrival times for the photons so generated that happen to impact the photomultiplier and are successfully detected. “True” here means the arrival times before time jitter is taken into account. (Time jitter is considered in a separate model.)  
         [0037]    Let N(t,t+h) denote the number of these arrival times in the interval (t,t+h] for t≧0 and h&gt;0. Photons in our apparatus are generated by quantum-mechanical events and are relatively rare (˜1 detected per 1 billion switching events per um of NMOS transistor width for 0.18 um bulk CMOS at 1.65V). Available single photon imaging cameras can have a detection deadtime of 10 μsec or less, so the fraction of potentially detectable photons that are lost due to deadtime is negligible except when imaging very wide devices (&gt;10 mm). Thus, it is reasonable to assume that the point process described by N(t,t+h) satisfies the postulates for a Non-Homogeneous Poisson Process (NHPP):  
           Pr{N ( t,t+h )=0}=1−λ( t ) h+o ( h ), and  (1)  
           Pr{N ( t,t+h )=1}=λ( t ) h+o ( h ) as h→0  (2)  
         [0038]    where λ(t)&gt;O is the intensity function of the process and o(h)/h→0 as h→0. For our purposes, o(h) is set to zero, and λ(t) represents the photon intensity at the detector. It is the modeling and analysis of the intensity function that is of interest, as it reflects the current flow through the imaged channel(s).  
         [0039]    The photon intensity λ(t) varies as a function of time. This is evident from the data in FIG. 5. FIG. 6 shows a parsimonious model of the intensity function that includes as parameters the values of interest: start time f, rise time r, height h, and fall time g. The simplest model having these four parameters is a scalene triangle. The triangular model (t) can be expressed in the following way:  
                 γ        (       s   ;   0     ,   r   ,   h   ,   g     )       =       h   rg          (       gs   ·     δ        (       r   -   s     ,   s     )         +       r        (     r   +   g   -   s     )       ·     δ        (       r   +   g   -   s     ,     s   -   r       )           )              
            γ        (       s   ;   f     ,   r   ,   h   ,   g     )       =     γ        (         s   -   f     ;   0     ,   r   ,   h   ,   g     )                 (   3   )                               
 
         [0040]    where δ(x,y) is a function that equals one when both arguments are greater than or equal to zero, and that equals zero otherwise. There is a background photon arrival rate evident in FIG. 5 (due to noise in the detector) which can modeled as a constant intensity α over a fixed interval of observation (a,b]. The combined intensity function is then:  
               λ        (   t   )       =       α   ·     δ        (       t   -   a     ,     b   -   t       )         +       ∑     i   =     1                 m                      γ        (       t   ;     f   i       ,     r   i     ,     h   i     ,     g   i       )                   (   4   )                               
 
         [0041]    where m is the number of switching events observed in the timing window by the set of pixels S. For example, in FIG. 5, there are four switching events. Some of the parameters may be assumed to be equal, e.g. in FIG. 5, it might be assumed that r 3 =r 1 , r 4 =r 2 , h 3 =h 1 , h 4 =h 2 , g 3 =g 1 , and g 4 =g 2  Such assumptions would provide a more parsimonious model for parameter extraction and simplify the confidence interval determination process.  
         [0042]    In general for a NHPP, the probability of having r photons detected within a given time interval (s,t] is:  
               Pr        [       N        (     s   ,   t     )       =   r     ]       =         exp        (     -     [       μ        (   t   )       -     μ        (   s   )         ]       )                [       μ        (   t   )       -     μ        (   s   )         ]     r       r   !          r     =     0   ,   1   ,              …               (   5   )                               
 
         [0043]    where  
         μ        (   s   )       =       ∫     -   ∞     s            λ        (   x   )                 x     .                               
 
         [0044]    is the cumulative intensity function of the NHPP. If a sample of arrival times x 1 , x 2 , . . . , x n  is observed over the interval (a,b] in a particular PICA imaging session, the likelihood function A for this sample is given by:  
               Λ        (       x   1     ,   …              ,       x   n     ;     λ        (   x   )           )       =       exp        (     -       ∫   a   b            λ        (   u   )               u           )              ∏     i   =   1     n                       λ        (     x   i     )       .                 (   6   )                               
 
         [0045]    Refer to explanation starting on p. 495 of Lawless, J. F.,  Statistical Models and Methods for Lifetime Data,  John Wiley &amp; Sons, New York, 1982, which is hereby incorporated by reference.  
         [0046]    However, where the data is reported in the preferred form, the sample of recorded arrival times x 1 , x 2 , . . . , x n , is reported as histogram, i.e., a sequence of photon counts c 1 ,c 2 , . . . , c k  over a sequence of k successive time bins covering a time interval (a,b]. (The preferred apparatus is designed to support time bins down to 2.5 ps wide.) For any partition of the time interval (a,b], say α=α 0 &lt;α 1 &lt; . . . &lt;α n =b, the likelihood of observing c 1  photons in (α 0 ,α 1 ], c 2  photons in (a 1 ,a 2 ], . . . , c n  photons in (α n−1 ,α n ], is given by:  
               Λ        (       c   1     ,   …              ,       c   n     ;     λ        (   s   )           )       =       exp        (     -     [       μ        (   b   )       -     μ        (   a   )         ]       )              ∏     i   =   1     n                         [       μ        (     a   i     )       -     μ        (     a     i   -   1       )         ]       c   i           c   i     !                   (   7   )                               
 
         [0047]    where  
         μ        (   s   )       =       ∫     -   ∞     s            λ        (   x   )                 x     .                               
 
         [0048]    Though in practical situations there is a known, finite time interval (a,b] outside of which λ(t)=0 and there is zero probability of observing an arrival time, equations (6) and (7) above hold true even if the time interval (a,b] is expanded to (−∞,∞), so long as ∫ a   b  λ(x)dx remains finite.  
         [0049]    Note that the likelihood equation is a function of the intensity function, which in turn is a function of the parameters α and ƒ i ,r i ,h i ,g i , i=1 . . . m. The MLE approach estimates these parameters by determining which parameter values maximize the likelihood function. A variety of numerical maximization techniques are known and may be used to identify these parameter values, typically by (equivalently) minimizing the negative of the logarithm of the likelihood function (NLL).  
         [0050]    A preferred numerical technique for minimizing the NLL is known as a Line Search Method. In this method a descent direction (in the multiple dimension space that contains the parameter values sought) is determined on each iteration so that the next iterate (vector of parameter values) is chosen to be a point in this descent direction that decreases the NLL. Many implementations of this method can utilize the Hessian (matrix of second derivatives with respect to the vector of parameters) of the NLL. A closed form for the Hessian can be computed for the NLL corresponding to Equation 7 above. Line search methods typically have good convergence properties even when the starting values given for the parameters are not close to the ultimate parameter values that minimize the NLL (the MLE parameter values).  
         [0051]    In any event, once the MLE parameters have been calculated, the intensity function can then be plotted as shown in FIG. 7. The parameter values may also be used to determine propagation delay, duty cycle, clock skew, and timing jitter. As just one example, if each peak in FIG. 7 is from a different transistor, and the time between the first and third peaks is solely due to signal propagation delay, the propagation delay can be determined as ƒ 3 −ƒ 1  (or if peak-to-peak delay is preferred, ƒ 3 +r 3 −ƒ 1 −r 1 ).  
         [0052]    Model for Photon Arrival Times (With Jitter)  
         [0053]    The above analysis has been performed without considering the effects of timing jitter. In fact, jitter is an important consideration since the data measurement is taken over a “wall-clock” time interval of several hours or more. Time measurements made on the picosecond scale, even though they are made in relative time, are nevertheless sensitive to jitter.  
         [0054]    The measured photon arrival times may be subject to jitter contributed by several sources including: the photon detection process, the timing measurement electronics, the trigger signal, the on-chip clocking signal, and switching fluctuations due to variations in the power supply voltage. For the true photon arrival times X 1 , X 2 , . . . , X r , assume the totality of time jitter affecting them is captured by the sequence of independent (of each other and of the X i &#39;s) random jitter times Z 1 , Z 2 , . . . Z r , with a common Gaussian distribution given by:  
                 Pr        {       Z   i     ≤   t     }       =         ∫     -   ∞     t            σ     -   1            ϕ        (     x   /   σ     )               x         =     Φ        (     t   /   σ     )           ,     
        where           (   8   )                 ϕ        (   x   )       =       1       2      π                       -     x   2       /   2       .               (   9   )                               
 
         [0055]    In this model what is actually recorded as arrival times are the times T i =X i +Z i , i=1, . . . , r.  
         [0056]    For simplicity, let us initially assume that the jittered arrival times are observed over the interval (−∞,∞). In particular, assume that an arrival time near zero that receives a negative jitter component can be observed even though it may be less than zero. Then by Theorem 1 (derived in the Appendix) the likelihood function for a sample of jittered arrival times t 1 ,t 2 , . . . ,t n , in the interval (−∞,∞) is given by:  
                 Λ        (       t   1     ,   …              ,       t   n     ;     λ        (   t   )         ,   σ     )       =       exp        (     -       ∫     -   ∞     ∞            λ        (   t   )               t           )              ∏     i   =   1     n                     λ        (       t   i     ,   σ     )             ,           (   10   )                               
 
         [0057]    where the “jittered” intensity function is given by the convolution of the original intensity function with the Gaussian pulse:  
               λ        (     s   ,   σ     )       =       ∫     -   ∞     ∞            λ        (     s   -   z     )            σ     -   1            ϕ        (     z   /   σ     )                 z     .                 (   11   )                               
 
         [0058]    The likelihood function in equation (10) is equivalent to the likelihood function that would result from directly assuming a NHPP model with intensity function λ(s, σ) with respect to the recorded photon arrival times T i , if it can be shown that the intensity function integrals are equal. The proof is as follows:  
                       ∫     s   =     -   ∞       ∞            λ        (     s   ,   σ     )               s         =       ∫     s   =     -   ∞       ∞            ∫     z   =     -   ∞       ∞            λ        (     s   -   z     )            σ     -   1            ϕ        (     z   /   σ     )               z                        s                       =       ∫     z   =     -   ∞       ∞            (       ∫     s   =     -   ∞       ∞            λ        (   s   )                          s         )          σ     -   1            ϕ        (     z   /   σ     )               z                     =       ∫     s   =     -   ∞       ∞            λ        (   s   )                          s                       (   12   )                               
 
         [0059]    where the limit interchange is allowed since  
         ∫     -   ∞     ∞            λ        (   s   )               s                             
 
         [0060]    is assumed to be finite.  
         [0061]    As before, the likelihood function can be restated for data reported in the histogram format:  
                 Λ        (       c   1     ,   …              ,       c   n     ;     λ        (   s   )         ,   σ     )       =       exp        (     -     [         μ   σ          (   b   )       -       μ   σ          (   a   )         ]       )              ∏     i   =   1     n              [         μ   σ          (     a   i     )       -       μ   σ          (     a     i   -   1       )         ]       c   i           c   i     !             ,           (   13   )                               
 
         [0062]    where  
           μ   σ          (   s   )       =       ∫     -   ∞     s            λ        (     x   ,   σ     )                 x     .                               
 
         [0063]    This cumulative photon intensity function can be stated in the following form (for a single peak, ignoring the constant background rate for the present):  
                 μ   σ          (       s   ;   0     ,   r   ,   h   ,   g     )       =              ∫     -   ∞     s            ∫     -   ∞     ∞            λ        (         s   -   z     ;   0     ,   r   ,   h   ,   g     )            σ     -   1            ϕ        (     z   σ     )               z                        x                   (   14   )                                  =                h     2      rg          gs                 σ                   ϕ        (     s   ϕ     )         -       (     r   +   g     )          (     s   -   r     )        σ                   ϕ        (       s   -   r     σ     )         +                                r        (     s   -     (     r   +   g     )       )          σ                   ϕ        (       s   -     (     r   +   g     )       σ     )         +                              g        (       s   2     +     σ   2       )            Φ        (     s   σ     )         -       (     r   +   g     )          (         (     s   -   r     )     2     +     σ   2       )          Φ        (       s   -   r     σ     )         +                            r        (         (     s   -     (     r   +   g     )       )     2     +     σ   2       )            Φ        (       s   -     (     r   +   g     )       σ     )                       (   15   )                               
 μ σ ( s;f,r,h,g )=μ σ ( s−f; 0 ,r,h,g )  (16)  
         [0064]    The complete cumulative photon intensity function with jitter, for multiple peaks and a constant background rate a over (a,b] is:  
                       μ   σ          (   s   )       =            α   (       σ                   ϕ        (       s   -   a     σ     )         -     σ                   ϕ        (       b   -   s     σ     )         +       (     b   -   s     )          Φ        (       s   -   b     σ     )         +                                      (     s   -   a     )          Φ        (       s   -   a     σ     )         )     +   ∑       i   =     1      …                 m                μ   σ          (       s   ;     f   i       ,     r   i     ,     h   i     ,     g   i       )                       (   17   )                               
 
         [0065]    Accounting for jitter adds one parameter to the model (standard deviation a), so for the triangular intensity function model with jitter, the parameter list is α, σ, ƒ i ,r i ,h i ,g i , i=1 . . . m. (As before, we assume that the number of peaks m is determinable from visual inspection or some other source). The parameter values that maximize the likelihood function are calculated as before, and the jittered intensity function is plotted in FIG. 8.  
         [0066]    Likelihood equation (13) applies, and a numerical maximization method may be employed to extract the parameters.  
         [0067]    Confidence Interval Determination  
         [0068]    In the following, we denote the set of parameters as a vector θ, and the set of MLE estimates of the parameter values as {circumflex over (θ)}. Once the MLE estimates {circumflex over (θ)} have been found, a confidence interval (or region) can be determined for any subset of the parameters θ 1 , where θ is expressible as (θ 1 ,θ 2 ). While various ways exist, a convenient approximation of the confidence interval can be calculated from the profile likelihood. Ignoring for the moment the likelihood function&#39;s dependence on the data, and focusing on the dependence on the parameters θ=(θ 1 ,θ 2 ), the profile likelihood of a given parameter subset θ 1  is defined as:  
               R        (     θ   1     )       =       max     θ   2            [       Λ        (     (       θ   1     ,     θ   2       )     )         Λ        (     θ   ^     )         ]               (   18   )                               
 
         [0069]    The approximate 100(1−ε)% confidence region for the given parameter subset θ 1  is the set of all values of θ 1  such that:  
           R (θ 1 )&gt; exp (−χ (1−ε,k)   2 /2)  (19)  
         [0070]    where χ (1−ε,k)   2  denotes the 100(1−ε)% quantile of a Chi-square distribution with k degrees of freedom, where k is the number of parameters in subset θ 1 .  
         [0071]    Another feature of MLE that is put to use in this approach is the determination of a confidence interval for a function of the parameters. Where g(θ)=β is one-to-one transformation of θ into a new parameter set βand g(θ) has continuous first partial derivatives with respect to θ, then the MLE of {circumflex over (β)}=g({circumflex over (θ)}). Moreover, the inverse relation g −1 [β]=θ holds and  
               R                   (     β   1     )       =       max     β   2            [       Λ        (       g     -   1            [     (       β   1     ,     β   2       )     ]       )         Λ        (     θ   ^     )         ]               (   20   )                               
 
         [0072]    for any subset of parameters β 1 , where β=(β 1 , β 2 ).  
         [0073]    As just one example, consider (for simplicity) modeling only the first two peaks in FIG. 7 with parameters θ=&lt;ƒ 1 ,r 1 ,h 1 ,g 1 ,f 2 , r 2 ,h 2 ,g 2 ,α,σ&gt; and seeking to estimate the peak-to-peak delay. Define the one-to-one transformation g as post multiplication by the full-rank matrix  
             A   =     [           -   1         1       1       0       0       0       0       0       0       0             -   1         0       1       0       0       0       0       0       0       0           0       0       0       1       0       0       0       0       0       0           0       0       0       0       1       0       0       0       0       0           1       0       0       0       0       0       0       0       0       0           1       0       0       0       0       1       0       0       0       0           0       0       0       0       0       0       1       0       0       0           0       0       0       0       0       0       0       1       0       0           0       0       0       0       0       0       0       0       1       0           0       0       0       0       0       0       0       0       0       1         ]             (   21   )                               
 
         [0074]    so that β=g(θ)=θ·A=(β 1 ,β 2 ), with  
         β 1 =&lt;ƒ 2   +r   2 −ƒ 1   −r   1 &gt; and β 2 =&lt;ƒ 1 ,ƒ 1   +r   1   ,h   1   ,g   1   ,r   2   ,h   2   ,g   2 ,α,σ&gt;  (22)  
         Then  
         R(β 1 )&gt;exp(−χ 1−ε,1   2 /2)  (23)  
         [0075]    defines an approximate 100(1−ε)% confidence interval for the peak-to-peak delay ƒ 2 +r 2 −ƒ 1 −r 1 .  
         [0076]    Extraction Method  
         [0077]    Turning now to FIG. 9, the preferred parameter extraction method is summarized in the form of a flowchart. The software executed by processing module  406  preferably performs most, but not necessarily all, of the actions in this flowchart. In block  902 , a parameterized intensity function model is chosen. A wide number of models may be chosen, but the preferred models include the triangular model over a constant background intensity, with and without jitter. In block  904 , the likelihood equations for the model are constructed, and a solution procedure is incorporated into the program code. In block  906 , the photon arrival data is obtained. Initial processing may be performed on the data, e.g., to determine an average number of background noise counts and to remove these from the data.  
         [0078]    In block  908 , the number of switching transition peaks is determined. There are a number of ways this may be accomplished. In one approach, the photon arrival data is displayed to a user, who visually determines the number of peaks in the data. In a second approach, the processing module  406  executes a peak-finding algorithm to determine the number of peaks in the data. In a third approach, the user determines an expected number of peaks from knowledge of the imaged area of the integrated circuit and knowledge of the test pattern. Other approaches are possible and may be used.  
         [0079]    In block  910 , extracts values for the model parameters from the photon arrival data. This is accomplished by using the solution routine created in block  904 . In block  912 , the desired circuit switching parameters are then determined using the model parameters and the model.  
         [0080]    Among the potential advantages that may be realized by the disclosed systems and methods are the potential for automation. The reliance on user interaction can be minimized, even to the point where the user merely powers-up the system. Further, the disclosed systems and methods can be applied to identify parameters of faint signals and/or signals that are acquired over a relatively short time.  
         [0081]    Numerous variations and modifications will become apparent to those skilled in the art once the above disclosure is fully appreciated. It is intended that the following claims be interpreted to embrace all such variations and modifications.  
       APPENDIX  
       [0082]    Theorem 1: Let t 1 ,t 2 , . . . ,t n  be a sample of photon arrival times observed over the interval (−∞,∞) from an underl;ying NHPP with intensity function λ(t) that has undergone N(0,σ) time jitter. If there exists a finite bound M such that  
         λ( t )≦ M  for all t∈(−∞,∞) and              ∫     -   ∞     ∞            λ        (   u   )                          u         ≤   M     ,                           
         [0083]    then the Likelihood function Λ for the photon arrival times is given by  
                 Λ        (       t   1     ,   …              ,       t   n     ;     λ        (   s   )         ,   σ     )       =       exp        (     -       ∫     -   ∞     ∞            λ        (   u   )               u           )              ∏     i   =   1     n          λ        (       t   i     ,   σ     )             ,                 where                   λ        (     s   ,   σ     )         =       ∫     -   ∞     ∞            λ        (     s   -   z     )            σ     -   1            ϕ        (     z   σ     )                 z     .                                     
 
         [0084]    Proof: For brievity of notation let  
         f        (   z   )       =       σ     -   1              ϕ        (     z   σ     )       .                             
 
         [0085]    Now, by definition of a Likelihood function for continuous random variables,  
         Λ        (       t   1     ,   …              ,       t   n     ;     λ        (   t   )         ,   σ     )       =       Lim     ɛ   →   0                     Pr          {                T   i     -     t   i            &lt;   ɛ     ,     i   =   1     ,   …              ,   n     }     ·         (     2                 ɛ     )       -   n       .                               
 
         [0086]    By the definition of the T i &#39;s, we have  
         Λ        (       t   1     ,   …              ,       t   n     ;     λ        (   t   )         ,   σ     )       =       Lim     ɛ   →   0                     Pr          {                X   i     +     Z   i     -     t   i            &lt;   ɛ     ,     i   =   1     ,   …              ,   n     }     ·         (     2                 ɛ     )       -   n       .                               
 
         [0087]    By the definition of condition probability, we have  
       Λ   =       Lim     ɛ   →   0                         ∫     -   ∞     ∞                     ⋯                     ∫     -   ∞     ∞          Pr        {                X   1     +     z   1     -     t   1            &lt;   ɛ     ,   …              ,                  X   n     +     z   n     -     t   n            &lt;   ɛ     |     Z   1       =     z   1       ,   …              ,       Z   n     =     z   n         }          f        (     z   1     )                     ⋯                   f        (     z   n     )                 z   1                     ⋯                          z   n       ·         (     2                 ɛ     )       -   n       .                                       
 
         [0088]    Using the limit interchange allowed by the uniform boundedness of λ(t) and f(z), we reach  
       Λ   =                  ∫     -   ∞     ∞                     ⋯                     ∫     -   ∞     ∞                       Lim     ɛ   →   0                     Pr          {                X   1     +     z   1     -     t   1            &lt;   ɛ     ,   …              ,                  X   n     +     z   n     -     t   n            &lt;   ɛ     |     Z   1       =     z   1       ,   …              ,       Z   n     =     z   n         }     ·       (     2                 ɛ     )       -   n              f        (     z   1     )                     ⋯                   f        (     z   n     )                 z   1                     ⋯                        z   n                                     
 
         [0089]    Then, by the definition of the Likelihood function for the underlying NHPP,  
       Λ   =                  ∫     -   ∞     ∞                     ⋯                     ∫     -   ∞     ∞                       λ        (       t   1     -     z   1       )                     ⋯                   λ        (       t   n     -     z   n       )                       exp        (     -       ∫     -   ∞     ∞            λ        (   u   )               u           )                       f        (     z   1     )                     ⋯                   f        (     z   n     )                 z   1                     ⋯                          z   n       .                                   
 
         [0090]    Finally, we rearrange  
       Λ   =       exp        (     -       ∫     -   ∞     ∞            λ        (   u   )               u           )              ∏     i   =   1     n            ∫     -   ∞     ∞            λ        (       t   i     -   z     )            f        (   z   )               z