Simulation of reduction of oxidized metal nanoparticles

I analyze theoretically the spatio-temporal kinetics of reduction of oxidized metal nanoparticles by hydrogen (or methane). The focus is on the experimentally observed formation of metal and oxide domains separated partly by pores. The interpretation of such multiphase processes in nanoparticles at the mean-field level is hardly possible primarily due to complex geometry, and accordingly I use the lattice Monte Carlo technique in order to tackle this problem. The main conclusions drawn from the corresponding generic simulations are as follows. (i) The patterns predicted are fairly sensitive to the metal-metal and metal-oxygen interactions. With decreasing the former interaction and increasing the latter interaction, there is transition from the formation of metal aggregates and voids to the formation of a metal film around the oxide core. (ii) During the initial phase of these kinetics, the extent of reduction can roughly be described by using the power law, and the corresponding exponent is about 0.3. (iii) With decreasing the hydrogen (or methane) pressure and/or increasing the oxide nanoparticle size, as expected, the kinetics are predicted to become longer. (iv) The dependence of the patterns on the presence of the support and/or Kirkendall void in an oxide nanoparticle is shown as well.


Introduction
Metal nanopartiles (MNPs) have long been in the center of heterogeneous catalysis [1] and are now in the center of nanoscience in general with numerous potential and already realized applications (reviewed e.g. in [2,3]). In reality, MNPs tend to deteriorate. In catalysis under relatively high temperatures, it often occurs via sintering (reviewed in [4][5][6]). In other applications, at lower temperatures, the deterioration of MNPs is frequently related to their oxidation, i.e., conversion to oxide NPs (ONPs; see e.g. recent reviews [7,8], experiments [9][10][11][12], kinetic models [13,14], and references therein). The metal properties of ONPs can be restored at least partly by reduction via reaction with hydrogen (or methane) as reviewed in [15] (for recent experiments, see also e.g. [16][17][18][19]). The available mean-field kinetic models of reduction of ONPs imply the formation of a core-shell NP structure (see e.g. [20][21][22][23]). In reality, this process can be more complex. For NiO NPs, for example, the experiments indicate that the reduction is accompanied by the formation of metal and oxide domains separated partly by pores [15,24]. The description of such processes at the mean-field level is hardly possible primarily due to complex geometry. The full-scale Monte Carlo (MC) simulations are also hardly possible due to the uncertainty with the choice of numerous parameter and large values of metalmetal and metal-oxygen interactions, and the corresponding models are now lacking. Herein, I present the first generic 2D lattice MC simulations illustrating the key features of phase separation during reduction of ONPs.
From the perspective of statistical physics, the process under consideration is a special case of phase separation under reactive conditions. Other examples of processes belonging to this class are available in heterogeneous catalysis, and the corresponding models are generic as well (see e.g. MC simulations [25,26]). The specifics of phase separation during reduction of ONPs and in heterogeneous catalytic reactions are, however, fairly different.

Model
The simulations are performed on a square L × L lattice. Each site can be occupied by one monomer, A (metal), B (oxygen), or C (mimic of hydrogen), or be vacant. An ONP is represented at t = 0 by a square c(2 × 2 ) l × l A-B array (see, e.g., Fig. 1). The other sites are initially (at t = 0 ) assumed to be vacant or occupied by C with probability p < 1 (this parameter is related to hydrogen pressure). The analysis includes diffusion of A, B, and C via jumps to nearest-neighbour (nn) vacant sites. These jumps are schematically represented as where R is A, B, or C, and Z is a nn vacant site. In addition, there is irreversible reaction between reactants located in nn sites, The driving force for the oxide formation is described by introducing the attractive nn interaction, AB < 0 , between nn A and B. The tendency of metal atoms to aggregate is reflected by introducing the attractive nn A-A interaction, AA < 0 . The other interactions ( AC , BB , BC , and CC ) are neglected.
In the framework of the transition state theory, the rate constant of a monomer jump to a nn vacant site is determined by the pre-exponential factor and the activation energy identified as usual with the difference of the monomer energies in the activated state, i.e., at the saddle points of the potential barriers, and in the ground state, i.e., near the bottom of potential wells [27] (for the other dynamics, see e.g. [28]). A monomer performing a jump interacts laterally (in the 2D case) with neighbours, which are in the ground state. Depending on the location of a jumping monomer, there are lateral interactions in the ground and activated states, X,i and * X,i , where X is the subscript characterizing the monomer type ( X ≡ A , B, or C), and i is the subscript characterizing the arrangement of neighbours. The difference of these energies determines the contribution of lateral interactions to the jump activation energy. To reduce the number of parameters, the lateral interaction in the activated states are here neglected, i.e., * X,i = 0 . The lateral interaction in the ground state is reduced to the nn interactions as described above. With this specification, the rate constants of a A, B, or C jump in one of the directions are represented as , and k • C are the maximal values of the rate constants, and n and m are the numbers of neighbours corresponding to a jumping A or B monomer and given i.
For MC simulations, one needs dimensionless probabilities ( p i ≤ 1 ) of possible events. To get such probabilities, the rate constants are usually normalized to a properly chosen rate constant which is larger than or equal to the rate constants of all the possible events. After such normalization, expressions (3) can be replaced by and p • C ≤ 1 are the jump probabilities in the absence of neighbours.
Reaction (2) is considered to be fast, i.e., to occur just after a diffusion jump of B (or C) if it has one or more nn C (or B). If after a diffusion jump of B (or C) it has more than one nn C (or B), the second reactant is chosen at random.
On the lattice boundary, the no-flux boundary condition are employed for A and B, i.e., the jumps of these monomers out of the lattice are not allowed. The C supply to the lattice is mimicked by using for this species the grand canonical distribution with the prescribed average C population, p (as at t = 0 ), on the border sites during trials of C jumps from the boundary sites to the interior of the lattice or from the interior to the boundary.

Algorithm of simulations
With the specification above, the algorithm of MC simulations is as follows: where 0 < ≤ 1 is a random number.
The jumps near and at the boundary sites are simulated with the specification described in the end of Sec. 2. On average, Δt = 1 corresponds to L 2 MC trials. In the simulations presented, as usual, Δt = 1 is identified with one MC step (MCS). To convert t into real time, it should be divided by the rate constant which was used for normalization of the jump rate constants. The simulations are, however, focused on the patterns arising during reduction of ONPs, and from this perspective the time units are not important. For this reason, the time is below given in MCS.
MC runs were performed on a lattice with L = 200 up to t = 10 7 MCS. The size of the A+B array at t = 0 was l = 100 or 70. The kinetics observed were characterized by calculating the extent of reduction, where n B (t) and n B (0) ≡ l 2 ∕2 are the current and initial B populations.

Results of simulations
In the model presented, the jump probabilities defined by (4) , p • C , and the ratios AA ∕k B T < 0 and AB ∕k B T < 0 . In reality, the absolute values of these ratios are large ( ≫ 1 ), and with their realistic values MC simulations are too slow. For this reason, the use of relatively small absolute values of these ratios is usually inevitable. With this reservation, the bulk of the MC simulations shown below was performed with AA ∕k B T = −4 and AB ∕k B T = −3 . To make the simulations faster, the diffusion jumps were run with p To illustrate the predicted reduction kinetics and the corresponding patterns, I present 6 MC runs with different sets of parameters or initial conditions. Run 1 was performed with L = 200 , l = 100 , p = 0.1 , AA ∕k B T = −4 , and AB ∕k B T = −3 . Initially, the A-B array was located in the center (Fig. 1). In Runs 2-6, one or two of these parameters or the initial conditions are changed in order to show the sensitivity of the results with respect to such changes. In particular, the initial arrangements were without vacant sites inside (as in Fig. 1) in Runs 1-5 and with an array of vacant sites in Run 6.
The reduction kinetics corresponding to Runs 1-5 are exhibited in Fig. 2. During the initial phase of these kinetics, the extent of reduction defined by (5) can roughly be described by using the power law, with ≃ 0.3 . This exponent is close to the Lifshitz-Slyozov exponent, LS = 1∕3 , for Ostwald ripening (see e.g. [29,30]). This coincidence is natural because the reduction occurs in parallel and is influenced by the growth of A aggregates. More specifically, the analysis of the kinetics shown in Fig. 2 indicates that = 0.32 for Run 1, 0.21 for Run 2, 0.51 for Run 3, 0.31 for Run 4, and 0.33 for Run 5. The Typical patterns observed during MC simulations with the first set of the parameters (Run 1) are shown in Fig. 3. The A aggregates and voids are seen to be rather small up to t = 10 4 MCS [Fig. 3a]. With increasing time, the A aggregates and voids grow and at t = 10 7 MCS [ Fig. 3d] their size becomes to be comparable with the initial size of the A-B array.
Regarding the variation of parameters, it is of interest to notice that the results of the simulations are fairly sensitive to the values of the ratios AA ∕k B T and AB ∕k B T . This is illustrated by using the second set of the parameters with AA ∕k B T = −3 and AB ∕k B T = −4 (Run 2, Fig. 4). With decreasing | AA ∕k B T| and increasing | AB ∕k B T| (compared to that in the first set of the parameters), the driving force for the A segregation becomes weaker, i.e., the relative stability of the A-B phase increases, and the model predicts formation of the A layer around the A-B core (Fig. 4), and the whole kinetics becomes slower compared to that  (Fig. 2). One of the manifestation of this A-layer-related slowdown is that is reduced down to 0.21.
The main goal of the work was to show the A segregation and void formation, and the first set of parameters (Run 1) was selected from this perspective. The other simulations presented below were performed with the same values of the ratios AA ∕k B T and AB ∕k B T as in the first set of parameters. In particular, the third set of the parameters (Run 3) was chosen as the first one except the p value which was reduced from 0.1 to 0.02. In reality, this physically corresponds to decrease of the hydrogen pressure. With this modification, as expected, the kinetics becomes slower than during Run 1 (Fig. 2). The patterns are, however, qualitatively similar to those observed in Run 1 (cf. Figure 3 and S1 in Supplementary  Information), i.e., there is no formation of the coherent A film in the external area of the array. For this reason, this slowdown of the kinetics does not reduce (as in the case of Run 2). In fact, increases up to 0.51, because B particles have more time to diffuse from the internal area to the external area and react there.
In the fourth set of the parameters (Run 4), l was reduced from 100 to 70. With this modification, the kinetics becomes shorter (see Fig. 2 and cf. Figures 3 and  S2 in Supplementary Information). The fifth set of the parameters (Run 5) was identical to the first one. The difference was in the initial conditions with one of the side of the A-B array contacting one of the lattice boundaries and prohibition of C supply at this boundary. These boundary conditions are aimed to simulate an oxide nanoparticle located at the support. In this case, the kinetics is somewhat slower than during Run 1 (Fig. 2), and the reduction occurs primarily in the upper part of the array (Fig. 5).
The sixth set of the parameters (Run 6) was identical to the first one as well, and the difference was also in the initial conditions. Here, the initial 100 × 100 A-B array had inside a circular region of vacant sites. The radius of this region was 30. These initial conditions are aimed to simulate an ONP with a void inside. Such ONPs are often observed in experiments after oxidation of MNPs and associated either with the specifics of diffusion of metal atoms and oxygen (Kirkendall effect [31,32]) and/or induction of tensile strain [33,34] during oxidation. With this setup, the initial phase of the kinetics and the corresponding patterns are similar to those predicted with the first set of the parameters (cf. Figures 3 and  S3 in Supplementary Information). After penetration of C to the central region, the kinetics becomes, however, faster, and there is difference in the patterns.

Conclusion
The MC simulations performed in the framework of the proposed generic 2D lattice model on the times up to 10 7 MCS illustrate the specifics of the formation of metal aggregates and voids during reduction of OMNs by hydrogen or methane. This duration of the MC runs is sufficient in order to reach appreciable extent of reduction. The conclusions drawn from the simulations are as follows: (i) As one could expect physically, the patterns observed in the simulations are fairly sensitive to the A-A (metal-metal) and A-B (metal-oxygen) interactions.
With decreasing the former interaction and increasing the latter interaction, the model predicts the transition from the formation of metal aggregates and voids to the formation of a metal film around the oxide core. This effect was not obvious in advance and accordingly can be classified as novel. (ii) During the initial phase of the kinetics under consideration, the extent of reduction can roughly be described by using the power law with ≃ 0.3 . This finding is novel as well. (iii) With decreasing p (hydrogen or methane pressure) and/or increasing the ONP size, the kinetics become longer. (iv) The patterns depend on the presence of the support and/or Kirkendall void in an ONM.
All these conclusions obtained in the framework of the 2D lattice model appear to be physically reasonable, and one can add that two of them, (iii) and (iv), are obvious and applicable in the 3D case as well. The other two, (i) and (ii), are, however, not trivial and merit additional more specific remarks. In particular, conclusion (i) is qualitative, could be physically expected, and likely holds in the 3D case. Conclusion (ii) is quantitative and physically supported by referring to the specifics of Ostwald ripening and the corresponding Lifshitz-Slyozov exponent which is applicable in the 2D and 3D cases. With this reservation, one cannot, however, exclude that in the 3D case the exponent for the ONP reduction will be somewhat different compared to the obtained one, ≃ 0.3 . In fact, the difference has been observed in the simulations presented ( = 0.21 for Run 2, and = 0.51 for Run 3).
Taken together, the results of the simulations are instructive. In particular, one can see that the qualitative interpretation of the specifics of the formation of metal aggregates and voids during reduction of ONPs by hydrogen or methane is possible at the generic level.
Finally, it is worth to notice that the model presented can be extended in various directions. One of the extensions is to choose the parameters so that the patterns become to be closer to those implied in the mean-field treatments mentioned in the Introduction. Another and perhaps more interesting extension is to perform simulations in the 3D case.
Author Contributions All the work was done by the author. 1 3