- Email: [email protected]
On obtaining surface time profiles from pulsed molecular beam scattering time-of-flight measurements
surface s c i e n c e ELSEVIER
Surface Science 365 (1996) 525-534
On obtaining surface time profiles from pulsed molecular beam scattering time-of-flight measurements. P . N . Brier, J.N. Fletcher, P.A. G o r r y * Department of Chemistry, Universityof Manchester, Manchester M13 9PL, UK
Received 27 January 1996;accepted for publication 28 March 1996
Abstract
The use of very short gas pulses to obtain the surface time response of a reactive system is a powerful technique. However, reconstructing the time profile from real data can be a non-trivial problem. In particular we consider the common situation where the residence time profile falls sharply at t = 0 and the flight time to the detector is of a similar magnitude to the residence time. In the presence of noise, deconvolution of the (Maxwell-Boltzmann) flight profile is not a simple matter. Three methods; model-based forward convolution, Optimum (Wiener) filtering/deconvolution and model-free forward convolution, are explored and an overall strategy is suggested. Keywords: Chlorine; Gallium arsenide; Models of surface-chemical reactions; Molecule-solid reactions; Molecule-solid scattering
and diffraction - inelastic; Surface chemical reaction
1. Introduction
Molecular beam scattering under U H V conditions is a powerful tool in the study of chemical processes on surfaces [ 1]. In addition to providing mass resolved angular and translational energy distributions the use of modulated or pulsed beams can provide real-time information on the residence time profiles of species on surfaces and hence provide direct information on the surface kinetics involved. The use of chopped molecular beams in Modulated Beam Relaxation Spectrometry (MBRS) has been well documented [ 2 - 5 ] . In its simplest form the scattered amplitude and phase
* Corresponding author. Fax: + 44 161 2754598; e-mail: [email protected]
delay of the products are measured with phase sensitive detection. Variation of the modulation frequency then allows the determination of the reaction probability and average residence time. A more modern approach is to measure the entire scattered profile with a transient digitiser and to utilise several frequency components from its fourier transform to determine the so-called product vector [ 1 - 5 ] . Polar plots of the product vector (impulse response transfer function) can then be compared with model functions from different kinetic schemes. The strengths and weaknesses of the MBRS technique have been discussed in detail by Olander [6]. In MBRS experiments the measured profile is a convolution of the impulse response of the system (to a ~-function of gas) and the chopping function profile. Since the important quantity is the impulse response, this must then be obtained by deconvolu-
0039-6028/96/$15.00 Copyright© 1996Elsevier ScienceB.V. All rights reserved PH S0039-6028 (96) 00718-2
526
P.N. Brier et al./Surface Science 365 (1996) 525-534
tion. As discussed below, deconvolution is not always a simple task, especially if high frequency components are required in the presence of noise. Another approach is the application of well separated gas pulses [5,7]. These have a flatter frequency response than chopped beams which, in particular, aids in the measurement of fast reactions. The principle advantage of gas pulses however, is that for very short pulses the measured time response is the impulse response itself and no deconvolution of a chopping profile is required [8]. Even with no instrumental broadening, the impulse response is itself a convolution of two time processes; the surface time profile arising from the reaction kinetics on the surface and the flight time from the surface to the detector. In this paper we examine several techniques for obtaining the surface time profile from the impulse response of the system arising from very short gas pulses.
2. Experimental The experimental data used here arise from our recent scattering studies of C12 on GaAs (100) [9,10]. An excellent review of molecular beam reactive scattering studies of surface chemistry on semi-conductors has recently been provided by Yu and DeLouise [8]. A fuller description of the experimental set-up is given elsewhere [ 10]. Briefly, short supersonic gas pulses (FWHM~35/~s, Mach N o e l 2 ) of C12/He are incident at 45 ° onto a heated GaAs (100) surface under UHV conditions. The resulting GaC1 product flux, measured along the surface normal, is then detected by an electron bombardment ioniser/quadrupole mass spectrometer/pulse counting system situated 500 mm away from the surface. Time-of-fright spectra (TOF) are recorded by a 100 MHz Stanford Research Systems SR430 multichannel scalar. Data are collected in a beamon/beam-off mode to provide background subtraction. The T O F spectra are corrected for the ion transit time in the detector and for the flight time from the pulsed beam to the surface.
3. Mathematical considerations Providing the width of the initial gas pulse approximates to a 3-function, the measured TOF distribution, 7-(0, is a convolution of a time-offlight distribution, Tf(t), and a surface time profile,
T,(t). T(t)= i T~(ts)Tt(t-ts)dts"
(l)
o
If the residence time on the surface is short compared to the flight time (q <
exp{ Tf(t) =
t4
,
(2)
where L = flight distance and c~= 2x/~kT/m.
3.1. Forward convolution The most robust method of finding T,(t) is to adopt a model form and then perform the convolu-
P.N. Brier et al./Surface Science 365 (1996) 525-534
tion integral Eq. (1) numerically. This has the advantage that the evaluation of the integral can be performed with high accuracy but the disadvantage that a number of model functions must be tried to find the best description for T~(t). Fig. 1 shows the experimental GaC1 T O F from a surface at 450°C. Also shown is the fit obtained from a two exponential form for T~(t). T~(t) = exp ( - t/z 1)+ B exp ( - t/z2).
527
recourse to model functions for T~(t). We now examine this approach. In general one has a measured distribution g(t) which is a convolution of the true signal, f(t) and a response function h(t) g(t)= i f(z)h(t-z) dr.
(4)
0
(3)
In Fourier space the convolution theorem gives
The parameters B, zl and r2 were found using a non-linear least squares procedure due to Marquardt [ 11,17]. The T~(t) is characterised by a fast and slow step for production of GaC1. It is likely that this reflects the desorption of GaC1 from two different sites with different AaaH values [ 10]. For the purposes of this paper we shall treat this as the "true" form for T~(t) when comparing methods.
G(f) = H(f). F(f) hence
G(f)
F ( f ) = H(f~)"
(5)
Unfortunately if H(f) contains zeros (or very small values) then non-zero values of G(f) may be greatly magnified, destroying any hope of recovering F(f) and hencef(t). This is often the case with real experimental data where most realistic response functions H ( f ) will contain some very small values, and random noise with its
3.2. The deconvolution problem Ideally we would like a procedure for deconvoluting Eq.(1) directly so that we do not need
~ " 1.0E -'I
>, if/ t-
a
0.5-
..Q
+
E Z
"
"+ +
.+
0.0-
1
0
5
10
15
20
time/ms Fig. 1. The T O F spectrum of GaCI from C12q- OaAs (100) at a surface temperature of 450°C. The channel width is 81.9 #s and the initial gas pulse is ~ 35 #s. Also shown is the forward convolution backfit (with two components) of Eq. (1) and the two-exponential model of Eq. (3) (zl =0.507 ms, z2=6.148 ms, B=0.1537).
528
P.N. Brier et al./Surface Science 365 (1996) 525-534
elements in # ( f ) with f-- f e u t . We chose fcut by examining the frequency components of G(f) and H ( f ) - care must be taken to ensure that truncation of the high frequency components does not distort the F(f) obtained. Fig. 2 shows the frequency components of the data in Fig. 1, the Maxwell-Boltzmann function, and the results of applying Eq. (6) with a cutoff filter with fc~t = 1600 Hz. We can see that the overall behaviour agrees well with the forward convolution fit but that considerable oscillations occur. Most significantly the behaviour close to t = 0 is very different. The simple filter deconvolution suggests a peak in GaC1 production at early time rather than decay from a maximum at t = 0. Such behaviour suggests a model with two sequen-
unrestricted power spectrum, ensures that there are non-zero elements in G(f) at all frequencies. This is indeed what happens with the Maxwell-Boltzmann expression where H ( f ) has very small values at high frequencies. 3.3. Filtering We can alleviate the problem by applying a filter 4)(f) to remove high-frequency components in G(f) so that they are not amplified during the deconvolution. F(f) =
G(f).~(f) = G(f). W(f). n(f)
(6)
The simplest filter is a simple cut-off where all
Frequency/ Hz 0
_
2'0t
_
1000
2000
3000
4000
5000
I
I
I
I
I
6000
Inverse
~Q. 1.5
Maxwell-Boltzmann x 10-6 E
o 1.00 o LL
TOF Data
~
j" J
~:=~
05-
S/'t
"
lI-Boltzmann
.,
J
7
0.01.0"o
E 0.5o t-
v
v m
"v VVV~v"v~-'~-V
0.0-
o
~"
^V^V"v^v^v' " ^ ^ v"v , ,
i'o time I ms
--
^^ v V
;5
Fig. 2. Upper panel: The frequency components of the TOE data of Fig. 1 and the Maxwell-Boltzmann function for T = 450°C. Also shown is the inverse of the Maxwdl-Boltzmann components, clearly showing the large magnification factor at high frequencies.
P.N. Brier et al./Surface Science365 (1996) 525 534 (in the least-squares sense) is given by [15,16]:
~ _
529
n(f)*
lD(f)[ 2
IN(012
W(f) =
"..~.~_ Frequency
Fig. 3. Typical power spectrum. ID(f)12is the power spectrum of the total measuredsignal, IN(f)12is the noisepower spectrum, and IG(f)12is the required signal power spectrum. The noise spectrum shownis for whitenoise with a flat frequencyresponse.
tial first order steps, as has been observed by Eldridge and Yu for oxygen on Si(100) [ 12], rather than two simultaneous ones. Unfortunately a careful study with synthetic data in the absence of noise reveals that the truncation of high frequency components makes it impossible to mimic the sharp fall from t = 0 of the forward convolution model. The earlier success of Eldridge and Yu [12] in deconvoluting the observed T O F was in part due to a fortuitous form of the surface time profile. The necessity to include high-frequency components in G(f) greatly hinders the use of the deconvolution method for processes involving a sharp fall from t = 0.
3.4. Optimum (Wiener)filtering Alternative iterative deconvolution methods are widely used. However, it can be shown that iterative deconvolution techniques, such as the Morrison [13] and van Cittert algorithms [14], can be expressed as a simple filter operation in Fourier space [ 15]. The question then arises as to what is the optimum filter to use in Eq. (6)? If we explicitly include noise, n(t), in our description we have a measured signal d(t) given by:
iN(f)l 2 . IH(f)12 + iV(f)12
(8)
where * represents the complex conjugate. The optimum filter thus involves not only the frequency spectrum of the response function (Maxwell-Boltzmann) but also the power spectra of the noise and F(f). Obviously, since F(f) is what we are seeking, we cannot make use of Eq. (8) directly - but we can construct some approximations to Eq. (8), based on various assumptions about N ( f ) and F(f), which will be better than a simple cut-off filter.
3.5. Approximate optimum filters 3.5.1. Constant noise power Fig. 3 shows a schematic view of a typical power spectrum. We can write Eq. (8) as
w(D=
1
[~(f)l 2
H(f) IG(f)l 2 + IN(f)l e
1 IO(f)12--lN(f)l 2 - H(f~) ID(f)] 2
(9)
In the absence of a separate determination of
N ( f ) the best assumption is that we have white noise with a flat noise spectrum. We can then use the high frequency end of the spectrum (where D ( f ) = N ( f ) ) to determine an average value for the noise power. N2av = (IN(f)] 2> for high frequencies
1 W ( f ) = H(f~)
ID(f)12-N2v
ID(f)l 2
(10)
Unfortunately this filter does not work well if Nay is small since at this limit it approaches the simple deconvolution form of Eq. (5) and suffers similar problems of noise amplification.
3.5.2. Constant noise:signal power ratio Returning to Eq. (8) we can replace the noise:signal term by a constant at all frequencies.
d(t) = g(t) + n(t) { hence D ( f ) = G ( f ) + N ( f ) }. (7) It can be shown that the optimal form for W(f)
_ ~ LN(f)[2 S n(t) 2 dt IF(f)[ 2 - Sf(t) 2 d~
(11)
530
P.N. Brier et al./Surface Science 365 (1996) 525-534
and for small signal components (or large noise)
H(f)* W(f) = iH(f)l z + ¢ .
(12)
The average noise:signal power ratio can often be estimated from the data or we can simply choose a small value for ~b, 10-3-10 -6 say. The filter in Eq. (12) is robust and fairly insensitive to the value of ¢ chosen. It utilises the most important (large value) frequency components in H(f) directly while automatically limiting the magnification from small values - since the filter tends to zero when H(f) is small i.e.
H(f) >>¢, W(f)--+1/H(f),
N]v
iFe(f)12 >>In(f)l 2, W(f)~O Eq.(13) has the advantages that it explicitly makes use of an approximate form for F(f) in constructing the filter and that the only subjective decision is the frequency above which to estimate the average noise power. In reality the difference between this filter and one for constant ~b is often very small. The result of applying Eq. (13) is also shown in Fig. 4. All filters suffer the similar problems at t = 0.
H(f) < ¢, W(f)--+H(f)/~ ~ O.
3.6. Model-freeforward convolution
Fig. 4 shows the form of Eq. (12) for two values of ¢ and the result of applying Eq. (12) to the data of Fig. 1 with ¢ = 5 x 10 -3. The result is an improvement on a simple cut-off filter but still suffers from the difficulty that it cannot reproduce the behaviour at t = 0 because of the truncation of the high frequency components.
A final approach is that of a model-free forward convolution. For n data points {di} we represent f(t) by n points {f/} that can be freely adjusted and use the forward convolution procedure of Eq. (1) to calculate the resulting T O F spectrum {gl) for comparison to the measured TOF data. The values of the {f~} points are adjusted to minimise the least-squares fit to the data. In matrix notation
3.5.3. Estimated F(f) The true optimum filter Eq. (8) contains F(f). Although we do not know this a priori we can use the results from a filter with constant ~b to obtain an approximation, F¢(f), to F(f) and hence construct an improved version of Eq. (12). In the absence of an independent determination of N ( f ) we again use an average value constructed from the high frequency part of D(f).
H(f)* F4,(f)=D(f ) iH(f)l z +gP and (13)
(14)
hki f .
and we minimise ~(2= ~
(gk--dk)2
k =1
(~k2
k=l
O'k2 i=1
= 2
~
hki(gk-- dk)
k =1
tTk2
2
~ fi
N Zv IH(f)l 2 + IFe(f)lz
~
i=1
~2
H(f)* W(f)=
gk =
,15,
using an iterative estimation of {f/} based on the steepest descent method for minimising Z2 [17]. ~X2
This filter is also robust at controlling magnification effects since for large signal components (or small noise) N2v - -
IF,(f)l 2
1
2,
W(f)~
- -
n(f)
f'~+1=f'r + Af'~ - -
0f,
(16)
where 2 is a scaling parameter used to control the step-size and aid convergence. This method requires an initial set of {f} values and, like all iterative methods, a good starting set greatly aids
531
P.N. Brier et aL/Surface Science365 (1996) 525 534
Frequency / Hz
,,-,, "10
1000
2000
3000
4000
5000
I
I
I
I
I
6000
0 = 5x10 -3
1.0-
09
E O
Z v
0.5-
n
0.0 1.0-
"o f,O
t~
E k--
0.5-
o cv
v ~ n
v
0.0-
T
0
5
I'0
I'5
--
V
2'0
time / ms Fig. 4. Upper panel: The frequency components of the filter W(f) of Eq. (12) for two different values of ~. They have been normalised to their peak value to facilitate comparison.
the convergence. We can use the results of either the best-model forward convolution procedure or the Wiener filter deconvolution to indicate the likely shape of f ( t ) . To test the convergence a number of starting sets {f/} are used. Fig. 5 shows the results of the best two-exponential model Eq.(3), the cut-off filter Eq.(6), the optimum Wiener filter Eq. (13), and the model-free forward convolution Eq. (16). In the latter case the starting set {f,.} was a single exponential curve with a time constant intermediate between the two values found for Eq.(3). The results of Eq.(16) have clearly converged to excellent agreement with the two exponential fit and show no evidence of a peak at t > 0 as suggested by the Wiener filter. We have sucessfully applied the model-free for-
ward convolution approach to the deconvolution of GaC1 T O F spectra over a range of surface temperatures and timescales, These all follow the same pattern of a two-exponential decrease. In order to test the algorithm further we have applied it to a range of synthetic data based on different kinetic models. The upper panel of Fig. 6 shows synthetic data generated by assuming GaCI is produced by two sequential first-order processes (as in the case of O2 on Si (100) [-12]) with the same time constants as before. 1 T~(t)
=
-
-
"El --'~2
[exp ( - t/z 1) - exp ( - t/r2)].
(17)
The model function is convoluted with the
532
P.N. Brier et al./Surface Science 365 (1996) 525-534
1.0"
I
t~ ~0.5 t-
v
v
I 0.0
o
)
"'s . . . . .
;
4
time / ms Fig. 5. A comparison of the behaviour off(t) at early time for all the methods. The curves are: the best forward convolution model of Eq. (3), the model-free forward convolution of Eq. ( 16 ) ( both solid ), the cut-off filter in Eq. (6) (dot-dash), and the Optimum Wiener filter (dashed). Only the forward convolution methods have the correct behaviour at t=0.
10
+
g
1.0-
0.5-
0.0 0
5
10
15
20
time/ms Fig. 6. Upper panel: synthetic data for two sequential first-order processes (z1=0.507 ms, %=6.148 ms), convoluted with a Maxwell-Boltzmann distribution at 450°C and Gaussian noise added. Lower panel: the model function of Eq. (17) (thick line) and the result of the forward convolution procedure (thin line).
P.N. Brier et al./Surface Science 365 (1996) 525-534 Maxwell-Boltzmann distribution at 450°C and Gaussian noise added. The lower panel shows the result of the forward convolution procedure using exactly the same starting guess used for Fig. 5. This initial guess, a single exponential decay, now has entirely the wrong shape since it has its maximum value at t = 0 instead of the required zero. The forward convolution model still converges to the correct shape. The success of the model-free forward convolution is at first surprising since n freely adjustable points would seem to be a recipe for an unstable solution. However, the behaviour of Eq.(16) is found to be very good providing the 2 value is kept small to inhibit large steps. The majority of deconvolution problems involve correction for instrument response where the convolution function is much narrower than the measured spectrum. In such a case a value at f ( t ) contributes only to a small number of data points close to t and noise can significantly affect the reconstruction o f f ( t ) locally. The opposite is the case here where the Maxwell-Boltzmann convolution function is very broad and each point in f ( t ) c o n t r i b u t e s to a very large range of experimental data points. This has the effect of greatly reducing the effect that noise has on determining t h e f ( t ) value and the n-point model-free set of points {.[i} smoothly converge to a stable solution.
4. Conclusions The task of obtaining the underlying shape of surface time profiles from the measured T O F spectra from pulsed molecular beam scattering is not a simple one. The three methods used, modelbased forward convolution, Wiener filter deconvolution, and model-free forward convolution all have their strengths and weaknesses. These can be summarised by: 4.1. Model-based forward convolution 1. Robust. 2. Accurate. 3. Requires many trial models.
4. Requires iterative optimisation parameters.
533 of model
4.2. Wiener filter deconvolution 1. 2. 3. 4.
Robust. Quick and simple. Insensitive to exact form of filter used. Poor for sharp features due to loss of high frequency components.
4.3. Model-free forward convolution 1. Can be unstable need to control step size. 2. Requires reasonable starting set of values. 3. Several starting sets required to test convergence. 4. G o o d noise rejection in final f(t). 5. Can reproduce sharp features. An approach that utilises all three techniques provides the best method for obtaining the surface time profiles. Initial use of the Wiener filter provides a rapid and robust method of obtaining a good indication of the overall behaviour. This is used to select a suitable set of starting values for the model-free forward convolution procedure. Finally this can be used to suggest suitable analytical models, usually with a physical interpretation behind them, for final fitting by forward convolution.
Acknowledgements This work was supported by the Engineering and Physical Science Research Council and the DRA, Electronics Division at Malvern.
References [ 1] M.Asscherand G.A. Somorjai, in: Atomic and Molecular Beam Methods, Vol.2, ed. G. Scoles (Oxford University Press, Oxford, 1992) p. 488. [2] J.A. Schwartz and R.J. Madix, Surf. Sci. 46 41974) 317. [3] M.P. D'Evelyn and R.J. Madix, Surf. Sci. Rep. 3 (1984) 413. 1-4] C.T. Foxon, M.R. Boudry and B.A. Joyce, Surf. Sci. 44 (1974) 69.
534
P.N. Brier et al./Surface Science 365 (1996) 525-534
[5] H.H. Sawin and R.P. Merrill, J. Vac. Sci. Technol. 19 (1981) 40. I-6] D.R. Olander, J. Coll. Interface Sci. 58 (1977) 169. [7] H.C. Chang and W.H. Weinberg, Appl. Surf. Sci. 3 (1979) 168. [8] M.L. Yu and L.A. DeLouise, Surf. Sci. Rep. 19 (1994) 285. I-9] P. Bond, P.N. Brier, J. Fletcher, P.A. Gorry and M.E. Pemble, Chem. Phys. Lett. 208 (1993) 269 and 211 (1993) 649. [10] P. Bond, P.N. Brier, J. Fletcher, P.A. Gorry and M.E. Pemble, in preparation [11] D.W. Marquardt, J. Soc. Ind. Appl. Math. 2 (1963) 431.
[12] B.N. Eldridge and M.L. Yu, Rev. Sci. Instrum. 58 (1987) 1014. [13] J.D. Morrison, J. Chem. Phys. 39 (1963) 200. [14] P.H. van Cittert, Z. Physik 69 (1931) 298. [15] C.C. Watson and D.V. Ellis, Nucl. Geophys. 5 (1991) 451. [16] N. Wiener, Extrapolation, Interpretation and Smoothing of Stationary Time Series (Technology Press of MIT and Wiley, New York, 1941) [ 17] P.R. Bevington and D.K. Robinson, Data Reduction and Error Analysis for the Physical Sciences, 2nd ed. (McGraw Hill, New York, 1992)
Surface Science 365 (1996) 525-534
On obtaining surface time profiles from pulsed molecular beam scattering time-of-flight measurements. P . N . Brier, J.N. Fletcher, P.A. G o r r y * Department of Chemistry, Universityof Manchester, Manchester M13 9PL, UK
Received 27 January 1996;accepted for publication 28 March 1996
Abstract
The use of very short gas pulses to obtain the surface time response of a reactive system is a powerful technique. However, reconstructing the time profile from real data can be a non-trivial problem. In particular we consider the common situation where the residence time profile falls sharply at t = 0 and the flight time to the detector is of a similar magnitude to the residence time. In the presence of noise, deconvolution of the (Maxwell-Boltzmann) flight profile is not a simple matter. Three methods; model-based forward convolution, Optimum (Wiener) filtering/deconvolution and model-free forward convolution, are explored and an overall strategy is suggested. Keywords: Chlorine; Gallium arsenide; Models of surface-chemical reactions; Molecule-solid reactions; Molecule-solid scattering
and diffraction - inelastic; Surface chemical reaction
1. Introduction
Molecular beam scattering under U H V conditions is a powerful tool in the study of chemical processes on surfaces [ 1]. In addition to providing mass resolved angular and translational energy distributions the use of modulated or pulsed beams can provide real-time information on the residence time profiles of species on surfaces and hence provide direct information on the surface kinetics involved. The use of chopped molecular beams in Modulated Beam Relaxation Spectrometry (MBRS) has been well documented [ 2 - 5 ] . In its simplest form the scattered amplitude and phase
* Corresponding author. Fax: + 44 161 2754598; e-mail: [email protected]
delay of the products are measured with phase sensitive detection. Variation of the modulation frequency then allows the determination of the reaction probability and average residence time. A more modern approach is to measure the entire scattered profile with a transient digitiser and to utilise several frequency components from its fourier transform to determine the so-called product vector [ 1 - 5 ] . Polar plots of the product vector (impulse response transfer function) can then be compared with model functions from different kinetic schemes. The strengths and weaknesses of the MBRS technique have been discussed in detail by Olander [6]. In MBRS experiments the measured profile is a convolution of the impulse response of the system (to a ~-function of gas) and the chopping function profile. Since the important quantity is the impulse response, this must then be obtained by deconvolu-
0039-6028/96/$15.00 Copyright© 1996Elsevier ScienceB.V. All rights reserved PH S0039-6028 (96) 00718-2
526
P.N. Brier et al./Surface Science 365 (1996) 525-534
tion. As discussed below, deconvolution is not always a simple task, especially if high frequency components are required in the presence of noise. Another approach is the application of well separated gas pulses [5,7]. These have a flatter frequency response than chopped beams which, in particular, aids in the measurement of fast reactions. The principle advantage of gas pulses however, is that for very short pulses the measured time response is the impulse response itself and no deconvolution of a chopping profile is required [8]. Even with no instrumental broadening, the impulse response is itself a convolution of two time processes; the surface time profile arising from the reaction kinetics on the surface and the flight time from the surface to the detector. In this paper we examine several techniques for obtaining the surface time profile from the impulse response of the system arising from very short gas pulses.
2. Experimental The experimental data used here arise from our recent scattering studies of C12 on GaAs (100) [9,10]. An excellent review of molecular beam reactive scattering studies of surface chemistry on semi-conductors has recently been provided by Yu and DeLouise [8]. A fuller description of the experimental set-up is given elsewhere [ 10]. Briefly, short supersonic gas pulses (FWHM~35/~s, Mach N o e l 2 ) of C12/He are incident at 45 ° onto a heated GaAs (100) surface under UHV conditions. The resulting GaC1 product flux, measured along the surface normal, is then detected by an electron bombardment ioniser/quadrupole mass spectrometer/pulse counting system situated 500 mm away from the surface. Time-of-fright spectra (TOF) are recorded by a 100 MHz Stanford Research Systems SR430 multichannel scalar. Data are collected in a beamon/beam-off mode to provide background subtraction. The T O F spectra are corrected for the ion transit time in the detector and for the flight time from the pulsed beam to the surface.
3. Mathematical considerations Providing the width of the initial gas pulse approximates to a 3-function, the measured TOF distribution, 7-(0, is a convolution of a time-offlight distribution, Tf(t), and a surface time profile,
T,(t). T(t)= i T~(ts)Tt(t-ts)dts"
(l)
o
If the residence time on the surface is short compared to the flight time (q <
exp{ Tf(t) =
t4
,
(2)
where L = flight distance and c~= 2x/~kT/m.
3.1. Forward convolution The most robust method of finding T,(t) is to adopt a model form and then perform the convolu-
P.N. Brier et al./Surface Science 365 (1996) 525-534
tion integral Eq. (1) numerically. This has the advantage that the evaluation of the integral can be performed with high accuracy but the disadvantage that a number of model functions must be tried to find the best description for T~(t). Fig. 1 shows the experimental GaC1 T O F from a surface at 450°C. Also shown is the fit obtained from a two exponential form for T~(t). T~(t) = exp ( - t/z 1)+ B exp ( - t/z2).
527
recourse to model functions for T~(t). We now examine this approach. In general one has a measured distribution g(t) which is a convolution of the true signal, f(t) and a response function h(t) g(t)= i f(z)h(t-z) dr.
(4)
0
(3)
In Fourier space the convolution theorem gives
The parameters B, zl and r2 were found using a non-linear least squares procedure due to Marquardt [ 11,17]. The T~(t) is characterised by a fast and slow step for production of GaC1. It is likely that this reflects the desorption of GaC1 from two different sites with different AaaH values [ 10]. For the purposes of this paper we shall treat this as the "true" form for T~(t) when comparing methods.
G(f) = H(f). F(f) hence
G(f)
F ( f ) = H(f~)"
(5)
Unfortunately if H(f) contains zeros (or very small values) then non-zero values of G(f) may be greatly magnified, destroying any hope of recovering F(f) and hencef(t). This is often the case with real experimental data where most realistic response functions H ( f ) will contain some very small values, and random noise with its
3.2. The deconvolution problem Ideally we would like a procedure for deconvoluting Eq.(1) directly so that we do not need
~ " 1.0E -'I
>, if/ t-
a
0.5-
..Q
+
E Z
"
"+ +
.+
0.0-
1
0
5
10
15
20
time/ms Fig. 1. The T O F spectrum of GaCI from C12q- OaAs (100) at a surface temperature of 450°C. The channel width is 81.9 #s and the initial gas pulse is ~ 35 #s. Also shown is the forward convolution backfit (with two components) of Eq. (1) and the two-exponential model of Eq. (3) (zl =0.507 ms, z2=6.148 ms, B=0.1537).
528
P.N. Brier et al./Surface Science 365 (1996) 525-534
elements in # ( f ) with f--
unrestricted power spectrum, ensures that there are non-zero elements in G(f) at all frequencies. This is indeed what happens with the Maxwell-Boltzmann expression where H ( f ) has very small values at high frequencies. 3.3. Filtering We can alleviate the problem by applying a filter 4)(f) to remove high-frequency components in G(f) so that they are not amplified during the deconvolution. F(f) =
G(f).~(f) = G(f). W(f). n(f)
(6)
The simplest filter is a simple cut-off where all
Frequency/ Hz 0
_
2'0t
_
1000
2000
3000
4000
5000
I
I
I
I
I
6000
Inverse
~Q. 1.5
Maxwell-Boltzmann x 10-6 E
o 1.00 o LL
TOF Data
~
j" J
~:=~
05-
S/'t
"
lI-Boltzmann
.,
J
7
0.01.0"o
E 0.5o t-
v
v m
"v VVV~v"v~-'~-V
0.0-
o
~"
^V^V"v^v^v' " ^ ^ v"v , ,
i'o time I ms
--
^^ v V
;5
Fig. 2. Upper panel: The frequency components of the TOE data of Fig. 1 and the Maxwell-Boltzmann function for T = 450°C. Also shown is the inverse of the Maxwdl-Boltzmann components, clearly showing the large magnification factor at high frequencies.
P.N. Brier et al./Surface Science365 (1996) 525 534 (in the least-squares sense) is given by [15,16]:
~ _
529
n(f)*
lD(f)[ 2
IN(012
W(f) =
"..~.~_ Frequency
Fig. 3. Typical power spectrum. ID(f)12is the power spectrum of the total measuredsignal, IN(f)12is the noisepower spectrum, and IG(f)12is the required signal power spectrum. The noise spectrum shownis for whitenoise with a flat frequencyresponse.
tial first order steps, as has been observed by Eldridge and Yu for oxygen on Si(100) [ 12], rather than two simultaneous ones. Unfortunately a careful study with synthetic data in the absence of noise reveals that the truncation of high frequency components makes it impossible to mimic the sharp fall from t = 0 of the forward convolution model. The earlier success of Eldridge and Yu [12] in deconvoluting the observed T O F was in part due to a fortuitous form of the surface time profile. The necessity to include high-frequency components in G(f) greatly hinders the use of the deconvolution method for processes involving a sharp fall from t = 0.
3.4. Optimum (Wiener)filtering Alternative iterative deconvolution methods are widely used. However, it can be shown that iterative deconvolution techniques, such as the Morrison [13] and van Cittert algorithms [14], can be expressed as a simple filter operation in Fourier space [ 15]. The question then arises as to what is the optimum filter to use in Eq. (6)? If we explicitly include noise, n(t), in our description we have a measured signal d(t) given by:
iN(f)l 2 . IH(f)12 + iV(f)12
(8)
where * represents the complex conjugate. The optimum filter thus involves not only the frequency spectrum of the response function (Maxwell-Boltzmann) but also the power spectra of the noise and F(f). Obviously, since F(f) is what we are seeking, we cannot make use of Eq. (8) directly - but we can construct some approximations to Eq. (8), based on various assumptions about N ( f ) and F(f), which will be better than a simple cut-off filter.
3.5. Approximate optimum filters 3.5.1. Constant noise power Fig. 3 shows a schematic view of a typical power spectrum. We can write Eq. (8) as
w(D=
1
[~(f)l 2
H(f) IG(f)l 2 + IN(f)l e
1 IO(f)12--lN(f)l 2 - H(f~) ID(f)] 2
(9)
In the absence of a separate determination of
N ( f ) the best assumption is that we have white noise with a flat noise spectrum. We can then use the high frequency end of the spectrum (where D ( f ) = N ( f ) ) to determine an average value for the noise power. N2av = (IN(f)] 2> for high frequencies
1 W ( f ) = H(f~)
ID(f)12-N2v
ID(f)l 2
(10)
Unfortunately this filter does not work well if Nay is small since at this limit it approaches the simple deconvolution form of Eq. (5) and suffers similar problems of noise amplification.
3.5.2. Constant noise:signal power ratio Returning to Eq. (8) we can replace the noise:signal term by a constant at all frequencies.
d(t) = g(t) + n(t) { hence D ( f ) = G ( f ) + N ( f ) }. (7) It can be shown that the optimal form for W(f)
_ ~ LN(f)[2 S n(t) 2 dt IF(f)[ 2 - Sf(t) 2 d~
(11)
530
P.N. Brier et al./Surface Science 365 (1996) 525-534
and for small signal components (or large noise)
H(f)* W(f) = iH(f)l z + ¢ .
(12)
The average noise:signal power ratio can often be estimated from the data or we can simply choose a small value for ~b, 10-3-10 -6 say. The filter in Eq. (12) is robust and fairly insensitive to the value of ¢ chosen. It utilises the most important (large value) frequency components in H(f) directly while automatically limiting the magnification from small values - since the filter tends to zero when H(f) is small i.e.
H(f) >>¢, W(f)--+1/H(f),
N]v
iFe(f)12 >>In(f)l 2, W(f)~O Eq.(13) has the advantages that it explicitly makes use of an approximate form for F(f) in constructing the filter and that the only subjective decision is the frequency above which to estimate the average noise power. In reality the difference between this filter and one for constant ~b is often very small. The result of applying Eq. (13) is also shown in Fig. 4. All filters suffer the similar problems at t = 0.
H(f) < ¢, W(f)--+H(f)/~ ~ O.
3.6. Model-freeforward convolution
Fig. 4 shows the form of Eq. (12) for two values of ¢ and the result of applying Eq. (12) to the data of Fig. 1 with ¢ = 5 x 10 -3. The result is an improvement on a simple cut-off filter but still suffers from the difficulty that it cannot reproduce the behaviour at t = 0 because of the truncation of the high frequency components.
A final approach is that of a model-free forward convolution. For n data points {di} we represent f(t) by n points {f/} that can be freely adjusted and use the forward convolution procedure of Eq. (1) to calculate the resulting T O F spectrum {gl) for comparison to the measured TOF data. The values of the {f~} points are adjusted to minimise the least-squares fit to the data. In matrix notation
3.5.3. Estimated F(f) The true optimum filter Eq. (8) contains F(f). Although we do not know this a priori we can use the results from a filter with constant ~b to obtain an approximation, F¢(f), to F(f) and hence construct an improved version of Eq. (12). In the absence of an independent determination of N ( f ) we again use an average value constructed from the high frequency part of D(f).
H(f)* F4,(f)=D(f ) iH(f)l z +gP and (13)
(14)
hki f .
and we minimise ~(2= ~
(gk--dk)2
k =1
(~k2
k=l
O'k2 i=1
= 2
~
hki(gk-- dk)
k =1
tTk2
2
~ fi
N Zv IH(f)l 2 + IFe(f)lz
~
i=1
~2
H(f)* W(f)=
gk =
,15,
using an iterative estimation of {f/} based on the steepest descent method for minimising Z2 [17]. ~X2
This filter is also robust at controlling magnification effects since for large signal components (or small noise) N2v - -
IF,(f)l 2
1
2,
W(f)~
- -
n(f)
f'~+1=f'r + Af'~ - -
0f,
(16)
where 2 is a scaling parameter used to control the step-size and aid convergence. This method requires an initial set of {f} values and, like all iterative methods, a good starting set greatly aids
531
P.N. Brier et aL/Surface Science365 (1996) 525 534
Frequency / Hz
,,-,, "10
1000
2000
3000
4000
5000
I
I
I
I
I
6000
0 = 5x10 -3
1.0-
09
E O
Z v
0.5-
n
0.0 1.0-
"o f,O
t~
E k--
0.5-
o cv
v ~ n
v
0.0-
T
0
5
I'0
I'5
--
V
2'0
time / ms Fig. 4. Upper panel: The frequency components of the filter W(f) of Eq. (12) for two different values of ~. They have been normalised to their peak value to facilitate comparison.
the convergence. We can use the results of either the best-model forward convolution procedure or the Wiener filter deconvolution to indicate the likely shape of f ( t ) . To test the convergence a number of starting sets {f/} are used. Fig. 5 shows the results of the best two-exponential model Eq.(3), the cut-off filter Eq.(6), the optimum Wiener filter Eq. (13), and the model-free forward convolution Eq. (16). In the latter case the starting set {f,.} was a single exponential curve with a time constant intermediate between the two values found for Eq.(3). The results of Eq.(16) have clearly converged to excellent agreement with the two exponential fit and show no evidence of a peak at t > 0 as suggested by the Wiener filter. We have sucessfully applied the model-free for-
ward convolution approach to the deconvolution of GaC1 T O F spectra over a range of surface temperatures and timescales, These all follow the same pattern of a two-exponential decrease. In order to test the algorithm further we have applied it to a range of synthetic data based on different kinetic models. The upper panel of Fig. 6 shows synthetic data generated by assuming GaCI is produced by two sequential first-order processes (as in the case of O2 on Si (100) [-12]) with the same time constants as before. 1 T~(t)
=
-
-
"El --'~2
[exp ( - t/z 1) - exp ( - t/r2)].
(17)
The model function is convoluted with the
532
P.N. Brier et al./Surface Science 365 (1996) 525-534
1.0"
I
t~ ~0.5 t-
v
v
I 0.0
o
)
"'s . . . . .
;
4
time / ms Fig. 5. A comparison of the behaviour off(t) at early time for all the methods. The curves are: the best forward convolution model of Eq. (3), the model-free forward convolution of Eq. ( 16 ) ( both solid ), the cut-off filter in Eq. (6) (dot-dash), and the Optimum Wiener filter (dashed). Only the forward convolution methods have the correct behaviour at t=0.
10
+
g
1.0-
0.5-
0.0 0
5
10
15
20
time/ms Fig. 6. Upper panel: synthetic data for two sequential first-order processes (z1=0.507 ms, %=6.148 ms), convoluted with a Maxwell-Boltzmann distribution at 450°C and Gaussian noise added. Lower panel: the model function of Eq. (17) (thick line) and the result of the forward convolution procedure (thin line).
P.N. Brier et al./Surface Science 365 (1996) 525-534 Maxwell-Boltzmann distribution at 450°C and Gaussian noise added. The lower panel shows the result of the forward convolution procedure using exactly the same starting guess used for Fig. 5. This initial guess, a single exponential decay, now has entirely the wrong shape since it has its maximum value at t = 0 instead of the required zero. The forward convolution model still converges to the correct shape. The success of the model-free forward convolution is at first surprising since n freely adjustable points would seem to be a recipe for an unstable solution. However, the behaviour of Eq.(16) is found to be very good providing the 2 value is kept small to inhibit large steps. The majority of deconvolution problems involve correction for instrument response where the convolution function is much narrower than the measured spectrum. In such a case a value at f ( t ) contributes only to a small number of data points close to t and noise can significantly affect the reconstruction o f f ( t ) locally. The opposite is the case here where the Maxwell-Boltzmann convolution function is very broad and each point in f ( t ) c o n t r i b u t e s to a very large range of experimental data points. This has the effect of greatly reducing the effect that noise has on determining t h e f ( t ) value and the n-point model-free set of points {.[i} smoothly converge to a stable solution.
4. Conclusions The task of obtaining the underlying shape of surface time profiles from the measured T O F spectra from pulsed molecular beam scattering is not a simple one. The three methods used, modelbased forward convolution, Wiener filter deconvolution, and model-free forward convolution all have their strengths and weaknesses. These can be summarised by: 4.1. Model-based forward convolution 1. Robust. 2. Accurate. 3. Requires many trial models.
4. Requires iterative optimisation parameters.
533 of model
4.2. Wiener filter deconvolution 1. 2. 3. 4.
Robust. Quick and simple. Insensitive to exact form of filter used. Poor for sharp features due to loss of high frequency components.
4.3. Model-free forward convolution 1. Can be unstable need to control step size. 2. Requires reasonable starting set of values. 3. Several starting sets required to test convergence. 4. G o o d noise rejection in final f(t). 5. Can reproduce sharp features. An approach that utilises all three techniques provides the best method for obtaining the surface time profiles. Initial use of the Wiener filter provides a rapid and robust method of obtaining a good indication of the overall behaviour. This is used to select a suitable set of starting values for the model-free forward convolution procedure. Finally this can be used to suggest suitable analytical models, usually with a physical interpretation behind them, for final fitting by forward convolution.
Acknowledgements This work was supported by the Engineering and Physical Science Research Council and the DRA, Electronics Division at Malvern.
References [ 1] M.Asscherand G.A. Somorjai, in: Atomic and Molecular Beam Methods, Vol.2, ed. G. Scoles (Oxford University Press, Oxford, 1992) p. 488. [2] J.A. Schwartz and R.J. Madix, Surf. Sci. 46 41974) 317. [3] M.P. D'Evelyn and R.J. Madix, Surf. Sci. Rep. 3 (1984) 413. 1-4] C.T. Foxon, M.R. Boudry and B.A. Joyce, Surf. Sci. 44 (1974) 69.
534
P.N. Brier et al./Surface Science 365 (1996) 525-534
[5] H.H. Sawin and R.P. Merrill, J. Vac. Sci. Technol. 19 (1981) 40. I-6] D.R. Olander, J. Coll. Interface Sci. 58 (1977) 169. [7] H.C. Chang and W.H. Weinberg, Appl. Surf. Sci. 3 (1979) 168. [8] M.L. Yu and L.A. DeLouise, Surf. Sci. Rep. 19 (1994) 285. I-9] P. Bond, P.N. Brier, J. Fletcher, P.A. Gorry and M.E. Pemble, Chem. Phys. Lett. 208 (1993) 269 and 211 (1993) 649. [10] P. Bond, P.N. Brier, J. Fletcher, P.A. Gorry and M.E. Pemble, in preparation [11] D.W. Marquardt, J. Soc. Ind. Appl. Math. 2 (1963) 431.
[12] B.N. Eldridge and M.L. Yu, Rev. Sci. Instrum. 58 (1987) 1014. [13] J.D. Morrison, J. Chem. Phys. 39 (1963) 200. [14] P.H. van Cittert, Z. Physik 69 (1931) 298. [15] C.C. Watson and D.V. Ellis, Nucl. Geophys. 5 (1991) 451. [16] N. Wiener, Extrapolation, Interpretation and Smoothing of Stationary Time Series (Technology Press of MIT and Wiley, New York, 1941) [ 17] P.R. Bevington and D.K. Robinson, Data Reduction and Error Analysis for the Physical Sciences, 2nd ed. (McGraw Hill, New York, 1992)