Theoretical Parameter-Free Analysis Model for Temperature-Programmed Desorption (TPD) Spectra.

This paper proposes a parameter-free mathematical model of analyzing either monosite or multisite temperature-programmed desorption (TPD) spectra. By linearizing the integral function difference, the desorption kinetic parameters, such as the desorption order n, the desorption activation energy E d, and the pre-exponential factor ν, can be extracted simultaneously with promising accuracy. A custom "ant" is further established in the model to explore the spectra by a "prediction-correction" loop, and the kinetics and the coverage distribution of the individual peak in the spectra can be solved sequentially. Meanwhile, eight cases on spectrum analysis, including but not limited to the spectrum with coverage-dependent kinetics, the spectrum affected by the noise, the practical spectrum, are demonstrated to fully understand the model's principle, process, and application. Moreover, the model optimization and resolution limitation are further discussed to stimulate the future potential of the innovative parameter-free model.

Introduction

Temperature-programmed desorption (TPD) is one of the powerful techniques developed for surface catalysis study.[1−3] The adsorbate interaction,[4] the molecule surface reaction,[5] and the standard entropies of adsorbed molecules[6,7] can be evaluated by analyzing their TPD kinetics. For instance, the interaction between CO2 and the dissociative water on the CaO films,[8] or the interaction between the water and the two-dimensional silica and aluminosilicate bilayers,[9] can be revealed by analyzing the TPD spectra. Besides, the TPD kinetics can help identify the surface reaction route for the redox reaction of H2 and H+ on the platinum/nitrogen-doped mesoporous carbon composites.[10] Usually, the adsorbate desorption rate is given by the Polanyi–Wigner equation[7,11]where θ is the transient coverage, t is the time, ν is the pre-exponential factor, n is the desorption kinetic order, Ed is the desorption activation energy, R is the ideal gas constant, and T is the temperature.

To extract the TPD kinetic parameters, Redhead proposed a method of determining the desorption activation energy based on the temperature corresponding to the maximum desorption rate. This method had some limitations, such as the activation energy and pre-exponential factor that were independent of the surface coverage, and the heating rate had to change by at least two orders of magnitude to achieve reasonable accuracy.[12] King overcame those limitations by establishing a complete lineshape analysis method, which, however, relied on a reasonable assumption of the pre-exponential factor.[11] Although Habenschaden and Küppers eliminated the guess of the pre-exponential factor, their method was strongly limited to the onset of the desorption at very low temperature, which required adequate measurement accuracy.[13] Miller et al. further classified those methods into two groups, namely, the integral and differential approaches. The integral approaches, such as the Redhead analysis, were thought to extract the coverage-independent kinetic parameters,[14] while the differential approaches, such as the lineshape-related analysis, could predict reliable coverage-dependent results.[15] In addition, other methods based on the statistical analysis[16] or the process variation[17,18] were also discussed. However, most methods or models assumed a preliminary pre-exponential factor,[19,20] or the parameter dependency,[11,21,22] and the experimental conditions were strictly limited to a specific process.[23] As a result, a variational method was proposed to determine the pre-exponential factor by minimizing the difference between the simulated TPD spectrum and the experimental results.[20,24,25] Alternatively, the pre-exponential factor was also approximated by the peak’s position and width.[26,27]

Nevertheless, few works made the simultaneous predictions of the kinetic order, the desorption activation energy, and the pre-exponential factor without assumptions. This paper innovatively establishes a parameter-free analysis model by linearizing the integral function difference between every two adjacent discretized nodes. Also, eight cases are demonstrated to help fully understand the model’s principle, process, and application. Besides, further discussions on model optimization and limitation are given to stimulate the future potential of the proposed model.

Results and Discussion

Theoretical Basis of the Parameter-Free Analysis Model

Under a linear temperature program, the Polanyi–Wigner equation giveswhere β = dT/dt is the heating rate.

Its integral form g(θ) over [θ, θ – Δθ] can be written aswhere Δθ is the change in the coverage over ΔT, while θ and (θ – Δθ) are the transient coverages of the two adjacent nodes in the corresponding segment, respectively, which would be explained later.

In case of zero/first/second-order desorption, g(θ, θ – Δθ) can be determined asHowever, the term on the right hand of eq does not have an analytical solution. Therefore, this term is usually replaced by several useful approximations.[28,29] Inspired by Ortega’s linear integral method,[30,31] when ΔT → 0 and Δθ → 0, the desorption activation energy Ed, the pre-exponential factor ν, and the experimental condition β in each segment can be considered as constant. Thus, eq can be rewritten asThen,The flow chart for analyzing a monosite TPD spectrum by the proposed parameter-free model is summarized in Figure a. A look at the TPD spectrum shows that the model first calculates the initial adsorbate coverage θ0 and the spectrum is discretized into segments by a user-defined ΔT as explained in Figure b. Then, g(θ, θ – Δθ) at every two adjacent nodes, which are marked yellow, are determined, and the corresponding natural logarithm results are plotted against the reciprocal of the temperature. The kinetic order n is judged by the goodness-of-linear fit of the different order desorption equations in eq . Accordingly, the desorption activation energy Ed and the pre-exponential factor ν are given by the slope and intercept of the best-fitting results, respectively.

Figure 1

Model description for analyzing a monosite TPD spectrum: (a) the prediction flow chart and (b) a user-defined ΔT-based discretization.

For a multisite TPD spectrum, the prediction flow chart is illuminated in Figure .

Figure 2

Flow chart for analyzing a multisite TPD spectrum by the proposed parameter-free model.

Taking a two-site TPD spectrum as an example, the model first reads the spectrum and calculates its initial coverage θ0. Because there are two desorption peaks in the spectrum, the discretized nodes cannot be treated as a whole. Alternatively, in addition to a user-defined ΔT, an “ant”, which consists of a certain number of continuous discretized nodes, is defined. Therefore, the ant has two parameters, namely, the ant’s length and its moving step in the form of the node number. Figure a schematically shows that a user-defined ant has a length of six discretized nodes. In the first cycle, the process is the same as the explanations in the aforementioned monosite TPD spectrum analysis, g(θ, θ – Δθ) at every two adjacent nodes of the six nodes are calculated sequentially. Then, their natural logarithm results are linearly fitted against the reciprocal of the temperature. Therefore, the kinetic parameters and the coefficients of determination (R2) are given. The ant continues to explore by moving two nodes as shown in Figure b. The updated six nodes are processed again in the same way to obtain the new results in the second cycle. Thereafter, the ant moves cycle by cycle to access all the nodes until the end, as shown in Figure c. As a result, the predicted kinetic parameters with respect to the increasing cycles are depicted, and the “terraces”, or the constant results are searched to identify the exact kinetic parameters for either peak.

Figure 3

Schematic diagram of a user-defined ant’s (a) length, (b) moving step, and (c) final cycle when analyzing a two-site TPD spectrum.

In the next step, the model introduces a “prediction–correction” loop to calculate the coverage θ of the previously identified peak. The θ–T curve is first initialized by the Runge–Kutta method based on the obtained kinetic parameters and a guessed coverage θ. Subsequently, the θr–T curve is obtained by subtracting the θ–T curve from the θ0–T curve. Its coverage θr is then calculated and checked if it is less than 0. If not, the previously guessed coverage θ is corrected by (θ + Δθ′). The loop does not terminate until the calculated θr is less than 0, and the peak’s coverage θ is finalized by subtracting Δθ′ from the coverage θ in the last loop. At last, the corresponding θr–T spectrum is revisited by the ant to extract its own kinetic parameters. Case II in the following section would demonstrate in detail how does the model handle a two-site TPD spectrum.

Case Studies for Model Basic Analysis

In this section, two cases are demonstrated to use the proposed model to analyze a monosite and a multisite artificial TPD spectrum, respectively.

Case I: Three typical monosite TPD spectra are firstly generated based on the Runge–Kutta method in Figure a, namely, the zero-order kinetic TPD spectrum (Ed = 30 kJ·mol–1, ν = 1 × 108 ML·s–1), the first-order kinetic TPD spectrum (Ed = 60 kJ·mol–1, ν = 1 × 1018 s–1), and the second-order kinetic TPD spectrum (Ed = 50 kJ·mol–1, ν = 1 × 1015 ML–1·s–1). The model discretizes these curves and calculates the difference of the g functions between any two adjacent nodes with three different kinetic order equations. The natural logarithms of the results are plotted against the reciprocal of the temperature in Figure b–d, respectively. Taking the artificial zero-order TPD spectrum as an example, Figure b clearly shows that the results of the g function in the form of the zero order have a perfect linear relationship, and the linear line’s slope and intercept are determined as −3605.04 and 16.10, respectively. Therefore, the input spectrum has a zero kinetic order and the Ed and ν are predicted as 29.97 kJ·mol–1 and 9.82 × 107 ML·s–1, respectively. All of the predicted kinetic parameters are compared with the input counterparts in Table , and the consistent results prove that the proposed model is reliable for extracting kinetic parameters from a monosite TPD spectrum.

Figure 4

Model analysis on three artificial monosite spectra: (a) the input monosite spectra and (b–d) the model prediction results in the form of zero-, first-, and second-order equations for each input spectrum.

Table 1

Comparison of the Input Kinetic Parameters and the Model Prediction Results for Case I

	input kinetic parameters			model prediction results
n	Ed (kJ·mol–1)	ν	n	Ed (kJ·mol–1)	ν
0	30	1 × 108 ML·s–1	0	29.97	9.82 × 107 ML·s–1
1	60	1 × 1018 s–1	1	59.73	8.24 × 1017 s–1
2	50	1 × 1015 ML–1·s–1	2	49.90	9.36 × 1014 ML–1·s–1

Case II: In Figure a, the main site desorption spectrum follows the first-order kinetics (Ed = 60 kJ·mol–1, ν = 1 × 1018 s–1), while the defect site desorption spectrum has second-order kinetics (Ed = 40 kJ·mol–1, ν = 5 × 1010 ML–1·s–1). Then, the two spectra are merged to produce a two-site TPD spectrum in Figure b.

Figure 5

Model analysis on a two-site spectrum: (a) the individual spectra simulated for the main and the defect sites (b) and then merged into a two-site spectrum for the model analysis; the comparison of the desorption activation energy, pre-exponential factor, and coefficient of determination in each cycle solved by the g functions in terms of the (c) zero-, (d) first-, and (e) second-order kinetics, respectively.

The two-site spectrum is first discretized into nodes by ΔT = 0.1 K, and a custom ant consists of 100 nodes; therefore, the ant itself covers a temperature range of 9.9 K. The ant’s step size is 20 nodes, which means that it will move by 2 K in each cycle. In this case, the start and end temperatures for the model analysis are 102 and 247.9 K, respectively. Therefore, the total number of loop cycles is 69. The natural logarithm of the g functions in terms of the zero-, first-, and second-order kinetics are calculated separately in Figure c–e, and then the results are linearly fitted against the reciprocal of the temperature to determine the desorption activation energy, the pre-exponential factor, and the coefficient of determination (R2).

By searching and identifying the terraces in the diagram, the model predicts that there is a second-order spectrum after 40 cycles in Figure e, and the predicted results are summarized in Table S1. Starting from 180 K, the desorption energies and the corresponding pre-exponential factors are almost constant, and the values of R2 are close to 1. Therefore, the average desorption energy between 180 and 223.9 K is approximately 40 kJ·mol–1, while the average pre-exponential factor is 4.98 × 1010 ML–1·s–1.

In the next step, the coverage θ of the above-mentioned second-order spectrum is initially guessed as 0.1 ML, and its coverage distribution θ–T (red solid line) is calculated by the Runge–Kutta method in Figure a. Then, the residual coverage distribution θr–T (red dotted line), is obtained by subtracting θ–T from the original coverage distribution θ0–T (blue solid line). Thereafter, as shown in Figure b–e, the loop begins by correcting the guessed coverage θ to (θ + Δθ′), where Δθ′ = 0.1 ML, and the loop does not terminate until the residual coverage after the subtraction at about 175 K is less than 0 in the last loop when θ = 0.5 ML. Therefore, the coverage θ of the above-mentioned second-order spectrum is 0.4 ML, which is the correct value in the previous loop. In short, a spectrum between 180 and 223.9 K is predicted to have a coverage of 0.4 ML with the second-order kinetics (Ed = 40 kJ·mol–1, ν = 4.98 × 1010 ML–1·s–1).

Figure 6

Prediction–correction loop of determining the coverage distribution for the second-order spectrum when the prediction coverages are (a) 0.1 ML, (b) 0.2 ML, (c) 0.3 ML, (d) 0.4 ML, and (e) 0.5 ML, respectively.

As a result, the residual coverage distribution θr–T is employed for the model analysis again. The left spectrum is roughly located between 154 and 174 K or between 27 and 32 cycles, and its kinetic parameters are extracted in a similar way; the results are compared in Figure . Only the first-order fitting results have the terrace-like behavior, and the values of R2 are close to 1, so the spectrum follows the first-order kinetics (Ed = 59.83 kJ·mol–1, ν = 9.27 × 1017 s–1).

Figure 7

Comparison of (a) the desorption activation energy, (b) the pre-exponential factor, and (c) the R2 extracted by the proposed model for the left spectrum between 27 and 32 cycles.

Case Studies for Model Advanced Analysis

In this section, another six cases, namely, Cases III–VIII, with either an artificial or a practical TPD spectrum, demonstrate the model’s potential of handling advanced analysis.

Case III: The heating rate is an important parameter in the TPD experiment. This case simulates the zero-order TPD spectra (Ed = 30 kJ·mol–1, ν = 1 × 108 ML·s–1) under three heating rates, namely, 0.8, 1.0, and 1.2 K·s–1, in Figure S1a. The model analysis results are demonstrated in Figure S1b–d. The results indicate that all the three spectra follow the zero-order desorption kinetics. Besides, the desorption energy Ed is approximately 29.97 kJ·mol–1, while the pre-exponential factors are between 9.79 × 107 and 9.84 × 107 ML·s–1. Therefore, the proposed model is proved to extract the kinetic parameters from the TPD spectra regardless of the heating rate.

Case IV: The initial surface coverage is another parameter in the TPD experiments. This case takes the second-order TPD spectrum (Ed = 50 kJ·mol–1, ν = 1 × 1015 ML–1·s–1) as an example, the initial surface coverage increases from 0.5 to 1.0 ML, then to 1.5 ML in Figure S2a. The model analysis results are compared in Figure S2b–d. All the spectra are proved to be of the second order, and the desorption energy Ed and the pre-exponential factor ν are evaluated to be 49.90 kJ·mol–1 and 9.36 × 1014 ML–1·s–1, respectively. In addition, the intercepts of the linear fitting lines in Figure S2b–d are constant, while that in Figure S1b–d vary from 15.92 to 16.32. It may indicate that the pre-exponential factor is more sensitive to the heating rate than to the initial coverage.

Case V: Different from the spectrum with the coverage-independent kinetic parameters in Case IV, this case employs three monosite TPD spectra with coverage-dependent parameters. The linear relation between the desorption energy and coverage is given by eq , while the compensation effect between the desorption energy and the pre-exponential factor is given by eq .The three first-order TPD spectra with coverage-dependent parameters are given in Figure S3a. The initial coverage increases Ed and ν simultaneously, from 57.5 kJ·mol–1 and 1.53 × 1017 s–1 when the coverage is 0.5 ML to 60.0 kJ·mol–1 and 1.00 × 1018 s–1 when the coverage is 1.0 ML, and then to 62.5 kJ·mol–1 and 6.52 × 1018 s–1 when the coverage is 1.5 ML, respectively. Those spectra are analyzed separately, and the results are shown in Figure S3b–d. First, all the spectra are proved to be of the first kinetic order. Second, the desorption energy and the pre-exponential factor are calculated and validated against the input parameters in Table . The consistent results agree that the proposed model is able to deal with the TPD spectrum with coverage-dependent kinetic parameters.

Table 2

Comparison of the Input Kinetic Parameters and the Model Prediction Results for Case V

	input kinetic parameters			model prediction results
n	Ed (kJ·mol–1)	ν (s–1)	n	Ed (kJ·mol–1)	ν (s–1)
1	57.5	1.53 × 1017	1	57.26	1.28 × 1017
1	60	1.00 × 1018	1	59.78	8.66 × 1017
1	62.5	6.52 × 1018	1	62.16	5.08 × 1018

Case VI: In this case, as shown in Figure a,b, the simulated TPD spectrum first follows the first order (Ed = 50 kJ·mol–1, ν = 1 × 1015 s–1) from coverage = 1.0 ML. Also, it later shifts to follow the second order (Ed = 35 kJ·mol–1, ν = 1 × 1010 ML–1·s–1) when the transient coverage is 0.8249 ML at 160 K. The input spectrum and its coverage distribution against the temperature (blue solid line) are depicted in the insets. Thereafter, the input spectrum is discretized by ΔT = 0.1 K. The ant in the proposed model consists of 100 nodes, and its moving step is 10 nodes. The terrace-searching results in Figure c clearly show that the spectrum between 47 and 50 cycles, corresponding to the temperature between 147.0 and 159.9 K, follows the first-order kinetics, while that between 61 and 64 cycles, corresponding to the temperature between 161.0 and 173.90 K, follows the second-order kinetics. Besides, the prediction kinetic parameters are quite consistent with the input values (black dotted line). Therefore, the ant in the proposed model has powerful flexibility in analyzing the complex TPD spectrum.

Figure 8

(a) Simulated TPD spectrum (blue) is generated by a first-order spectrum (black) starting from coverage = 1.0 ML and later shifting to a second-order spectrum (red) when (b) the transient coverage is 0.8249 ML at 160 K, and (c) the terrace-searching results by the ant in the proposed model are demonstrated in terms of the kinetic desorption energy, the pre-exponential factor, and coefficient of determination.

Case VII: The noise is very common in the TPD measurement. In this case, an artificial second-order TPD spectrum (Ed = 50 kJ·mol–1, ν = 1 × 1013 ML–1·s–1) (red) is first generated in Figure . Then, three new spectra (blue) in Figure a–c are produced by adding 5, 10, and 15% relative errors to the previous spectrum, respectively. Meanwhile, another spectrum (blue) in Figure d is simulated with a fixed error range of ±0.001 ML·s–1, which is randomly added to the spectrum. Without data smooth operation, the noise-affected spectra are discretized by ΔT = 0.1 K. Then, the 5% relative noise affected spectrum is examined by the proposed model when the ant’s length has 500 nodes and its moving step is 10 nodes. For the 10 and 15% relative noise affected spectra, the ant’s length is optimized to 530 nodes. Since the fluctuation brought by the noise definitely causes the values of the nodes to be unstable, the R2 results accordingly decrease. After 90 cycles, or 190.0 K, the greatest R2 results are reached, and the prediction kinetic parameters for the 5% relative noise affected spectrum in each cycle are summarized in Table S2. In contrast, the main peak in Figure d is much less affected by the fixed noise, and the values of R2 are over 0.84 after 81 cycles at 181.0 K, and the accompanying prediction parameters are summarized in Table S3. All the prediction results are compared with the input parameters in Table . As a result, the proposed model is able to extract the kinetic parameters from a TPD spectrum with either a relative or a fixed noise under present simulation conditions.

Figure 9

Original simulated TPD spectrum (red solid line) and the reproduced spectrum with (a) 5%, (b) 10%, (c) 15% relative noise and (d) the fixed noise within ±0.001 ML·s–1 (blue solid line).

Table 3

Comparison of the Input Kinetic Parameters and the Model Prediction Results for Case VII

		input kinetic parameters			model prediction results
noise type	n	Ed (kJ·mol–1)	ν (ML–1·s–1)	n	Ed (kJ·mol–1)	ν (ML–1·s–1)
5%-relative	2	50	1 × 1013	2	50.26	0.95 × 1013
10%-relative				2	49.34	1.22 × 1013
15%-relative				2	48.04	1.09 × 1013
fixed (±0.001 ML·s–1)				2	50.38	1.24 × 1013

Case VIII: In the last case, the experimental TPD spectrum of ammonia adsorbed at 293 K (initial surface coverage = 100%) by Abello et al.[32] is reproduced to further test the proposed model. The raw data (Table S4) are reproduced in Figure a and then discretized by ΔT = 1 K. The ant’s length is optimized to have 50 nodes, and the ant moves by 1 node in each cycle. The model only predicts the proper desorption energy distribution under the zero kinetic order, and the results are depicted in Figure b. A quasi-terrace appears with a relatively steady R2 between 81 and 92 cycles, which correspond to the temperature between 376 and 436 K. Therefore, the average desorption energy is approximately 53.23 kJ·mol–1, which is consistent with the reported results of 29.26–66.04 kJ·mol–1.[32−34]

Figure 10

(a) Reproduced experimental TPD spectrum of ammonia adsorbed at 293 K (initial surface coverage = 100%) by Abello et al. (reuse with permission from ref (32) Copyright 1995 Published by Elsevier B.V.) and (b) the prediction desorption energy distribution with respect to the calculation cycle under zero-order kinetics.

Further Discussion

Discretization Parameter ΔT

The proposed parameter-free analysis model is established based on the linear integral function difference, so the discretization parameter ΔT directly affects the prediction accuracy. Therefore, taking a first-order kinetic TPD spectrum (Ed = 60 kJ·mol–1, ν = 1 × 1018·s–1) as an example, the kinetic parameters are predicted when ΔT is 0.01, 0.1, and 1 K, respectively. The results are compared in Table . The finer ΔT generally improves the prediction accuracy of the desorption energy, while it has a limited impact on the pre-exponential factor. On the other hand, the finer ΔT does not only require the accuracy of the experimental measurement but also requires more computing resources for model prediction. Therefore, it is better to choose a reasonable discretization parameter by model optimization.

Table 4

Comparison of Model Prediction Kinetic Parameters under Different ΔT

	n	Ed (kJ·mol–1)	ν (s–1)
input	1	60	1 × 1018
output (ΔT = 0.01 K)	1	60.15	1.12 × 1018
output (ΔT = 0.1 K)	1	59.73	8.24 × 1017
output (ΔT = 1 K)	1	59.66	8.98 × 1017

Noninteger Kinetic Order

Although the proposed model only shows how to predict an integer kinetic order of a TPD spectrum in the case studies, in theory, it can extract a noninteger kinetic order. As a result, the integral g function difference in eq can be rewritten in a more general form

Peak-Distinguished Resolution

Peak-distinguished resolution in a multisite TPD spectrum is of great interest for evaluating the proposed model. It is worth stressing again that the ant’s length and its moving step size in the model have a direct and important effect on the peak-distinguished resolution. For instance, the two-site spectrum in Case II has a peak–peak overlap region between 155 and 178 K in Figure S4, where the ant cannot distinguish one from the other. Therefore, it explains why the extracted kinetic parameters fluctuate a lot in this temperature range in Figure . As soon as the ant’s entire body has left the overlap region after about 40 cycles, it continues to follow the defect-site spectrum. Therefore, the prediction results are stable and form the corresponding terraces. Also, the spectra overlap regions are located between 147 and 160 K in Cases VI and VIII, between 161 and 174 K in Case VI, and between 376 and 419 K in Case VIII. Therefore, the custom ants in these cases usually have a shorter length with fewer nodes to reduce the errors caused by fluctuations. Briefly, the peak-distinguished resolution strongly depends on whether the ant’s movement stays within a nonoverlap region or not. On the other hand, the kinetic parameters are predicted by the linear fit over all the nodes inside the custom ant. Therefore, more nodes can produce more reliable results. Besides, the larger step size can accelerate the ant’s movement to save the computing time, while the smaller step size can generate a more significant terrace. To sum up, it is quite difficult to define the standard parameters of the proposed model, so it is recommended to carry out optimization work in advance to figure out the most suitable parameters for a specific case.

Conclusions

A parameter-free analysis model for simultaneously extracting the kinetic order, the desorption activation energy, and the pre-exponential factor of a TPD spectrum is established. The innovative model is feasible in predicting the coverage-independent or -dependent kinetic parameters under different heating rates. In addition to handling a monosite TPD spectrum, the model can also separate a multisite spectrum through a coverage prediction–correction loop method. Besides, eight case studies on the TPD spectrum prediction are demonstrated, including but not limited to the spectrum with coverage-dependent desorption energy, the spectrum with noise, and the practical spectrum, to help fully understand the principle, the process, and the future potential of the proposed model. Moreover, the discretization grid on the prediction accuracy, the noninteger kinetic order prediction, and the peak-distinguished resolution in a multisite TPD spectrum analysis are further discussed.