Analysis of the Diffusion in a Multilayer Structure under a Constant Heating Rate: The Calculation of Activation Energy from the In Situ Neutron Reflectometry Measurement

Understanding mass transport in micro- and nanostructures is of paramount importance in improving the performance and reliability of the micro- and nanostructures. In this work, we solve the diffusion problem in a multilayer structure with periodic conditions under a constant heating rate via a Fourier series. Analytical relation is established between the coefficients of eigenfunctions and the intensity of X-ray or neutron Bragg peak. The logarithm of temporal variation of the intensity of X-ray or neutron Bragg peak is a linear function of the nominal diffusion time, with the nominal diffusion time being dependent on the heating rate. This linear relation is validated by experimental data. The Taylor series expansion of the linear relation to the first order of the diffusion time yields an approximately linear relation between the logarithm of temporal variation of the intensity of X-ray or neutron peak and the diffusion time for small diffusion times, which can be likely used to calculate the activation energy for the diffusion in a multilayer structure. The validation of such an approach is subjected to the fact that the characteristic time for heat conduction is much less than the characteristic time for the ramp heating as well as the characteristic time for diffusion.


■ INTRODUCTION
Understanding atomic diffusion in solids 1 is of practical importance in controlling the synthesis of micro-and nanostructures and the performance and stability of nanoelectronic devices and systems. For example, lithium-ion diffusion is essential for the operation and integrity of a lithium-ion battery. 2−4 In general, atomic diffusion is a thermally activated process, and the Arrhenius relation can be used to describe the temperature dependence of the diffusion coefficient of atomic diffusion as 1 = D D e Q RT 0 / (1) where D is the diffusion coefficient, D 0 is a prefactor, Q is the activation energy, R is the gas constant, and T is the absolute temperature.
There are various techniques available to measure the diffusion coefficient of atomic diffusion in solids. The most common technique is the tracer method, 5−10 which tracks the spatial distribution of isotopic elements at different times under the isothermal condition with different techniques, including secondary ion mass spectrometry (SIMS). 11 −20 Using nuclear magnetic resonance (NMR) 21−24 or quasielastic neutron scattering (QENS) 25−27 makes it possible to analyze atomic diffusion in crystalline materials and in battery research. Diffusion coefficients of atomic diffusion have also been measured by low-angle X-ray diffraction (XRD) (X-ray reflectometry, XRR) in solids, 28,29 XRD, 30,31 energy-dispersive X-ray spectroscopy, 28,32,33 and low-angle neutron diffraction (neutron reflectometry, NR). 34−54 Most studies reported in the literature have been based on the isothermal condition, which require multiple measurements at different temperatures in order to determine the activation energy for atomic diffusion in solids. For example, Mizoguchi and Murata 29 conducted low-angle XRD (XRR) analysis of compositionally modulated amorphous Co−Zr films isothermally annealed consecutively at different temperatures in a range of 120−180°C and calculated the activation energy for the interdiffusion in the modulated amorphous Co−Zr films from the variation of the intensity of the X-ray peak.
Diffusion coefficients and activation energies for selfdiffusion have been obtained with NR 37−54 in a variety of materials. In nonoxide high-temperature superhard ceramic coatings, nitrogen self-diffusion was measured in amorphous Si 3 N 4 35,37 and in silicon carbo-nitrides. 38,39 For crystalline metal oxides, lithium self-diffusion was examined in LiNbO 3 . 40,41 Self-diffusion was measured in metals such as in nanocrystalline copper, 42 iron, 43,44 and Fe-based compounds. 37 In semiconductors, self-diffusion was measured in crystalline silicon 45,46 and germanium 47 and in amorphous silicon 48 and germanium. 49 Extrinsic diffusion (so-called "Fremddiffusion") was also suitably measured with NR 50−54 for lithium permeation through amorphous silicon, 50−52 amorphous lithium silicide, 53 and nanocrystalline chromium 54 thin layers. The diffusivities determined by NR were found to be in agreement with the corresponding ones determined from other established techniques such as SIMS. 40,41,46,47,52 The diffusivities and the activation energies for atomic diffusion were generally measured under isothermal-diffusion annealing conditions. The multi-isothermal annealing for the determination of the activation energy of atomic motion is timeconsuming and can last more than several hours. 49 To substantially reduce the time needed in the measurement of the temperature dependence of atomic diffusion, Huger et al. 55 introduced the concept of increasing temperature at a constant rate in-situ NR measurements for the determination of the temperature dependence of atomic diffusion in solids. They have demonstrated the feasibility of this method in the study of atomic diffusion in amorphous 73 Ge/ nat Ge multilayers. 55 In the heart of this concept is the curve fitting of the intensity of Bragg peaks with a fitting parameter of the activation energy, Q, as 55 where I(t) and I 0 are the intensities of the Bragg peak at time t and initial time (t = 0), respectively, d is the bilayer thickness, T 0 is the onset temperature, and λ is the heating rate (the increasing rate of temperature). It needs to be pointed out that eq 2 is a semi-empirical relation and the integration presented in eq 2 complexes the fitting process. Developing a simple polynomial formulation would make it readily to determine the activation energy from the temporal evolution of the intensity of the Bragg peak. Also, Huger et al. 55 and Schmidt et al. 37 did not provide detailed information on the derivation of eq 2 theoretically as well as any numerical calculation of the concentration evolution during the diffusion annealing of a multilayer structure.
Considering the important role of atomic diffusion in the synthesis of nanostructured materials and in the structural integrity of multilayer structures, we analyze one-dimensional diffusion problem in a multilayer structure and introduce the concept of nominal characteristic time. The multilayer structure consists of a stack of a bilayer structure, as shown in Figure 1, and experiences uniform diffusion heating at a constant heating rate. The effect of the heating rate on the spatiotemporal evolution of the concentration of the diffusive species is illustrated. A polynomial equation is proposed for the determination of the activation energy of the diffusive species for the in-situ measurement of atomic diffusion by X-ray and neutron reflectometry. The lower and upper bounds of the heating rate are discussed. ■ RESULTS Figure 1 shows a multilayer structure consisting of a stack of a bilayer structure. The thickness of each individual layer of the bilayer structure is a 1 for the first layer and a 2 for the second layer, as shown in the figure. Practically, both the layers in the bilayer structure likely contain diffusive species of different amounts (concentrations) before the start of diffusion, with one layer having a higher concentration of diffusive species and the other having a lower concentration of diffusive species. Without loss of generality, we assume that, at the onset of diffusion, there is uniform concentration, c 0 , of the diffusive species present in the first layer, and no diffusive species is present in the second layer. Note that the analysis method presented here can be readily extended to the case with uniform concentration, c 0 , of the diffusive species presented in the second layer and no diffusive species in the second layer. The multilayer structure experiences uniform diffusion heating under a constant heating rate of λ, and the diffusivity of the diffusive species according to eq 1 can be expressed as Note that eq 3 is based on the fact that the diffusivity of the diffusive species is independent of local concentration of the diffusive species and there is no phase change during the heating.
Considering the periodic characteristic of the multilayer structure, the diffusion problem in the multilayer structure can be simplified to as the diffusion problem in a symmetric trilayer structure (see Figure 1) as follows. For one-dimensional diffusion in the trilayer, the diffusion equation is with c being the concentration of the diffusive species and D being given in eq 3. The initial condition is The symmetric characteristic of the problem yields the boundary conditions as Note that the initial condition of (5) can also be used for the outer layers with uniform distribution of the diffusive species, c s0 , by introducing an auxiliary variable of c̃(=c − c s0 ).
Using the method of separation of variables, the concentration of the diffusive species can be expressed as Figure 1. Multilayer structure consisting of a stack of a bilayer structure.
in which C n (x) and Φ n (x) satisfy the following differential equations The solutions of Φ n (x) and C n (x), which satisfy the boundary conditions of (6), are with B n being the coefficients to be determined from the initial condition. The eigenvalues of χ n are = + = n a a n 2 1, 2, 3, ...
Using the initial condition of (5), the coefficients B n are found to be c n na a a n 2 sin 1, 2, 3, ...
Thus, the spatiotemporal evolution of the concentration of the diffusive species is in which 1/2 is included to satisfy the condition of c(x, 0) = 0 for a 1 /2 < |x| < (a 1 + a 2 )/2. According to Mizoguchi and Murata 29 and Schmidt et al., 37 the principles of using X-ray and NR in the characterization of the diffusion coefficient of the diffusive species is that the intensities of the X-ray and NR Bragg peaks are proportional to the square of the coefficient of cos χ n x in eq 13. Thus, we have which gives It is evident that eq 15 is slightly different from eq 2 with the contribution of a constant term to the peak intensity of X-rays or neutrons.
Generally, it would be relatively easy to curve-fit experimental results with an explicit formulation. The integration in eq 15 increases the complexity in the curve fitting, and a simple polynomial relation between I(t) and Q is desirable. Expanding the exponential term in eq 15 in a Taylor series, eq 15 can be expressed as To the order of t 4 , the temporal variation of ln[I(t)/I(0)] is found to be Huger et al. 55 reported Q = 2.14 eV (207.4 kJ/mol) and D 0 = 1.25 × 10 −5 m 2 /s for the self-diffusion of Ge in a multilayer structure with 2 min in the counting time for a single reflectivity pattern. The multilayer structure consisted of 10 bilayers of 73 Ge (165 Å)/ nat Ge (165 Å), which was annealed at a ramping rate of 1°C/min in a temperature range of 320−500°C . Using their data and eq 13, the spatial distribution of the diffusive species (Ge in this case) is shown in Figure 2 for a 1 = a 2 = 165 Å at different instants. At the diffusion time of t = 1 s, there is limited diffusion occurring at the position of x = a 1 /2 with the concentration of the diffusive species at c 0 in the most portion of the region [0, a 1 /2) and at 0 in the most portion of the region (a 1 /2, a 2 /2]. Increasing the diffusion time leads to that more diffusive species migrate across the cross-section of x  Figure 3a shows temporal evolution of the concentration of the diffusive species at the center of the trilayer (x = 0) for five different heating rates of 0, 1, 2.5, 5, and 10 K/min over a period of 12 ks (3.33 h). Without the increase in temperature (λ = 0), there is no observable change in the concentration of the diffusive species over a period of 12 ks. Increasing the heating rate causes the increase in the temperature of the trilayer and accelerates the diffusion of the diffusive species. The time for the onset of the migration of the diffusive species away from the center of the trilayer decreases with increasing the heating rate due to the increase in the diffusivity. Figure 3b presents the temporal evolution of the concentration of the diffusive species at the center of the trilayer (x = 0) with the nominal diffusion time. It is interesting to notice that there is no observable difference of the correlation between the variation of the concentration of the diffusive species at the center of the trilayer and the nominal diffusion time. Such a result suggests that the effect of the constant heating rate on the diffusion of the diffusive species can likely be uniquely described by using the nominal diffusion time given in eq 18.
The variation of the nominal diffusion time with the diffusion time is presented in Figure 4. It is evident that the nominal diffusion time is dependent on the heating rate, as expected. For isothermal heating, there is no significant change in the logarithm of the nominal diffusion time for the diffusion time in the range of 2−12 ks. Note that the nominal diffusion time is proportional to the diffusion time for isothermal heating. The rapid change in the logarithm of the nominal diffusion time occurs for the diffusion time in the range of 0−1 ks, and the range of the diffusion time for the rapid change in the logarithm of the nominal diffusion increases with the increase of the heating rate. Also, the larger the heating rate, the larger the nominal diffusion time at the same diffusion time.
To analyze the variation of the intensities of X-ray and neutron Bragg peaks with the diffusion time, we limit the analysis to the case of n = 1. For n = 1, eq 16 yields and eq 17 gives for the Taylor series expansions of eq 19 to the orders of t 2 and t 4 , respectively. Figure 5a shows the variation of the intensity of X-ray or neutron Bragg peak    Figure 5b for a small diffusion time. For the diffusion time in the range of 0−1 ks, the results from eq 20 are nearly overlapped with the corresponding ones from eq 19; for the diffusion time in the range of 0−3 ks, the results from eq 21 are nearly overlapped with the corresponding ones from eq 19. These trends imply that eqs 20 and 21 can be likely used to approximately calculate the activation energy and the prefactor for small diffusion times (small range of temperatures).

■ DISCUSSION
According to the results shown in Figure 3b, the correlation between the concentration of the diffusive species and the nominal diffusive time defined in eq 18 is independent of the heating rate for the heating rate in the range of 0−10 K/min. Such a result implies that the variation of the intensity of the Xray or neutron Bragg peak (ln[I(t)/I(0)]) with the nominal diffusion time is likely also independent of the heating rate for the heating rate in the range of 0−10 K/min.
Comparing eq 18 with eq 15, we obtain The logarithm of the temporal variation of the intensity of the X-ray or neutron Bragg peak (ln[I(τ)/I(0)]) is a linear function of the nominal diffusion time. Eq 22 reveals that the logarithm of the temporal variation of the intensity of the X-ray or neutron Bragg peak for a multilayer structure under continuous heating with a constant heating rate is a linear function of the nominal diffusion time if the diffusion mechanism remains the same for the temperature range under investigation. Note that such a linear relation is not discussed in the work of Huger et al. 55 Huger et al. 55 also reported the temporal evolution of the integrated and normalized intensity of the Bragg peak for the self-diffusion of Ge in the multilayer structure with 5 min in the counting time for a single reflectivity pattern. Figure 6a shows the variation of the integrated and normalized intensity of the Bragg peak of their experimental data with the nominal diffusion time of eq 18 for three different activation energies of 2.0, 2.14, and 2.29 eV. It is evident that the logarithm of the integrated and normalized intensity of the Bragg peak of their experimental data is nearly a linear function of the nominal diffusion time except the last experimental data of I/I 0 = 0.051. Figure 6b presents linear-regression fitting of the data shown in Figure 6a for I/I 0 larger than 0.051. For the comparison, the fitting curves are also included in the figure. The square of the goodness of fit (R) is 0.99 for all the three fitting, which validates the relation of eq 22.
The curve fitting yields 2π 2 D 0 /a 1 2 to be 0.918 × 10 11 , 0.932 × 10 12 , and 1.11 × 10 13 s −1 for the activation energies of 2.0, 2.14, and 2.29 eV, respectively. Thus, the prefactor D 0 is found to be 1.27 × 10 −5 , 1.28 × 10 −4 , and 1.53 × 10 −3 m 2 /s for the activation energies of 2.0, 2.14, and 2.29 eV, respectively. Both  the activation energies and the corresponding prefactors fall in the ranges of the activation energy and the prefactor given by Huger et al. 55 As discussed above, eq 20 can be used to approximately calculate the activation energy and the prefactor for small diffusion times (small range of temperatures). Figure 7 shows the variation of the integrated and normalized intensity of the Bragg peak of the experimental data 55 with the annealing time less than 3500 s. It is evident that the integrated and normalized intensity of the Bragg peak is nearly a linear function of the annealing time. According to eq 20 and the experimental data shown in Figure 7, the term of t 2 in eq 20 is thus negligible. This result suggests that the coefficient of the term of t 2 is much less than the corresponding one for the term of t in agreement with the condition for the Taylor series expansion. Such a linear relation between the integrated and normalized intensity of the Bragg peak of the experimental data and the annealing time implies that one can use the linear relation to determine the activation energy from the ramp heating for a small annealing time as follows.
For a small annealing time, the linear relation between the integrated and normalized intensity of the Bragg peak and the annealing time can be expressed as For different starting temperatures of the annealing processes of the multilayer structure, the coefficient of the term of t is an exponential function of the starting temperature. Defining χ as the coefficient of the term of t, we have with T 1 and T 2 being the starting temperatures of two different annealing processes. From the temporal variation of the integrated and normalized intensity of the Bragg peak for small annealing times, we can find the coefficients of the term of t for the two different annealing processes and then calculate the activation energy of the diffusion process in the multilayer structure from eq 25. Using linear-regression fitting to fit the experimental data in Figure 7, we obtain a slope of −8.6 × 10 −6 s −1 , i.e., χ(593 K) = −8.6 × 10 −6 s −1 , which is comparable to −6.5 × 10 −6 s −1 obtained from eq 24 with D 0 = 19.95 × 10 −4 m 2 /s and Q = 2.20 eV. Note that the numerical values of D 0 = 19.95 × 10 −4 m 2 /s and Q = 2.20 eV fall in the ranges of D 0 = 1.55 −1.3 +18.4 × 10 −4 m 2 /s and Q = 2.29 −0.15 +0.17 eV given by Huger et al. 55 It is worth mentioning that the experimental data of I/I 0 presented in Figures 6 and 7 were obtained by subtracting the background of all the measured Bragg peaks, i.e., the experimental data of I/I 0 presented in Figures 6 and 7 do not contain any information of the first term on the right-hand side of eqs 22 and 23. According to eq 24, the information of the first term on the right-hand side of eqs 22 and 23 has no effects on the calculation of the activation energy of the diffusion process in the multilayer structure. Thus, the slope obtained from Figure 7 can be used in the calculation of the activation energy of the diffusion process from eq 24.
It should be pointed out that the above analysis is based on the fact that the multilayer experiences uniform change in temperature during the diffusion annealing under a constant heating rate. There are actually two processes involved in the diffusion annealing of the multilayer structure: one is the heat conduction and the other is the mass transport. For the heat conduction, the characteristic time, t h , is calculated as with L being the thickness of the multilayer structure and k c , ρ, and C p the effective thermal conductivity, effective density, and specific heat at constant pressure of the multilayer, respectively. For the diffusion, the characteristic time, t d , is calculated as with D eff being the effective diffusion coefficient of the diffusive species in the multilayer structure.
To have a uniform distribution of temperature under the diffusion annealing for the in situ measurement of the activation energy for the diffusion of the diffusive species, it requires (28) in which f t h and t d are the lower and upper bounds of the characteristic time for diffusion annealing in the in-situ measurement. Eq 28 provides the condition of determining the heating rate applicable to the in situ measurement of the activation energy for the diffusion of the diffusive species in a multilayer structure. For the multilayer structure used by Huger et al., 55 there are 10 bilayers of 73 Ge (165 Å)/ nat Ge (165 Å) and L = 0.33 μm. For germanium, k c = 58 W/m·K, C p = 320 J/kg·K, ρ = 5323 kg/m 3 and D = 1.20 × 10 −19 m 2 /s as calculated from Q = 2.14 eV (207.4 kJ/mol), D 0 = 1.25 × 10 −5 m 2 /s, and T f = 793 K. We have t h = 3.2 × 10 −9 s and t d = 908,010 s. For the heating rate of 1−10 K/min, the corresponding heating time per degree is 60 to 6 s, which satisfies eq 28. Thus, the use of uniform temperature in the multilayer under the heating condition used in the analysis is reasonable.

■ SUMMARY
Mass transport plays an important role in the formation and growth of new phases and new compounds, which determine the structural integrity of multilayer structures used in a variety of optoelectronic devices and solid-state energy storage. 2−4 One of the challenges in understanding mass transport on micro-and nanoscales is the determination of the activation energy responsible for the rate process controlling the mass transport. The works presented by Mizoguchi and Murata 29 and extended by Huger et al. 55 point to the feasibility of determining the activation energy for atomic diffusion in multilayer structures by using a constant heating rate during the annealing process if there are no changes in the activation energy for diffusion.
We have analytically solved the diffusion problem in a multilayer structure with periodic conditions under a constant heating rate. Correlating the coefficients of eigenfunctions with the intensity of the X-ray or neutron Bragg peak establishes the analytical expression for the determination of the activation energy from the annealing under a constant heating rate. Introducing a nominal diffusion time of eq 18, the logarithm of the temporal variation of the intensity of the X-ray peak or Bragg peak for a multilayer structure under continuous heating with a constant heating rate is found to be a linear function of the nominal diffusion time if the diffusion mechanism remains the same for the temperature range under investigation. The experimental results given by Huger et al. 55 support the linear relationship used to determine the activation energy for the diffusion in multilayer structures. From this relationship, a simple method is proposed for the determination of the activation energy of the diffusion in multilayer structures under a small annealing time.
It should be pointed out that the method presented by Mizoguchi and Murata 29 and extended by Huger et al. 55 is based on the uniform distribution of temperature in the multilayer structure. This requires that the characteristic time for heat conduction is much less than the characteristic time for the ramp heating as well as the characteristic time for diffusion, as shown in eq 28.