Method of using control models for data compression

BACKGROUND OF THE INVENTION

1. Field of the Invention

This invention relates generally to semiconductor fabrication technology, and, more particularly, to efficiently storing data relating to semiconductor fabrication processes.

2. Description of the Related Art

There is a constant drive within the semiconductor industry to increase the quality, reliability and throughput of integrated circuit devices, e.g., microprocessors, memory devices, and the like. This drive is fueled by consumer demands for higher quality computers and electronic devices that operate more reliably. These demands have resulted in a continual improvement in the manufacture of semiconductor devices, e.g., transistors, as well as in the manufacture of integrated circuit devices incorporating such transistors. Additionally, reducing defects in the manufacture of the components of a typical transistor also lowers the overall cost per transistor as well as the cost of integrated circuit devices incorporating such transistors.

The technologies underlying semiconductor processing tools have attracted increased attention over the last several years, resulting in substantial refinements. However, despite the advances made in this area, many of the processing tools that are currently commercially available suffer certain deficiencies. In particular, such tools often lack advanced process data monitoring capabilities, such as the ability to provide historical parametric data in a user-friendly format, as well as event logging, real-time graphical display of both current processing parameters and the processing parameters of the entire run, and remote, i.e., local site and worldwide, monitoring. These deficiencies can engender nonoptimal control of critical processing parameters, such as throughput accuracy, stability and repeatability, processing temperatures, mechanical tool parameters, and the like. This variability manifests itself as within-run disparities, run-to-run disparities and tool-to-tool disparities that can propagate into deviations in product quality and performance.

An additional set of problems is posed by the sheer amount of processing data available to monitor. For example, in monitoring an etching process using optical emission spectroscopy (OES), problems are posed by the sheer number of OES frequencies or wavelengths available to monitor. The monitoring and/or control of each of the many processing steps in a wafer fabrication facility typically generates a large amount of processing data. For example, a data file for each wafer monitored, in an etching process using optical emission spectroscopy (OES), may be as large as 2-3 megabytes (MB), and each etcher can typically process about 500-700 wafers per day. Conventional storage methods would require over a gigabyte (GB) of storage space per etcher per day and over 365 GB per etcher per year. Further, the raw OES data generated in such monitoring is typically “noisy” and unenhanced.

The present invention is directed to overcoming, or at least reducing the effects of, one or more of the problems set forth above.

SUMMARY OF THE INVENTION

In one aspect of the present invention, a method is provided for compressing data. The method comprises collecting processing data representative of a process. The method further comprises modeling the collected processing data using a control model. The method also comprises applying the control model to compress the collected processing data.

In another aspect of the present invention, a computer-readable, program storage device is provided, encoded with instructions that, when executed by a computer, perform a method, the method comprising collecting processing data representative of a process. The method further comprises modeling the collected processing data using a control model. The method also comprises applying the control model to compress the collected processing data.

In yet another aspect of the present invention, a computer programmed to perform a method is provided, the method comprising collecting processing data representative of a process. The method further comprises modeling the collected processing data using a control model. The method also comprises applying the control model to compress the collected processing data.

BRIEF DESCRIPTION OF THE DRAWINGS

The invention may be understood by reference to the following description taken in conjunction with the accompanying drawings, in which the leftmost significant digit(s) in the reference numerals denote(s) the first figure in which the respective reference numerals appear, and in which:

FIGS. 1-7

schematically illustrate a flow diagram for various embodiments of a method according to the present invention;

FIGS. 8-14

schematically illustrate a flow diagram for various alternative embodiments of a method according to the present invention;

FIGS. 15-21

schematically illustrate a flow diagram for yet other various embodiments of a method according to the present invention;

FIGS. 22 and 23

schematically illustrate first and second Principal Components for respective data sets;

FIG. 24

schematically illustrates OES spectrometer counts plotted against wavelengths;

FIG. 25

schematically illustrates a time trace of OES spectrometer counts at a particular wavelength;

FIG. 26

schematically illustrates representative mean-scaled spectrometer counts for OES traces of a contact hole etch plotted against wavelengths and time;

FIG. 27

schematically illustrates a time trace of Scores for the second Principal Component used to determine an etch endpoint;

FIGS. 28 and 29

schematically illustrate geometrically Principal Components Analysis for respective data sets;

FIG. 30

schematically illustrates a time trace of OES spectrometer counts at a particular wavelength and a reconstructed time trace of the OES spectrometer counts at the particular wavelength;

FIG. 31

schematically illustrates a method in accordance with the present invention;

FIG. 32

schematically illustrates workpieces being processed using a processing tool, using a plurality of control input signals, in accordance with the present invention; and

FIGS. 33-36

schematically illustrate geometrically polynomial least-squares fitting, in accordance with the present invention.

While the invention is susceptible to various modifications and alternative forms, specific embodiments thereof have been shown by way of example in the drawings and are herein described in detail. It should be understood, however, that the description herein of specific embodiments is not intended to limit the invention to the particular forms disclosed, but on the contrary, the intention is to cover all modifications, equivalents, and alternatives falling within the spirit and scope of the invention as defined by the appended claims.

DETAILED DESCRIPTION OF SPECIFIC EMBODIMENTS

Illustrative embodiments of the invention are described below. In the interest of clarity, not all features of an actual implementation are described in this specification. It will of course be appreciated that in the development of any such actual embodiment, numerous implementation-specific decisions must be made to achieve the developers' specific goals, such as compliance with system-related and business-related constraints, which will vary from one implementation to another. Moreover, it will be appreciated that such a development effort might be complex and time-consuming, but would nevertheless be a routine undertaking for those of ordinary skill in the art having the benefit of this disclosure.

In various illustrative embodiments, control systems may be provided that use modeling of data in their functionality. Control models may deal with large volumes of data. This large volume of data may be fit by the control model so that a much smaller data set of control model parameters is created. When the much smaller data set of control model parameters is put into the control model, the original much larger data set is accurately recreated (to within a specified tolerance and/or approximation). These control models are already created for control purposes, so the resulting control model parameters may be stored to sufficiently capture, in the compressed control model parameter data set, the original voluminous data set with much smaller storage requirements. For example, fault detection models, in particular, deal with large volumes of data. This large volume of data may be fit by the fault detection model so that a much smaller data set of fault detection model parameters is created. When the much smaller data set of fault detection model parameters is put into the fault detection model, the original much larger data set is accurately recreated (to within a specified tolerance and/or approximation).

In various alternative illustrative embodiments, control systems that use modeling of data in their functionality may be built. In these various alternative illustrative embodiments, building the control model comprises fitting the collected processing data using at least one of polynomial curve fitting, least-squares fitting, polynomial least-squares fitting, non-polynomial least-squares fitting, weighted least-squares fitting, weighted polynomial least-squares fitting, weighted non-polynomial least-squares fitting, Partial Least Squares (PLS), and Principal Components Analysis (PCA).

In various illustrative embodiments, samples may be collected for N+1 data points (x

i

,y

i

), where i=1, 2, . . . , N, N+1, and a polynomial of degree N,

P_{N} (x) = a_{0} + a_{1} x + a_{2} x^{2} + \dots + a_{k} x^{k} + \dots + a_{N} x^{N} = \sum_{k = 0}^{N} a_{k} x^{k},

may be fit to the N+1 data points (x

i

,y

i

). For example, 100 time data points (N=99) may be taken relating the pressure p, the flowrate f, and/or the temperature T, of a process gas to the thickness t of a process layer being deposited in a deposition processing step, resulting in respective sets of N+1 data points (p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

). For each of these sets of N+1 data points (p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

), the conventional raw storage requirements would be on the order of 2(N+1). By way of comparison, data compression, according to various illustrative embodiments of the present invention, using polynomial interpolation, for example, would reduce the storage requirements, for each of these sets of N+1 data points (p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

), to the storage requirements for the respective coefficients a

k

, for k=0, 1, . . . , N, a storage savings of a factor of 2 (equivalent to a data compression ratio of 2:1) for each of these sets of N+1 data points (p

i,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

). Polynomial interpolation is described, for example, in

Numerical Methods for Scientists and Engineers

, by R. W. Hamming, Dover Publications, New York, 1986, at pages 230-235. The requirement that the polynomial P

N

(x) pass through the N+1 data points (x

i

,y

i

) is

y_{i} = P_{N} (x_{i}) = \sum_{k = 0}^{N} a_{k} x_{i}^{k}, for i = 1, 2, \dots, N, N + 1,

a set of N+1 conditions. These N+1 conditions then completely determine the N+1 coefficients a

k

, for k=0, 1, . . . , N.

The determinant of the coefficients of the unknown coefficients a

k

is the

Vandermonde determinant : V_{N + 1} = &LeftBracketingBar; \begin{matrix} 1 & x_{1} & x_{1}^{2} & \dots & x_{1}^{N} \\ 1 & x_{2} & x_{2}^{2} & \dots & x_{2}^{N} \\ ⋮ & ⋮ & ⋮ & ⋰ & ⋮ \\ 1 & x_{N} & x_{N}^{2} & \dots & x_{N}^{N} \\ 1 & x_{N + 1} & x_{N + 1}^{2} & \dots & x_{N + 1}^{N} \end{matrix} &RightBracketingBar; = &LeftBracketingBar; x_{i}^{k} &RightBracketingBar;, where i = 1, 2, \dots, N + 1,

and k=0, 1, . . . , N. The Vandermonde determinant V

N+1

, considered as a function of the variables x

i

, V

N+1

=V

N+1

(x

1

x

2

, . . . , x

N

,x

N+1

), is clearly a polynomial in the variables x

i

, as may be seen by expanding out the determinant, and a count of the exponents shows that the degree of the polynomial is

0 + 1 + 2 + 3 + \dots + k + \dots + N = \sum_{k = 0}^{N} k = \frac{N (N + 1)}{2}

(for example, the diagonal term of the Vandermonde determinant V

N+1

is 1·x

2

·x

3

2

. . . x

N

N−1

·x

N+1

N

).

Now, if x

N+1

=x

j

, for j=1, 2, . . . , N, then the Vandermonde determinant V

N+1

=0, since any determinant with two identical rows vanishes, so the Vandermonde determinant V

N+1

must have the factors (x

N+1

−x

j

), for j=1, 2, . . . , N, corresponding to the N factors

\prod_{j = 1}^{N} (x_{N + 1} - x_{j}) .

Similarly, if x

N

=x

j

, for j=1, 2, . . . , N−1, then the Vandermonde determinant V

N+1

=0, so the Vandermonde determinant V

N+1

must also have the factors (x

N−x

j

), for j=1, 2, . . . , N−1, corresponding to the N−1 factors

\prod_{j = 1}^{N - 1} (x_{N} - x_{j}) .

Generally, if x

m

=x

j

, for j<m, where m=2, . . . , N, N+1, then the Vandermonde determinant V

N+1

=0, so the Vandermonde determinant V

N+1

must have all the factors (x

m−x

j

), for j<m, where m=2, . . . , N, N+1, corresponding to the factors

\prod_{m > j = 1}^{N + 1} (x_{m} - x_{j}) .

Altogether, this represents a polynomial of degree

N + (N - 1) + \dots + 2 + 1 = \sum_{k = 1}^{N} k = \frac{N (N + 1)}{2},

since, when m=N+1, for example, j may take on any of N values, j=1, 2, . . . , N, and when m=N, j may take on any of N−1 values, j=1, 2, . . . , N−1, and so forth (for example, when m=3, j may take only two values, j=1, 2, and when m=2, j may take only one value, j=1), which means that all the factors have been accounted for and all that remains is to find any multiplicative constant by which these two representations for the Vandermonde determinant V

N+1

might differ. As noted above, the diagonal term of the Vandermonde determinant V

N+1

is 1·x

2

·x

3

2

. . . x

N

N−1

·x

N+1

, and this may be compared to the term from the left-hand sides of the product of factors

\prod_{m > j = 1}^{N + 1} (x_{m} - x_{j}) = \prod_{j = 1}^{N} (x_{N + 1} - x_{j}) \prod_{j = 1}^{N - 1} (x_{N} - x_{j}) \dots \prod_{j = 1}^{2} (x_{3} - x_{j}) \prod_{j = 1}^{1} (x_{2} - x_{j}), x_{N + 1}^{N} \cdot x_{N}^{N - 1} \dots x_{3}^{2} \cdot x_{2},

which is identical, so the multiplicative constant is unity and the Vandermonde determinant

V_{N + 1} is V_{N + 1} (x_{1}, x_{2}, \dots x_{N + 1}) = &LeftBracketingBar; x_{i}^{k} &RightBracketingBar; = \prod_{m > j = 1}^{N + 1} (x_{m} - x_{j}) .

This factorization of the Vandermonde determinant V

N+1

shows that if x

i

≠x

j

, for i≠j, then the Vandermonde determinant V

N+1

cannot be zero, which means that it is always possible to solve for the unknown coefficients a

k

, since the Vandermonde determinant V

N+1

is the determinant of the coefficients of the unknown coefficients a

k

. Solving for the unknown coefficients a

k

, using determinants, for example, substituting the results into the polynomial of degree

N, P_{N} (x) = \sum_{k = 0}^{N} a_{k} x^{k},

and rearranging suitably gives the determinant equation

&LeftBracketingBar; \begin{matrix} y & 1 & x & x^{2} & \dots & x^{N} \\ y_{1} & 1 & x_{1} & x_{1}^{2} & \dots & x_{1}^{N} \\ y_{2} & 1 & x_{2} & x_{2}^{2} & \dots & x_{2}^{N} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋰ & ⋮ \\ y_{N} & 1 & x_{N} & x_{N}^{2} & \dots & x_{N}^{N} \\ y_{N + 1} & 1 & x_{N + 1} & x_{N + 1}^{2} & \dots & x_{N + 1}^{N} \end{matrix} &RightBracketingBar; = 0,

which is the solution to the polynomial fit. This may be seen directly as follows. Expanding this determinant by the elements of the top row, this is clearly a polynomial of degree N. The coefficient of the element y in the first row in the expansion of this determinant by the elements of the top row is none other than the Vandermonde determinant V

N+1

. In other words, the cofactor of the element y in the first row is, in fact, the Vandermonde determinant V

N+1

. Indeed, the cofactor of the nth element in the first row, where n=2, . . . , N+2, is the product of the coefficient an 2 in the polynomial expansion

y = P_{N} (x) = \sum_{k = 0}^{N} a_{k} x^{k}

with the Vandermonde determinant V

N+1

. Furthermore, if x and y take on any of the sample values x

i

and y

i

, for i=1, 2, . . . , N, N+1, then two rows of the determinant would be the same and the determinant must then vanish. Thus, the requirement that the polynomial y=P

N

(x) pass through the N+1 data points (x

i

,y

i

),

y = P_{N} (x) = \sum_{k = 0}^{N} a_{k} x_{i}^{k}, for i = 1, 2, \dots, N, N + 1,

is satisfied.

For example, a quadratic curve may be found that goes through the sample data set (−1,a), (0,b), and (1,c). The three equations are P

2

(−1)=a=

0

−a

1

+a

2

, P

2

(0)=b=a

0

, and P

2

(1)=c=a

0

+a

1

+a

2

, which imply that b=a

0

, c−a=2a

1

, and c+a−2b=2a

2

, so that

y (x) = P_{2} (x) = b + \frac{c - a}{2} x + \frac{c + a - 2 b}{2} x^{2},

which is also the result of expanding

&LeftBracketingBar; \begin{matrix} y & 1 & x & x^{2} \\ a & 1 & - 1 & 1 \\ b & 1 & 0 & 0 \\ c & 1 & 1 & 1 \end{matrix} &RightBracketingBar; = 0 = y &LeftBracketingBar; \begin{matrix} 1 & - 1 & 1 \\ 1 & 0 & 0 \\ 1 & 1 & 1 \end{matrix} &RightBracketingBar; - 1 &LeftBracketingBar; \begin{matrix} a & - 1 & 1 \\ b & 0 & 0 \\ c & 1 & 1 \end{matrix} &RightBracketingBar; + x &LeftBracketingBar; \begin{matrix} a & 1 & 1 \\ b & 1 & 0 \\ c & 1 & 1 \end{matrix} &RightBracketingBar; - x^{2} &LeftBracketingBar; \begin{matrix} a & 1 & - 1 \\ b & 1 & 0 \\ c & 1 & 1 \end{matrix} &RightBracketingBar;,

the coefficient of y being the respective Vandermonde determinant V

3

=2.

Similarly, a quartic curve may be found that goes through the sample data set (−2, a), (−1,b), (0,c), (1,b), and (2,a). The five equations are P

4

(−2)=a=a

0

−2a

1

+4a

2

−8a

3

+16a

4

, P

4

(−1)=b=a

0

−a

1

+a

2

−a

3

+a

4

, P

4

(0)=c=a

0

, P

4

(1)=b=a

0

+a

1

+a

2

+a

4

, and P

4

(2)=a=a

0

+2a

1

+4a

2

+8 a

3

+16a

4

, which imply that c=a

0

, 0=a

1

=a

3

(which also follows from the symmetry of the data set), (a−c)−16(b−c)=−12a

2

, and (a−c)−4(b−c)=12a

4

, so that

y (x) = P_{4} (x) = c - \frac{a - 16 b + 15 c}{12} x^{2} + \frac{a - 4 b + 3 c}{12} x^{4} .

In various alternative illustrative embodiments, samples may be collected for M data points (x

i

,y

i

), where i=1, 2, . . . , M, and a first degree polynomial (a straight line),

P_{1} (x) = a_{0} + a_{1} x = \sum_{k = 0}^{1} a_{k} x^{k},

P_{1} (x) = a_{0} + a_{1} x = \sum_{k = 0}^{1} a_{k} x^{k},

fit (in a least-squares sense) to the M data points (x

i

y

i

), only requires two parameters, a

0

and a

1

, which would reduce the storage requirements, for each of these sets of M data points (p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

), to the storage requirements for the respective coefficients a

k

, for k=0, 1, a storage savings of a factor of 2M/2 (equivalent to a data compression ratio of M:1) for each of these sets of M data points (p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

). Least-squares fitting is described, for example, in

Numerical Methods for Scientists and Engineers

, by R. W. Hamming, Dover Publications, New York, 1986, at pages 427-443.

The least-squares criterion may be used in situations where there is much more data available than parameters so that exact matching (to within round-off) is out of the question. Polynomials are most commonly used in least-squares matching, although any linear family of suitable functions may work as well. Suppose some quantity x is being measured by making M measurements x

i

, for i=1, 2, . . . , M, and suppose that the measurements x

i

, are related to the “true” quantity x by the relation x

i

=x+ε

i

, for i=1, 2, . . . , M, where the residuals ε

i

are regarded as noise. The principle of least-squares states that the best estimate ξ of the true value x is the number that minimizes the sum of the squares of the deviations of the data from their estimate

f (ξ) = \sum_{i = 1}^{M} ε_{i}^{2} = \sum_{i = 1}^{M} {(x_{i} - ξ)}^{2},

which is equivalent to the assumption that the average x

a

, where

x_{a} = \frac{1}{M} \sum_{i = 1}^{M} x_{i},

is the best estimate ξ of the true value x. This equivalence may be shown as follows. First, the principle of least-squares leads to the average x

a

. Regarding

f (ξ) = \sum_{i = 1}^{M} ε_{i}^{2} = \sum_{i = 1}^{M} {(x_{i} - ξ)}^{2}

as a function of the best estimate ξ, minimization with respect to the best estimate ξ may proceed by differentiation:

\frac{ⅆ f (ξ)}{ⅆ ξ} = - 2 \sum_{i = 1}^{M} (x_{i} - ξ) = 0,

which implies that

\sum_{i = 1}^{M} x_{i} - \sum_{i = 1}^{M} ξ = 0 = \sum_{i = 1}^{M} x_{i} - M ξ,

so that

ξ = \frac{1}{M} \sum_{i = 1}^{M} x_{i} = x_{a},

or, in other words, that the choice x

a

=ξ minimizes the sum of the squares of the residuals ε

i

. Noting also that

\frac{ⅆ^{2} f (ξ)}{ⅆ ξ^{2}} = 2 \sum_{i = 1}^{M} 1 = 2 M > 0,

the criterion for a minimum is established.

Conversely, if the average x

a

is picked as the best choice x

a

=ξ it can be shown that this choice, indeed, minimizes the sum of the squares of the residuals ε

i

. Set

f (x_{a}) = \sum_{i = 1}^{M} {(x_{i} - x_{a})}^{2} = \sum_{i = 1}^{M} x_{i}^{2} - 2 x_{a} \sum_{i = 1}^{M} x_{i} + \sum_{i = 1}^{M} x_{a}^{2} = \sum_{i = 1}^{M} x_{i}^{2} - 2 x_{a} M x_{a} + M x_{a}^{2} = \sum_{i = 1}^{M} x_{i}^{2} - M x_{a}^{2} .

If any other value x

b

is picked, then, plugging that other value x

b

into f(x) gives

f (x_{b}) = \sum_{i = 1}^{M} {(x_{i} - x_{b})}^{2} = \sum_{i = 1}^{M} x_{i}^{2} - 2 x_{b} \sum_{i = 1}^{M} x_{i} + \sum_{i = 1}^{M} x_{b}^{2} = \sum_{i = 1}^{M} x_{i}^{2} - 2 x_{b} M x_{a} + M x_{b}^{2} .

Subtracting f(x

a

) from f(x

b

) gives f(x

b

)−f(x

a

)=M[x

a

2

−2x

a

x

b

+x

b

2

]=M(x

a

−x

b

)

2

≧0, so that f(x

b

)≧f(x

a

), with equality if, and only if, x

b

=x

a

. In other words, the average x

a

, indeed, minimizes the sum of the squares of the residuals ε

i

. Thus, it has been shown that the principle of least-squares and the choice of the average as the best estimate are equivalent.

There may be other choices besides the least-squares choice. Again, suppose some quantity x is being measured by making M measurements x

i

, for i=1, 2, . . . , M, and suppose that the measurements x

i

, are related to the “true” quantity x by the relation x

i

=x+ε

i

, for i=1, 2, . . . , M, where the residuals ε

i

are regarded as noise. An alternative to the least-squares choice may be that another estimate χ of the true value x is the number that minimizes the sum of the absolute values of the deviations of the data from their estimate

f (χ) = \sum_{i = 1}^{M} &LeftBracketingBar; ε_{i} &RightBracketingBar; = \sum_{i = 1}^{M} &LeftBracketingBar; x_{i} - χ &RightBracketingBar;,

which is equivalent to the assumption that the median or middle value x

m

of the M measurements x

i

, for i=1, 2, . . . , M (if M is even, then average the two middle values), is the other estimate χ of the true value x. Suppose that there are an odd number M=2k+1 of measurements x

i

, for i=1, 2, . . . , M, and choose the median or middle value x

m

as the estimate χ of the true value x that minimizes the sum of the absolute values of the residuals ε

i

. Any upward shift in this value x

m

would increase the k terms |x

i

−x| that have x

i

below x

m

, and would decrease the k terms |x

i

−x| that have x

i

above x

m

, each by the same amount. However, the upward shift in this value x

m

would also increase the term |x

m

−x| and, thus, increase the sum of the absolute values of all the residuals ε

i

. Yet another choice, instead of minimizing the sum of the squares of the residuals ε

i

, would be to choose to minimize the maximum deviation, which leads to

\frac{x_{\min} + x_{\max}}{2} = x_{midrange},

the midrange estimate of the best value.

Returning to the various alternative illustrative embodiments in which samples may be collected for M data points (x

i

,y

i

), where i=1, 2, . . . , M, and a first degree polynomial (a straight line),

P_{1} (x) = a_{0} + a_{1} x = \sum_{k = 0}^{1} a_{k} x^{k},

may be fit (in a least-squares sense) to the M data points (x

i

, y

i

), there are two parameters, a

0

and a

1

, and a function F(a

0

,a

1

) that needs to be minimized as follows. The function F(a

0

,a

1

) is given by

F (a_{0}, a_{1}) = \sum_{i = 1}^{M} ε_{i}^{2} = \sum_{i = 1}^{M} {[P_{1} (x_{i}) - y_{i}]}^{2} = \sum_{i = 1}^{M} {[a_{0} + a_{1} x_{i} - y_{i}]}^{2}

and setting the partial derivatives of F(a

0

,a

1

) with respect to a

0

and a

1

equal to zero gives

\frac{\partial F (a_{0}, a_{1})}{\partial a_{0}} = 2 \sum_{i = 1}^{M} [a_{0} + a_{1} x_{i} - y_{i}] = 0

and

\frac{\partial F (a_{0}, a_{1})}{\partial a_{0}} = 2 \sum_{i = 1}^{M} [a_{0} + a_{1} x_{i} - y_{i}] x_{i} = 0,

respectively. Simplifying and rearranging gives

a_{0} M + a_{1} \sum_{i = 1}^{M} x_{i} = \sum_{i = 1}^{M} y_{i} and a_{0} \sum_{i = 1}^{M} x_{i} + a_{1} \sum_{i = 1}^{M} x_{i}^{2} = \sum_{i = 1}^{M} x_{i} y_{i},

respectively, where there are two equations for the two unknown parameters a

0

and a

1

, readily yielding a solution.

In various other alternative illustrative embodiments, samples may be collected for M data points (x

i

,y

i

), where i=1, 2, . . . , M, and a polynomial of degree N,

P_{N} (x) = a_{0} + a_{1} x + a_{2} x^{2} + \dots + a_{k} x^{k} + \dots + a_{N} x^{N} = \sum_{k = 0}^{N} a_{k} x^{k},

may be fit (in a least-squares sense) to the M data points (x

i

,y

i

). For example, 100 time data points (M=100) may be taken relating the pressure p, the flowrate f, and/or the temperature T, of a process gas to the thickness t of a process layer being deposited in a deposition processing step, resulting in the M data points (p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

), respectively. For each of these sets of M data points (p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

), the conventional raw storage requirements would be on the order of 2M. By way of comparison, data compression, according to various illustrative embodiments of the present invention, using a polynomial of degree N that is fit (in a least-squares sense) to the M data points (x

i

,y

i

), only requires N+1 parameters, a

k

, for k=0, 1, . . . , N, which would reduce the storage requirements, for each of these sets of M data points (p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

), to the storage requirements for the respective coefficients a

k

, for k=0, 1, . . . , N, a storage savings of a factor of 2M/(N+1), equivalent to a data compression ratio of 2M:(N+1), for each of these sets of M data points (p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

). In one illustrative embodiment, the degree N of the polynomial is at least 10 times smaller than M. For semiconductor manufacturing processes, data sufficient for fault detection purposes ranges from about 100 to about 10000 samples for the value of M. Requirements for polynomial order N are in the range of 1 to about 10.

The function F(a

0

,a

1

, . . . , a

N

) may be minimized as follows. The function F(a

0

,a

1

, . . . , a

N

) is given by

F (a_{0}, a_{1}, \dots, a_{N}) = \sum_{i = 1}^{M} ε_{i}^{2} = \sum_{i = 1}^{M} {[P_{N} (x_{i}) - y_{i}]}^{2}

and setting the partial derivatives of F(a

0

,a

1

, . . . , a

N

) with respect to a

j

, for j=0, 1, . . . , N, equal to zero gives

\begin{matrix} \frac{\partial F (a_{0}, a_{1}, \dots, a_{N})}{\partial a_{j}} = 2 \sum_{i = 1}^{M} [P_{N} (x_{i}) - y_{i}] x_{i}^{j} \\ = 2 \sum_{i = 1}^{M} [\sum_{k = 0}^{N} a_{k} x_{i}^{k} - y_{i}] x_{i}^{j} = 0, \end{matrix}

for j=0, 1, . . . , N, sine (x

i

)

j

is the coefficient of a

j

in the polynomial

P_{N} (x_{i}) = \sum_{k = 0}^{N} a_{k} x_{i}^{k} .

Simplifying and rearranging gives

\begin{matrix} \sum_{i = 1}^{M} [\sum_{k = 0}^{N} a_{k} x_{i}^{k}] x_{i}^{j} = \sum_{k = 0}^{N} a_{k} [\sum_{i = 1}^{M} x_{i}^{k + j}] \\ \equiv \sum_{k = 0}^{N} a_{k} S_{k + j} = \sum_{i = 1}^{M} x_{i}^{j} y_{i} \equiv T_{j}, \end{matrix}

for j=0, 1, . . . , N, where

\sum_{i = 1}^{M} x_{i}^{k + j} \equiv S_{k + j} and \sum_{i = 1}^{M} x_{i}^{j} y_{i} \equiv T_{j},

respectively. There are N+1 equations

\sum_{k = 0}^{N} a_{k} S_{k + j} = T_{j},

for j=0, 1, . . . , N, also known as the normal equations, for the N+1 unknown parameters a

k

, for k=0, 1, . . . , N, readily yielding a solution, provided that the determinant of the normal equations is not zero. This may be demonstrated by showing that the homogeneous equations

\sum_{k = 0}^{N} a_{k} S_{k + j} = 0

only have the trivial solution a

k

=0, for k=0, 1, . . . , N, which may be shown as follows. Multiply the jth homogeneous equation by a

j

and sum over all j, from j=0 to j=N,

\sum_{j = 0}^{N} a_{j} \sum_{k = 0}^{N} a_{k} S_{k + j} = \sum_{j = 0}^{N} a_{j} \sum_{k = 0}^{N} a_{k} \sum_{i = 1}^{M} x_{i}^{k} x_{i}^{j} = \sum_{i = 1}^{M} (\sum_{k = 0}^{N} a_{k} x_{i}^{k}) (\sum_{j = 0}^{N} a_{j} x_{i}^{j}) = \sum_{i = 1}^{M} {(P_{N} (x_{i}))}^{2} = 0,

which would imply that P

N

(x

i

)≡0, and, hence, that a

k

=0, for k=0, 1, . . . , N, the trivial solution. Therefore, the determinant of the normal equations is not zero, and the normal equations may be solved for the N+1 parameters a

k

, for k=0, 1, . . . , N, the coefficients of the least-squares polynomial of degree N,

P_{N} (x) = \sum_{k = 0}^{N} a_{k} x^{k},

that may be fit to the M data points (x

i

,y

i

).

Finding the least-squares polynomial of degree N,

P_{N} (x) = \sum_{k = 0}^{N} a_{k} x^{k},

that may be fit to the M data points (x

i

,y

i

) may not be easy when the degree N of the least-squares polynomial is very large. The N+1 normal equations

\sum_{k = 0}^{N} a_{k} S_{k + j} = T_{j},

for j=0, 1, . . . , N, for the N+1 unknown parameters a

k

, for k=0, 1, . . . , N, may not be easy to solve, for example, when the degree N of the least-squares polynomial is much greater than about 10. This may be demonstrated as follows. Suppose that the M data points (x

i

,y

i

) are more or less uniformly distributed in the interval 0≦x≦1, so that

S_{k + j} = \sum_{i = 0}^{M} x_{i}^{k + j} \approx M \int_{0}^{i} x^{k + j} ⅆ x = \frac{M}{k + j + 1} .

The resulting determinant for the normal equations is then approximately given by

&LeftBracketingBar; S_{k + j} &RightBracketingBar; \approx &LeftBracketingBar; \frac{M}{k + j + 1} &RightBracketingBar; \approx M^{N + 1} &LeftBracketingBar; \frac{1}{k + j + 1} &RightBracketingBar; = M^{N + 1} H_{N + 1},

for j,k=0, 1, . . . , N, where H

N

, for j,k=0, 1, . . . , N−1, is the Hilbert determinant of order N, which has the value

H_{N} = \frac{{[0! 1! 2! 3! \dots (N - 1)!]}^{3}}{N! (N + 1)! (N + 2)! \dots (2 N - 1)!}

that approaches zero very rapidly. For example,

H_{1} = \frac{{[0!]}^{3}}{1!} = &LeftBracketingBar; 1 &RightBracketingBar; = 1, H_{2} = \frac{{[0! 1!]}^{3}}{2! 3!} = &LeftBracketingBar; \begin{matrix} 1 & \frac{1}{2} \\ \frac{1}{2} & \frac{1}{3} \end{matrix} &RightBracketingBar; = \frac{1}{3} - \frac{1}{4} = \frac{1}{12}, and

H_{3} = \frac{{[0! 1! 2!]}^{3}}{3! 4! 5!} = &LeftBracketingBar; \begin{matrix} 1 & \frac{1}{2} & \frac{1}{3} \\ \frac{1}{2} & \frac{1}{3} & \frac{1}{4} \\ \frac{1}{3} & \frac{1}{4} & \frac{1}{5} \end{matrix} &RightBracketingBar;,

where

\begin{matrix} &LeftBracketingBar; \begin{matrix} 1 & \frac{1}{2} & \frac{1}{3} \\ \frac{1}{2} & \frac{1}{3} & \frac{1}{4} \\ \frac{1}{3} & \frac{1}{4} & \frac{1}{5} \end{matrix} &RightBracketingBar; = \frac{1}{3} &LeftBracketingBar; \begin{matrix} \frac{1}{2} & \frac{1}{3} \\ \frac{1}{3} & \frac{1}{4} \end{matrix} &RightBracketingBar; - \frac{1}{4} &LeftBracketingBar; \begin{matrix} 1 & \frac{1}{2} \\ \frac{1}{3} & \frac{1}{4} \end{matrix} &RightBracketingBar; + \frac{1}{5} &LeftBracketingBar; \begin{matrix} 1 & \frac{1}{2} \\ \frac{1}{2} & \frac{1}{3} \end{matrix} &RightBracketingBar; \\ = \frac{1}{3} \frac{1}{72} - \frac{1}{4} \frac{1}{12} + \frac{1}{5} \frac{1}{12} = \frac{1}{216} - \frac{1}{240} = \frac{1}{2160} . \end{matrix}

This suggests that the system of normal equations is ill-conditioned and, hence, difficult to solve when the degree N of the least-squares polynomial is very large. Sets of orthogonal polynomials tend to be better behaved.

In various other alternative illustrative embodiments, samples may be collected for M data points (x

i

,y

i

), where i=1, 2, . . . , M, and a linearly independent set of N+1 functions f

i

(x), for j=0, 1, 2, . . . , N,

\begin{matrix} y (x) = a_{0} f_{0} (x) + a_{1} f_{1} (x) + \dots + a_{j} f_{j} (x) + \dots + a_{N} f_{N} (x) \\ = \sum_{j = 0}^{N} a_{j} f_{j} (x), \end{matrix}

may be fit (in a non-polynomial least-squares sense) to the M data points (x

i

,y

i

). For example, 100 time data points (M=100) may be taken relating the pressure p, the flowrate f, and/or the temperature T, of a process gas to the thickness t of a process layer being deposited in a deposition processing step, resulting in the M data points (p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

), respectively. For each of these sets of M data points (p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

), the conventional raw storage requirements would be on the order of 2M. By way of comparison, data compression, according to various illustrative embodiments of the present invention, using a linearly independent set of N+1 functions f

j

(x) as a basis for a representation that is fit (in a non-polynomial least-squares sense) to the M data points (x

i

,y

i

), only requires N+1 parameters, a

k

, for k=0, 1, . . . , N, which would reduce the storage requirements, for each of these sets of M data points (p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

), to the storage requirements for the respective coefficients a

k

, for k=0, 1, . . . , N, a storage savings of a factor of 2M/(N+1), equivalent to a data compression ratio of 2M:(N+1), for each of these sets of M data points (p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

). In one illustrative embodiment, the number N+1 of the linearly independent set of basis functions f

j

(x) is at least 10 times smaller than M. For semiconductor manufacturing processes, data sufficient for fault detection purposes ranges from about 100 to about 10000 samples for the value of M Requirements for polynomial order N are in the range of 1 to about 10.

The function F(a

0

,a

1

, . . . , a

N

) may be minimized as follows. The function F(a

0

,a

1

, . . . , a

N

) is given by

F (a_{0}, a_{1}, \dots, a_{N}) = \sum_{i = 1}^{M} ε_{i}^{2} = \sum_{i = 1}^{M} {[y (x_{i}) - y_{i}]}^{2}

and setting the partial derivatives of F(a

0

,a

1

, . . . , a

N

) with respect to a

j

, for j=0, 1, . . . , N, equal to zero gives

\begin{matrix} \frac{\partial F (a_{0}, a_{1}, \dots, a_{N})}{\partial a_{j}} = 2 \sum_{i = 1}^{M} [y (x_{i}) - y_{i}] x_{i}^{j} \\ = 2 \sum_{i = 1}^{M} [\sum_{k = 0}^{N} a_{k} f_{k} (x_{i}) - y_{i}] f_{j} (x_{i}) = 0, \end{matrix}

for j 0, 1, . . . , N, since f

j

(x

i

) is the coefficient of a

j

in the representation

y (x_{i}) = \sum_{k = 0}^{N} a_{k} f_{k} (x_{i}) .

Simplifying gives

\begin{matrix} \sum_{i = 1}^{M} [\sum_{k = 0}^{N} a_{k} f_{k} (x_{i})] f_{j} (x_{i}) = \sum_{k = 0}^{N} a_{k} [\sum_{i = 1}^{M} f_{k} (x_{i}) f_{j} (x_{i})] \\ \equiv \sum_{k = 0}^{N} a_{k} S_{k, j} = \sum_{i = 1}^{M} f_{j} (x_{i}) y_{i} \equiv T_{j}, \end{matrix}

for j=0, 1, . . . , N, where

\sum_{i = 1}^{M} f_{k} (x_{i}) f_{j} (x_{i}) \equiv S_{k, j} and \sum_{i = 1}^{M} f_{j} (x_{i}) y_{i} \equiv T_{j},

respectively. There are N+1 equations

\sum_{k = 0}^{N} a_{k} S_{k, j} = T_{j},

\sum_{k = 0}^{N} a_{k} S_{k, j} = 0

only have the trivial solution a

k

=0, for k=0, 1, . . . , N, which may be shown as follows. Multiply the jth homogeneous equation by a

j

and sum over all j,

\begin{matrix} \sum_{j = 0}^{N} a_{j} \sum_{k = 0}^{N} a_{k} S_{k, j} = \sum_{j = 0}^{N} a_{j} \sum_{k = 0}^{N} a_{k} \sum_{i = 1}^{M} f_{k} (x_{i}) f_{j} (x_{i}) \\ = \sum_{i = 1}^{M} (\sum_{k = 0}^{N} a_{k} f_{k} (x_{i})) (\sum_{j = 0}^{N} a_{j} f_{j} (x_{i})), \end{matrix}

but

\sum_{i = 1}^{M} (\sum_{k = 0}^{N} a_{k} f_{k} (x_{i})) (\sum_{j = 0}^{N} a_{j} f_{j} (x_{i})) = \sum_{i = 1}^{M} {(y (x_{i}))}^{2} = 0,

which would imply that y(x

i

)≡0, and, hence, that a

k

=0, for k=0, 1, . . . , N, the trivial solution. Therefore, the determinant of the normal equations is not zero, and the normal equations may be solved for the N+1 parameters a

k

, for k 0, 1, . . . , N, the coefficients of the non-polynomial least-squares representation

y (x) = \sum_{j = 0}^{N} a_{j} f_{j} (x)

that may be fit to the M data points (x

i

,y

i

), using the linearly independent set of N+1 functions f

j

(x) as the basis for the non-polynomial least-squares representation

y (x) = \sum_{j = 0}^{N} a_{j} f_{j} (x) .

If the data points (x

i

,y

i

) are not equally reliable for all M, it may be desirable to weight the data by using non-negative weighting factors w

i

. The function F(a

0

,a

1

, . . . , a

N

) may be minimized as follows. The function F(a

0

,a

1

, . . . , a

N

) is given by

F (a_{0}, a_{1}, \dots, a_{N}) = \sum_{i = 1}^{M} w_{i} ε_{i}^{2} = \sum_{i = 1}^{M} {w_{i} [y (x_{i}) - y_{i}]}^{2}

and setting the partial derivatives of F(a

0

,a

1

, . . . , a

N

) with respect to a

j

, for j=0, 1, . . . , N, equal to zero gives

\begin{matrix} \frac{\partial F (a_{0}, a_{1}, \dots, a_{N})}{\partial a_{j}} = 2 \sum_{i = 1}^{M} w_{i} [y (x_{i}) - y_{i}] x_{i}^{j} \\ = 2 \sum_{i = 1}^{M} w_{i} [\sum_{k = 0}^{N} a_{k} f_{k} (x_{i}) - y_{i}] f_{j} (x_{i}) = 0, \end{matrix}

for j=0, 1, . . . , N, since f

j

(x

i

) is the coefficient of a

j

in the representation

y (x_{i}) = \sum_{j = 0}^{N} a_{k} f_{k} (x_{i}) .

Simplifying gives

\begin{matrix} \sum_{i = 1}^{M} w_{i} [\sum_{k = 0}^{N} a_{k} f_{k} (x_{i})] f_{j} (x_{i}) = \sum_{k = 0}^{N} a_{k} [\sum_{i = 1}^{M} w_{i} f_{k} (x_{i}) f_{j} (x_{i})] \\ = \sum_{i = 1}^{M} w_{i} f_{j} (x_{i}) y_{i}, or \\ \sum_{k = 0}^{N} a_{k} [\sum_{i = 1}^{M} w_{i} f_{k} (x_{i}) f_{j} (x_{i})] \\ \equiv \sum_{k = 0}^{N} a_{k} S_{k, j} = \sum_{i = 1}^{M} w_{i} f_{j} (x_{i}) y_{i} \equiv T_{j} \end{matrix}

for j=0, 1, . . . , N, where

\sum_{i = 1}^{M} w_{i} f_{k} (x_{i}) f_{j} (x_{i}) \equiv S_{k, j} and \sum_{i = 1}^{M} w_{i} f_{j} (x_{i}) y_{i} \equiv T_{j},

respectively. There are N+1 equations

\sum_{k = 0}^{N} a_{k} S_{k, j} = T_{j},

for j=0, 1, . . . , N, also known as the normal equations, including the non-negative weighting factors w

i

, for the N+1 unknown parameters a

k

, for k=0, 1, . . . , N, readily yielding a solution, provided that the determinant of the normal equations is not zero. This may be demonstrated by showing that the homogeneous equations

\sum_{k = 0}^{N} a_{k} S_{k, j} = 0

only have the trivial solution a

k

=0, for k=0, 1, . . . , N, which may be shown as follows. Multiply the jth homogeneous equation by a

j

and sum over all j,

\begin{matrix} \sum_{j = 0}^{N} a_{j} \sum_{k = 0}^{N} a_{k} S_{k, j} = \sum_{j = 0}^{N} a_{j} \sum_{k = 0}^{N} a_{k} \sum_{i = 1}^{M} w_{i} f_{k} (x_{i}) f_{j} (x_{i}) \\ = \sum_{i = 1}^{M} w_{i} (\sum_{k = 0}^{N} a_{k} f_{k} (x_{i})) (\sum_{j = 0}^{N} a_{j} f_{j} (x_{i})), \end{matrix}

but

\sum_{i = 1}^{M} w_{i} (\sum_{k = 0}^{N} a_{k} f_{k} (x_{i})) (\sum_{j = 0}^{N} a_{j} f_{j} (x_{i})) = \sum_{i = 1}^{M} w_{i} {(y (x_{i}))}^{2} = 0,

which would imply that y(x

i

)≡0, and, hence, that a

k

=0, for k=0, 1, . . . , N, the trivial solution. Therefore, the determinant of the normal equations is not zero, and the normal equations, including the non-negative weighting factors w

i

, may be solved for the N+1 parameters a

k

, for k=0, 1, . . . , N, the coefficients of the non-polynomial least-squares representation

y (x) = \sum_{j = 0}^{N} a_{j} f_{j} (x)

that may be fit to the M data points (x

i

,y

i

), using the linearly independent set of N+1 functions f

j

(x) as the basis for the non-polynomial least-squares representation

y (x) = \sum_{j = 0}^{N} a_{j} f_{j} (x),

and including the non-negative weighting factors w

i

.

As discussed above, in various alternative illustrative embodiments, control systems that use modeling of data in their functionality may be built. In various of these alternative illustrative embodiments, building the control model comprises fitting the collected processing data using Principal Components Analysis (PCA). This may be shown, for illustrative purposes only, as follows, by applying Principal Components Analysis (PCA) in process control (such as in etch endpoint determination), but, of course, the present invention is not limited to data compression in etch processing, nor to control models built by fitting the collected processing data using Principal Components Analysis (PCA).

Using optical emission spectroscopy (OES) as an analytical tool in process control (such as in etch endpoint determination) affords the opportunity for on-line measurements in real time. A calibration set of OES spectra bounding the acceptable process space within which a particular property (such as the etch endpoint) is to be controlled may be obtained by conventional means. Applying a multivariant statistical method such as Principal Components Analysis (PCA) to the calibration set of OES spectra affords a method of identifying the most important characteristics (Principal Components and respective Loadings and corresponding Scores) of the set of OES spectra that govern the controlled property, and are inherently related to the process. Control then is effected by using only up to four of such characteristics (Principal Components and respective Loadings and corresponding Scores), which can be determined quickly and simply from the measured spectra, as the control criteria to be applied to the process as a whole. The result is a very effective way of controlling a complex process using at most four criteria (the first through fourth Principal Components and respective Loadings and corresponding Scores) objectively determined from a calibration set, which can be applied in real time and virtually continuously, resulting in a well-controlled process that is (ideally) invariant.

In particular, it has been found that the second Principal Component contains a very robust, high signal-to-noise indicator for etch endpoint determination. It has also been found that the first four Principal Components similarly may be useful as indicators for etch endpoint determination as well as for data compression of OES data. In various illustrative embodiments, PCA may be applied to the OES data, either the whole spectrum or at least a portion of the whole spectrum. If the engineer and/or controller knows that only a portion of the OES data contains useful information, PCA may be applied only to that portion, for example.

In one illustrative embodiment of a method according to the present invention, as shown in

FIGS. 1-7

, described in more detail below, archived data sets of OES wavelengths (or frequencies), from wafers that had previously been plasma etched, may be processed and Loadings for the first through fourth Principal Components determined from the archived OES data sets may be used as model Loadings to calculate approximate Scores corresponding to newly acquired OES data. These approximate Scores, along with the mean values for each wavelength, may then be stored as compressed OES data.

In another illustrative embodiment of a method according to the present invention, as shown in

FIGS. 8-14

, described in more detail below, archived data sets of OES wavelengths (or frequencies), from wafers that had previously been plasma etched, may be processed and Loadings for the first through fourth Principal Components determined from the archived OES data sets may be used as model Loadings to calculate approximate Scores corresponding to newly acquired OES data. These approximate Scores may be used as an etch endpoint indicator to determine an endpoint for an etch process.

In yet another illustrative embodiment of a method according to the present invention, as shown in

FIGS. 15-21

, described in more detail below, archived data sets of OES wavelengths (or frequencies), from wafers that had previously been plasma etched, may be processed and Loadings for the first through fourth Principal Components, determined from the archived OES data sets, may be used as model Loadings to calculate approximate Scores corresponding to newly acquired OES data. These approximate Scores, along with model mean values also determined from the archived OES data sets, may then be stored as more compressed OES data. These approximate Scores may also be used as an etch endpoint indicator to determine an endpoint for an etch process.

Any of these three embodiments may be applied in real-time etch processing. Alternatively, either of the last two illustrative embodiments may be used as an identification technique when using batch etch processing, with archived data being applied statistically, to determine an etch endpoint for the batch.

Various embodiments of the present invention are applicable to any plasma etching process affording a characteristic data set whose quality may be said to define the “success” of the plasma etching process and whose identity may be monitored by suitable spectroscopic techniques such as OES. The nature of the plasma etching process itself is not critical, nor is the specific identity of the workpieces (such as semiconducting silicon wafers) whose spectra are being obtained. However, an “ideal” or “target” characteristic data set should be able to be defined, having a known and/or determinable OES spectrum, and variations in the plasma etching process away from the target characteristic data set should be able to be controlled (i.e., reduced) using identifiable independent process variables, e.g., etch endpoints, plasma energy and/or temperature, plasma chamber pressure, etchant concentrations, flow rates, and the like. The ultimate goal of process control is to maintain a properly operating plasma etching process at status quo. When a plasma etching process has been deemed to have reached a proper operating condition, process control should be able to maintain it by proper adjustment of plasma etching process parameters such as etch endpoints, temperature, pressure, flow rate, residence time in the plasma chamber, and the like. At proper operating conditions, the plasma etching process affords the “ideal” or “target” characteristics.

One feature of illustrative embodiments of our control process is that the spectrum of the stream, hereafter referred to as the test stream, may be obtained continuously (or in a near-continuous manner) on-line, and compared to the spectrum of the “target” stream, hereafter referred to as the standard stream. The difference in the two spectra may then be used to adjust one or more of the process variables so as to produce a test stream more nearly identical with the standard stream. However, as envisioned in practice, the complex spectral data will be reduced and/or compressed to no more than 4 numerical values that define the coordinates of the spectrum in the Principal Component or Factor space of the process subject to control. Typically small, incremental adjustments will be made so that the test stream approaches the standard stream while minimizing oscillation, particularly oscillations that will tend to lose process control rather than exercise control. Adjustments may be made according to one or more suitable algorithms based on process modeling, process experience, and/or artificial intelligence feedback, for example.

In principle, more than one process variable may be subject to control, although it is apparent that as the number of process variables under control increases so does the complexity of the control process. Similarly, more than one test stream and more than one standard stream may be sampled, either simultaneously or concurrently, again with an increase in complexity. In its simplest form, where there is one test stream and one process variable under control, one may analogize the foregoing to the use of a thermocouple in a reaction chamber to generate a heating voltage based on the sensed temperature, and to use the difference between the generated heating voltage and a set point heating voltage to send power to the reaction chamber in proportion to the difference between the actual temperature and the desired temperature which, presumably, has been predetermined to be the optimum temperature. Since the result of a given change in a process variable can be determined quickly, this new approach opens up the possibility of controlling the process by an automated “trial and error” feedback system, since unfavorable changes can be detected quickly. Illustrative embodiments of the present invention may operate as null-detectors with feedback from the set point of operation, where the feedback signal represents the deviation of the total composition of the stream from a target composition.

In one illustrative embodiment, OES spectra are taken of characteristic data sets of plasma etching processes of various grades and quality, spanning the maximum range of values typical of the particular plasma etching processes. Such spectra then are representative of the entire range of plasma etching processes and are often referred to as calibration samples. Note that because the characteristic data sets are representative of those formed in the plasma etching processes, the characteristic data sets constitute a subset of those that define the boundaries of representative processes. It will be recognized that there is no subset that is unique, that many different subsets may be used to define the boundaries, and that the specific samples selected are not critical.

Subsequently, the spectra of the calibration samples are subjected to the well-known statistical technique of Principal Component Analysis (PCA) to afford a small number of Principal Components (or Factors) that largely determine the spectrum of any sample. The Principal Components, which represent the major contributions to the spectral changes, are obtained from the calibration samples by PCA or Partial Least Squares (PLS). Thereafter, any new sample may be assigned various contributions of these Principal Components that would reproduce the spectrum of the new sample. The amount of each Principal Component required is called a Score, and time traces of these Scores, tracking how various of the Scores are changing with time, are used to detect deviations from the “target” spectrum.

In mathematical terms, a set of m time samples of an OES for a workpiece (such as a semiconductor wafer having various process layers formed thereon) taken at n channels or wavelengths or frequencies may be arranged as a rectangular n×m matrix X. In other words, the rectangular n×m matrix X may be comprised of 1 to n rows (each row corresponding to a separate OES channel or wavelength or frequency time sample) and 1 to m columns (each column corresponding to a separate OES spectrum time sample). The values of the rectangular n×m matrix X may be counts representing the intensity of the OES spectrum, or ratios of spectral intensities (normalized to a reference intensity), or logarithms of such ratios, for example. The rectangular n×m matrix X may have rank r, where r≦min{m,n} is the maximum number of independent variables in the matrix X. The use of PCA, for example, generates a set of Principal Components P (whose “Loadings,” or components, represent the contributions of the various spectral components) as an eigenmatrix (a matrix whose columns are eigenvectors) of the equation ((X−M)(X−M)

T

)P=Λ

2

P, where M is a rectangular n×m matrix of the mean values of the columns of X (the m columns of M are each the column mean vector μ

n×1

of X

n×m

), Λ

2

is an n×n diagonal matrix of the squares of the eigenvalues λ

i

, i=1,2, . . . , r, of the mean-scaled matrix X−M, and a Scores matrix, T, with X−M=PT

T

and (X−M)

T

=(PT

T

)

T

=(T

T

)

T

P

T

=TP

T

, so that ((X−M)(X−M)

T

)P=((PT

T

)(TP

T

))P and ((PT

T

)(TP

T

))P=(P(T

T

T)P

T

)P=P(T

T

T)=Λ

2

2P. The rectangular n×m matrix X, also denoted X

n×m

, may have elements x

ij

, where i=1,2, . . . , n, and j=1,2, . . . , m, and the rectangular m×n matrix X

T

, the transpose of the rectangular n×m matrix X, also denoted (X

T

)

m×n

, may have elements x

ji

, where i=1,2, . . . , n, and j=1,2, . . . , m. The n×n matrix (X−M)(X−M)

T

is (m−1) times the covariance matrix S

n×n

, having elements s

ij

, where i=1,2, . . . , n, and j=1,2, . . . , n, defined so that:

s_{ij} = \frac{m \sum_{k = 1}^{m} x_{ik} x_{jk} - \sum_{k = 1}^{m} x_{ik} \sum_{k = 1}^{m} x_{jk}}{m (m - 1)},

corresponding to the rectangular n×m matrix X

n×m

.

For the purposes of the process control envisioned in this application, it has been found that at most 4 Principal Components are needed to accommodate the data for a large range of plasma etching processes from a variety of plasma etching chambers. The spectrum of the standard sample is then expressed in terms of time traces of the Scores of the 4 Principal Components used, the spectrum of the test sample is similarly expressed, and the differences in the time traces of the Scores are used to control the process variables. Thus, no direct correlations between the sample spectrum and plasma etch endpoints need be known. In fact, the nature of the sample itself need not be known, as long as there is a standard, the OES spectrum of the standard is known, a set of at most 4 Principal Components is identified for the class of test stream samples, and one can establish how to use the 4 Principal Components to control the process variables (as discussed more fully below).

The spectrum of any sample may be expressed as a 2-dimensional representation of the intensity of emission at a particular wavelength vs. the wavelength. That is, one axis represents intensity, the other wavelength. The foregoing characterization of a spectrum is intended to incorporate various transformations that are mathematically covariant; e.g., instead of emission one might use absorption and/or transmission, and either may be expressed as a percentage or logarithmically. Whatever the details, each spectrum may be viewed as a vector. The group of spectra arising from a group of samples similarly corresponds to a group of vectors. If the number of samples is N, there are at most N distinct spectra. If, in fact, none of the spectra can be expressed as a linear combination of the other spectra, then the set of spectra define an N-dimensional spectrum space. However, in the cases of interest here, where a particular stream in an invariant plasma etching process is being sampled, it has been observed that, in fact, any particular spectrum may be accurately represented as a linear combination of a small number, M=4, of other spectra—their “Principal Components” that are referred to as “working” spectra. These “working” spectra may be viewed as the new basis set, i.e., linearly independent vectors that define the 4-dimensional spectrum space in which the samples reside. The spectrum of any other sample is then a linear combination of the “working” spectra (or is at least projectable onto the 4-dimensional spectrum space spanned by the “working” spectra). Our experience demonstrates that the samples typically reside in, at most, a 4-dimensional spectrum space, and that this 4-dimensional model suffices as a practical matter.

Statistical methods are available to determine the set of “working” spectra appropriate for any sample set, and the method of PCA is the one most favored in the practice of the present invention, although other methods, e.g., partial least squares, non-linear partial least squares, (or, with less confidence, multiple linear regression), also may be utilized. The “working” spectra, or the linearly independent vectors defining the sample spectrum space, are called Principal Components or Factors. Thus the spectrum of any sample is a linear combination of the Factors. The fractional contribution of any Factor is called the Score for the respective Factor. Hence, the spectrum of any sample completely defines a set of Scores that greatly reduces the apparent complexity of comparing different spectra. In fact, it has been found that for many processes of interest in semiconductor processing, a very small number of Factors, no more than 4, suffice to define accurately the sample spectrum space for the purpose of process control. This means that the process of characterizing the difference between a test sample and the standard sample comes down to the difference between only 4 numbers—the Scores of the respective 4 Factors for the sample and “target.” It is significant to note that the small number of Scores embodies a great deal of information about the samples and the process, and that only 4 numbers are adequate to control the process within quite close tolerances. By using the null approach of illustrative embodiments of the present invention, the use of Scores is simplified to teaching small shifts (and restoring them to zero) rather than drawing conclusions and/or correlations from the absolute values of the Scores.

Although other methods may exist, four methods for computing Principal Components are as follows:

1. eigenanalysis (EIG);

2. singular value decomposition (SVD);

3. nonlinear iterative partial least squares (NIPALS); and

4. power method.

Each of the first two methods, EIG and SVD, simultaneously calculates all possible Principal Components, whereas the NIPALS method allows for calculation of one Principal Component at a time. However, the power method, described more fully below, is an iterative approach to finding eigenvalues and eigenvectors, and also allows for calculation of one Principal Component at a time. There are as many Principal Components as there are channels (or wavelengths or frequencies). The power method may efficiently use computing time.

For example, consider the 3×2 matrix A, its transpose, the 2×3 matrix A

T

, their 2×2 matrix product A

T

A, and their 3×3 matrix product AA

T

:

A = (\begin{matrix} 1 & 1 \\ 1 & 0 \\ 1 & - 1 \end{matrix})

A^{T} = (\begin{matrix} 1 & 1 & 1 \\ 1 & 0 & - 1 \end{matrix})

A^{T} A = (\begin{matrix} 1 & 1 & 1 \\ 1 & 0 & - 1 \end{matrix}) (\begin{matrix} 1 & 1 \\ 1 & 0 \\ 1 & - 1 \end{matrix}) = (\begin{matrix} 3 & 0 \\ 0 & 2 \end{matrix})

{AA}^{T} = (\begin{matrix} 1 & 1 \\ 1 & 0 \\ 1 & - 1 \end{matrix}) (\begin{matrix} 1 & 1 & 1 \\ 1 & 0 & - 1 \end{matrix}) = (\begin{matrix} 2 & 1 & 0 \\ 1 & 1 & 1 \\ 0 & 1 & 2 \end{matrix}) .

EIG reveals that the eigenvalues λ of the matrix product A

T

A are 3 and 2. The eigenvectors of the matrix product A

T

A are solutions t of the equation (A

T

A)

t

=λ

t

, and may be seen by inspection to be

t

1

T

=(1,0) and

t

2

T

=(0,1), belonging to the eigenvalues λ

1

=3 and λ

2

=2, respectively.

The power method, for example, may be used to determine the eigenvalues λ and eigenvectors p of the matrix product AA

T

, where the eigenvalues λ and the eigenvectors p are solutions p of the equation (AA

T

)

p

=λ

p

. A trial eigenvector

p

T

=(1,1,1) may be used:

({AA}^{T}) \underline{p} = (\begin{matrix} 2 & 1 & 0 \\ 1 & 1 & 1 \\ 0 & 1 & 2 \end{matrix}) (\begin{matrix} 1 \\ 1 \\ 1 \end{matrix}) = (\begin{matrix} 3 \\ 3 \\ 3 \end{matrix}) = 3 (\begin{matrix} 1 \\ 1 \\ 1 \end{matrix}) = λ_{1} {\underline{p}}_{1} .

This indicates that the trial eigenvector

p

T

=(1,1,1) happened to correspond to the eigenvector

p

1

T

=(1,1,1) belonging to the eigenvalue λ

1

=3. The power method then proceeds by subtracting the outer product matrix

p

1

p

1

T

from the matrix product AA

T

to form a residual matrix R

1

:

R_{1} = (\begin{matrix} 2 & 1 & 0 \\ 1 & 1 & 1 \\ 0 & 1 & 2 \end{matrix}) - (\begin{matrix} 1 \\ 1 \\ 1 \end{matrix}) (\begin{matrix} 1 & 1 & 1 \end{matrix}) = (\begin{matrix} 2 & 1 & 0 \\ 1 & 1 & 1 \\ 0 & 1 & 2 \end{matrix}) - (\begin{matrix} 1 & 1 & 1 \\ 1 & 1 & 1 \\ 1 & 1 & 1 \end{matrix}) = (\begin{matrix} 1 & 0 & - 1 \\ 0 & 0 & 0 \\ - 1 & 0 & 1 \end{matrix}) .

Another trial eigenvector

p

T

=(1,0,−1) may be used:

({AA}^{T} - {\underline{p}}_{1} {\underline{p}}_{1}^{T}) \underline{p} = R_{1} \underline{p} = (\begin{matrix} 1 & 0 & - 1 \\ 0 & 0 & 0 \\ - 1 & 0 & 1 \end{matrix}) (\begin{matrix} 1 \\ 0 \\ - 1 \end{matrix}) = (\begin{matrix} 2 \\ 0 \\ - 2 \end{matrix}) = 2 (\begin{matrix} 1 \\ 0 \\ - 1 \end{matrix}) = λ_{2} {\underline{p}}_{2} .

This indicates that the trial eigenvector

p

T

=(1,0,−1) happened to correspond to the eigenvector

p

2

T

=(1,0,−1) belonging to the eigenvalue λ

2

=2. The power method then proceeds by subtracting the outer product matrix

p

2

p

2

T

from the residual matrix R

1

to form a second residual matrix R

2

:

R_{2} = (\begin{matrix} 1 & 0 & - 1 \\ 0 & 0 & 0 \\ - 1 & 0 & 1 \end{matrix}) - (\begin{matrix} 1 \\ 0 \\ - 1 \end{matrix}) (\begin{matrix} 1 & 0 & - 1 \end{matrix}) = (\begin{matrix} 1 & 0 & - 1 \\ 0 & 0 & 0 \\ - 1 & 0 & 1 \end{matrix}) - (\begin{matrix} 1 & 0 & - 1 \\ 0 & 0 & 0 \\ - 1 & 0 & 1 \end{matrix}) = (\begin{matrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix}) .

The fact that the second residual matrix R

2

vanishes indicates that the eigenvalue λ

3

=0 and that the eigenvector

p

3

is completely arbitrary. The eigenvector

p

3

may be conveniently chosen to be orthogonal to the eigenvectors

p

1

T=(1,1,1) and

p

2

T=(1,0,−1), so that the eigenvector

p

3

T=(1,−2,1). Indeed, one may readily verify that:

({AA}^{T}) {\underline{p}}_{3} = (\begin{matrix} 2 & 1 & 0 \\ 1 & 1 & 1 \\ 0 & 1 & 2 \end{matrix}) (\begin{matrix} 1 \\ 2 \\ - 1 \end{matrix}) = (\begin{matrix} 0 \\ 0 \\ 0 \end{matrix}) = 0 (\begin{matrix} 1 \\ 2 \\ - 1 \end{matrix}) = λ_{3} {\underline{p}}_{3} .

Similarly, SVD of A shows that A=PT

T

, where P is the Principal Component matrix and T is the Scores matrix:

A = (\begin{matrix} 1 / \sqrt{3} & 1 / \sqrt{2} & 1 / \sqrt{6} \\ 1 / \sqrt{3} & 0 & - 2 / \sqrt{6} \\ 1 / \sqrt{3} & - 1 / \sqrt{2} & 1 / \sqrt{6} \end{matrix}) (\begin{matrix} \sqrt{3} & 0 \\ 0 & \sqrt{2} \\ 0 & 0 \end{matrix}) (\begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix}) = (\begin{matrix} 1 / \sqrt{3} & 1 / \sqrt{2} & 1 / \sqrt{6} \\ 1 / \sqrt{3} & 0 & - 2 / \sqrt{6} \\ 1 / \sqrt{3} & - 1 / \sqrt{2} & 1 / \sqrt{6} \end{matrix}) (\begin{matrix} \sqrt{3} & 0 \\ 0 & \sqrt{2} \\ 0 & 0 \end{matrix}) .

SVD confirms that the singular values of A are 3 and 2, the positive square roots of the eigenvalues λ

1

=3 and λ

2

=2 of the matrix product A

T

A. Note that the columns of the Principal Component matrix P are the orthonormalized eigenvectors of the matrix product AA

T

.

Likewise, SVD of A

T

shows that A

T

=TP

T

:

A^{T} = (\begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix}) (\begin{matrix} \sqrt{3} & 0 & 0 \\ 0 & \sqrt{2} & 0 \end{matrix}) (\begin{matrix} \frac{1}{\sqrt{3}} & \frac{1}{\sqrt{3}} & \frac{1}{\sqrt{3}} \\ \frac{1}{\sqrt{2}} & 0 & \frac{- 1}{\sqrt{2}} \\ \frac{1}{\sqrt{6}} & \frac{- 2}{\sqrt{6}} & \frac{1}{\sqrt{6}} \end{matrix}) = A^{T} = (\begin{matrix} \sqrt{3} & 0 & 0 \\ 0 & \sqrt{2} & 0 \end{matrix}) (\begin{matrix} \frac{1}{\sqrt{3}} & \frac{1}{\sqrt{3}} & \frac{1}{\sqrt{3}} \\ \frac{1}{\sqrt{2}} & 0 & \frac{- 1}{\sqrt{2}} \\ \frac{1}{\sqrt{6}} & \frac{- 2}{\sqrt{6}} & \frac{1}{\sqrt{6}} \end{matrix}) = {TP}^{T} .

SVD confirms that the (non-zero) singular values of AT are 3 and 2, the positive square roots of the eigenvalues λ

1

=3 and λ

2

=2 of the matrix product AA

T

. Note that the columns of the Principal Component matrix P (the rows of the Principal Component matrix P

T

) are the orthonormalized eigenvectors of the matrix product AA

T

. Also note that the non-zero elements of the Scores matrix T are the positive square roots 3 and 2 of the (non-zero) eigenvalues λ

1

=3 and λ

2

=2 of both of the matrix products A

T

A and AA

T

.

Taking another example, consider the 4×3 matrix B, its transpose, the 3×4 matrix B

T

, their 3×3 matrix product B

T

B, and their 4×4 matrix product BB

T

:

B = (\begin{matrix} 1 & 1 & 0 \\ 1 & 0 & 1 \\ 1 & 0 & - 1 \\ 1 & - 1 & 0 \end{matrix})

B^{T} = (\begin{matrix} 1 & 1 & 1 & 1 \\ 1 & 0 & 0 & - 1 \\ 0 & 1 & - 1 & 0 \end{matrix})

B^{T} B = (\begin{matrix} 1 & 1 & 1 & 1 \\ 1 & 0 & 0 & - 1 \\ 0 & 1 & - 1 & 0 \end{matrix}) (\begin{matrix} 1 & 1 & 0 \\ 1 & 0 & 1 \\ 1 & 0 & - 1 \\ 1 & - 1 & 0 \end{matrix}) = (\begin{matrix} 4 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 2 \end{matrix})

{BB}^{T} = (\begin{matrix} 1 & 1 & 0 \\ 1 & 0 & 1 \\ 1 & 0 & - 1 \\ 1 & - 1 & 0 \end{matrix}) (\begin{matrix} 1 & 1 & 1 & 1 \\ 1 & 0 & 0 & - 1 \\ 0 & 1 & - 1 & 0 \end{matrix}) = (\begin{matrix} 2 & 1 & 1 & 0 \\ 1 & 2 & 0 & 1 \\ 1 & 0 & 2 & 1 \\ 0 & 1 & 1 & 2 \end{matrix}) .

EIG reveals that the eigenvalues of the matrix product B

T

B are 4, 2 and 2. The eigenvectors of the matrix product B

T

B are solutions t of the equation (B

T

B)

t

=λ

t

, and may be seen by inspection to be

t

1

T

=(1,0,0),

t

2

T

=(0,1,0), and

t

3

T

=(0,0,1), belonging to the eigenvalues λ

1

=4, λ

2

=2, and λ

3

=2, respectively.

The power method, for example, may be used to determine the eigenvalues λ and eigenvectors

p

of the matrix product BB

T

, where the eigenvalues λ and the eigenvectors

p

are solutions

p

of the equation (BB

T

)

p

=λ

p

. A trial eigenvector

p

T

=(1,1,1,1) may be used:

({BB}^{T}) \underline{p} = (\begin{matrix} 2 & 1 & 1 & 0 \\ 1 & 2 & 0 & 1 \\ 1 & 0 & 2 & 1 \\ 0 & 1 & 1 & 2 \end{matrix}) (\begin{matrix} 1 \\ 1 \\ 1 \\ 1 \end{matrix}) = (\begin{matrix} 4 \\ 4 \\ 4 \\ 4 \end{matrix}) = 4 (\begin{matrix} 1 \\ 1 \\ 1 \\ 1 \end{matrix}) = λ_{1} {\underline{p}}_{1} .

This indicates that the trial eigenvector

p

T

=(1,1,1,1) happened to correspond to the eigenvector

p

1

T

=(1,1,1,1) belonging to the eigenvalue λ

1

=4. The power method then proceeds by subtracting the outer product matrix

p

1

p

1

T

from the matrix product BB

T

to form a residual matrix R

1

:

R_{1} = (\begin{matrix} 2 & 1 & 1 & 0 \\ 1 & 2 & 0 & 1 \\ 1 & 0 & 2 & 1 \\ 0 & 1 & 1 & 2 \end{matrix}) - (\begin{matrix} 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 \\ 1 & 1 & 1 & 1 \end{matrix}) = (\begin{matrix} 1 & 0 & 0 & - 1 \\ 0 & 1 & - 1 & 0 \\ 0 & - 1 & 1 & 0 \\ - 1 & 0 & 0 & 1 \end{matrix}) .

Another trial eigenvector

p

T

=(1,0,0,−1) may be used:

({BB}^{T} - {\underline{p}}_{1} {\underline{p}}_{1}^{T}) \underline{p} = R_{1} \underline{p} = (\begin{matrix} 1 & 0 & 0 & - 1 \\ 0 & 1 & - 1 & 0 \\ 0 & - 1 & 1 & 0 \\ - 1 & 0 & 0 & 1 \end{matrix}) (\begin{matrix} 1 \\ 0 \\ 0 \\ - 1 \end{matrix}) = (\begin{matrix} 2 \\ 0 \\ 0 \\ - 2 \end{matrix}) = 2 (\begin{matrix} 1 \\ 0 \\ 0 \\ - 1 \end{matrix}) = λ_{2} {\underline{p}}_{2} .

This indicates that the trial eigenvector

p

T

=(1,0,0,−1) happened to correspond to the eigenvector

p

2

T

=(1,0,0,−1) belonging to the eigenvalue λ

2

=2. The power method then proceeds by subtracting the outer product matrix

p

2

p

2

T

from the residual matrix BB

T

to form a second residual matrix R

2

:

R_{2} = (\begin{matrix} 1 & 0 & 0 & - 1 \\ 0 & 1 & - 1 & 0 \\ 0 & - 1 & 1 & 0 \\ - 1 & 0 & 0 & 1 \end{matrix}) - (\begin{matrix} 1 & 0 & 0 & - 1 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ - 1 & 0 & 0 & 1 \end{matrix}) = (\begin{matrix} 0 & 0 & 0 & 0 \\ 0 & 1 & - 1 & 0 \\ 0 & - 1 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{matrix}) .

Another trial eigenvector

p

T

=(0,1,−1,0) may be used:

({BB}^{T} - {\underline{p}}_{2} {\underline{p}}_{2}^{T}) \underline{p} = R_{2} \underline{p} = (\begin{matrix} 0 & 0 & 0 & 0 \\ 0 & 1 & - 1 & 0 \\ 0 & - 1 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{matrix}) (\begin{matrix} 0 \\ 1 \\ - 1 \\ 0 \end{matrix}) = (\begin{matrix} 0 \\ 2 \\ - 2 \\ 0 \end{matrix}) = 2 (\begin{matrix} 0 \\ 1 \\ - 1 \\ 0 \end{matrix}) = λ_{3} {\underline{p}}_{3} .

This indicates that the trial eigenvector

p

T

=(0,1,−1,0) happened to correspond to the eigenvector

p

3

T

=(0,1,−1,0) belonging to the eigenvalue λ

3

=2. The power method then proceeds by subtracting the outer product matrix

p

3

p

3

T

from the second residual matrix R

2

to form a third residual matrix R

3

:

R_{3} = (\begin{matrix} 0 & 0 & 0 & 0 \\ 0 & 1 & - 1 & 0 \\ 0 & - 1 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{matrix}) - (\begin{matrix} 0 & 0 & 0 & 0 \\ 0 & 1 & - 1 & 0 \\ 0 & - 1 & 1 & 0 \\ 0 & 0 & 0 & 0 \end{matrix}) = (\begin{matrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{matrix}) .

The fact that the third residual matrix R

3

vanishes indicates that the eigenvalue λ

4

=0 and that the eigenvector

p

4

is completely arbitrary. The eigenvector

p

4

may be conveniently chosen to be orthogonal to the eigenvectors

p

1

T

=(1,1,1,1),

p

2

T

=(1,0,0,−1), and

p

3

T

=(0,1,−1,0), so that the eigenvector

p

4

T

=(1,−1,−1,1). Indeed, one may readily verify that:

({BB}^{T}) {\underline{p}}_{4} = (\begin{matrix} 2 & 1 & 1 & 0 \\ 1 & 2 & 0 & 1 \\ 1 & 0 & 2 & 1 \\ 0 & 1 & 1 & 2 \end{matrix}) (\begin{matrix} 1 \\ - 1 \\ - 1 \\ 1 \end{matrix}) = (\begin{matrix} 0 \\ 0 \\ 0 \\ 0 \end{matrix}) = 0 (\begin{matrix} 1 \\ - 1 \\ - 1 \\ 1 \end{matrix}) = λ_{4} {\underline{p}}_{4} .

In this case, since the eigenvalues λ

2

=2 and λ

3

=2 are equal, and, hence, degenerate, the eigenvectors

p

2

T

=(1,0,0,−1) and

p

3

T

=(0,1,−1,0) belonging to the degenerate eigenvalues λ

2

=2=λ

3

may be conveniently chosen to be orthonormal. A Gram-Schmidt orthonormalization procedure may be used, for example.

Similarly, SVD of B shows that B=PT

T

, where P is the Principal Component matrix and T is the Scores matrix:

B = (\begin{matrix} 1 / 2 & 1 / \sqrt{2} & 0 & 1 / 2 \\ 1 / 2 & 0 & 1 / \sqrt{2} & - 1 / 2 \\ 1 / 2 & 0 & - 1 / \sqrt{2} & - 1 / 2 \\ 1 / 2 & - 1 / \sqrt{2} & 0 & 1 / 2 \end{matrix}) (\begin{matrix} 2 & 0 & 0 \\ 0 & \sqrt{2} & 0 \\ 0 & 0 & \sqrt{2} \\ 0 & 0 & 0 \end{matrix}) (\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}) = B = (\begin{matrix} 1 / 2 & 1 / \sqrt{2} & 0 & 1 / 2 \\ 1 / 2 & 0 & 1 / \sqrt{2} & - 1 / 2 \\ 1 / 2 & 0 & - 1 / \sqrt{2} & - 1 / 2 \\ 1 / 2 & - 1 / \sqrt{2} & 0 & 1 / 2 \end{matrix}) (\begin{matrix} 2 & 0 & 0 \\ 0 & \sqrt{2} & 0 \\ 0 & 0 & \sqrt{2} \\ 0 & 0 & 0 \end{matrix}) = {PT}^{T} .

SVD confirms that the singular values of B are 2, 2 and 2, the positive square roots of the eigenvalues λ

1

=4, λ

2

=2 and λ

3

=2 of the matrix product B

T

B.

Likewise, SVD of B

T

shows that:

B^{T} = (\begin{matrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{matrix}) (\begin{matrix} 2 & 0 & 0 & 0 \\ 0 & \sqrt{2} & 0 & 0 \\ 0 & 0 & \sqrt{2} & 0 \end{matrix}) (\begin{matrix} 1 / 2 & 1 / 2 & 1 / 2 & 1 / 2 \\ 1 / \sqrt{2} & 0 & 0 & - 1 / \sqrt{2} \\ 0 & 1 / \sqrt{2} & - 1 / \sqrt{2} & 0 \\ 1 / 2 & - 1 / 2 & - 1 / 2 & 1 / 2 \end{matrix}) = B^{T} = (\begin{matrix} 2 & 0 & 0 & 0 \\ 0 & \sqrt{2} & 0 & 0 \\ 0 & 0 & \sqrt{2} & 0 \end{matrix}) (\begin{matrix} 1 / 2 & 1 / 2 & 1 / 2 & 1 / 2 \\ 1 / \sqrt{2} & 0 & 0 & - 1 / \sqrt{2} \\ 0 & 1 / \sqrt{2} & - 1 / \sqrt{2} & 0 \\ 1 / 2 & - 1 / 2 & - 1 / 2 & 1 / 2 \end{matrix}) = {TP}^{T} .

SVD confirms that the (non-zero) singular values of B

T

are 2, 2 and 2, the positive square roots of the eigenvalues λ

1

=4, λ

2

=2 and λ

3

3

=2 of the matrix product AA

T

. Note that the columns of the Principal Component matrix P (the rows of the Principal Component matrix P

T

) are the orthonormalized eigenvectors of the matrix product BB

T

. Also note that the non-zero elements of the Scores matrix T are the positive square roots 2, 2, and 2 of the (non-zero) eigenvalues λ

1

=4, λ

2

=2 and λ

3

=2 of both of the matrix products B

T

B and BB

T

.

The matrices A and B discussed above have been used for the sake of simplifying the presentation of PCA and the power method, and are much smaller than the data matrices encountered in illustrative embodiments of the present invention. For example, in one illustrative embodiment, for each wafer, 8 scans of OES data over 495 wavelengths may be taken during an etching step, with about a 13 second interval between scans. In this illustrative embodiment, 18 wafers may be run and corresponding OES data collected. The data may be organized as a set of 18×495 matrices X

s

=(X

ij

)

s

, where s=1,2, . . . , 8, for each of the different scans, and X

ij

is the intensity of the ith wafer run at the jth wavelength. Putting all 8 of the 18×495 matrices X

s

=(X

ij

)

s

, for s=1,2, . . . , 8, next to each other produces the overall OES data matrix X, an 18×3960 matrix X=[X

1

,X

2

, . . . , X

8

]=[(X

ij

)

1

,(X

ij

)

2

, . . . , (X

ij

)

8

,]. Each row in X represents the OES data from 8 scans over 495 wavelengths for a run. Brute force modeling using all 8 scans and all 495 wavelengths would entail using 3960 input variables to predict the etching behavior of 18 sample wafers, an ill-conditioned regression problem. Techniques such as PCA and/or partial least squares (PLS, also known as projection to latent structures) reduce the complexity in such cases by revealing the hierarchical ordering of the data based on levels of decreasing variability. In PCA, this involves finding successive Principal Components. In PLS techniques such as NIPALS, this involves finding successive latent vectors.

As shown in

FIG. 22

, a scatterplot

2200

of data points 2210 may be plotted in an n-dimensional variable space (n=3 in FIG.

22

). The mean vector

2220

may lie at the center of a p-dimensional Principal Component ellipsoid

2230

(p=2 in FIG.

22

). The mean vector

2220

may be determined by taking the average of the columns of the overall OES data matrix X. The Principal Component ellipsoid

2230

may have a first Principal Component

2240

(major axis in FIG.

22

), with a length equal to the largest eigenvalue of the mean-scaled OES data matrix X−M, and a second Principal Component

2250

(minor axis in FIG.

22

), with a length equal to the next largest eigenvalue of the mean-scaled OES data matrix X−M.

For example, the 3×4 matrix B

T

given above may be taken as the overall OES data matrix X (again for the sake of simplicity), corresponding to 4 scans taken at 3 wavelengths. As shown in

FIG. 23

, a scatterplot

2300

of data points

2310

may be plotted in a 3-dimensional variable space. The mean vector

2320

μ

may lie at the center of a 2-dimensional Principal Component ellipsoid

2330

(really a circle, a degenerate ellipsoid). The mean vector

2320

μ

a may be determined by taking the average of the columns of the overall OES 3×4 data matrix B

T

. The Principal Component ellipsoid

2330

may have a first Principal Component

2340

(“major” axis in

FIG. 23

) and a second Principal Component

2350

(“minor” axis in FIG.

23

). Here, the eigenvalues of the mean-scaled OES data matrix B

T

−M are equal and degenerate, so the lengths of the “major” and “minor” axes in

FIG. 23

are equal. As shown in

FIG. 23

, the mean vector

2320

μ

is given by:

\underline{μ} = \frac{1}{4} [(\begin{matrix} 1 \\ 1 \\ 0 \end{matrix}) + (\begin{matrix} 1 \\ 0 \\ 1 \end{matrix}) + (\begin{matrix} 1 \\ 0 \\ - 1 \end{matrix}) + (\begin{matrix} 1 \\ - 1 \\ 0 \end{matrix})] = (\begin{matrix} 1 \\ 0 \\ 0 \end{matrix}),

and the matrix M has the mean vector

2320

μ

for all 4 columns.

In another illustrative embodiment, 5500 samples of each wafer may be taken on wavelengths between about 240-1100 nm at a high sample rate of about one per second. For example, 5551 sampling points/spectrum/second (corresponding to 1 scan per wafer per second taken at 5551 wavelengths, or 7 scans per wafer per second taken at 793 wavelengths, or 13 scans per wafer per second taken at 427 wavelengths, or 61 scans per wafer per second taken at 91 wavelengths) may be collected in real time, during etching of a contact hole using an Applied Materials AMAT 5300 Centura etching chamber, to produce high resolution and broad band OES spectra.

As shown in

FIG. 24

, a representative OES spectrum

2400

of a contact hole etch is illustrated. Wavelengths, measured in nanometers (nm) are plotted along the horizontal axis against spectrometer counts plotted along the vertical axis.

As shown in

FIG. 25

, a representative OES trace

2500

of a contact hole etch is illustrated. Time, measured in seconds (sec) is plotted along the horizontal axis against spectrometer counts plotted along the vertical axis. As shown in

FIG. 25

, by about 40 seconds into the etching process, as indicated by dashed line

2510

, the OES trace

2500

of spectrometer counts “settles down” to a range of values less than or about 300, for example.

As shown in

FIG. 26

, representative OES traces

2600

of a contact hole etch are illustrated. Wavelengths, measured in nanometers (nm) are plotted along a first axis, time, measured in seconds (sec) is plotted along a second axis, and mean-scaled OES spectrometer counts, for example, are plotted along a third (vertical) axis. As shown in

FIG. 26

, over the course of about 150 seconds of etching, three clusters of wavelengths

2610

,

2620

and

2630

, respectively, show variations in the respective mean-scaled OES spectrometer counts. In one illustrative embodiment, any one of the three clusters of wavelengths

2610

,

2620

and

2630

may be used, either taken alone or taken in any combination with any one (or both) of the others, as an indicator variable signaling an etch endpoint. In an alternative illustrative embodiment, only the two clusters of wavelengths

2620

and

2630

having absolute values of mean-scaled OES spectrometer counts that exceed a preselected threshold absolute mean-scaled OES spectrometer count value (for example, about 200, as shown in

FIG. 26

) may be used, either taken alone or taken together, as an indicator variable signaling an etch endpoint. In yet another alternative illustrative embodiment, only one cluster of wavelengths

2630

having an absolute value of mean-scaled OES spectrometer counts that exceeds a preselected threshold absolute mean-scaled OES spectrometer count value (for example, about 300, as shown in

FIG. 26

) may be used as an indicator variable signaling an etch endpoint.

As shown in

FIG. 27

, a representative Scores time trace

2700

corresponding to the second Principal Component during a contact hole etch is illustrated. Time, measured in seconds (sec) is plotted along the horizontal axis against Scores (in arbitrary units) plotted along the vertical axis. As shown in

FIG. 27

, the Scores time trace

2700

corresponding to the second Principal Component during a contact hole etch may start at a relatively high value initially, decrease with time, pass through a minimum value, and then begin increasing before leveling off. It has been found that the inflection point (indicated by dashed line

2710

, and approximately where the second derivative of the Scores time trace

2700

with respect to time vanishes) is a robust indicator for the etch endpoint.

Principal Components Analysis (PCA) may be illustrated geometrically. For example, the 3×2 matrix C (similar to the 3×2 matrix A given above):

C = (\begin{matrix} 1 & - 1 \\ 1 & 0 \\ 1 & 1 \end{matrix})

may be taken as the overall OES data matrix X (again for the sake of simplicity), corresponding to 2 scans taken at 3 wavelengths. As shown in

FIG. 28

, a scatterplot

2800

of OES data points

2810

and

2820

, with coordinates (1,1,1) and (−1,0,1), respectively, may be plotted in a 3-dimensional variable space where the variables are respective spectrometer counts for each of the 3 wavelengths. The mean vector

2830

μ

may lie at the center of a 1-dimensional Principal Component ellipsoid

2840

(really a line, a very degenerate ellipsoid). The mean vector

2830

μ

may be determined by taking the average of the columns of the overall OES 3×2 matrix C. The Principal Component ellipsoid

2840

may have a first Principal Component

2850

(the “major” axis in

FIG. 28

, with length 5, lying along a first Principal Component axis

2860

) and no second or third Principal Component lying along second or third Principal Component axes

2870

and

2880

, respectively. Here, two of the eigenvalues of the mean-scaled OES data matrix C−M are equal to zero, so the lengths of the “minor” axes in

FIG. 28

are both equal to zero. As shown in

FIG. 28

, the mean vector

2830

μ

is given by:

\underline{μ} = \frac{1}{2} [(\begin{matrix} 1 \\ 1 \\ 1 \end{matrix}) + (\begin{matrix} - 1 \\ 0 \\ 1 \end{matrix})] = (\begin{matrix} 0 \\ 1 / 2 \\ 1 \end{matrix}),

and the matrix M has the mean vector

2830

μ

for both columns. As shown in

FIG. 28

, PCA is nothing more than a principal axis rotation of the original variable axes (here, the OES spectrometer counts for 3 wavelengths) about the endpoint of the mean vector

2830

μ

, with coordinates (0,1/2,1) with respect to the original coordinate axes and coordinates [0,0,0] with respect to the new Principal Component axes

2860

,

2870

and

2880

. The Loadings are merely the direction cosines of the new Principal Component axes

2860

,

2870

and

2880

with respect to the original variable axes. The Scores are simply the coordinates of the OES data points

2810

and

2820

, [5

0.5

/2,0,0] and [−5

0.5

/2,0,0], respectively, referred to the new Principal Component axes

2860

,

2870

and

2880

.

The mean-scaled 3×2 OES data matrix C−M, its transpose, the 2×3 matrix (C−M)

T

, their 2×2 matrix product (C−M)

T

(C−M), and their 3×3 matrix product (C−M) (C−M)

T

are given by:

C - M = (\begin{matrix} 1 & - 1 \\ 1 & 0 \\ 1 & 1 \end{matrix}) - (\begin{matrix} 0 & 0 \\ 1 / 2 & 1 / 2 \\ 1 & 1 \end{matrix}) = (\begin{matrix} 1 & - 1 \\ 1 / 2 & 1 / 2 \\ 0 & 0 \end{matrix})

{(C - M)}^{T} = (\begin{matrix} 1 & 1 / 2 & 0 \\ - 1 & - 1 / 2 & 0 \end{matrix})

{(C - M)}^{T} (C - M) = (\begin{matrix} 1 & 1 / 2 & 0 \\ - 1 & - 1 / 2 & 0 \end{matrix}) (\begin{matrix} 1 & - 1 \\ 1 / 2 & 1 / 2 \\ 0 & 0 \end{matrix}) = (\begin{matrix} 5 / 4 & - 5 / 4 \\ - 5 / 4 & 5 / 4 \end{matrix})

(C - M) {(C - M)}^{T} = (\begin{matrix} 1 & - 1 \\ 1 / 2 & 1 / 2 \\ 0 & 0 \end{matrix}) (\begin{matrix} 1 & 1 / 2 & 0 \\ - 1 & - 1 / 2 & 0 \end{matrix}) = (\begin{matrix} 2 & 1 & 0 \\ 1 & 1 / 2 & 0 \\ 0 & 0 & 0 \end{matrix})

The 3×3 matrix (C−M)(C−M)

T

is the covariance matrix S

3×3

, having elements s

ij

, where i=1,2,3, and j=1,2,3, defined so that:

s_{ij} = \frac{2 \sum_{k = 1}^{2} c_{ik} c_{jk} - \sum_{k = 1}^{2} c_{ik} \sum_{k = 1}^{2} c_{jk}}{2 (2 - 1)},

corresponding to the rectangular 3×2 matrix C

3×2

.

EIG reveals that the eigenvalues λ of the matrix product (C−M)

T

(C−M) are 5/2 and 0, for example, by finding solutions to the secular equation:

&LeftBracketingBar; \begin{matrix} 5 / 4 - λ & - 5 / 4 \\ - 5 / 4 & 5 / 4 - λ \end{matrix} &RightBracketingBar; = 0.

The eigenvectors of the matrix product (C−M)

T

(C−M) are solutions t of the equation (C−M)

T

(C−M)

t

=λ

t

, which may be rewritten as ((C−M)

T

(C−M)−λ)

t

=0. For the eigenvalue λ

1

=5/2, the eigenvector

t

1

may be seen by

(\begin{matrix} 5 / 4 - λ & - 5 / 4 \\ - 5 / 4 & 5 / 4 - λ \end{matrix}) \underline{t} = (\begin{matrix} - 5 / 4 & - 5 / 4 \\ - 5 / 4 & - 5 / 4 \end{matrix}) \underline{t} = 0

to be

t

1

T

=(1,−1). For the eigenvalue λ

1

=0, the eigenvector

t

2

may be seen by

(\begin{matrix} 5 / 4 - λ & - 5 / 4 \\ - 5 / 4 & 5 / 4 - λ \end{matrix}) \underline{t} = (\begin{matrix} 5 / 4 & - 5 / 4 \\ - 5 / 4 & 5 / 4 \end{matrix}) \underline{t} = 0

The power method, for example, may be used to determine the eigenvalues λ and eigenvectors

p

of the matrix product (C−M)(C−M)

T

, where the eigenvalues λ and the eigenvectors

p

are solutions

p

of the equation ((C−M)(C−M)

T

)p=λ

p

. A trial eigenvector

p

T

=(1,1,1) may be used:

((C - M) {(C - M)}^{T}) \underline{p} = (\begin{matrix} 2 & 1 & 0 \\ 1 & 1 / 2 & 0 \\ 0 & 0 & 0 \end{matrix}) (\begin{matrix} 1 \\ 1 \\ 1 \end{matrix}) = (\begin{matrix} 3 \\ 3 / 2 \\ 0 \end{matrix}) = 3 (\begin{matrix} 1 \\ 1 / 2 \\ 0 \end{matrix}) = 3 \underline{q} ((C - M) {(C - M)}^{T}) \underline{q} = (\begin{matrix} 2 & 1 & 0 \\ 1 & 1 / 2 & 0 \\ 0 & 0 & 0 \end{matrix}) (\begin{matrix} 1 \\ 1 / 2 \\ 0 \end{matrix}) = (\begin{matrix} 5 / 2 \\ 5 / 4 \\ 0 \end{matrix}) = 5 / 2 (\begin{matrix} 1 \\ 1 / 2 \\ 0 \end{matrix}) = λ_{1} {\underline{p}}_{1} .

This illustrates that the trial eigenvector

p

T

=(1,1,1) gets replaced by the improved trial eigenvector

q

T

=(1,1/2,0) that happened to correspond to the eigenvector

p

1

T

=(1,1/2,0) belonging to the eigenvalue λ

1

=5/2. The power method then proceeds by subtracting the outer product matrix

p

1

p

1

T

from the matrix product (C−M)(C−M)

T

to form a residual matrix R

1

:

R_{1} = (\begin{matrix} 2 & 1 & 0 \\ 1 & 1 / 2 & 0 \\ 0 & 0 & 0 \end{matrix}) - (\begin{matrix} 1 \\ 1 / 2 \\ 0 \end{matrix}) (\begin{matrix} 1 & 1 / 2 & 0 \end{matrix}) = (\begin{matrix} 2 & 1 & 0 \\ 1 & 1 / 2 & 0 \\ 0 & 0 & 0 \end{matrix}) - (\begin{matrix} 1 & 1 / 2 & 0 \\ 1 / 2 & 1 / 4 & 0 \\ 0 & 0 & 0 \end{matrix})

R_{1} = (\begin{matrix} 1 & 1 / 2 & 0 \\ 1 / 2 & 1 / 4 & 0 \\ 0 & 0 & 0 \end{matrix}) .

Another trial eigenvector

p

T=(−1,2,0), orthogonal to the eigenvector

p

1

T=(1,1/2,0) may be used:

((C - M) {(C - M)}^{T} - {\underline{p}}_{1} {\underline{p}}_{1}^{T}) \underline{p} = R_{1} \underline{p} = (\begin{matrix} 1 & 1 / 2 & 0 \\ 1 / 2 & 1 / 4 & 0 \\ 0 & 0 & 0 \end{matrix}) (\begin{matrix} - 1 \\ 2 \\ 0 \end{matrix}) = (\begin{matrix} 0 \\ 0 \\ 0 \end{matrix}) = 0 (\begin{matrix} - 1 \\ 2 \\ 0 \end{matrix}) = λ_{2} {\underline{p}}_{2} .

This indicates that the trial eigenvector

p

T

=(−1,2,0) happened to correspond to the eigenvector

p

2

T=(−1,2,0) belonging to the eigenvalue x

2

=0. The power method then proceeds by subtracting the outer product matrix

p

2

p

2

T

from the residual matrix R

1

to form a second residual matrix R

2

:

R_{2} = (\begin{matrix} 1 & 1 / 2 & 0 \\ 1 / 2 & 1 / 4 & 0 \\ 0 & 0 & 0 \end{matrix}) - (\begin{matrix} - 1 \\ 2 \\ 0 \end{matrix}) (\begin{matrix} - 1 & 2 & 0 \end{matrix}) = (\begin{matrix} 1 & 1 / 2 & 0 \\ 1 / 2 & 1 / 4 & 0 \\ 0 & 0 & 0 \end{matrix}) - (\begin{matrix} 1 & - 2 & 0 \\ - 2 & 4 & 0 \\ 0 & 0 & 0 \end{matrix})

R_{2} = (\begin{matrix} 0 & 5 / 2 & 0 \\ 5 / 2 & - 15 / 4 & 0 \\ 0 & 0 & 0 \end{matrix}) .

Another trial elgenvector

p

T

=(0,0,1), orthogonal to the eigenvectors

p

1

T

=(1,1/2,0) and

p

2

T

(−1,2,0) may be used:

((C - M) {(C - M)}^{T} - {\underline{p}}_{1} {\underline{p}}_{1}^{T} - {\underline{p}}_{2} {\underline{p}}_{2}^{T}) \underline{p} = R_{2} \underline{p} = (\begin{matrix} 1 & 5 / 2 & 0 \\ 5 / 2 & - 15 / 4 & 0 \\ 0 & 0 & 0 \end{matrix}) (\begin{matrix} 0 \\ 0 \\ 1 \end{matrix}) = (\begin{matrix} 0 \\ 0 \\ 0 \end{matrix})

((C - M) {(C - M)}^{T} - {\underline{p}}_{1} {\underline{p}}_{1}^{T} - {\underline{p}}_{2} {\underline{p}}_{2}^{T}) \underline{p} = R_{2} \underline{p} = 0 (\begin{matrix} 0 \\ 0 \\ 1 \end{matrix}) = λ_{3} {\underline{p}}_{3} .

This indicates that the trial eigenvector

p

T

=(0,0,1) happened to correspond to the eigenvector

p

3

T

=(0,0,1) belonging to the eigenvalue λ

3

=0. Indeed, one may readily verify that:

((C - M) {(C - M)}^{T}) {\underline{p}}_{3} = (\begin{matrix} 2 & 1 & 0 \\ 1 & 1 / 2 & 0 \\ 0 & 0 & 0 \end{matrix}) (\begin{matrix} 0 \\ 0 \\ 1 \end{matrix}) = (\begin{matrix} 0 \\ 0 \\ 0 \end{matrix}) = 0 (\begin{matrix} 0 \\ 0 \\ 1 \end{matrix}) = λ_{3} {\underline{p}}_{3} .

Similarly, SVD of C−M shows that C−M=PT

T

, where P is the Principal Component matrix (whose columns are orthonormalized eigenvectors proportional to

p

1

,

p

2

and

p

3

, and whose elements are the Loadings, the direction cosines of the new Principal Component axes

2860

,

2870

and

2880

related to the original variable axes) and T is the Scores matrix (whose rows are the coordinates of the OES data points

2810

and

2820

, referred to the new Principal Component axes

2860

,

2870

and

2880

):

C - M = (\begin{matrix} 2 / \sqrt{5} & - 1 / \sqrt{5} & 0 \\ 1 / \sqrt{5} & 2 / \sqrt{5} & 0 \\ 0 & 0 & 1 \end{matrix}) (\begin{matrix} \sqrt{5} / \sqrt{2} & 0 \\ 0 & 0 \\ 0 & 0 \end{matrix}) (\begin{matrix} 1 / \sqrt{2} & - 1 / \sqrt{2} \\ 1 / \sqrt{2} & 1 / \sqrt{2} \end{matrix})

C - M = (\begin{matrix} 2 / \sqrt{5} & - 1 / \sqrt{5} & 0 \\ 1 / \sqrt{5} & 2 / \sqrt{5} & 0 \\ 0 & 0 & 1 \end{matrix}) (\begin{matrix} \sqrt{5} / 2 & - \sqrt{5} / 2 \\ 0 & 0 \\ 0 & 0 \end{matrix}) = {PT}^{T} .

The transpose of the Scores matrix (T

T

) is given by the product of the matrix of eigenvalues of C−M with a matrix whose rows are orthonormalized eigenvectors proportional to

t

1

and

t

2

. As shown in

FIG. 28

, the direction cosine (Loading) of the first Principal Component axis

2860

with respect to the wavelength

1

counts axis is given by cos Θ

11

=2/{square root over (5)}, and the direction cosine (Loading) of the first Principal Component axis

2860

with respect to the wavelength

2

counts axis is given by cos Θ

21

=1/{square root over (5)}. Similarly, the direction cosine (Loading) of the first Principal Component axis

2860

with respect to the wavelength

3

counts axis is given by cos Θ

12

=cos(π/2)=0. Similarly, the direction cosine (Loading) of the second Principal Component axis

2870

with respect to the wavelength

1

counts axis is given by cos Θ

12

=−1/{square root over (5)}, the direction cosine (Loading) of the second Principal Component axis

2870

with respect to the wavelength

2

counts axis is given by cos Θ

22

=2/{square root over (5)}, and the direction cosine (Loading) of the second Principal Component axis

2870

with respect to the wavelength

3

counts axis is given by cos Θ

32

=cos(π/2)=0. Lastly, the direction cosine (Loading) of the third Principal Component axis

2880

with respect to the wavelength

1

counts axis is given by cos Θ

13

=cos(π/2)=0, the direction cosine (Loading) of the third Principal Component axis

2880

with respect to the wavelength

2

counts axis is given by cos Θ

23

=cos (π/2)=0, and the direction cosine (Loading) of the third Principal Component axis

2880

with respect to the wavelength

3

counts axis is given by cos Θ

33

=cos(0)=1.

SVD confirms that the singular values of C−M are 5/2 and 0, the non-negative square roots of the eigenvalues λ

1

=5/2 and λ

2

=0 of the matrix product (C−M)

T

(C−M). Note that the columns of the Principal Component matrix P are the orthonormalized eigenvectors of the matrix product (C−M)(C−M)

T

.

Taking another example, a 3×4 matrix D (identical to the 3×4 matrix B

T

given above):

D = (\begin{matrix} 1 & 1 & 1 & 1 \\ 1 & 0 & 0 & - 1 \\ 0 & 1 & - 1 & 0 \end{matrix})

may be taken as the overall OES data matrix X (again for the sake of simplicity), corresponding to 4 scans taken at 3 wavelengths. As shown in

FIG. 29

, a scatterplot

2900

of OES data points with coordinates (1,1,0), (1,0,1), (1,0,−1) and (1,−1,0), respectively, may be plotted in a 3-dimensional variable space where the variables are respective spectrometer counts for each of the 3 wavelengths. The mean vector

2920

μ

may lie at the center of a 2-dimensional Principal Component ellipsoid

2930

(really a circle, a somewhat degenerate ellipsoid). The mean vector

2920

μ

may be determined by taking the average of the columns of the overall OES 3×4 matrix D. The Principal Component ellipsoid

2930

may have a first Principal Component

2940

(the “major” axis in

FIG. 29

, with length

2

, lying along a first Principal Component axis

2950

), a second Principal Component

2960

(the “minor” axis in

FIG. 29

, also with length

2

, lying along a second Principal Component axis

2970

), and no third Principal Component lying along a third Principal Component axis

2980

. Here, two of the eigenvalues of the mean-scaled OES data matrix D−M are equal, so the lengths of the “major” and “minor” axes of the Principal Component ellipsoid

2930

in

FIG. 29

are both equal, and the remaining eigenvalue is equal to zero, so the length of the other “minor” axis of the Principal Component ellipsoid

2930

in

FIG. 29

is equal to zero. As shown in

FIG. 29

, the mean vector

2920

μ

is given by:

\underline{μ} = \frac{1}{4} [(\begin{matrix} 1 \\ 1 \\ 0 \end{matrix}) + (\begin{matrix} 1 \\ 0 \\ 1 \end{matrix}) + (\begin{matrix} 1 \\ 0 \\ - 1 \end{matrix}) + (\begin{matrix} 1 \\ - 1 \\ 0 \end{matrix})] = (\begin{matrix} 1 \\ 0 \\ 0 \end{matrix})

and the matrix M has the mean vector

29201

μ

for all 4 columns. As shown in

FIG. 29

, PCA is nothing more than a principal axis rotation of the original variable axes (here, the OES spectrometer counts for 3 wavelengths) about the endpoint of the mean vector

2920

μ

, with coordinates (1,0,0) with respect to the original coordinate axes and coordinates [0,0,0] with respect to the new Principal Component axes

2950

,

2970

and

2980

. The Loadings are merely the direction cosines of the new Principal Component axes

2950

,

2970

and

2980

with respect to the original variable axes. The Scores are simply the coordinates of the OES data points, [1,0,0], [0,1,0], [0,−1,0] and [−1,0,0], respectively, referred to the new Principal Component axes

2950

,

2970

and

2980

.

The 3×3 matrix product (D−M)(D−M) is given by:

(D - M) {(D - M)}^{T} = (\begin{matrix} 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & - 1 \\ 0 & 1 & - 1 & 0 \end{matrix}) (\begin{matrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 0 & 0 & - 1 \\ 0 & - 1 & 0 \end{matrix}) = (\begin{matrix} 0 & 0 & 0 \\ 0 & 2 & 0 \\ 0 & 0 & 2 \end{matrix}) .

The 3×3 matrix (D−M)(D−M)

T

is 3 times the covariance matrix S

3×3

, having elements s

ij

, where i=1,2,3, and j=1,2,3, defined so that:

s_{ij} = \frac{4 \sum_{k = 1}^{4} d_{ik} d_{jk} - \sum_{k = 1}^{4} d_{ik} \sum_{k = 1}^{4} d_{jk}}{4 (4 - 1)},

corresponding to the rectangular 3×4 matrix D

3×4

.

EIG reveals that the eigenvalues of the matrix product (D−M)(D−M)

T

are 0, 2 and 2. The eigenvectors of the matrix product (D−M)(D−M)

T

are solutions p of the equation ((D−M)(D−M)

T

)

p

=λ

p

, and may be seen by inspection to be

p

1

T

=(0,1,0),

p

2

T

=(0,0,1), and

p

3

T

=(1,0,0), belonging to the eigenvalues λ

1

=2, λ

2

=2, and λ

3

=0, respectively (following the convention of placing the largest eigenvalue first).

As may be seen in

FIG. 29

, the direction cosine (Loading) of the first Principal Component axis

2950

with respect to the wavelength

1

counts axis is given by cos Θ

11

=cos(π/2)=0, the direction cosine (Loading) of the first Principal Component axis

2970

with respect to the wavelength

2

counts axis is given by cos Θ

21

=cos(π/2)=1, and the direction cosine (Loading) of the first Principal Component axis

2860

with respect to the wavelength

3

counts axis is given by cos Θ

31

=cos(π/2)=0. Similarly, the direction cosine (Loading) of the second Principal Component axis

2970

with respect to the wavelength

1

counts axis is given by cos Θ

12

=cost(π/2)=0, the direction cosine (Loading) of the second Principal Component axis

2970

with respect to the wavelength

2

counts axis is given by cos Θ

22

=cos(π/2)=0, and the direction cosine (Loading) of the second Principal Component axis

2970

with respect to the wavelength

3

counts axis is given by cos Θ

32

=cos(π/2)=1. Lastly, the direction cosine (Loading) of the third Principal Component axis

2980

with respect to the wavelength

1

counts axis is given by cos Θ

13

=cos(π/2)=1, the direction cosine (Loading) of the third Principal Component axis

2980

with respect to the wavelength

2

counts axis is given by cos Θ

23

=cos(π/2)=0, and the direction cosine (Loading) of the third Principal Component axis

2980

with respect to the wavelength

3

counts axis is given by cos Θ

33

=cos(π/2)=0.

The transpose of the Scores matrix T

T

may be obtained simply by multiplying the mean-scaled OES data matrix D−M on the left by the transpose of the Principal Component matrix P, whose columns are

p

1

,

p

2

,

p

3

, the orthonormalized eigenvectors of the matrix product (D−M)(D−M)

T

:

T^{T} = P^{T} (D - M) = (\begin{matrix} 0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{matrix}) (\begin{matrix} 0 & 0 & 0 & 0 \\ 1 & 0 & 0 & - 1 \\ 0 & 1 & - 1 & 0 \end{matrix}) = (\begin{matrix} 1 & 0 & 0 & - 1 \\ 0 & 1 & - 1 & 0 \\ 0 & 0 & 0 & 0 \end{matrix}) .

The columns of the transpose of the Scores matrix T

T

(or, equivalently, the rows of the Scores matrix T) are, indeed, the coordinates of the OES data points, [1,0,0], [0,1,0], [0,−1,0] and [−1,0,0], respectively, referred to the new Principal Component axes

2950

,

2970

and

2980

.

It has been found that the second Principal Component contains a very robust, high signal-to-noise indicator for etch endpoint determination. The overall mean-scaled OES rectangular n×m data matrix X

nm

−M

nm

may be decomposed into a portion corresponding to the first and second Principal Components and respective Loadings and Scores, and a residual portion:

X - M = P_{PC} T_{PC}^{T} + P_{res} T_{res}^{T} = (\begin{matrix} p_{1} & p_{2} \end{matrix}) (\begin{matrix} t_{1}^{T} \\ t_{2}^{T} \end{matrix}) + P_{res} T_{res}^{T}

where P

PC

is an n×2 matrix X−M=p

1

t

1

T

+p

2

t

2

T

=x

PC

+X

res

, whose columns are the first and second Principal Components, T

PC

is an m×2 Scores matrix for the first and second Principal Components, T

PC

T

is a 2×m Scores matrix transpose, P

PC

T

PC

T

=x

PC

is an n×m matrix, P

res

is an n×(m−2) matrix whose columns are the residual Principal Components, T

res

is an m×(m−2) Scores matrix for the residual Principal Components, T

res

is an (m−2)×m Scores matrix transpose, and P

res

T

res

T

=X

res

is an n×m matrix. The kth column of X−M, x

k

, k=1,2, . . . , m, an n×1 matrix, similarly decomposes into x

k

=(x

PC

)

k

+(x

res

)

k

, where (x

PC

)

k

=P

PC

P

PC

T

x

k

is the projection of x

k

into the Principal Component subspace (PCS) spanned by the first and second Principal Components, (x

res

)

k

=(I

n×n

−P

PC

P

PC

T

)x

k

is the projection of x

k

into the residual subspace orthogonal to the PCS spanned by the first and second Principal Components, (x

PC

)

k

T

=x

k

T

P

PC

P

PC

T

is the projection of x

k

T

into the Principal Component subspace (PCS) spanned by the first and second Principal Components, (x

res

)

k

T

=x

k

T

(I

n×n

−P

PC

P

PC

T

) is the projection of x

k

T

into the residual subspace orthogonal to the PCS spanned by the first and second Principal Components, and I

n×n

is the n×n identity matrix.

Using the projection of x

k

into the PCS spanned by the first and second Principal Components and the projection of x

k

into the residual subspace orthogonal to the PCS spanned by the first and second Principal Components, there are two tests of the amount of variance in the data that are accounted for by the first and second Principal Components. One test is the Q-test, also known as the Squared Prediction Error (SPE) test, where SPE=∥(I

n×n

−P

PC

P

PC

T

)x

k

∥

2

=∥(x

res

)

k

∥

2

≦δ

α

2

. Another test is the Hotelling T

2

test, a multivariate generalization of the well-known Student's t-test, where T

2

=x

k

T

P

PC

Λ

−2

P

PC

T

x

k

=x

k

T

P

PC

P

PC

T

Λ

−2

P

PC

P

PC

T

x

k

=(x

PC

)

k

T

Λ

−2

(x

PC

)

k

≦χ

α

2

, where Λ

2

is 2×2 diagonal matrix of the squares of the eigenvalues λ

i

, i=1,2, belonging to the first and second Principal Components of the overall mean-scaled OES rectangular n×m data matrix X

nm

−M

nm

. Both the SPE test and the Hotelling T

2

test can be used to monitor the etching process, for example.

It has also been found that the first through fourth Principal Components are similarly useful for containing high signal-to-noise indicators for etch endpoint determination as well as being useful for OES data compression. The overall mean-scaled OES rectangular n×m data matrix X

nm

−M

nm

may be decomposed into a portion corresponding to the first through fourth Principal Components and respective Loadings and Scores, and a residual portion:

X - M = P_{PC} T_{PC}^{T} + P_{res} T_{res}^{T} = (\begin{matrix} p_{1} & p_{2} & p_{3} & p_{4} \end{matrix}) (\begin{matrix} t_{1}^{T} \\ t_{2}^{T} \\ t_{3}^{T} \\ t_{4}^{T} \end{matrix}) + P_{res} T_{res}^{T}

whose columns are the first through fourth Principal Components, T

PC

is an m×4 Scores matrix for the first through fourth Principal Components, T

PC

T

is a 4×m Scores matrix transpose, P

PC

T

PC

T

=x

PC

is an n×m matrix, P

res

is an n×(m−4) matrix whose columns are the residual Principal Components, T

res

is an m×(m−4) Scores matrix for the residual Principal Components, T

res

T

is a (m−4)×m Scores matrix transpose, and P

res

T

res

T

=X

res

is an n×m matrix. The kth column of X−M, x

k

, k=1,2, . . . , m, an n×1 matrix, similarly decomposes into x

k

=(x

PC

)

k

+(x

res

)

k

, where (x

PC

)

k

=P

PC

P

PC

T

x

k

is the projection of x

k

into the Principal Component subspace (PCS) spanned by the first through fourth Principal Components, (x

res

)

k

=(I

n×n

−P

PC

P

PC

T

)x

k

is the projection of x

k

into the residual subspace orthogonal to the PCS spanned by the first through fourth Principal Components, (x

PC

)

k

T

=x

k

T

P

PC

P

PC

T

is the projection of x

k

T

into the Principal Component subspace (PCS) spanned by the first through fourth Principal Components, (x

res

)

k

T

=x

k

T

(I

n×n

−P

PC

P

PC

T

) is the projection of x

k

T

into the residual subspace orthogonal to the PCS spanned by the first through fourth Principal Components, and I

n×n

is the n×n identity matrix.

Using the projection of x

k

into the PCS spanned by the first through fourth Principal Components and the projection of x

k

into the residual subspace orthogonal to the PCS spanned by the first through fourth Principal Components, there are again two tests of the amount of variance in the data that are accounted for by the first through fourth Principal Components. One test is the Q-test, also known as the Squared Prediction Error (SPE) test, where SPE=∥(I

n×n

−P

PC

P

PC

T

)x

k

∥

2

=∥(x

res

)

k

∥

2

≦δ

α

2

. Another test is the Hotelling T

2

test, a multivariate generalization of the well-known Student's t-test, where T

2

=x

k

T

P

PC

Λ

−2

P

PC

T

x

k

=x

k

T

P

PC

P

PC

T

Λ

−2

P

PC

P

PC

T

x

k

=(x

PC

)

k

T

Λ

−2

(x

PC

)

k

≦χ

α

2

, where Λ

2

is a 4×4 diagonal matrix of the squares of the eigenvalues λ

i

, i=1,2,3,4 belonging to the first through fourth Principal Components of the overall mean-scaled OES rectangular n×m data matrix X

nm

−M

nm

.

More generally, the overall mean-scaled OES rectangular n×m data matrix X

nm

−M

nm

of rank r, where r≦min{m,n}, may be decomposed into a portion corresponding to the first through rth Principal Components and respective Scores, and a residual portion:

X - M = P_{PC} T_{PC}^{T} + P_{res} T_{res}^{T} = (\begin{matrix} p_{1} & p_{2} & \dots & p_{r - 1} & p_{r} \end{matrix}) (\begin{matrix} t_{1}^{T} \\ t_{2}^{T} \\ ⋮ \\ t_{r - 1}^{T} \\ t_{r}^{T} \end{matrix}) + P_{res} T_{res}^{T}

which expands out to X−M=p

1

t

1

+p

2

t

2

T

+ . . . +p

r−1

t

r−1

T

+p

r

t

r

T

+P

res

T

res

T

=x

PC

+X

res

where P

PC

is an n×r matrix whose columns are the first through rth Principal Components, T

PC

is an m×r Scores matrix for the first through rth Principal Components, T

PC

T

is an r×m Scores matrix transpose, P

PC

T

PC

T

=x

PC

is an n×m matrix, P

res

is an n×(m−r) matrix whose columns are the residual Principal Components (if m=r, P

res

=0), T

res

is an m×(m−r) Scores matrix for the residual Principal Components (if m=r, T

res

=0), T

res

is an (m−r)×m Scores matrix transpose (if m=r, T

res

T

=0), and T

res

P

res

T

=X

res

is an n×m matrix (if m=r, X

res

=0). The kth column of X−M, x

k

, k=1,2, . . . , m, an n×1 matrix, similarly decomposes into x

k

=(x

PC

)

k

+(x

res

)

k

, where (x

PC

)

k

=P

PC

P

PC

T

x

k

is the projection of x

k

into the Principal Component subspace (PCS) spanned by the first through rth Principal Components, (x

res

)

k

=(I

n×n

−P

PC

P

PC

T

)x

k

is the projection of x

k

into the residual subspace orthogonal to the PCS spanned by the first through rth Principal Components, (x

PC

)

k

T

=x

k

T

P

PC

P

PC

T

is the projection of x

k

T

into the Principal Component subspace (PCS) spanned by the first through rth Principal Components, (x

res

)

k

T

=x

k

T

(I

n×n

−P

PC

P

PC

T

) is the projection of x

k

T

into the residual subspace orthogonal to the PCS spanned by the first through rth Principal Components, and I

n×n

is the n×n identity matrix.

Using the projection of x

k

into the PCS spanned by the first through rth Principal Components and the projection of x

k

into the residual subspace orthogonal to the PCS spanned by the first through rth Principal Components, there are likewise two tests of the amount of variance in the data that are accounted for by the first through rth Principal Components. One test is the Q-test, also known as the Squared Prediction Error (SPE) test, where SPE=∥(I

n×n

−P

PC

P

PC

T

)x

k

∥

2

=∥(x

res

)

k

∥

2

≦δ

α

2

. Another test is the Hotelling T

2

test, a multivariate generalization of the well-known Student's t-test, where T

2

=x

k

T

P

PC

Λ

−2

P

PC

T

x

k

=x

k

T

P

PC

P

PC

T

Λ

−2

P

PC

P

PC

T

x

k

=(x

PC

)

k

T

Λ

−2

(x

PC

)

k

≦χ

α

2

, where Λ

2

is an r×r diagonal matrix of the squares of the eigenvalues λ

i

, i=1,2, . . . , r belonging to the first through rth Principal Components of the overall mean-scaled OES rectangular n×m data matrix X

nm

−M

nm

of rank r, where r≦min{m,n}.

In one illustrative embodiment of a method according to the present invention, as shown in

FIGS. 1-7

, archived data sets (Y

n×m

) of OES wavelengths (or frequencies), from wafers that had previously been plasma etched, may be processed and the weighted linear combination of the intensity data, representative of the archived OES wavelengths (or frequencies) collected over time during the plasma etch, defined by the first through pth Principal Components, may be used to compress newly acquired OES data. The rectangular n×m matrix Y (Y

n×m

) may have rank r, where r≦min{m,n} is the maximum number of independent variables in the matrix Y. Here p≦r; in various illustrative embodiments, p is in a range of 1-4; in various alternative illustrative embodiments, p=2. The first through pth Principal Components may be determined from the archived data sets (Y

n×m

) of OES wavelengths (or frequencies), for example, by any of the techniques discussed above.

For example, archived data sets (Y

n×m

) of OES wavelengths (or frequencies), from wafers that had previously been plasma etched, may be processed and Loadings (Q

n×4

) for the first through fourth Principal Components determined from the archived OES data sets (Y

n×m

) may be used as model Loadings (Q

n×4

) to calculate approximate Scores (T

m×4

) corresponding to newly acquired OES data (X

n×m

). These approximate Scores (T

m×4

), along with the mean values for each wavelength (M

n×m

), effectively the column mean vector (

μ

n×1

) of the newly acquired OES data (X

n×m

), may then be stored as compressed OES data.

As shown in

FIG. 1

, a workpiece

100

, such as a semiconducting substrate or wafer, having one or more process layers and/or semiconductor devices such as an MOS transistor disposed thereon, for example, is delivered to an etching preprocessing step j

105

, where j may have any value from j=1 to j N−1. The total number N of processing steps, such as masking, etching, depositing material and the like, used to form the finished workpiece

100

, may range from N=1 to about any finite value.

As shown in

FIG. 2

, the workpiece

100

is sent from the etching preprocessing step j

105

to an etching step j+1

110

. In the etching step j+1

110

, the workpiece

100

is etched to remove selected portions from one or more process layers formed in any of the previous processing steps (such as etching preprocessing step j

105

, where j may have any value from j=1 to j=N−1). As shown in

FIG. 2

, if there is further processing to do on the workpiece

100

(if j<N−1), then the workpiece

100

may be sent from the etching step j+1

110

and delivered to a postetching processing step j+2

115

for further postetch processing, and then sent on from the postetching processing step j+2

115

. Alternatively, the etching step j+1

110

may be the final step in the processing of the workpiece

100

. In the etching step j+1

110

, OES spectra are measured in situ by an OES spectrometer (not shown), producing raw OES data

120

(X

n×m

) indicative of the state of the workpiece

100

during the etching.

In one illustrative embodiment, about 5500 samples of each wafer may be taken on wavelengths between about 240-1100 nm at a high sample rate of about one per second. For example, 5551 sampling points/spectrum/second (corresponding to 1 scan per wafer per second taken at 5551 wavelengths) may be collected in real time, during etching of a contact hole using an Applied Materials AMAT 5300 Centura etching chamber, to produce high resolution and broad band OES spectra.

As shown in

FIG. 3

, the raw OES data

120

(X

n×m

) is sent from the etching step j+1

110

and delivered to a mean-scaling step

125

, producing a means matrix (M

n×m

), whose m columns are each the column mean vector (bend) of the raw OES data

120

(X

n×m

), and mean-scaled OES data (X

n×m

−M

n×m

). In the mean-scaling step

125

, in various illustrative embodiments, the mean values are treated as part of a model built from the archived data sets (Y

n×m

) of OES wavelengths (or frequencies), from wafers that had previously been plasma etched. In other words, a means matrix (N

n×m

) previously determined from the archived data sets (Y

n×m

) of OES wavelengths (or frequencies), from wafers that had previously been plasma etched, is used to generate alternative mean-scaled OES data (X

n×m

−N

n×m

). In various alternative illustrative embodiments, the mean values for each wafer and/or mean value for each wavelength, for example, are determined as discussed above, and are used to generate the mean-scaled OES data (X

n×m

−M

n×m

).

As shown in

FIG. 4

, the means matrix (M

n×m

) and the mean-scaled OES data (X

n×m

−M

n×m

) 130 are sent from the mean scaling step

125

to a Scores calculating step

135

, producing approximate Scores (T

m×p

). In the Scores calculating step

135

, in various illustrative embodiments, the mean-scaled OES data (X

n×m

−M

n×m

) are multiplied on the left by the transpose of the Principal Component (Loadings) matrix Q

n×p

with columns

q

1

,

q

2

, . . . ,

q

p

, that are the first p orthonormalized eigenvectors of the matrix product (Y−N)(Y−N)

T

: (T

T

)

p×m

=(Q

T

)

p×n

(X−M)

n×m

, producing the transpose of the Scores matrix T

T

of the approximate Scores (T

m×p

). The columns of the transpose of the Scores matrix T

T

, or, equivalently, the rows of the approximate Scores matrix (T

m×p

), are the coordinates of the OES data points referred to new approximate Principal Component axes.

In the Scores calculating step

135

, in various illustrative embodiments, the approximate Scores (T

m×p

) are calculated using the Loadings (Q

n×p

) derived from the model built from the archived mean-scaled data sets (Y

n×m

−N

n×m

) of OES wavelengths (or frequencies), from wafers that had previously been plasma etched. In other words, the Loadings (Q

n×p

) previously determined from the archived mean-scaled data sets (Y

n×m

−N

n×m

) of OES wavelengths (or frequencies), from wafers that had previously been plasma etched, are used to generate the approximate Scores (T

m×p

) corresponding to the mean-scaled OES data (X

n×m

−M

n×m

) derived from the raw OES data

120

(X

n×m

).

The Loadings (Q

n×p

) are defined by the first through pth Principal Components. Here p≦r; in various illustrative embodiments, p is in a range of 1-4; in various alternative illustrative embodiments, p=2. The first through pth Principal Components may be determined off-line from the archived data sets (Y

n×m

) of OES wavelengths (or frequencies), for example, by any of the techniques discussed above. The values of the rectangular n×m matrix Y (Y

n×m

) for the archived data sets may be counts representing the intensity of the archived OES spectrum, or ratios of spectral intensities (normalized to a reference intensity), or logarithms of such ratios, for example. The rectangular n×m matrix Y (Y

n×m

) for the archived data sets may have rank r, where r≦min{m,n} is the maximum number of independent variables in the matrix Y. The use of PCA, for example, generates a set of Principal Component Loadings (Q

n×p

), representing contributing spectral components, an eigenmatrix (a matrix whose columns are eigenvectors) of the equation ((Y−N)(Y−N)

T

)Q=Λ

2

P, where N is a rectangular n×m matrix of the mean values of the columns of Y (the m columns of N are each the column mean vector

μ

n×1

of Y

n×m

), Λ

2

is an n×n diagonal matrix of the squares of the eigenvalues λ

i

, i=1,2, . . . , r, of the mean-scaled matrix Y−N, and a Scores matrix, U, with Y−N=QU

T

and (Y−N)

T

=(QU

T

)

T

=(U

T

)

T

Q

T

=UQ

T

, so that ((Y−N)(Y−N)

T

)Q=((QU

T

)(UQ

T

))Q and ((QU

T

)(UQ

T

))Q=(Q(U

T

U)Q

T

)Q=Q(U

T

U)=Λ

2

Q. The rectangular n×m matrix Y, also denoted Y

n×m

, may have elements y

ij

, where i=1,2, . . . , n, and j=1,2, . . . , m, and the rectangular m×n matrix Y

T

, the transpose of the rectangular n×m matrix Y, also denoted (Y

T

)

m×n

, may have elements y

ji

, where i=1,2, . . . , n, and j=1,2, . . . , m. The n×n matrix (Y−N)(Y−N)

T

is (m−1) times the covariance matrix S

n×n

, having elements s

ij

, where m i=1,2, . . . , n, and j=1,2, . . . , n, defined so that:

s_{ij} = \frac{m \sum_{k = 1}^{m} y_{ik} y_{jk} - \sum_{k = 1}^{m} y_{ik} \sum_{k = 1}^{m} y_{jk}}{m (m - 1)},

corresponding to the rectangular n×m matrix Y

n×m

.

As shown in

FIG. 5

, a feedback control signal

140

may be sent from the Scores calculating step

135

to the etching step j+1

110

to adjust the processing performed in the etching step j+1

110

. For example, based on the determination of the approximate Scores (T

m×p

) calculated using the Loadings (Q

n×p

) derived from the model built from the archived mean-scaled data sets (Y

n×m

−N

n×m

) of OES wavelengths (or frequencies), from wafers that had previously been plasma etched, the feedback control signal

140

may be used to signal the etch endpoint.

As shown in

FIG. 6

, the means matrix (M

n×m

) and the approximate Scores (T

m×p

)

145

are sent from the Scores calculating step

135

and delivered to a save compressed PCA data step

150

. In the save compressed PCA data step

150

, the means matrix (M

n×m

) and the approximate Scores (T

m×p

)

145

are saved and/or stored to be used in reconstructing {circumflex over (X)}

n×m

, the decompressed approximation to the raw OES data

120

(X

n×m

). The decompressed approximation {circumflex over (X)}

n×m

to the raw OES data

120

(X

n×m

) may be reconstructed from the means matrix (M

n×m

) and the approximate Scores (T

m×p

)

145

as follows:

In one illustrative embodiment, n=5551, m=100, and p=4, so that the raw OES data

120

(X

5551×100

) requires a storage volume of 5551×100, and generates a means matrix (M

5551×100

) that only requires a storage volume of 5551×1, since all the 100 columns of the means matrix (M

5551×100

) are identical (each of the 5551 rows of each column being the mean value for that wavelength or channel over the 100 scans). The Loadings (Q

5551×4

) are determined off-line from archived data sets (Y

5551×100

) of OES wavelengths (or frequencies), for example, by any of the techniques discussed above, and need not be separately stored with each wafer OES data set, so the storage volume of 5551×4 for the Loadings (Q

5551×4

) does not have to be taken into account in determining an effective compression ratio for the OES wafer data. The approximate Scores (T

100×4

) only require a storage volume of 100×4. Therefore, the effective compression ratio for the OES wafer data in this illustrative embodiment is about (5551×100)/(5551×1) or about 100 to 1 (100:1). More precisely, the compression ratio in this illustrative embodiment is about (5551×100)/(5551×1+100×4) or about 93.3 to 1 (93.3:1).

In another illustrative embodiment, n=5551, m=100, and p=4, so that the raw OES data

120

(X

5551×100

) requires a storage volume of 5551×100, and generates a means matrix (M

5551×100

) that only requires a storage volume of 793×1, since all the 100 columns of the means matrix (M

5551×100

) are identical and the rows are arranged into 793 sets of 7 identical rows (each of the 5551 rows of each column being the mean value for that wavelength or channel over the 100 scans, averaged over the respective 7 wavelengths). The Loadings (Q

5551×4

) are determined off-line from archived data sets (Y

5551×100

) of OES wavelengths (or frequencies), for example, by any of the techniques discussed above, and need not be separately stored with each wafer OES data set, so the storage volume of 5551×4 for the Loadings (Q

5551×4

) does not have to be taken into account in determining an effective compression ratio for the OES wafer data. The approximate Scores (T

100×4

) only require a storage volume of 100×4. Therefore, the effective compression ratio for the OES wafer data in this illustrative embodiment is about (5551×100)/(793×1) or about 700 to 1 (100:1). More precisely, the compression ratio in this illustrative embodiment is about (5551×100)/(793×1+100×4) or about 465 to 1 (465:1).

In yet another illustrative embodiment, n=5551, m=100, and p=4, so that the raw OES data

120

(X

5551×100

) requires a storage volume of 5551×100, and generates a means matrix (M

5551×100

) that only requires a storage volume of 5551×1, since all the 100 columns of the means matrix (M

5551×100

) are identical (each of the 5551 rows of each column being the mean value for that wavelength or channel over the 100 scans). The Loadings (Q

5551×4

) are determined off-line from archived data sets (Y

5551×100

) of OES wavelengths (or frequencies), for example, by any of the techniques discussed above, and require a storage volume of 5551×4, in this alternative embodiment. The approximate Scores (T

100×4

) only require a storage volume of 1 00×4. Therefore, the effective compression ratio for the OES wafer data in this illustrative embodiment is about (5551×100)/(5551×5) or about 20 to 1 (20:1). More precisely, the compression ratio in this illustrative embodiment is about (5551×100)/(5551 ×5+100×4) or about 19.7 to 1 (19.7:1).

In still yet another illustrative embodiment, n=5551, m=100, and p=4, so that the raw OES data

120

(X

5551×100

) requires a storage volume of 5551×100, and generates a means matrix (M

5551×100

) that only requires a storage volume of less than or equal to about 5551×1, since all the 100 columns of the means matrix (M

5551×100

) are identical and means thresholding is used to decimate the rows (each of the 5551 rows of each column being the mean value for that wavelength or channel over the 100 scans, if that mean value is greater than or equal to a specified threshold value, such as a specified threshold value in a range of about 30-50, or zero, otherwise). The Loadings (Q

5551×4

) are determined off-line from archived data sets (Y

5551×100

) of OES wavelengths (or frequencies), for example, by any of the techniques discussed above, and need not be separately stored with each wafer OES data set, so the storage volume of 5551×4 for the Loadings (Q

5551×4

) does not have to be taken into account in determining an effective compression ratio for the OES wafer data. The approximate Scores (T

100×4

) only require a storage volume of 100×4. Therefore, the effective compression ratio for the OES wafer data in this illustrative embodiment is less than or equal to about (5551×100)/(5551×1) or about 100 to 1 (100:1).

As shown in

FIG. 30

, a representative OES trace

3000

of a contact hole etch is illustrated. Time, measured in seconds (sec) is plotted along the horizontal axis against spectrometer counts plotted along the vertical axis. As shown in

FIG. 30

, by about 40 seconds into the etching process, as indicated by dashed line

3010

, the OES trace

3000

of spectrometer counts “settles down” to a range of values less than or about 300, for example. A representative reconstructed OES trace

3020

(corresponding to {circumflex over (X)}

n×m

), for times to the right of the dashed line

3010

(greater than or equal to about 40 seconds, for example), is schematically illustrated and compared with the corresponding noisy raw OES trace

3030

(corresponding to X

n×m

), also for times to the right of the dashed line

3010

. The reconstructed OES trace

3020

(corresponding to {circumflex over (X)}

n×m

) is much smoother and less noisy than the raw OES trace

3030

(corresponding to X

n×m

).

As shown in

FIG. 7

, in addition to, and/or instead of, the feedback control signal

140

, a feedback control signal

155

may be sent from the save compressed PCA data step

150

to the etching step j+1

110

to adjust the processing performed in the etching step j+1

110

. For example, based on the determination of the approximate Scores (T

m×p

) calculated using the Loadings (Q

n×p

) derived from the model built from the archived mean-scaled data sets (Y

n×m

−N

n×m

) of OES wavelengths (or frequencies), from wafers that had previously been plasma etched, the feedback control signal

155

may be used to signal the etch endpoint.

In another illustrative embodiment of a method according to the present invention, as shown in

FIGS. 8-14

, archived data sets (Y

n×m

) of OES wavelengths (or frequencies), from wafers that had previously been plasma etched, may be processed and Loadings (Q

n×4

) for the first through fourth Principal Components determined from the archived OES data sets (Y

n×m

) may be used as model Loadings (Q

n×4

) to calculate approximate Scores (T

m×4

) corresponding to newly acquired OES data (X

n×m

). These approximate Scores (T

m×4

) may be used as an etch endpoint indicator to determine an endpoint for an etch process.

As shown in

FIG. 8

, a workpiece

800

, such as a semiconducting substrate or wafer, having one or more process layers and/or semiconductor devices such as an MOS transistor disposed thereon, for example, is delivered to an etching preprocessing step j

805

, where j may have any value from j=1 to j=N−1. The total number N of processing steps, such as masking, etching, depositing material and the like, used to form the finished workpiece

800

, may range from N=1 to about any finite value.

As shown in

FIG. 9

, the workpiece

800

is sent from the etching preprocessing step j

805

to an etching step j+1

810

. In the etching step j+1

810

, the workpiece

800

is etched to remove selected portions from one or more process layers formed in any of the previous processing steps (such as etching preprocessing step j

805

, where j may have any value from j=1 to j=N−1). As shown in

FIG. 9

, if there is further processing to do on the workpiece

800

(if j<N−1), then the workpiece

800

may be sent from the etching step j+1

810

and delivered to a postetching processing step j+2

815

for further postetch processing, and then sent on from the postetching processing step j+2

815

. Alternatively, the etching step j+1

810

may be the final step in the processing of the workpiece

800

. In the etching step j+1

810

, OES spectra are measured in situ by an OES spectrometer (not shown), producing raw OES data

820

(X

n×m

) indicative of the state of the workpiece

800

during the etching.

In one illustrative embodiment, about 5500 samples of each wafer may be taken on wavelengths between about 240-1100 nm at a high sample rate of about one per second. For example, 5551 sampling points/spectrum/second (corresponding to 1 scan per wafer per second taken at 5551 wavelengths) may be collected in real time, during etching of a contact hole using an Applied Materials AMAT 5300 Centura etching chamber, to produce high resolution and broad band OES spectra.

As shown in

FIG. 10

, the raw OES data

820

(X

n×m

) is sent from the etching step j+1

810

and delivered to a mean-scaling step

825

, producing a means matrix (M

n×m

), whose m columns are each the column mean vector (

μ

n×1

) of the raw OES data

820

(X

n×m

), and mean-scaled OES data (X

n×m

−M

n×m

). In the mean-scaling step

825

, in various illustrative embodiments, the mean values are treated as part of a model built from the archived data sets (Y

n×m

) of OES wavelengths (or frequencies), from wafers that had previously been plasma etched. In other words, a means matrix (N

n×m

) previously determined from the archived data sets (Y

n×m

) of OES wavelengths (or frequencies), from wafers that had previously been plasma etched, is used to generate alternative mean-scaled OES data (X

n×m

−N

n×m

). In various alternative illustrative embodiments, the mean values for each wafer and/or mean value for each wavelength, for example, are determined as discussed above, and are used to generate the mean-scaled OES data (X

n×m

−M

n×m

).

As shown in

FIG. 11

, the means matrix (M

n×m

) and the mean-scaled OES data (X

n×m

−M

n×m

)

830

are sent from the mean scaling step

825

to a Scores calculating step

835

, producing approximate Scores (T

m×p

). In the Scores calculating step

835

, in various illustrative embodiments, the mean-scaled OES data (X

n×m

−M

n×m

) are multiplied on the left by the transpose of the Principal Component (Loadings) matrix Q

n×p

, with columns

q

1

,

q

2

, . . . ,

q

p

, that are the first p orthonormalized eigenvectors of the matrix product (Y−N)(Y−N)

T

: (T

T

)

p×m

=(Q

T

)

p×n

(X−M)

n×m

, producing the transpose of the Scores matrix T

T

of the approximate Scores (T

m×p

). The columns of the transpose of the Scores matrix T

T

, or, equivalently, the rows of the approximate Scores matrix (T

m×p

), are the coordinates of the OES data points referred to new approximate Principal Component axes.

In the Scores calculating step

835

, in various illustrative embodiments, the approximate Scores (T

m×p

) are calculated using the Loadings (Q

n×p

) derived from the model built from the archived mean-scaled data sets (Y

n×m

−N

n×m

) of OES wavelengths (or frequencies), from wafers that had previously been plasma etched. In other words, the Loadings (Q

n×p

), previously determined from the archived mean-scaled data sets (Y

n×m

−N

n×m

) of OES wavelengths (or frequencies), from wafers that had previously been plasma etched, are used to generate the approximate Scores (T

m×p

) corresponding to the mean-scaled OES data (X

n×n

−M

n×m

) derived from the raw OES data

820

(X

n×m

).

The Loadings (Q

n×p

) are defined by the first through pth Principal Components. Here p≦r; in various illustrative embodiments, p is in a range of 1-4; in various alternative illustrative embodiments, p=2. The first through pth Principal Components may be determined off-line from the archived data sets (Y

n×m

) of OES wavelengths (or frequencies), for example, by any of the techniques discussed above. The values of the rectangular n×m matrix Y (Y

n×m

) for the archived data sets may be counts representing the intensity of the archived OES spectrum, or ratios of spectral intensities (normalized to a reference intensity), or logarithms of such ratios, for example. The rectangular n×m matrix Y (Y

n×m

) for the archived data sets may have rank r, where r≦min{m,n} is the maximum number of independent variables in the matrix Y. The use of PCA, for example, generates a set of Principal Component Loadings (Q

n×p

), representing contributing spectral components, an eigenmatrix (a matrix whose columns are eigenvectors) of the equation ((Y−N)(Y−N)

T

)Q=Λ

2

P, where N is a rectangular n×m matrix of the mean values of the columns of Y (the m columns of N are each the column mean vector

μ

n×1

of Y

n×m

), Λ

2

is an n×n diagonal matrix of the squares of the eigenvalues λ

i

, i=1,2, . . . , r, of the mean-scaled matrix Y−N, and a Scores matrix, U, with Y−N=QU

T

and (Y−N)

T

=(QU

T

)

T

=(U

T

)

T

Q

T

=UQ

T

, so that ((Y−N)(Y−N)

T

)Q=((QU

T

)(UQ

T

))Q and ((QU

T

)(UQ

T

))Q=(Q(U

T

U)Q

T

)Q=Q(U

T

U)=Λ

2

Q. The rectangular n×m matrix Y, also denoted Y

n×m

, may have elements y

ij

, where i=1,2, . . . , n, and j=1,2, . . . , m, and the rectangular m×n matrix Y the transpose of the rectangular n×m matrix Y, also denoted (Y

T

)

m×n

, may have elements y

ji

, where i=1,2, . . . , n, and j=1,2, . . . , m. The n×n matrix (Y−N)(Y−N)

T

is (m−1) times the covariance matrix S

n×n

, having elements s

ij

, where i=1,2, . . . , n, and j=1,2, . . . , n, defined so that:

s_{ij} = \frac{m \sum_{k = 1}^{m} y_{ik} y_{jk} - \sum_{k = 1}^{m} y_{ik} \sum_{k = 1}^{m} y_{jk}}{m (m - 1)},

corresponding to the rectangular n×m matrix Y

n×m

.

As shown in

FIG. 12

, a feedback control signal

840

may be sent from the Scores calculating step

835

to the etching step j+1

810

to adjust the processing performed in the etching step j+1

810

. For example, based on the determination of the approximate Scores (T

m×p

) calculated using the Loadings (Q

n×p

) derived from the model built from the archived mean-scaled data sets (Y

n×m

−N

n×m

) of OES wavelengths (or frequencies), from wafers that had previously been plasma etched, the feedback control signal

840

may be used to signal the etch endpoint.

As shown in

FIG. 13

, the approximate Scores (T

m×p

)

845

are sent from the Scores calculating step

835

and delivered to a use Scores as etch indicator step

850

. In the use Scores as etch indicator step

850

, the approximate Scores (T

m×p

)

845

are used as an etch indicator. For example, as shown in

FIG. 27

, a representative Scores time trace

2700

corresponding to the second Principal Component during a contact hole etch is illustrated. Time, measured in seconds (sec) is plotted along the horizontal axis against Scores (in arbitrary units) plotted along the vertical axis. As shown in

FIG. 27

, the Scores time trace

2700

corresponding to the second Principal Component during a contact hole etch may start at a relatively high value initially, decrease with time, pass through a minimum value, and then begin increasing before leveling off. It has been found that the inflection point (indicated by dashed line

2710

, and approximately where the second derivative of the Scores time trace

2700

with respect to time vanishes) is a robust indicator for the etch endpoint.

As shown in

FIG. 14

, in addition to, and/or instead of, the feedback control signal

840

, a feedback control signal

855

may be sent from the use Scores as etch indicator step

850

to the etching step j+1

810

to adjust the processing performed in the etching step j+1

810

. For example, based on the determination of the approximate Scores (T

m×p

) calculated using the Loadings (Q

n×p

) derived from the model built from the archived mean-scaled data sets (Y

n×m

−N

n×m

) of OES wavelengths (or frequencies), from wafers that had previously been plasma etched, the feedback control signal

855

may be used to signal the etch endpoint.

In yet another illustrative embodiment of a method according to the present invention, as shown in

FIGS. 15-21

, archived data sets (Y

n×m

) of OES wavelengths (or frequencies), from wafers that had previously been plasma etched, may be processed and Loadings (Q

n×4

) for the first through fourth Principal Components determined from the archived OES data sets (Y

n×m

) may be used as model Loadings (Q

n×4

) to calculate approximate Scores (T

m×4

) corresponding to newly acquired OES data (X

n×m

). These approximate Scores (T

m×4

), with or without the mean values for each wavelength (N

n×m

), effectively the column mean vector (

μ

n×n

) of the archived OES data (Y

n×m

), may then be stored as compressed OES data. These approximate Scores (T

m×4

) may also be used as an etch endpoint indicator to determine an endpoint for an etch process.

As shown in

FIG. 15

, a workpiece

1500

, such as a semiconducting substrate or wafer, having one or more process layers and/or semiconductor devices such as an MOS transistor disposed thereon, for example, is delivered to an etching preprocessing step j

1505

, where j may have any value from j=1 to j=N−1. The total number N of processing steps, such as masking, etching, depositing material and the like, used to form the finished workpiece

1500

, may range from N=1 to about any finite value.

As shown in

FIG. 16

, the workpiece

1500

is sent from the etching preprocessing step j

1505

to an etching step j+1

1510

. In the etching step j+1

1510

, the workpiece

1500

is etched to remove selected portions from one or more process layers formed in any of the previous processing steps (such as etching preprocessing step j

1505

, where j may have any value from j=1 to j=N−1). As shown in

FIG. 16

, if there is further processing to do on the workpiece

1500

(if j<N−1), then the workpiece

1500

may be sent from the etching step j+1

1510

and delivered to a postetching processing step j+2

1515

for further postetch processing, and then sent on from the postetching processing step j+2

1515

. Alternatively, the etching step j+1

1510

may be the final step in the processing of the workpiece

1500

. In the etching step j+1

1510

, OES spectra are measured in situ by an OES spectrometer (not shown), producing raw OES data

1520

(X

n×m

) indicative of the state of the workpiece

1500

during the etching.

In one illustrative embodiment, about 5500 samples of each wafer may be taken on wavelengths between about 240-1100 nm at a high sample rate of about one per second. For example, 5551 sampling points/spectrum/second (corresponding to 1 scan per wafer per second taken at 5551 wavelengths) may be collected in real time, during etching of a contact hole using an Applied Materials AMAT 5300 Centura etching chamber, to produce high resolution and broad band OES spectra.

As shown in

FIG. 17

, the raw OES data

1520

(X

n×m

) is sent from the etching step j+1

1510

and delivered to a Scores calculating step

1525

, where a means matrix (N

n×m

), whose m columns are each the column mean vector (

μ

n×1

) of the archived OES data (Y

n×m

), is used to produce alternative mean-scaled OES data (X

n×m

−N

n×m

). In the Scores calculating step

1525

, in various illustrative embodiments, the mean values are treated as part of a model built from the archived data sets (Y

n×m

) of OES wavelengths (or frequencies), from wafers that had previously been plasma etched. In other words, a means matrix (N

n×m

) previously determined from the archived data sets (Y

n×m

) of OES wavelengths (or frequencies), from wafers that had previously been plasma etched, is used to generate the alternative mean-scaled OES data (X

n×m

−N

n×m

).

As shown in

FIG. 17

, the alternative mean-scaled OES data (X

n×m

−N

n×m

) are used to produce alternative approximate Scores (T

m×p

). In the Scores calculating step

1525

, in various illustrative embodiments, the mean-scaled OES data (X

n×m

−N

n×m

) are multiplied on the left by the transpose of the Principal Component (Loadings) matrix Q

n×p

, with columns

q

1

,

q

2

, . . . ,

q

p

, that are the first p orthonormalized eigenvectors of the matrix product (Y−N)(Y−N)

T

: (T

T

)

p×m

=(Q

T

)

p×n

(X−N)

n×m

, producing the transpose of the Scores matrix T

T

of the alternative approximate Scores (T

m×p

). The columns of the transpose of the Scores matrix T

T

or, equivalently, the rows of the alternative approximate Scores matrix (T

m×p

), are the alternative coordinates of the OES data points referred to new approximate Principal Component axes.

In the Scores calculating step

1525

, in various illustrative embodiments, the alternative approximate Scores (T

m×p

) are calculated using the Loadings (Q

n×p

) derived from the model built from the archived mean-scaled data sets (Y

n×m

−N

n×m

) of OES wavelengths (or frequencies), from wafers that had previously been plasma etched. In other words, the Loadings (Q

n×p

), previously determined from the archived mean-scaled data sets (Y

n×m

−N

n×m

) of OES wavelengths (or frequencies), from wafers that had previously been plasma etched, are used to generate the alternative approximate Scores (T

m×p

) corresponding to the mean-scaled OES data (X

n×m

−N

n×m

) derived from the raw OES data

1520

(X

n×m

). The Loadings (Q

n×p

) are defined by the first through pth Principal Components. Here p≦r; in various illustrative embodiments, p is in a range of 1-4; in various alternative illustrative embodiments, p=2. The first through pth Principal Components may be determined off-line from the archived data sets (Y

n×m

) of OES wavelengths (or frequencies), for example, by any of the techniques discussed above. The values of the rectangular n×m matrix Y (Y

n×m

) for the archived data sets may be counts representing the intensity of the archived OES spectrum, or ratios of spectral intensities (normalized to a reference intensity), or logarithms of such ratios, for example. The rectangular n×m matrix Y (Y

n×m

) for the archived data sets may have rank r, where r≦min{m,n} is the maximum number of independent variables in the matrix Y. The use of PCA, for example, generates a set of Principal Component Loadings (Q

n×p

), representing contributing spectral components, an eigenmatrix (a matrix whose columns are eigenvectors) of the equation ((Y−N)(Y−N)

T

)Q=Λ

2

P, where N is a rectangular n×m matrix of the mean values of the columns of Y (the m columns of N are each the column mean vector

μ

n×1

of Y

n×m

), Λ

2

is an n×n diagonal matrix of the squares of the eigenvalues λ

i

, i=1,2, . . . , r, of the mean-scaled matrix Y−N, and a Scores matrix, U, with Y−N=QU

T

and (Y−N)

T

=(QU

T

)

T

=(U

T

)

T

Q

T

=UQ

T

, so that ((Y−N)(Y−N)

T

)Q=((QU

T

)(UQ

T

))Q and ((QU

T

)(UQ

T

))Q=(Q(U

T

U)Q

T

)Q=Q(U

T

U)=Λ

2

Q. The rectangular n×m matrix Y, also denoted Y

n×m

, may have elements y

ij

, where i=1,2, . . . , n, and j=1,2, . . . , m, and the rectangular m×n matrix Y

T

, the transpose of the rectangular n×m matrix Y, also denoted (Y

T

)

m×n

, may have elements y

ji

, where i=1,2, . . . , n, and j=1,2, . . . , m. The n×n matrix (Y−N)(Y−N)

T

is (m−1) times the covariance matrix S

n×n

, having elements s

ij

, where i=1,2, . . . , n, and j=1,2, . . . , n, defined so that:

s_{ij} = \frac{m \sum_{k = 1}^{m} y_{ik} y_{jk} - \sum_{k = 1}^{m} y_{ik} \sum_{k = 1}^{m} y_{jk}}{m (m - 1)},

corresponding to the rectangular n×m matrix Y

n×m

.

As shown in

FIG. 18

, a feedback control signal

1530

may be sent from the Scores calculating step

1535

to the etching step j+1

1510

to adjust the processing performed in the etching step j+1

1510

. For example, based on the determination of the alternative approximate Scores (T

m×p

) calculated using the Loadings (Q

n×p

) derived from the model built from the archived mean-scaled data sets (Y

n×m

−N

n×m

) of OES wavelengths (or frequencies), from wafers that had previously been plasma etched, the feedback control signal

1530

may be used to signal the etch endpoint.

As shown in

FIG. 19

, the alternative approximate Scores (T

m×p

)

1535

are sent from the Scores calculating step

1525

and delivered to a save compressed PCA data step

1540

. In the save compressed PCA data step

1540

, the alternative approximate Scores (T

m×p

)

1535

are saved and/or stored to be used in reconstructing {circumflex over (X)}

n×m

, the decompressed alternative approximation to the raw OES data

1520

(X

n×m

). The decompressed alternative approximation {circumflex over (X)}

n×m

to the raw OES data

1520

(X

n×m

) may be reconstructed from the means matrix (N

n×m

) and the alternative approximate Scores (T

m×p

)

1545

as follows: {circumflex over (X)}

n×m

=Q

n×p

(T

T

)

p×m

+N

n×m

.

In one illustrative embodiment, n=5551, m=100, and p=4, so that the raw OES data

1520

(X

5551×100

) requires a storage volume of 5551×100. The means matrix (M

5551×100

) is determined off-line from archived data sets (Y

5551×100

) of OES wavelengths (or frequencies), and need not be separately stored with each wafer OES data set. Thus, the storage volume of 5551×1 for the means matrix (N

5551×100

), where all the 100 columns of the means matrix (N

5551×100

) are identical (each of the 5551 rows of each column being the mean value for that wavelength or channel over the 100 scans), does not have to be taken into account in determining an effective compression ratio for the OES wafer data. The Loadings (Q

5551×4

) are also determined off-line from archived data sets (Y

5551×100

) of OES wavelengths (or frequencies), for example, by any of the techniques discussed above, and also need not be separately stored with each wafer OES data set, so the storage volume of 5551×4 for the Loadings (Q

5551×4

) also does not have to be taken into account in determining an effective compression ratio for the OES wafer data. The approximate Scores (T

100×4

) only require a storage volume of 100×4. Therefore, the effective compression ratio for the OES wafer data in this illustrative embodiment is about (5551×100)/(100×4) or about 1387.75 to 1 (1387.75:1).

As shown in

FIG. 30

, a representative OES trace

3000

of a contact hole etch is illustrated. Time, measured in seconds (sec) is plotted along the horizontal axis against spectrometer counts plotted along the vertical axis. As shown in

FIG. 30

, by about 40 seconds into the etching process, as indicated by dashed line

3010

, the OES trace

3000

of spectrometer counts “settles down” to a range of values less than or about 300, for example. A representative reconstructed OES trace

3020

(corresponding to {circumflex over (X)}

n×m

), for times to the right of the dashed line

3010

(greater than or equal to about 40 seconds, for example), is schematically illustrated and compared with the corresponding noisy raw OES trace

3030

(corresponding to X

n×m

), also for times to the right of the dashed line

3010

. The reconstructed OES trace

3020

(corresponding to {circumflex over (X)}

n×m

) is much smoother and less noisy than the raw OES trace

3030

(corresponding to X

n×m

).

As shown in

FIG. 20

, the alternative approximate Scores (T

m×p

)

1545

are sent from the save compressed PCA data step

1540

and delivered to a use Scores as etch indicator step

1550

. In the use Scores as etch indicator step

1550

, the alternative approximate Scores (T

m×p

)

1545

are used as an etch indicator. For example, as shown in

FIG. 27

, a representative Scores time trace

2700

corresponding to the second Principal Component during a contact hole etch is illustrated. Time, measured in seconds (sec) is plotted along the horizontal axis against Scores (in arbitrary units) plotted along the vertical axis. As shown in

FIG. 27

, the Scores time trace

2700

corresponding to the second Principal Component during a contact hole etch may start at a relatively high value initially, decrease with time, pass through a minimum value, and then begin increasing before leveling off. It has been found that the inflection point (indicated by dashed line

2710

, and approximately where the second derivative of the Scores time trace

2700

with respect to time vanishes) is a robust indicator for the etch endpoint.

As shown in

FIG. 21

, in addition to, and/or instead of, the feedback control signal

1530

, a feedback control signal

1555

may be sent from the use Scores as etch indicator step

1550

to the etching step j+1

1510

to adjust the processing performed in the etching step j+1

1510

. For example, based on the determination of the alternative approximate Scores (T

m×p

) calculated using the Loadings (Q

n×p

) derived from the model built from the archived mean-scaled data sets (Y

n×m

−N

n×m

) of OES wavelengths (or frequencies), from wafers that had previously been plasma etched, the feedback control signal

1555

may be used to signal the etch endpoint.

FIG. 31

illustrates one particular embodiment of a method

3100

practiced in accordance with the present invention.

FIG. 32

illustrates one particular apparatus

3200

with which the method

3100

may be practiced. For the sake of clarity, and to further an understanding of the invention, the method

3100

shall be disclosed in the context of the apparatus

3200

. However, the invention is not so limited and admits wide variation, as is discussed further below.

Referring now to both

FIGS. 31 and 32

, a batch or lot of workpieces or wafers

3205

is being processed through a processing tool

3210

. The processing tool

3210

may be any processing tool known to the art, provided it includes the requisite control capabilities. The processing tool

3210

includes a processing tool controller

3215

for this purpose. The nature and function of the processing tool controller

3215

will be implementation specific. For instance, for etch processing, an etch processing tool controller

3215

may control etch control input parameters such as etch recipe control input parameters and etch endpoint control parameters, and the like. Four workpieces

3205

are shown in

FIG. 32

, but the lot of workpieces or wafers, i.e., the “wafer lot,” may be any practicable number of wafers from one to any finite number.

The method

3100

begins, as set forth in box

3120

, by measuring processing data or parameters (such as OES spectral data for etch processing) characteristic of the processing performed on the workpiece

3205

in the processing tool

3210

. The nature, identity, and measurement of characteristic parameters will be largely implementation specific and even tool specific. For instance, capabilities for monitoring process parameters vary, to some degree, from tool to tool. Greater sensing capabilities may permit wider latitude in the characteristic parameters that are identified and measured and the manner in which this is done. Conversely, lesser sensing capabilities may restrict this latitude.

Turning to

FIG. 32

, in this particular embodiment, the processing data, such as one or more characteristic parameters, is measured and/or monitored by one or more tool sensors (not shown). The output(s) of these tool sensors are transmitted to a computer system

3230

over a line

3220

. The computer system

3230

analyzes these sensor outputs to identify the processing data (characteristic parameters).

Returning, to

FIG. 31

, once the processing data is identified and measured, the method

3100

proceeds by modeling the processing data using a control model, as set forth in box

3130

. The computer system

3230

in

FIG. 32

is, in this particular embodiment, programmed to model the processing data using a control model (such as a PCA control model for an etching process). The manner in which this control modeling occurs will be implementation specific.

In the embodiment of

FIG. 32

, a database

3235

stores a plurality of control models (such as PCA control models) and/or archived control data sets (such as PCA control data sets) that might potentially be applied, depending upon the processing data measured and/or which characteristic parameter is identified. This particular embodiment, therefore, requires some a priori knowledge of the processing data and/or the characteristic parameters that might be measured. The computer system

3230

then extracts an appropriate control model (such as a PCA control model) from the database

3235

of potential control models to apply to the processing data measured and/or the identified characteristic parameters. If the database

3235

does not include an appropriate control model, then processing data and/or the characteristic parameter may be ignored, or the computer system

3230

may attempt to develop one, if so programmed. The database

3235

may be stored on any kind of computer-readable, program storage medium, such as an optical disk

3240

, a floppy disk

3245

, or a hard disk drive (not shown) of the computer system

3230

. The database

3235

may also be stored on a separate computer system (not shown) that interfaces with the computer system

3230

.

Modeling of the processing data measured and/or the identified characteristic parameter may be implemented differently in alternative embodiments. For instance, the computer system

3230

may be programmed using some form of artificial intelligence to analyze the sensor outputs and controller inputs to develop a control model (such as a PCA control model) on-the-fly in a real-time implementation. This approach might be a useful adjunct to the embodiment illustrated in

FIG. 32

, and discussed above, where processing data and/or characteristic parameters are measured and identified for which the database

3235

has no appropriate control model.

The method

3100

of

FIG. 31

then proceeds by applying the control model (such as the PCA control model) to compress the processing data (such as OES spectral data for etch processing) and/or control the processing (such as determine an etch endpoint for etch processing), as set forth in box

3140

. The processing data compressed using a control model according to any of the various illustrative embodiments of the present invention may be stored on any kind of computer-readable, program storage medium, such as an optical disk

3240

, a floppy disk

3245

, or a hard disk drive (not shown) of the computer system

3230

, and/or together with the database

3235

. The processing data compressed using a control model according to any of the various illustrative embodiments of the present invention may also be stored on a separate computer system (not shown) that interfaces with the computer system

3230

.

Depending on the implementation, applying the control model may yield either a new value for the processing control parameter (such as the etch endpoint control parameter for etch processing) or a correction and/or update to the existing processing control parameter. The new processing control parameter is then formulated from the value yielded by the control model and is transmitted to the processing tool controller

3215

over the line

3220

. The processing tool controller

3215

then controls subsequent process operations in accordance with the new processing control inputs.

Some alternative embodiments may employ a form of feedback to improve the control modeling of processing data and/or characteristic parameters. The implementation of this feedback is dependent on several disparate facts, including the tool's sensing capabilities and economics. One technique for doing this would be to monitor at least one effect of the control model's implementation and update the control model based on the effect(s) monitored. The update may also depend on the control model. For instance, a linear control model may require a different update than would a non-linear control model, all other factors being the same.

As is evident from the discussion above, some features of the present invention are implemented in software. For instance, the acts set forth in the boxes

3120

-

3140

in

FIG. 31

are, in the illustrated embodiment, software-implemented, in whole or in part. Thus, some features of the present invention are implemented as instructions encoded on a computer-readable, program storage medium. The program storage medium may be of any type suitable to the particular implementation. However, the program storage medium will typically be magnetic, such as the floppy disk

3245

or the computer

3230

hard disk drive (not shown), or optical, such as the optical disk

3240

. When these instructions are executed by a computer, they perform the disclosed functions. The computer may be a desktop computer, such as the computer

3230

. However, the computer might alternatively be a processor embedded in the processing tool

3210

. The computer might also be a laptop, a workstation, or a mainframe in various other embodiments. The scope of the invention is not limited by the type or nature of the program storage medium or computer with which embodiments of the invention might be implemented.

Thus, some portions of the detailed descriptions herein are, or may be, presented in terms of algorithms, functions, techniques, and/or processes. These terms enable those skilled in the art most effectively to convey the substance of their work to others skilled in the art. These terms are here, and are generally, conceived to be a self-consistent sequence of steps leading to a desired result. The steps are those requiring physical manipulations of physical quantities. Usually, though not necessarily, these quantities take the form of electromagnetic signals capable of being stored, transferred, combined, compared, and otherwise manipulated.

It has proven convenient at times, principally for reasons of common usage, to refer to these signals as bits, values, elements, symbols, characters, terms, numbers, and the like. All of these and similar terms are to be associated with the appropriate physical quantities and are merely convenient labels applied to these quantities and actions. Unless specifically stated otherwise, or as may be apparent from the discussion, terms such as “processing,” “computing,” “calculating,” “determining,” “displaying,” and the like, used herein refer to the action(s) and processes of a computer system, or similar electronic and/or mechanical computing device, that manipulates and transforms data, represented as physical (electromagnetic) quantities within the computer system's registers and/or memories, into other data similarly represented as physical quantities within the computer system's memories and/or registers and/or other such information storage, transmission and/or display devices.

Any of these illustrative embodiments may be applied in real-time etch processing. Alternatively, either of the illustrative embodiments shown in

FIGS. 8-14

and

15

-

21

may be used as an identification technique when using batch etch processing, with archived data being applied statistically, to determine an etch endpoint for the batch.

In various illustrative embodiments, a process engineer may be provided with advanced process data monitoring capabilities, such as the ability to provide historical parametric data in a user-friendly format, as well as event logging, real-time graphical display of both current processing parameters and the processing parameters of the entire run, and remote, i.e., local site and worldwide, monitoring. These capabilities may engender more optimal control of critical processing parameters, such as throughput accuracy, stability and repeatability, processing temperatures, mechanical tool parameters, and the like. This more optimal control of critical processing parameters reduces this variability. This reduction in variability manifests itself as fewer within-run disparities, fewer run-to-run disparities and fewer tool-to-tool disparities. This reduction in the number of these disparities that can propagate means fewer deviations in product quality and performance. In such an illustrative embodiment of a method of manufacturing according to the present invention, a monitoring and diagnostics system may be provided that monitors this variability and optimizes control of critical parameters.

An etch endpoint determination signal as in any of the embodiments disclosed above may have a high signal-to-noise ratio and may be reproducible over the variations of the incoming wafers and the state of the processing chamber, for example, whether or not the internal hardware in the processing chamber is worn or new, or whether or not the processing chamber is in a “clean” or a “dirty” condition. Further, in particular applications, for example, the etching of contact and/or via holes, an etch endpoint determination signal as in any of the embodiments disclosed above may have a high enough signal-to-noise ratio to overcome the inherently very low signal-to-noise ratio that may arise simply by virtue of the small percentage (1% or so) of surface area being etched.

In various illustrative embodiments, the etch endpoint signal becomes very stable, and may throughput may be improved by reducing the main etch time from approximately 145 seconds, for example, to approximately 90-100 seconds, depending on the thickness of the oxide. In the absence of an etch endpoint determination signal as in any of the embodiments disclosed above, a longer etch time is conventionally needed to insure that all the material to be etched away has been adequately etched away, even in vias and contact holes with high aspect ratios. The presence of a robust etch endpoint determination signal as in any of the embodiments disclosed above thus allows for a shorter etch time, and, hence, increased throughput, compared to conventional etching processes.

Thus, embodiments of the present invention fill a need in present day and future technology for optimizing selection of wavelengths to monitor for endpoint determination or detection during etching. Similarly, embodiments of the present invention fill a need in present day and future technology for being able to determine an etch endpoint expeditiously, robustly, reliably and reproducibly under a variety of different conditions, even in real-time processing.

Further, it should be understood that the present invention is applicable to any plasma etching system, including reactive ion etching (RIE), high-density, inductively coupled plasma (ICP) systems, electron cyclotron resonance (ECR) systems, radio frequency induction systems, and the like.

Data compression of OES spectra as in any of the embodiments disclosed above may solve the set of problems is posed by the sheer number of OES frequencies or wavelengths available to monitor. The monitoring typically generates a large amount of data. For example, a data file for each wafer monitored may be as large as 2-3 megabytes (MB), and each etcher can typically process about 500-700 wafers per day. Conventional storage methods would require over a gigabytes (GB) of storage space per etcher per day and over 365 GB per etcher per year. Further, the raw OES data generated in such monitoring is typically “noisy” and unenhanced. Compression ratios of 100:1, as in various of the illustrative embodiments disclosed above, would only require tens of MB of storage per etcher per day and only about 4 GB per etcher per year. In addition, data compression and reconstruction of OES spectra as in any of the embodiments disclosed above may smooth and enhance the otherwise noisy and unenhanced raw OES data generated in etch monitoring.

Data compression of OES spectra as in any of the embodiments disclosed above may feature high compression ratios for the raw OES data ranging from about 20:1 to about 400:1, provide for efficient noise removal from the raw OES data and preserve relevant features of raw OES data. Data compression of OES spectra as in any of the embodiments disclosed above may also have fast computation characteristics, such as fast compression and fast reconstruction, so as to be suitable for robust on-line, real-time implementation. Data compression of OES spectra as in any of the embodiments disclosed above may further be compatible with real-time, PCA-based fault detection and classification (FDC), improve the accuracy and efficiency of etch processing, simplify manufacturing, lower overall costs and increase throughput.

As discussed above, in various illustrative embodiments, control systems may be provided that use modeling of data in their functionality. Control models may deal with large volumes of data. This large volume of data may be fit by the control model so that a much smaller data set of control model parameters is created. When the much smaller data set of control model parameters is put into the control model, the original much larger data set is accurately recreated (to within a specified tolerance and/or approximation). These control models are already created for control purposes, so the resulting control model parameters may be stored to sufficiently capture, in the compressed control model parameter data set, the original voluminous data set with much smaller storage requirements. For example, fault detection models, in particular, deal with large volumes of data. This large volume of data may be fit by the fault detection model so that a much smaller data set of fault detection model parameters is created. When the much smaller data set of fault detection model parameters is put into the fault detection model, the original much larger data set is accurately recreated (to within a specified tolerance and/or approximation).

As discussed above, in various alternative illustrative embodiments, control systems that use modeling of data in their functionality may be built. In these various alternative illustrative embodiments, building the control model comprises fitting the collected processing data using at least one of polynomial curve fitting, least-squares fitting, polynomial least-squares fitting, non-polynomial least-squares fitting, weighted least-squares fitting, weighted polynomial least-squares fitting, weighted non-polynomial least-squares fitting, Partial Least Squares (PLS), and Principal Components Analysis (PCA).

As discussed above, in various illustrative embodiments, samples may be collected for N+1 data points (x

i

,y

i

), where i=1, 2, . . . , N, N+1, and a polynomial of degree N,

P_{N} (x) = a_{0} + a_{1} x + a_{2} x^{2} + \dots + a_{k} x^{k} + \dots + a_{N} x^{N} = \sum_{k = 0}^{N} a_{k} x^{k},

may be fit to the N+1 data points (x

i

,y

i

). For example, 100 time data points (N=99) may be taken relating the pressure p, the flowrate f, and/or the temperature T, of a process gas to the thickness t of a process layer being deposited in a deposition processing step, resulting in respective sets of N+1 data points (p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

). For each of these sets of N+1 data points (p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

), the conventional raw storage requirements would be on the order of 2(N+1). By way of comparison, data compression, according to various illustrative embodiments of the present invention, using polynomial interpolation, for example, would reduce the storage requirements, for each of these sets of N+1 data points (p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

), to the storage requirements for the respective coefficients a

k

, for k=0, 1, . . . , N, a storage savings of a factor of 2 (equivalent to a data compression ratio of 2:1) for each of these sets of N+1 data points (p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

). Polynomial interpolation is described, for example, in

Numerical Methods for Scientists and Engineers

, by R. W. Hamming, Dover Publications, New York, 1986, at pages 230-235. The requirement that the polynomial P

N

(X) pass through the N+1 data points (x

i

,y

i

) is

y_{i} = P_{N} (x_{i}) = \sum_{k = 0}^{N} a_{k} x_{i}^{k}, for i = 1, 2, \dots, N, N + 1,

Vandermonde determinant : V_{N + 1} = &LeftBracketingBar; \begin{matrix} 1 & x_{1} & x_{1}^{2} & \dots & x_{1}^{N} \\ 1 & x_{2} & x_{2}^{2} & \dots & x_{2}^{N} \\ ⋮ & ⋮ & ⋮ & ⋰ & ⋮ \\ 1 & x_{N} & x_{N}^{2} & \dots & x_{N}^{N} \\ 1 & x_{N + 1} & x_{N + 1}^{2} & \dots & x_{N + 1}^{N} \end{matrix} &RightBracketingBar; = &LeftBracketingBar; x_{i}^{k} &RightBracketingBar;, where i = 1, 2, \dots, N + 1,

and k=0, 1, . . . , N. The Vandermonde determinant V

N+1

, considered as a function of the variables x

i

, V

N+1

=V

N+1

(x

1

,x

2

, . . . , x

N

,x

N+1

), is clearly a polynomial in the variables x

i

, as may be seen by expanding out the determinant, and a count of the exponents shows that the degree of the polynomial is

0 + 1 + 2 + 3 + \dots + k + \dots + N = \sum_{k = 0}^{N} k = \frac{N (N + 1)}{2}

\prod_{j = 1}^{N} (x_{N + 1} - x_{j}) .

Similarly, if x

N

=x

j

, for j=1, 2, . . . , N−1, then the Vandermonde determinant V

N+1

=0, so the Vandermonde determinant V

N+1

must also have the factors (x

N

−x

j

), for j=1, 2, . . . , N−1, corresponding to the N−1 factors

\prod_{j = 1}^{N} (x_{N} - x_{j}) .

Generally, if x

m

=x

j

, for j<m, where m=2, . . . , N, N+1, then the Vandermonde determinant V

N+1

=0, so the Vandermonde determinant V

N+1

must have all the factors (x

m

−x

j

), for j<m, where m=2, . . . , N, N+1, corresponding to the factors

\prod_{m > j = 1}^{N + 1} (x_{m} - x_{j}) .

Altogether, this represents a polynomial of degree

N + (N - 1) + \dots + k + \dots + 2 + 1 = \sum_{k = 1}^{N} k = \frac{N (N + 1)}{2},

since, when m=N+1, for example, j may take on any of N values, j=1, 2, . . . , N, and when m=N, j may take on any of N−1 values, j=1, 2, . . . , N−1, and so forth (for example, when m=3, j may take only two values, j=1, 2, and when m=2, j may take only one value, j=1), which means that all the factors have been accounted for and all that remains is to find any multiplicative constant by which these two representations for the Vandermonde determinant V

N+1

might differ. As noted above, the diagonal term of the Vandermonde determinant V

N+1

is 1·x

2

·x

3

. . . x

N

N−1

·x

N+1

N

, and this may be compared to the term from the left-hand sides of the product of factors

\prod_{m > j = 1}^{N + 1} (x_{m} - x_{j}) = \prod_{j = 1}^{N} (x_{N + 1} - x_{j}) \prod_{j = 1}^{N - 1} (x_{N} - x_{j}) \dots \prod_{j = 1}^{2} (x_{3} - x_{j}) \prod_{j = 1}^{1} (x_{2} - x_{j}), x_{N + 1}^{N} \cdot x_{N}^{N - 1} \dots x_{3}^{2} \cdot x_{2},

which is identical, so the multiplicative constant is unity and the Vandermonde determinant

V_{N + 1} is V_{N + 1} (x_{1}, x_{2}, \dots, x_{N + 1}) = &LeftBracketingBar; x_{i}^{k} &RightBracketingBar; = \prod_{m > j = 1}^{N + 1} (x_{m} - x_{j}) .

P_{N} (x) = \sum_{k = 0}^{N} a_{k} x^{k},

and rearranging suitably gives the determinant equation

&LeftBracketingBar; \begin{matrix} y & 1 & x & x^{2} & \dots & x^{N} \\ y_{1} & 1 & x_{1} & x_{1}^{2} & \dots & x_{1}^{N} \\ y_{2} & 1 & x_{2} & x_{2}^{2} & \dots & x_{2}^{N} \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋰ & ⋮ \\ y_{N} & 1 & x_{N} & x_{N}^{2} & \dots & x_{N}^{N} \\ y_{N + 1} & 1 & x_{N + 1} & x_{N + 1}^{2} & \dots & x_{N + 1}^{N} \end{matrix} &RightBracketingBar; = 0,

y = P_{N} (x) = \sum_{k = 0}^{N} a_{k} x^{k}

y_{i} = P_{N} (x_{i}) = \sum_{k = 0}^{N} a_{k} x_{i}^{k}, for i = 1, 2, \dots, N, N + 1,

is satisfied.

For example, a quadratic curve may be found that goes through the sample data set (−1,a), (0,b), and (1,c). The three equations are P

2

(−1)=a=a

0

−a

1

+a

2

, P

2

(0)=b=a

0

, and P

2

(1)=c=a

0

+a

1

+a

2

which imply that b=a

0

, c−a=2a

1

, and c+a−2b=2a

2

, so that

y (x) = P_{2} (x) = b + \frac{c - a}{2} x + \frac{c + a - 2 b}{2} x^{2},

which is also the result of expanding

&LeftBracketingBar; \begin{matrix} y & 1 & x & x^{2} \\ a & 1 & - 1 & 1 \\ b & 1 & 0 & 0 \\ c & 1 & 1 & 1 \end{matrix} &RightBracketingBar; = 0 = y &LeftBracketingBar; \begin{matrix} 1 & - 1 & 1 \\ 1 & 0 & 0 \\ 1 & 1 & 1 \end{matrix} &RightBracketingBar; - 1 &LeftBracketingBar; \begin{matrix} a & - 1 & 1 \\ b & 0 & 0 \\ c & 1 & 1 \end{matrix} &RightBracketingBar; + x &LeftBracketingBar; \begin{matrix} a & 1 & 1 \\ b & 1 & 0 \\ c & 1 & 1 \end{matrix} &RightBracketingBar; - x^{2} &LeftBracketingBar; \begin{matrix} a & 1 & - 1 \\ b & 1 & 0 \\ c & 1 & 1 \end{matrix} &RightBracketingBar;,

the coefficient of y being the respective Vandermonde determinant V

3

=2.

Similarly, a quartic curve may be found that goes through the sample data set (−2,a), (−1,b), (0,c), (1,b), and (2,a). The five equations are P

4

(−2)=a=a

0

−2a

1

+4a

2

−8a

3

+16a

4

, P

4

(−1)=b=a

0

−a

1

+a

2

−a

3

+a

4

, P

4

(0)=c=a

0

, P

4

(1)=b=a

0

+a

1

+a

2

+a

3

+a

4

, and P

4

(2)=a=a

0

+2a

1

+4a

2

+8a

3

+16a

4

, which imply that c=a

0

, 0=a

1

=a

3

(which also follows from the symmetry of the data set), (a−c)−16(b−c)=−12a

2

, and (a−c)−4(b−c)=12a

4

, so that

y (x) = P_{4} (x) = c - \frac{a - 16 b + 15 c}{12} x^{2} + \frac{a - 4 b + 3 c}{12} x^{4} .

As discussed above, in various alternative illustrative embodiments, samples may be collected for M data points (x

i

,y

i

), where i=1, 2, . . . , M, and a first degree polynomial (a straight line),

P_{1} (x) = a_{0} + a_{1} x = \sum_{k = 0}^{1} a_{k} x^{k},

may be fit (in a least-squares sense) to the M data points (x

1

,y

i

). For example, 100 time data points (M=100) may be taken relating the pressure p, the flowrate f, and/or the temperature T, of a process gas to the thickness t of a process layer being deposited in a deposition processing step, resulting in the M data points (p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

), respectively. For each of these sets of M data points (p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

), the conventional raw storage requirements would be on the order of 2M. By way of comparison, data compression, according to various illustrative embodiments of the present invention, using a first degree polynomial (a straight line),

P_{1} (x) = a_{0} + a_{1} x = \sum_{k = 0}^{1} a_{k} x^{k},

fit (in a least-squares sense) to the M data points (x

i

,y

i

), only requires two parameters, a

0

and a

1

, which would reduce the storage requirements, for each of these sets of M data points (p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

), to the storage requirements for the respective coefficients a

k

, for k=0, 1, a storage savings of a factor of 2M/2 (equivalent to a data compression ratio of M:1) for each of these sets of M data points (p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

). Least-squares fitting is described, for example, in

Numerical Methods for Scientists and Engineers

, by R. W. Hamming, Dover Publications, New York, 1986, at pages 427-443.

The least-squares criterion may be used in situations where there is much more data available than parameters so that exact matching (to within round-off) is out of the question. Polynomials are most commonly used in least-squares matching, although any linear family of suitable functions may work as well. Suppose some quantity x is being measured by making M measurements x

i

, for i=1, 2, . . . , M, and suppose that the measurements x

i

, are related to the “true” quantity x by the relation x

i

=x+ε

i

, for i=1, 2, . . . , M, where the residuals ε

i

are regarded as noise. The principle of least-squares states that the best estimate ξ of the true value x is the number that minimizes the sum of the squares of the deviations of the data from their estimate

f (ξ) = \sum_{i = 1}^{M} ε_{i}^{2} = \sum_{i = 1}^{M} {(x_{i} - ξ)}^{2},

which is equivalent to the assumption that the average x

a

, where

x_{a} = \frac{1}{M} \sum_{i = 1}^{M} x_{i},

is the best estimate ξ of the true value x. This equivalence may be shown as follows. First, the principle of least-squares leads to the average x

a

. Regarding

f (ξ) = \sum_{i = 1}^{M} ε_{i}^{2} = \sum_{i = 1}^{M} {(x_{i} - ξ)}^{2}

as a function of the best estimate ξ, minimization with respect to the best estimate ξ may proceed by differentiation:

\frac{ⅆ f (ξ)}{ⅆ ξ} = - 2 \sum_{i = 1}^{M} (x_{i} - ξ) = 0,

which implies that

\sum_{i = 1}^{M} x_{i} - \sum_{i = 1}^{M} ξ = 0 = \sum_{i = 1}^{M} x_{i} - M ξ,

so that

ξ = \frac{1}{M} \sum_{i = 1}^{M} x_{i} = x_{a},

or, in other words, that the choice x

a

=ξ minimizes the sum of the squares of the residuals ε

i

. Noting also that

\frac{ⅆ^{2} f (ξ)}{ⅆ ξ^{2}} = 2 \sum_{i = 1}^{M} 1 = 2 M > 0,

the criterion for a minimum is established.

Conversely, if the average x

a

is picked as the best choice x

a

=ξ, it can be shown that this choice, indeed, minimizes the sum of the squares of the residuals ε

i

. Set

\begin{matrix} f (x_{a}) = \sum_{i = 1}^{M} {(x_{i} - x_{a})}^{2} \\ = \sum_{i = 1}^{M} x_{i}^{2} - 2 x_{a} \sum_{i = 1}^{M} x_{i} + \sum_{i = 1}^{M} x_{a}^{2} \\ = \sum_{i = 1}^{M} x_{i}^{2} - 2 x_{a} M x_{a} + M x_{a}^{2} \\ = \sum_{i = 1}^{M} x_{i}^{2} - M x_{a}^{2} . \end{matrix}

If any other value x

b

is picked, then, plugging that other value x

b

into f(x) gives

\begin{matrix} f (x_{b}) = \sum_{i = 1}^{M} {(x_{i} - x_{b})}^{2} \\ = \sum_{i = 1}^{M} x_{i}^{2} - 2 x_{b} \sum_{i = 1}^{M} x_{i} + \sum_{i = 1}^{M} x_{b}^{2} \\ = \sum_{i = 1}^{M} x_{i}^{2} - 2 x_{b} M x_{a} + M x_{b}^{2} . \end{matrix}

f (χ) = \sum_{i = 1}^{M} &LeftBracketingBar; ε_{i} &RightBracketingBar; = \sum_{i = 1}^{M} &LeftBracketingBar; x_{i} - χ &RightBracketingBar;,

\frac{x_{\min} + x_{\max}}{2} = x_{midrange},

P_{1} (x) = a_{0} + a_{1} x = \sum_{k = 0}^{1} a_{k} x^{k},

may be fit (in a least-squares sense) to the M data points (x

i

,y

i

), there are two parameters, a

0

and a

1

, and a function F(a

0

,a

1

) that needs to be minimized as follows. The function F(a

0

,a

1

) is given by

F (a_{0}, a_{1}) = \sum_{i = 1}^{M} ε_{i}^{2} = \sum_{i = 1}^{M} {[P_{1} (x_{i}) - y_{i}]}^{2} = \sum_{i = 1}^{M} {[a_{0} + a_{1} x_{i} - y_{i}]}^{2}

and setting the partial derivatives of F(a

0

,a

1

) with respect to a

0

and a

1

equal to zero gives

\frac{\partial F (a_{0}, a_{1})}{\partial a_{0}} = 2 \sum_{i = 1}^{M} [a_{0} + a_{1} x_{i} - y_{i}] = 0 and

\frac{\partial F (a_{0}, a_{1})}{\partial a_{1}} = 2 \sum_{i = 1}^{M} [a_{0} + a_{1} x_{i} - y_{i}] x_{i} = 0,

respectively. Simplifying and rearranging gives

a_{0} M + a_{1} \sum_{i = 1}^{M} x_{i} = \sum_{i = 1}^{M} y_{i} and a_{0} \sum_{i = 1}^{M} x_{i} + a_{1} \sum_{i = 1}^{M} x_{i}^{2} = \sum_{i = 1}^{M} x_{i} y_{i},

respectively, where there are two equations for the two unknown parameters a

0

and a

1

, readily yielding a solution.

As shown in

FIG. 33

, for example, a first degree polynomial (a straight line),

P_{1} (x) = a_{0} + a_{1} x = \sum_{k = 0}^{1} a_{k} x^{k},

may be fit (in a least-squares sense) to the M=5 data points (1,0), (2,2), (3,2), (4,5), and (5,4). The residuals ε

i

, for i=1, 2, . . . , 5, are schematically illustrated in FIG.

33

. The equations for the two parameters a

0

and a

1

are 5a

0

+15a

1

=13 and 15a

0

+55a

1

=50, respectively, so that, upon multiplying the first equation by 3 and then subtracting that from the second equation, thereby eliminating a

0

, the solution for the parameter a

1

becomes a

1

=11/10, and this implies that the solution for the parameter a

0

becomes a

0

=−7/10. The first degree polynomial (the straight line) that provides the best fit, in the least-squares sense, is

P_{1} (x) = - \frac{7}{10} + \frac{11}{10} x = \frac{1}{10} (- 7 + 11 x),

as shown in FIG.

33

.

As shown in

FIG. 34

, for example, a first degree polynomial (a straight line),

P_{1} (x) = a_{0} + a_{1} x = \sum_{k = 0}^{1} a_{k} x^{k},

may be fit (in a least-squares sense) to the M=7 data points (−3,4), (−2,4), (−1,2), (0,2), (1,1), (2,0), and (3,0). The residuals ε

i

, for i=1, 2, . . . , 7, are schematically illustrated in FIG.

34

. The equations for the two parameters a

0

and a

1

are

a_{0} M + a_{1} \sum_{i = 1}^{M} x_{i} = \sum_{i = 1}^{M} y_{i} = 7 a_{0} + a_{1} (- 3 - 2 - 1 + 0 + 1 + 2 + 3) = (4 + 4 + 2 + 2 + 1 + 0 + 0) and

a_{0} \sum_{i = 1}^{M} x_{i} + a_{1} \sum_{i = 1}^{M} x_{i}^{2} = \sum_{i = 1}^{M} x_{i} y_{i} = a_{1} (9 + 4 + 1 + 0 + 1 + 4 + 9) = (- 12 - 8 - 2 + 0 + 1 + 0 + 0),

respectively, which give 7a

0

=13 and 28a

1

=−1, respectively. In other words, a

0

=13/7 and a

1

=−3/4, so that, the first degree polynomial (the straight line) that provides the best fit, in the least-squares sense, is

P_{1} (x) = \frac{13}{7} - \frac{3}{4} x,

as shown in FIG.

34

.

As discussed above, in various other alternative illustrative embodiments, samples may be collected for M data points (x

i

,y

i

), where i=1, 2, . . . , M, and a polynomial of degree N,

P_{N} (x) = a_{0} + a_{1} x + a_{2} x^{2} + \dots + a_{k} x^{k} + \dots + a_{N} x^{N} = \sum_{k = 0}^{N} a_{k} x^{k},

\begin{matrix} F (a_{0}, a_{1}, \dots, a_{N}) = \sum_{i = 1}^{M} ε_{i}^{2} \\ = \sum_{i = 1}^{M} {[P_{N} (x_{i}) - y_{i}]}^{2} \end{matrix}

and setting the partial derivatives of F(a

0

,a

1

, . . . , a

N

) with respect to a

j

, for j=0, 1, . . . , N, equal to zero gives

\begin{matrix} \frac{\partial F (a_{0}, a_{1}, \dots, a_{N})}{\partial a_{j}} = 2 \sum_{i = 1}^{M} [P_{N} (x_{i}) - y_{i}] x_{i}^{j} \\ = 2 \sum_{i = 1}^{M} [\sum_{k = 0}^{N} a_{k} x_{i}^{k} - y_{i}] x_{i}^{j} \\ = 0, \end{matrix}

for j=0, 1, . . . , N, since (x

i

)

j

is the coefficient of a

j

in the polynomial

P_{N} (x_{i}) = \sum_{k = 0}^{N} a_{k} x_{i}^{k} .

Simplifying and rearranging gives

\begin{matrix} \sum_{i = 1}^{M} [\sum_{k = 0}^{N} a_{k} x_{i}^{k}] x_{i}^{j} = \sum_{k = 0}^{N} a_{k} [{\sum x_{i}^{k + j}}_{i = 1}^{M}] \\ \equiv \sum_{k = 0}^{N} a_{k} S_{k + j} \\ = \sum_{i = 1}^{M} x_{i}^{j} y_{i} \\ \equiv T_{j}, \end{matrix}

for j=0, 1, . . . , N, where

\sum_{i = 1}^{M} x_{i}^{k + j} \equiv S_{k + j} and \sum_{i = 1}^{M} x_{i}^{j} y_{i} \equiv T_{j},

respectively. There are N+1 equations

\sum_{k = 0}^{N} a_{k} S_{k + j} = T_{j},

\sum_{k = 0}^{N} a_{k} S_{k + j} = 0

only have the trivial solution a

k

=0, for k 0, 1, . . . , N, which may be shown as follows. Multiply the jth homogeneous equation by a

j

and sum over all j, from j=0 to

\begin{matrix} j = N, \\ \sum_{j = 0}^{N} a_{j} \sum_{k = 0}^{N} a_{k} S_{k + j} = \sum_{j = 0}^{N} a_{j} \sum_{k = 0}^{N} a_{k} \sum_{i = 1}^{M} x_{i}^{k} x_{i}^{j} \\ = \sum_{i = 1}^{M} (\sum_{k = 0}^{N} a_{k} x_{i}^{k}) (\sum_{j = 0}^{N} a_{j} x_{i}^{j}) \\ = \sum_{i = 1}^{M} {(P_{N} (x_{i}))}^{2} \\ = 0, \end{matrix}

P_{N} (x) = \sum_{k = 0}^{N} a_{k} x^{k},

that may be fit to the M data points (x

i

,y

i

).

Finding the least-squares polynomial of degree N,

P_{N} (x) = \sum_{k = 0}^{N} a_{k} x^{4},

that may be fit to the M data points (x

i

,y

i

) may not be easy when the degree N of the least-squares polynomial is very large. The N+1 normal equations

\sum_{k = 0}^{N} a_{k} S_{k + j} = T_{j},

S_{k + j} = \sum_{i = 0}^{M} x_{i}^{k + j} \approx M \int_{0}^{i} x^{k + j} ⅆ x = \frac{M}{k + j + 1} .

The resulting determinant for the normal equations is then approximately given by

&LeftBracketingBar; S_{k + j} &RightBracketingBar; \approx &LeftBracketingBar; \frac{M}{k + j + 1} &RightBracketingBar; \approx M^{N + 1} &LeftBracketingBar; \frac{1}{k + j + 1} &RightBracketingBar; = M^{N + 1} H_{N + 1},

for j,k=0, 1, . . . , N, where H

N

, for j,k=0, 1, . . . , N−1, is the Hilbert determinant of order N, which has the value

H_{N} = \frac{{[0! 1! 2! 3! \dots (N - 1)!]}^{3}}{N! (N + 1)! (N + 2)! \dots (2 N - 1)!}

that approaches zero very rapidly. For example,

H_{1} = \frac{{[0!]}^{3}}{1!} = &LeftBracketingBar; 1 &RightBracketingBar; = 1, H_{2} = \frac{{[0! 1!]}^{3}}{2! 3!} = &LeftBracketingBar; \begin{matrix} 1 & \frac{1}{2} \\ \frac{1}{2} & \frac{1}{3} \end{matrix} &RightBracketingBar; = \frac{1}{3} - \frac{1}{4} = \frac{1}{12}, and

H_{3} = \frac{{[0! 1! 2!]}^{3}}{3! 4! 5!} = &LeftBracketingBar; \begin{matrix} 1 & \frac{1}{2} & \frac{1}{3} \\ \frac{1}{2} & \frac{1}{3} & \frac{1}{4} \\ \frac{1}{3} & \frac{1}{4} & \frac{1}{5} \end{matrix} &RightBracketingBar;,

where

&LeftBracketingBar; \begin{matrix} 1 & \frac{1}{2} & \frac{1}{3} \\ \frac{1}{2} & \frac{1}{3} & \frac{1}{4} \\ \frac{1}{3} & \frac{1}{4} & \frac{1}{5} \end{matrix} &RightBracketingBar; = \frac{1}{3} &LeftBracketingBar; \begin{matrix} \frac{1}{2} & \frac{1}{3} \\ \frac{1}{3} & \frac{1}{4} \end{matrix} &RightBracketingBar; - \frac{1}{4} &LeftBracketingBar; \begin{matrix} 1 & \frac{1}{2} \\ \frac{1}{3} & \frac{1}{4} \end{matrix} &RightBracketingBar; + \frac{1}{5} &LeftBracketingBar; \begin{matrix} 1 & \frac{1}{2} \\ \frac{1}{2} & \frac{1}{3} \end{matrix} &RightBracketingBar; = \frac{1}{3} \frac{1}{72} - \frac{1}{4} \frac{1}{12} + \frac{1}{5} \frac{1}{12} = \frac{1}{216} - \frac{1}{240} = \frac{1}{2160} .

P_{2} (x) = a_{0} + a_{1} x + a_{2} x^{2} = \sum_{k = 0}^{2} a_{k} x^{k},

may be fit (in a least-squares sense) to the M=7 data points (−3,4), (−2,2), (−1,3), (0,0), (1,−1), (2,−2), and (3,−5). The residuals ε

i

, for i =1, 2, . . . , 7, are schematically illustrated in FIG.

35

. The three normal equations are

\sum_{k = 0}^{2} a_{k} S_{k + j} = T_{j},

for j=0, 1, 2, where

\sum_{i = 1}^{7} x_{i}^{k + j} \equiv S_{k + j} and \sum_{i = 1}^{7} x_{i}^{j} y_{i} \equiv T_{j},

respectively, for the three parameters a

0

, and a

1

and a

2

. This gives

\sum_{k = 0}^{2} a_{k} S_{k} = T_{0}, \sum_{k = 0}^{2} a_{k} S_{k + 1} = T_{1}, and \sum_{k = 0}^{2} a_{k} S_{k + 2} = T_{2},

where

S_{0} = \sum_{i = 1}^{7} x_{i}^{0} = 7;

S_{1} = \sum_{i = 1}^{7} x_{i} = (- 3 - 2 - 1 + 1 + 2 + 3) = 0;

S_{2} = \sum_{i = 1}^{7} x_{i}^{2} = (9 + 4 + 1 + 1 + 4 + 9) = 28;

S_{3} = \sum_{i = 1}^{7} x_{i}^{3} = (- 27 - 8 - 1 + 0 + 1 + 8 + 27) = 0;

S_{4} = \sum_{i = 1}^{7} x_{i}^{4} = (81 + 16 + 1 + 0 + 1 + 16 + 81) = 196;

T_{0} = \sum_{i = 1}^{7} y_{i} = (4 + 2 + 3 + 0 - 1 - 2 - 5) = 1;

T_{1} = \sum_{i = 1}^{7} x_{i} y_{i} = (12 - 4 - 3 + 0 - 1 - 4 - 15) = - 39; and

T_{2} = \sum_{i = 1}^{7} x_{i}^{2} y_{i} = (36 + 8 + 3 + 0 - 1 - 8 - 45) = - 7,

so that the normal equations becomes

\sum_{k = 0}^{2} a_{k} S_{k} = T_{0} = 1 = 7 a_{0} + 0 a_{1} + 28 a_{2} = 7 a_{0} + 28 a_{2},

\sum_{k = 0}^{2} a_{k} S_{k + 1} = T_{1} = - 39 = 0 a_{0} + 28 a_{1} + 0 a_{2}, and \sum_{k = 0}^{2} a_{k} S_{k + 2} = T_{2} = - 7 = 28 a_{0} + 0 a_{1} + 196 a_{2} = 28 a_{0} + 196 a_{2},

respectively, which imply (upon multiplying the first normal equation by 7 and then subtracting that from the third normal equation) that −14=−21a

0

, that 28a

1

=−39 (from the second normal equation), and (upon multiplying the first normal equation by 4 and then subtracting that from the third normal equation) that −11=84a

2

, giving 3a

0

=2, 28a

1

=−39, and 84a

2

=−11, respectively. In other words, a

0

=2/3, a

1

=−39/28, and a

2

=−11/84, so that, the second degree polynomial (the quadratic) that provides the best fit, in the least-squares sense, is

P_{2} (x) = \frac{2}{3} - \frac{39}{28} x - \frac{11}{84} x^{2} = \frac{1}{84} (56 - 117 x - 11 x^{2}),

as shown in FIG.

35

.

As shown in

FIG. 36

, for example, a second degree polynomial (a quadratic),

P_{2} (x) = a_{0} + a_{1} x + a_{2} x^{2} = \sum_{k = 0}^{2} a_{k} x^{k},

may be fit (in a least-squares sense) to the M=6 data points (0,4), (1,7), (2,10), (3,13), (4,16), and (5,19). The residuals ε

i

, for i=1, 2, . . . , 6, are schematically illustrated in FIG.

36

. The three normal equations are

\sum_{k = 0}^{2} a_{k} S_{k + j} = T_{j},

for j=0, 1, 2, where

\sum_{i = 1}^{6} x_{i}^{k + j} \equiv S_{k + j} and \sum_{i = 1}^{6} x_{i}^{j} y_{i} \equiv T_{j},

respectively, for the three parameters a

0

, and a

1

and a

2

. This gives

\sum_{k = 0}^{2} a_{k} S_{k} = T_{0}, \sum_{k = 0}^{2} a_{k} S_{k + 1} = T_{1}, and \sum_{k = 0}^{2} a_{k} S_{k + 2} = T_{2},

where

S_{0} = \sum_{i = 1}^{6} x_{i}^{0} = 6;

S_{1} = \sum_{i = 1}^{6} x_{i} = (0 + 1 + 2 + 3 + 4 + 5) = 15;

S_{2} = \sum_{i = 1}^{6} x_{i}^{2} = (1 + 4 + 9 + 16 + 25) = 55;

S_{3} = \sum_{i = 1}^{6} x_{i}^{3} = (0 + 1 + 8 + 27 + 64 + 125) = 225;

S_{4} = \sum_{i = 1}^{6} x_{i}^{4} = (0 + 1 + 16 + 81 + 256 + 625) = 979;

T_{0} = \sum_{i = 1}^{6} y_{i} = (4 + 7 + 10 + 13 + 16 + 19) = 69;

T_{1} = \sum_{i = 1}^{6} x_{i} y_{i} = (0 + 7 + 20 + 39 + 64 + 95) = 225; and

T_{2} = \sum_{i = 1}^{6} x_{i}^{2} y_{i} = (0 + 7 + 40 + 117 + 256 + 475) = 895,

so that the normal equations becomes

\sum_{k = 0}^{2} a_{k} S_{k} = T_{0} = 69 = 6 a_{0} + 15 a_{1} + 55 a_{2}, \sum_{k = 0}^{2} a_{k} S_{k + 1} = T_{1} = 225 = 15 a_{0} + 55 a_{1} + 225 a_{2}, and

\sum_{k = 0}^{2} a_{k} S_{k + 2} = T_{2} = 895 = 55 a_{0} + 225 a_{1} + 979 a_{2},

respectively, which imply (upon multiplying the second normal equation by 4 and then subtracting that from the first normal equation multiplied by 10) that −10=−70a

1

−350a

2

, and (upon multiplying the second normal equation by 11 and then subtracting that from the third normal equation multiplied by 3) that 210=70a

1

+66a

2

. However, adding these last two results together shows that 0=a

2

. Furthermore, 3=a

1

. Therefore, using the fact that 3=a

1

and 0=a

2

, the normal equations become

\sum_{k = 0}^{2} a_{k} S_{k} = T_{0} = 69 = 6 a_{0} + 45, \sum_{k = 0}^{2} a_{k} S_{k + 1} = T_{1} = 225 = 15 a_{0} + 165, and

\sum_{k = 0}^{2} a_{k} S_{k + 2} = T_{2} = 895 = 55 a_{0} + 675,

respectively, which all imply that 4=a

0

. In other words, a

0

=4, a

1

=3, and a

2

=0, so that, the second degree polynomial (the quadratic) that provides the best fit, in the least-squares sense, is P

2

(x)=4+3x+0x

2

=4+3x, which is really just a straight line, as shown in FIG.

36

. The residuals ε

i

, for i=1, 2, . . . , 6, all vanish identically in this case, as schematically illustrated in FIG.

36

.

As discussed above, in various other alternative illustrative embodiments, samples may be collected for M data points (x

i

,y

i

), where i=1, 2, . . . , M, and a linearly independent set of N+, functions f

j

(x), for j=0, 1, 2, . . . , N,

y (x) = a_{0} f_{0} (x) + a_{1} f_{1} (x) + \dots + a_{j} f_{j} (x) + \dots + a_{N} f_{N} (x) = \sum_{j = 0}^{N} a_{j} f_{j} (x),

may be fit (in a non-polynomial least-squares sense) to the M data points (x

i

,y

i

). For example, 100 time data points (M=100) may be taken relating the pressure p, the flowrate f, and/or the temperature T, of a process gas to the thickness t of a process layer being deposited in a deposition processing step, resulting in the M data points (p

i

,t

i

) (f

i

,t

i

), and/or (T

i

,t

i

), respectively. For each of these sets of M data points (p

i

,t

i

) (f

i

,t

i

), and/or (T

i

,t

i

), the conventional raw storage requirements would be on the order of 2M. By way of comparison, data compression, according to various illustrative embodiments of the present invention, using a linearly independent set of N+1 functions f

j

(x) as a basis for a representation that is fit (in a non-polynomial least-squares sense) to the M data points (x

i

,y

i

), only requires N+1 parameters, a

k

, for k=0, 1, . . . , N, which would reduce the storage requirements, for each of these sets of M data points (p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

), to the storage requirements for the respective coefficients a

k

, for k=0, 1, . . . , N, a storage savings of a factor of 2M/(N+1), equivalent to a data compression ratio of 2M:(N+1), for each of these sets of M data points p

i

,t

i

), (f

i

,t

i

), and/or (T

i

,t

i

). In one illustrative embodiment, the number N+1 of the linearly independent set of basis functions f

j

(x) is at least 10 times smaller than M. For semiconductor manufacturing processes, data sufficient for fault detection purposes ranges from about 100 to about 10000 samples for the value of M. Requirements for polynomial order N are in the range of 1 to about 10.

The function F(a

0

,a

1

, . . . , a

N

) may be minimized as follows. The function F(a

0

,a

1

, . . . , a

N

) is given by

F (a_{0}, a_{1}, \dots, a_{N}) = \sum_{i = 1}^{M} ε_{i}^{2} = \sum_{i = 1}^{M} {[y (x_{i}) - y_{i}]}^{2}

and setting the partial derivatives of F(a

0

,a

1

, . . . , a

N

) with respect to a

j

, for j=0, 1, . . . , N, equal to zero gives

\frac{\partial F (a_{0}, a_{1}, \dots, a_{N})}{\partial a_{j}} = 2 \sum_{i = 1}^{M} [y (x_{i}) - y_{i}] x_{i}^{j} = 2 \sum_{i = 1}^{M} [\sum_{k = 0}^{N} a_{k} f_{k} (x_{i}) - y_{i}] f_{j} (x_{i}) = 0,

for j=0, 1, . . . , N, since f

j

(x

i

) is the coefficient of a

j

in the representation

y (x_{i}) = \sum_{k = 0}^{N} a_{k} f_{k} (x_{i}) .

Simplifying gives

\sum_{i = 1}^{M} [\sum_{k = 0}^{N} a_{k} f_{k} (x_{i})] f_{j} (x_{i}) = \sum_{k = 0}^{N} a_{k} [\sum_{i = 1}^{M} f_{k} (x_{i}) f_{j} (x_{i})] \equiv \sum_{k = 0}^{N} a_{k} S_{k, j} = \sum_{i = 1}^{M} f_{j} (x_{i}) y_{i} \equiv T_{j},

for j=0, 1, . . . , N, where

\sum_{i = 1}^{M} f_{k} (x_{i}) f_{j} (x_{i}) \equiv S_{k, j} and \sum_{i = 1}^{M} f_{j} (x_{i}) y_{i} \equiv T_{j},

respectively. There are N+1 equations

\sum_{k = 0}^{N} a_{k} S_{k, j} = T_{j},

\sum_{k = 0}^{N} a_{k} S_{k, j} = 0

only have the trivial solution a

k

=0, for k=0, 1, . . . , N, which may be shown as follows. Multiply the jth homogeneous equation by a

j

and sum over all j,

\sum_{j = 0}^{N} a_{j} \sum_{k = 0}^{N} a_{k} S_{k, j} = \sum_{j = 0}^{N} a_{j} \sum_{k = 0}^{N} a_{k} \sum_{i = 1}^{M} f_{k} (x_{i}) f_{j} (x_{i}) = \sum_{i = 1}^{M} (\sum_{k = 0}^{N} a_{k} f_{k} (x_{i})) (\sum_{j = 0}^{N} a_{j} f_{j} (x_{i})),

but

\sum_{i = 1}^{M} (\sum_{k = 0}^{N} a_{k} f_{k} (x_{i})) (\sum_{j = 0}^{N} a_{j} f_{j} (x_{i})) = \sum_{i = 1}^{M} {(y (x_{i}))}^{2} = 0,

which would imply that y(x

i

)≡0, and, hence, that a

k

=0, for k=0, 1, . . . , N, the trivial solution. Therefore, the determinant of the normal equations is not zero, and the normal equations may be solved for the N+1 parameters a

k

, for k=0, 1, . . . , N, the coefficients of the non-polynomial least-squares representation

y (x) = \sum_{j = 0}^{N} a_{j} f_{j} (x)

that may be fit to the M data points (x

i

,y

i

), using the linearly independent set of N+1 functions f

j

(x) as the basis for the non-polynomial least-squares representation

y (x) = \sum_{j = 0}^{N} a_{j} f_{j} (x) .

F (a_{0}, a_{1}, \dots, a_{N}) = \sum_{i = 1}^{M} w_{i} ε_{i}^{2} = \sum_{i = 1}^{M} {w_{i} [y (x_{i}) - y_{i}]}^{2}

and setting the partial derivatives of F(a

0

,a

1

, . . . , a

N

) with respect to a

j

, for j=0, 1, . . . , N, equal to zero gives

\frac{\partial F (a_{0}, a_{1}, \dots, a_{N})}{\partial a_{j}} = 2 \sum_{i = 1}^{M} w_{i} [y (x_{i}) - y_{i}] x_{i}^{j} = 2 \sum_{i = 1}^{M} w_{i} [\sum_{k = 0}^{N} a_{k} f_{k} (x_{i}) - y_{i}] f_{j} (x_{i}) = 0,

for j=0, 1, . . . , N, since f

j

(x

i

) is the coefficient of a

j

in the representation

y (x_{i}) = \sum_{k = 0}^{N} a_{k} f_{k} (x_{i}) .

Simplifying gives

\sum_{i = 1}^{M} w_{i} [\sum_{k = 0}^{N} a_{k} f_{k} (x_{i})] f_{j} (x_{i}) = \sum_{k = 0}^{N} a_{k} [\sum_{i = 1}^{M} w_{i} f_{k} (x_{i}) f_{j} (x_{i})] = \sum_{i = 1}^{M} w_{i} f_{j} (x_{i}) y_{i},

\sum_{k = 0}^{N} a_{k} [\sum_{i = 1}^{M} w_{i} f_{k} (x_{i}) f_{j} (x_{i})] \equiv \sum_{k = 0}^{N} a_{k} S_{k, j} = \sum_{i = 1}^{M} w_{i} f_{j} (x_{i}) y_{i} \equiv T_{j}

for j = 0, 1, \dots, N,

where

\sum_{i = 1}^{M} w_{i} f_{k} (x_{i}) f_{j} (x_{i}) \equiv S_{k, j} and \sum_{i = 1}^{M} w_{i} f_{j} (x_{i}) y_{i} \equiv T_{j},

respectively. There are N+1 equations

\sum_{k = 0}^{N} a_{k} S_{k, j} = T_{j}, for j = 0, 1, \dots, N,

also known as the normal equations, including the non-negative weighting factors w

i

, for the N+1 unknown parameters a

k

, for k=0, 1, . . . , N, readily yielding a solution, provided that the determinant of the normal equations is not zero. This may be demonstrated by showing that the homogeneous equations

\sum_{k = 0}^{N} a_{k} S_{k, j} = 0

only have the trivial solution a

k

=0, for k=0, 1, . . . , N, which may be shown as follows. Multiply the jth homogeneous equation by a

j

and sum over all j,

\sum_{j = 0}^{N} a_{j} \sum_{k = 0}^{N} a_{k} S_{k, j} = \sum_{j = 0}^{N} a_{j} \sum_{k = 0}^{N} a_{k} \sum_{i = 1}^{M} w_{i} f_{k} (x_{i}) f_{j} (x_{i}) = \sum_{i = 1}^{M} w_{i} (\sum_{k = 0}^{N} a_{k} f_{k} (x_{i})) (\sum_{j = 0}^{N} a_{j} f_{j} (x_{i})),

but

\sum_{i = 1}^{M} w_{i} (\sum_{k = 0}^{N} a_{k} f_{k} (x_{i})) (\sum_{j = 0}^{N} a_{j} f_{j} (x_{i})) = \sum_{i = 1}^{M} {w_{i} (y (x_{i}))}^{2} = 0,

which would imply that y(x

i

)≡0, and, hence, that a

k

=0, for k=0, 1, . . . , N, the trivial solution. Therefore, the determinant of the normal equations is not zero, and the normal equations, including the non-negative weighting factors w

i

, may be solved for the N+1 parameters a

k

, for k=0, 1, . . . , N, the coefficients of the non-polynomial least-squares representation

y (x) = \sum_{j = 0}^{N} a_{j} f_{j} (x)

that may be fit to the M data points (x

i

,y

i

), using the linearly independent set of N+1 functions f

j

(x) as the basis for the non-polynomial least-squares representation

y (x) = \sum_{j = 0}^{N} a_{j} f_{j} (x),

and including the non-negative weighting factors w

i

.

The particular embodiments disclosed above are illustrative only, as the invention may be modified and practiced in different but equivalent manners apparent to those skilled in the art having the benefit of the teachings herein. Furthermore, no limitations are intended to the details of construction or design herein shown, other than as described in the claims below. It is therefore evident that the particular embodiments disclosed above may be altered or modified and all such variations are considered within the scope and spirit of the invention. Accordingly, the protection sought herein is as set forth in the claims below.

Number	Name	Date	Kind
5864773	Barna et al.	Jan 1999	A
5877954	Klimasauskas et al.	Mar 1999	A
6356857	Qin et al.	Mar 2002	B1
6597959	Backa et al.	Jul 2003	B1

Method of using control models for data compression

Information

Patent Number

Date Filed

Date Issued

Inventors

Original Assignees

Examiners

Agents

CPC

US Classifications

Field of Search

US

International Classifications

Term Extension

Abstract

Description

Claims

US Referenced Citations (4)

Non-Patent Literature Citations (1)