Diffusion-limited reactions on a two-dimensional lattice with binary disorder
aa r X i v : . [ c ond - m a t . d i s - nn ] N ov Diffusion-limited reactions on a two-dimensional lattice with binary disorder
Andrea Wolff, ∗ Ingo Lohmar, † and Joachim Krug ‡ Institute for Theoretical Physics, University of Cologne, Z¨ulpicher Strasse 77, 50937 K¨oln, Germany
Yechiel Frank § and Ofer Biham ¶ Racah Institute of Physics, The Hebrew University, Jerusalem 91904, Israel (Dated: September 2, 2018)Reaction-diffusion systems where transition rates exhibit quenched disorder are common in physical andchemical systems. We study pair reactions on a periodic two-dimensional lattice, including continuous deposi-tion and spontaneous desorption of particles. Hopping and desorption are taken to be thermally activated pro-cesses. The activation energies are drawn from a binary distribution of well depths, corresponding to ‘shallow’and ‘deep’ sites. This is the simplest non-trivial distribution, which we use to examine and explain fundamentalfeatures of the system. We simulate the system using kinetic Monte Carlo methods and provide a thorough un-derstanding of our findings. We show that the combination of shallow and deep sites broadens the temperaturewindow in which the reaction is efficient, compared to either homogeneous system. We also examine the roleof spatial correlations, including systems where one type of site is arranged in a cluster or a sublattice. Finally,we show that a simple rate equation model reproduces simulation results with very good accuracy.
PACS numbers: 98.38.Bn, 68.43.-h, 98.38.Cp
I. INTRODUCTION
Reaction-diffusion systems are successful models to de-scribe a large variety of phenomena in physics, chemistry, andbiology [1, 2]. They may involve one or more reactant speciesthat diffuse and react with each other on a surface or in thebulk. In particular, surfaces often catalyze chemical reactionsbetween adsorbed atoms and molecules. The densities of theadsorbed chemical species and their reaction rates depend onparameters of the surface and on the temperature. Microscop-ically, one can describe the diffusion of particles on the sur-face as a random walk between adsorption sites. In homoge-neous systems all the adsorption sites are identical. However,most systems are heterogeneous, involving different types ofadsorption sites with a broad distribution of binding energies.The present study is motivated by a specific example of animportant surface process, namely the formation of molecularhydrogen on dust grains in the interstellar medium [3–8]. Hy-drogen atoms impinge from the gas phase onto a grain, anddiffuse on its surface. They may either desorb thermally fromthe surface, or encounter each other and form a molecule. Thisdefines a reaction-diffusion system in a spatially confined re-gion. In this article we are concerned with steady-state sys-tems, when the hydrogen recombination efficiency is definedas the fraction of impinging particles that end up (and eventu-ally desorb) in molecular form. This efficiency plays an im-portant role in the evolution of interstellar clouds. Typically,there is a narrow window of temperatures in which recombi-nation is efficient. At lower temperatures, the atoms are not ∗ Electronic address: [email protected] † Electronic address: [email protected]; now at: Racah Institute of Physics,The Hebrew University, Jerusalem 91904, Israel ‡ Electronic address: [email protected] § Electronic address: [email protected] ¶ Electronic address: [email protected] sufficiently mobile to react, whereas at higher temperaturesthey desorb too quickly.Assuming that all rates are spatially homogeneous, the sys-tem is well understood analytically. A zero-dimensional mas-ter equation for the particle number distribution [9, 10], to-gether with a proper definition and calculation of the reactionrate coefficient in terms of a first-passage problem [11, 12],suffices to accurately describe the many-particle system [13].It is very important, however, to consider disorder in the localrates of hopping and desorption of the particles. As we alludedto earlier, this is not only of theoretical interest. In fact, in theastrophysical context, the disordered case is much more real-istic, and it is long known that disorder potentially enhancesthe efficiency dramatically [5]. However, the combination of aconfined two-dimensional region, rate disorder and the many-particle reaction-diffusion dynamics makes this problem no-toriously hard to tackle analytically. Kinetic Monte Carlo(KMC) methods can be used to simulate such systems [e.g.14, 15], and algorithms are still subject to improvement [16,which also compares related approaches]. They remain com-putationally expensive, however, and a systematic understand-ing of the effects of disorder is still missing.Here we start such an analysis for the simplest form of ratedisorder, where each lattice position corresponds to either astandard (‘shallow’) site, or to a strong-binding (‘deep’) site,with enhanced binding energy. While we strive to keep thisa theoretical self-contained work, our models and questionsare motivated by applications and should easily translate topractice. This is one reason why we have chosen thermallyactivated rates throughout, and present most results in termsof temperature, and on scales relevant to the astrophysicalproblem just described. In the latter context, our work is rel-evant to systems combining physisorption (shallow sites) andchemisorption (deep sites) [17] , aside from features particularto specific material systems. Using such a discrete distributionturns out to be conceptually different from the case of con-tinuous distributions of binding strength [14], in which welldepths drawn from tails of the distribution may significantlyaffect the temperature window of high efficiency.Our goal in this paper is to provide a thorough under-standing of all relevant mechanisms of the described reaction-diffusion system. Most importantly, if we start from homo-geneous systems of either standard or deep sites, their tem-perature windows of high efficiency will typically be sepa-rated by a gap. It is a natural question whether a mixture ofthe two types of sites still exhibits two separated peaks, orwhether (and under what conditions) the efficiency is high forin-between temperatures.Our findings are relevant for other systems as diverse ascatalysts [18], exciton trapping in photosynthesis [19], excitontransport in semiconducting nanosystems [20], and diffusion-limited reactions on biomembranes [21]. The generalizationof our results to these and other related contexts should bestraightforward.The paper is organized as follows. In Sec. II we define thesystem and our notation. The following Sec. III provides aqualitative picture which identifies three temperature regimesand describes the relevant processes in each. In Sec. IV wegive a systematic account of extensive KMC simulations anddiscuss the observed behavior in detail. This includes thestudy of spatial correlations in the quenched disorder. Sec-tion V presents a simple yet accurate rate equation model, andwe explain the difference to the homogeneous case. We derivean expression for the efficiency in the most interesting regime.Finally, we present our conclusions in Sec. VI.
II. MODEL AND DEFINITIONS
We consider a system of a single particle species on a two-dimensional square lattice of S sites with periodic boundaryconditions. Each lattice site is characterized by a binding en-ergy, which can take one of two values — we call this a bi-nary lattice . The number of sites of either type is denotedby S i ( i = , S = S + S . Particles impinge onto thelattice at a homogeneous rate f per site. If a site is alreadyoccupied, the impinging particle is rejected. In the contextof surface chemistry this is known as Langmuir-Hinshelwood(LH) rejection [22].Particles explore the lattice by hopping to neighboring siteswith an (undirected) rate a , and they can desorb from a sitewith rate W . Both rates depend on the binding energy at theparticle position. If two particles meet on one site, they forma dimer and leave the system immediately. The key quantityof such a system is the efficiency h , defined as the ratio be-tween the number of particles that react and the total numberof impinging particles, when the system is in a steady state.In view of possible applications, we choose rates to be ther-mally activated by a system temperature T . The activationenergy for desorption is denoted E W i , which we identify withthe binding energy at the particle position. Similarly, hoppingfrom a type- i site has an activation energy E a i . All rates sharethe attempt frequency n , so that, e.g., W i = n exp ( − E W i / T ) —here and in the following energies are measured in tempera-ture units. We want to ensure detailed balance . The simplest E W i − E a i = const E a E a E W E W ∆ E FIG. 1: One-dimensional cut through the energy landscape of ourmodel. way to achieve this is by choosing W / a = W / a , or equiva-lently, E W − E a = E W − E a , and we will employ this choicethroughout. The number of sites visited by a single particlebefore desorption becomes then independent of disorder.To establish a connection to surface chemistry problems,we think of type-1 sites as standard or ‘shallow’ adsorptionsites, and of type-2 sites as strong-binding or ‘deep’ sites, with E W > E W . A one-dimensional cut through such an energylandscape is sketched in Fig. 1. III. QUALITATIVE DISCUSSIONA. Homogeneous systems
For a homogeneous lattice, the dependence of the efficiencyon the temperature, h ( T ) , is known [9, 10]. At low tempera-tures the particles are nearly immobile, thus they do not meetother particles. Therefore, the lattice is highly occupied, in-coming particles are mostly rejected, and the efficiency is low.For higher temperatures, hopping processes are activated andthe particles begin to explore the lattice. This leads to morefrequent encounters, so the efficiency rises. When the tem-perature is increased even further, the particles tend to desorbbefore encountering each other, the lattice coverage becomessmall, and the efficiency decreases.For the homogeneous system, the corresponding tempera-ture bounds have been obtained using rate equations [10]: T up = E W − E a ln ( n / f ) (1)is the temperature above which the kinetics becomes second-(instead of first-) order, whence desorption ends the typicalparticle residence and the efficiency is low, and T low = E a ln ( n / f ) (2)is the temperature below which particles arrive faster than theyhop, leading to dominant LH rejection and low efficiency. Theaverage of these two bounds reads T max = E W ln ( n / f ) (3)and corresponds to the temperature of maximum efficiency.If the binding energy E W is increased, the efficiency max-imum is shifted towards higher temperatures, and the shift isdirectly proportional to the change in binding energy. For abinding energy difference D E = E W − E W of two (otherwiseequal) lattices of either type-1 or type-2 sites, the relation be-tween the temperatures of maximal efficiency is given by T max2 = T max1 · (cid:18) + D EE W (cid:19) , (4)whereas the peak width T up i − T low i is the same. B. Binary systems
Now we consider the binary lattice introduced in Sec. II,with binding energies E W and E W . To each site, we randomlyassign a binding energy. There is a typical length for a parti-cle to find a strong-binding site. This length obviously short-ens when there are more and more of these sites on the lat-tice. At low temperatures around the efficiency maximum ofthe type-1 sites, particles can only diffuse on and desorb fromthese shallow sites, while particles landing on or hopping ontostrong-binding sites cannot leave by hopping or desorption,since the binding energy is too high. Recombinations eithertake place on the shallow sites, or by hopping to an occupiedneighboring strong-binding site. For very high temperaturesaround the efficiency maximum of the deep wells, the particlesdiffuse on, desorb from and recombine on those, while on theshallow sites, they desorb too quickly to allow any other pro-cesses. But in the intermediate temperature regime — rightof the shallow peak, left of the strong-binding peak — some-thing different happens. Here, the temperature is too low fordynamics on deep wells, so particles encountering a deep wellare stuck. On the other hand, particles on shallow sites tend todesorb rather quickly, and thus do not recombine on such sites.But if they find a deep well before desorbing, they are trappeduntil another adatom shares their fate and they recombine.The simple random walk with traps has been studied exten-sively [e.g. 19, 23]. To leading order, the average number ofsteps a random walker performs before trapping is given by h n i ≈ p S S ln S , (5)where S is the number of deep wells and S is the total numberof sites on the lattice. This leads to a trapping length ℓ trap = p h n i . (6)On the other hand, the typical radius of the area a walker ex-plores on standard sites before desorption is the random walklength [12] ℓ rw = r a W . (7)Trapping now competes with desorption from shallow sites;the former only depends on the number of traps S , while the latter is a function of temperature. As long as the randomwalk length ℓ rw is larger than the trapping length ℓ trap , the par-ticles are — on average — trapped before they can leave thelattice. For a given number of traps S this implies a high ef-ficiency approximately up to the temperature T eq where bothlengths become equal. If this temperature lies above the inter-mediate temperature range where both pure systems have poorefficiency, we can expect a high efficiency throughout, hencea full ‘bridging’ of the gap. Since the efficiency is high overthis whole temperature range then, we call this an efficiency plateau . We will calculate the value of the efficiency on sucha plateau in a rate equation model in Sec. V C.Summing up, we can divide the temperature axis into threeregions. The lowest temperatures where only particles onshallow sites are mobile, the intermediate regime where theparticles behave like random walkers on a lattice with traps,and the high temperatures where particles become mobile onstrong-binding sites. We now check this qualitative picturewith KMC simulations. IV. KINETIC MONTE CARLO SIMULATIONSA. Setup
In order to test our predictions, we carried out extensivekinetic Monte Carlo simulations. The standard algorithm pro-ceeds as follows [cf. 24, for a review]. We keep track of thefull microscopic dynamics of continuous-time random walk-ers [25] with standard exponential waiting time distributions.In each simulation step, the current system configuration de-termines the list of possible elementary processes and theirrates. By comparing a random number with the normalizedpartial sums of these rates we find the process to execute next.The simulation time is then advanced according to the totalsum of rates and the configuration is updated.For a given realization, we wait for the system to reach thesteady state before we measure the efficiency over 10 im-pingements. We use a square lattice of S = ×
100 sites.We choose the other model parameters inspired by an exem-plary system in the astrophysical application, to show the rel-evance of our work in this field, and since the correspondingsystem is known to exhibit interesting kinetic regimes. Theflux of hydrogen atoms per unit surface area depends on gasdensity and temperature. The flux per surface site is givenby the ratio between the flux per unit area and the densityof adsorption sites, hence it depends on the surface morphol-ogy. More precisely, the flux per site is given by f = r v / ( s ) ,where r is the density of hydrogen atoms in the gas phase, v is their average thermal velocity and s is the density of ad-sorption sites on the surface. To obtain typical values we use r =
10 cm − , v = . × cm / s (which corresponds to agas temperature of 100 K) and s = × cm − which is themeasured density of adsorption sites on the amorphous car-bon sample studied in Ref. 26. This results in a flux per site of f = . × − s − . For the attempt frequency we choose thestandard value of 10 s − which is commonly used through-out surface science. With each site we associate either thestandard binding energy E W =
658 K, as found for hydro-gen atoms on amorphous carbon [27], or an enhanced energy E W = E W + D E with D E = E a =
511 K or E a = E a + D E ,respectively.In each case, we determine the efficiency as a function ofthe temperature T , as well as of the relative frequency ofstrong-binding sites S / S . We do this for up to four differ-ent ways of distributing the binding strengths. For dynam-ics with nearest-neighbor hopping of the particles, we either randomly assign to each site a binding energy with proba-bilities p and p = − p , respectively, or we arrange thestrong-binding sites in a regular sublattice , or we concentrateall strong-binding sites in a single square cluster . In the caseof random assignment, S is then binomially distributed withparameter p . To eliminate the fluctuations in S , we averagethe efficiency over 20 realizations. In the following discussion,we can therefore identify S i / S with its average p i . For compar-ison with the rate equation model to be introduced in Sec. V,we also implement another kind of dynamics ( ‘longhop’ case ),namely hopping from any site to any other site of the lattice.This switches off any spatial correlations between the latticesites and thus is best suited for comparison with an effectivezero-dimensional model.Binary disorder models very similar to random assignmentand the clustered case have been simulated before [14]. Theauthors were predominantly concerned with showing thatsuch models can exhibit efficient reaction over a broader rangeof temperatures than homogeneous systems. Here we extendthese findings to a systematic picture for the effect of the deep-site fraction p and the energy gap D E . More importantly, weprovide detailed explanations and analytic results which ex-plain all notable features of the simulation outcome in termsof microscopic physical processes. B. Results
Figure 2 shows the results of our simulations. For each D E ,we simulated systems with 1, 4, 25 and 50% of strong-bindingsites. The random distribution is probably the most interestingregarding applications. Following the series of Figs. for each D E , we observe that the intermediate temperature regime isbridged in each case. This is in accordance with the analyticprediction of Sec. III B, since already for moderate deep-sitefraction, the trapping length ℓ trap is smaller than the randomwalk length ℓ rw for all intermediate temperatures. The obser-vation holds at least up to D E = D E that we have considered; beyond thisenergy scale one enters the regime of chemisorption, which isnot our focus in this work. The bigger the difference of thebinding energies, the more strong-binding sites are needed toform a genuine plateau, where the efficiency does not dependon the temperature. This complies with the ideas of Sec. III;when the deep-site peak is shifted to higher temperatures, T eq has to increase to warrant formation of a plateau. This isachieved by increasing the deep-site fraction. The variance ofthe efficiency between different realizations of random land- scapes was found to be negligible throughout. We also exam-ined the longhop case on such landscapes, and found that theefficiency varies just as much. Since this cannot be affectedby any spatial correlations, we conclude that this variation isalways due to fluctuations in the number of deep sites S only.The arrangement of strong-binding sites in a sublattice per-forms slightly better, compared to random assignment. Thisis not astonishing since the sublattice optimizes the distancebetween the traps. In the random case, small clusters of strong-binding sites can occur, in which a single trap is less efficient.An alternative picture is that the capture zones of individualtraps typically have an overlap, which is minimized in the sub-lattice case.On a lattice with a single square cluster of strong-bindingsites, there is no bridging effect for any energy differenceor frequency of strong-binding sites. For high frequenciesof either shallow or strong-binding sites, only one restrictedpeak emerges, while for intermediate frequencies of strongand shallow sites two nearly separated peaks appear. The ef-ficiency does not drop to zero in the intermediate temperatureregime, because an exchange between shallow and deep sitestakes place along the boundary of the cluster. However, sincethe boundary length scales as √ S , the fraction of boundarysites decreases with increasing S , and correspondingly the sup-pression of the efficiency in this regime becomes even morepronounced for larger systems. We checked this for a systemof 500 ×
500 sites (not shown). This is in contrast to the well-mixed case, where a finite fraction of sites are boundary sites(see below).For the longhop case, we first verified that results on a sub-lattice and a cluster landscape coincide, ensuring the correct-ness of the algorithm. The efficiency for this kind of dynamicsoutperforms even the sublattice results for nearest-neighborhopping. This is because in the sublattice case, there is stillthe necessity for a particle to actually travel to a trap insteadof having a non-zero probability to reach a trap on every step.A further analysis of this model is provided in Sec. V.In addition to our qualitative explanations, we numericallyexamine the dependence of the plateau efficiency value on thenumber of deep wells. First we note that for our choice ofparameters, the efficiency value at T max1 ≈
14 K always corre-sponds to the plateau value. From the results for D E =
750 Kshown in Fig. 3 we infer that the way of distributing the strong-binding sites is of crucial importance. In the case of a sin-gle square cluster of deep wells, the efficiency decreases lin-early as 1 − S / S , while for the random distribution the ef-ficiency first decreases more slowly (for less than 50% ofstrong-binding sites) and faster to the end (more than 50%).We propose that this effect is related to the border length be-tween shallow and deep sites, and use this connection to de-rive an empirical formula for the plateau efficiency. For ran-domly distributed deep wells, we calculate the border length L as function of S / S (cf. Fig. 4). We find a shallow site next toa deep site with probability ( S / S )( − S / S ) . Since the orien-tation of the pair does not matter, we gain an additional factorof 2. Furthermore we have 2 S possibilities to place such a pairof sites on a square lattice with S sites and periodic boundaryconditions. So we find the following expression for the border
10 12 14 16 18 20 22 240 . . . . . . T / K η Hom. systems10 12 14 16 18 20 22 240 . . . . . . T / K η . . . . . . T / K η . . . . . . T / K η . . . . . . T / K η
50% 10 15 20 25 30 350 . . . . . . T / K η Hom. systems10 15 20 25 30 350 . . . . . . T / K η . . . . . . T / K η . . . . . . T / K η . . . . . . T / K η
50% 10 20 30 40 500 . . . . . . T / K η Hom. systems10 20 30 40 500 . . . . . . T / K η . . . . . . T / K η . . . . . . T / K η . . . . . . T / K η FIG. 2: (Color online) Efficiency versus temperature for various fractions of deep sites. Left column D E =
250 K, middle D E =
750 K, right D E = D E =
750 K: sublattice(orange dashed line, crosses), and cluster (black line, squares). Vertical green line at T eq . The first row shows the results for homogeneoussystems of only standard or only deep sites, respectively. length between shallow and deep sites L = S · S S (cid:18) − S S (cid:19) . (8)Fitting the efficiency difference D h = h random − h cluster to a multiple of this border length yields D h = C · L , (9)with C = ( . ± . ) × − or h random ≈ (cid:18) − S S (cid:19) · (cid:18) + ( . ± . ) S S (cid:19) (10) . . . . . . . . . . . . . . . S / S η T =
14 K, ∆ E =
750 K
FIG. 3: (Color online) Efficiency as function of the deep-site frac-tion for T =
14 K and D E =
750 K, for clustered deep sites (black,squares, h cluster ) and randomly assigned energies (blue, diamonds, h random ). . . . . . . . . . . . . . . . S / SL · − T =
14 K, ∆ E =
750 K 0 . . . . . . . . ∆ η FIG. 4: (Color online) Border length L (green dashed line, left axis)and efficiency difference D h (dark red, right axis), as function ofdeep-site fraction, for T =
14 K and D E =
750 K. Vertical axis scal-ing taken from data fit. for the empirical plateau efficiency value. The quality of thefit
D h (cid:181) L for KMC results underlines the role of the borderlength, and this corroborates our picture that the dominant re-action process on the plateau is by hopping from standard todeep sites. Further insight into the origin of Eq. (10) will beprovided below in Sec. V C. V. RATE EQUATION MODEL
Rate equations have been used previously to study reactionson finite surfaces with different types of sites. Using surfaceswith varying roughness, where binding energies at a site aregiven by a vertical bond strength plus an additional lateralbond strength per in-layer neighbor of the landscape, KMCsimulations were performed [15]. A rate equation model wasthen used to check that such a rough landscape model is con-sistent with surfaces deemed astrophysically relevant and ex-amined in the laboratory. To this end, the rate equations with standard energy parameters were time-integrated to predictthe results of TPD experiments.Here we apply a rate equation model to quantitatively repro-duce our KMC findings as well as to further our qualitative un-derstanding of the system’s behavior. From the definition ofSec. II we derive a set of rate equations for the total number N i of particles on sites of type i . Terms for desorption and influxare easily written down, whereas the reaction terms are moresubtle, partly since in general, the reaction is not an elemen-tary process with given rate. For the time being, we denotethe appropriate rate coefficients as A i , and refer to Sec. V Afor details. The rate equations then take the form [15, 28]N. t. = f ( S − N ) − W N − A N ( S − N ) − A N N − A N + A N ( S − N ) − A N N , N. t. = f ( S − N ) − W N − A N ( S − N ) − A N N − A N + A N ( S − N ) − A N N . (11)Here the first two contributions cater for the impingement fluxwith rejection and the desorption of particles. For clarity weseparated the remaining terms. The next two terms describeleaving to a site of the opposite type (either to an empty orto an occupied site). Then we account for reactions insideone population due to hops between sites of the same type, re-moving two atoms. The remaining two contributions describegaining a particle by a hop from the other site type, and finally,losing one particle due to the reaction with a particle comingfrom the other population.It is tempting to substitute the ‘internal’ reaction term2 A i N i by 2 A i N i ( N i − ) , since it should really depend on thenumber of pairs . This is not adequate: In the rate equationtreatment the N i are continuous and can drop below unity, suchthat the reaction term (which we are ultimately interested in)could then become negative. The assumption that the reactionrate can be written as above is at the heart of the rate equationapproach (“mass action law”). Equations (11) are easily de-rived from the full master equation using this assumption inthe forms h N i ( N i − ) i ≈ N i and h N N i ≈ N N (where theexpectation is over the joint probability distribution P ( N , N ) and the r.h.s. N i ’s are already the mean values as above).The reaction terms also provide the recombination rate ofthe process. Adding up all terms proportional to the A i inN. / t. = N. / t. + N. / t., mere hopping terms (not leading to a re-action) cancel. Using that the reaction consumes two particles,we obtain the rate at which particles are removed by the reac-tion as 2 R = A N + A N + ( A + A ) N N , (12)which can be simplified to 2 R = ( A N + A N )( N + N ) .Relating this to the particle influx f ( S + S ) = f S gives the efficiency h = R / ( f S ) . A. The reaction rate coefficient
The homogeneous system was treated analytically by rateequations [10, 27, 29], the master equation [9, 10], and mo-ment equations [30, 31]. For these methods just as for stochas-tic or numerical methods based on these approaches, the re-action rate coefficient is a crucial quantity, typically approxi-mated as A ≈ a / S [32]. We have argued elsewhere that thisneglects the nature of two-dimensional diffusion (“back diffu-sion”) as well as the fundamental first-passage problem, thecompetition between a meeting (hence reaction) of particlesand the prior desorption of a reactant. Hence we put someeffort into a proper definition and evaluation of A [11, 12],and we claimed that these results should be applied in all men-tioned frameworks, including the rate equation treatment [13].Here, we return to the choice A i = a i / S , since the situa-tion is different. In rate equations such as Eqs. (11), there isno way to genuinely incorporate any spatial structure. Thisholds true for all zero-dimensional approaches, e.g., the mas-ter equation as well. However, in the homogeneous systemsstudied before, this neglect only concerns the spatial correla-tions in the particle residence probability, with well-studiedeffects [10, 11, 33]. In the heterogeneous system with its sep-arated populations, this approach additionally neglects sitetype correlations.For consistency, we are then forced to assume that a particlecan reach any other site by a single hop. In particular, it hopsto a site of type i with probability S i / S , and it meets a particleon an i -site with probability N i / S . The conventional choice A i ≈ a i / S thus arises naturally if we use rate equations to de-scribe a system with site disorder, and we adopt this choice inthe following. For a system with quenched spatial structureand nearest-neighbor hops only, this description correspondsmost closely to the well-mixed case. B. Comparison with KMC simulations
The rate equations (11) are exactly solvable at steady stateby finding the real positive root of a third-order polynomial.However, the results are cumbersome and less than illumi-nating. We therefore directly opted for a numerical solverthroughout.We find that the rate equations for the binary system repro-duce the outcome of extensive KMC (longhop) simulationsfor a wide parameter range of practical relevance to excellentaccuracy (see Fig. 5). As noted in prior work [11], however,since we present our results as functions of temperature andparameters are thermally activated, we typically have rathersteep rises or declines, when even factors of two or three inthe efficiency need not appear substantial. This hardly ex-plains the overall accuracy, especially on plateaus and mod-erate peaks for h considerably smaller than unity.Results on the validity of rate equations to describe themodel in the homogeneous case have shown that confinementto a finite surface renders the discreteness of particles andfluctuations in the particle number important [10, 11, 33–36].Consequently, the mean-field approach of rate equations con- siderably overestimates the recombination efficiency in smallsystems. We do not see such effects for several reasons. Weare interested in the behavior of the system with a substantialnumber of particles, when the effects of discreteness and offluctuations in this particle number are strongly reduced. Fur-ther, the confinement of particles to a finite surface is alsofar less important than for the homogeneous system, becausethe majority of these particles is trapped in deep wells in theregimes of most interest, anyway. Finally, our system can-not be considered small, and we cannot preclude completelythat differences might be more pronounced for smaller systemsizes or different activation energies. C. Plateau efficiency
A key question of this work concerns the bridging betweenthe two efficiency peaks corresponding to homogeneous sys-tems of one type of sites. We have found a convincing quali-tative picture before in Sec. III, and we observe an efficiencyplateau between the two (virtual) peaks for a wide range ofconditions. What is the value of the efficiency along thisplateau?We need to further simplify our rate equation model to ar-rive at a simple analytic answer. Several efforts to derivethis from first principles were not met with success, hence westart from some observations: Figure 5 shows that whenevera plateau emerges in the efficiency, practically all recombina-tions are due to hops between the two types of sites. Onlyfor very low concentrations of deep wells and when the ef-ficiency is close to unity, the lower temperature end of theplateau also includes a substantial contribution from recom-bination on standard sites. On the high-temperature end, anysizable contribution from reactions on deep sites already re-sults in an efficiency peak atop the plateau value anyway.For the plateau efficiency, we can therefore capture theessence of the model accounting only for reactions between the two populations. Type 1 denotes the standard sites, sowe clearly have A ≫ A . We use this to neglect all termsproportional to A , since they are small compared to their A counterparts, but keep all flux and desorption terms. Our rea-soning is to retain as many terms as possible, to remove thosefor recombination inside the N i populations, and since we canneglect the reaction A N N ≪ A N N , we have to leave outthe corresponding A hopping term for consistency as well.This leads to the simplified steady-state equations0 = f ( S − N ) − W N − A N S , = f ( S − N ) − W N + A N S − A N N , (13)which yield an efficiency h p = A N N f S = f A S S ( V + A S ) S ( V + A S )[ V ( V + A S ) + f A S ] , (14)where V i = W i + f . We could now evaluate this at a tempera-ture right on the plateau. It will turn out, however, that we canmake two more assumptions for this case.
10 15 20 250 . . . . . . T / K η Hom. systems10 15 20 250 . . . . . . T / K η . . . . . . T / K η . . . . . . T / K η . . . . . . T / K η
50% 5 10 15 20 25 30 35 400 . . . . . . T / K η Hom. systems5 10 15 20 25 30 35 400 . . . . . . T / K η . . . . . . T / K η . . . . . . T / K η . . . . . . T / K η
50% 10 20 30 40 50 600 . . . . . . T / K η Hom. systems10 20 30 40 50 600 . . . . . . T / K η . . . . . . T / K η . . . . . . T / K η . . . . . . T / K η FIG. 5: (Color online) Efficiency versus temperature for various fractions of deep sites. Left column D E =
250 K, middle D E =
750 K, right D E = A i = a i / S (solid), blue lines contributions by reaction on the 1- and 2-sites (dashed), and by switching between the types (dot-dashed). Dotted greenline the results of the plateau model, and green horizontal line the simple (16). Vertical green line at T eq . The first row shows the results forhomogeneous systems of only standard or only deep sites, respectively. First, we also neglect desorption from the 2-sites, so V = f ,and using A = a / S , Eq. (14) reduces to h p = ( S / S )( S / S )( + V / a )( V / a + S / S )( + V / a + S / S ) . (15)Second, on the plateau and for a reasonable deep-site fraction S / S , we have V / a ≪ S / S <
1. This yields h p ≈ S / S + , (16)which no longer depends on any energy scales, and which wefind to be in excellent agreement with both KMC and full rateequation results (Fig. 5): Whenever a plateau forms (i.e., if D E is large enough to separate the homogeneous-system peaks,and if there are enough deep wells if D E is fairly large), theabove expression is valid.To check the validity of these approximations, we recallfrom Sec. IV B that the peak temperature T max1 for standard-site parameters was found to always belong to the plateau.It is large enough not to lie on the low-temperature rise tothe standard-site peak, yet minimal so as not to depend onthe peak separation governed by D E . At T = T max1 , V = W + f = f and V = W + f = f [( f / n ) D E / E W + ] . Reason-ably, f / n ≪
1, while the smallest interesting D E ∼ E W − E a ,such that the ratio D E / E W is not excessively smaller thanunity. This justifies the approximation V ≈ f , immediatelyeliminating D E from the game, as suggested by Fig. 5. Wenow check the order of V / a = f / a = ( f / n ) ( E W − E a ) / E W .The exponent is about 0 .
22 for amorphous carbon, and withthe corresponding standard flux we have V / a ≈ . × − (cf. Sec. IV A). This is negligible compared to any interestingdeep-well fraction S / S , which completes the argument forEq. (16). (We checked that this holds at least equally well forstandard olivine parameters [27].)Knowing what terms can be neglected, this result is also eas-ily derived from further simplified rate equations. We ratherprovide an intuitive explanation. We consider the system inthe steady state, so all particles entering the system also haveto leave. They enter by impingement to any site, and leaveonly from 2-sites, by LH rejection or by reaction with an in-coming 1-particle. Consequently, particles from 1-sites arriveat a rate f S / S at each 2-site. This implies a rate 2 f S / S · N of particles to leave the system, as the reaction takes away two atoms. Alternatively, particles leave by LH rejection (rate-wise, this is merely a separate desorption process) at a rate f · N . The efficiency is the fraction of impinging particles thatreact; in the steady state, this is just the rate at which particlesleave due to reaction, normalized by the total rate to leave (byreaction or by LH rejection). This yields h p = S S + S , (17)which coincides with Eq. (16). We note that Eq. (17) can berewritten as h p = − S / S − S / ( S ) ≈ (cid:18) − S S (cid:19) (cid:18) + S S (cid:19) (18)for S / S ≪
1, which is precisely of the form of the empiri- cal relation (10). The coefficient inside the second bracket inEq. (10) deviates from 1 / S / S ∈ [ , ] , whereas Eq. (18) isstrictly valid only when S / S is small. VI. CONCLUSIONS
We have studied diffusion-limited reactions of particles ona two-dimensional lattice which consists of shallow and deepsites, using KMC simulations and rate equations. In the casewhen the two types of sites are randomly mixed, we foundthat the temperature range in which the reaction is efficientdramatically broadens compared to a homogeneous systemthat includes only shallow or only deep sites. The rate equa-tions are found to provide a good description of the systemand are in perfect agreement with the KMC results. We havealso studied a system in which the deep sites are clustered to-gether. In this case the hopping between shallow and deepsites is suppressed. As a result, the recombination efficiencyis dramatically reduced in comparison with the case in whichthe shallow and deep sites are randomly mixed.We expect that the qualitative features observed for the bi-nary distribution will hold for a broader class of models withdifferent distributions of binding energies. The results pre-sented in this paper are also relevant in the context of molec-ular hydrogen formation in the interstellar medium. Morespecifically, high abundances of molecular hydrogen are ob-served in photon-dominated regions [37]. In these regions,the grain temperatures are too high to form molecular hydro-gen from weakly adsorbed hydrogen atoms. It was proposedthat strong-binding sites in conjunction with the weak-bindingsites enable the efficient formation of molecular hydrogen un-der these conditions [38, 39]. Our work provides a quantita-tive basis for this mechanism.
Acknowledgments
This work was supported by Deutsche Forschungsgemein-schaft within SFB/TR-12
Symmetries and Universality inMesoscopic Systems and the Bonn-Cologne Graduate Schoolof Physics and Astronomy, and by the US-Israel BinationalScience Foundation. JK acknowledges the kind support andhospitality of the Hebrew University through the Lady DavisFellowship Trust. [1] D. ben-Avraham and S. Havlin,
Diffusion and Reactions inFractals and Disordered Systems (Cambridge University Press,2000).[2] B. A. Grzybowski, K. J. M. Bishop, C. J. Campbell, M. Fi-alkowski, and S. K. Smoukov, Soft Matter , 114 (2005).[3] R. J. Gould and E. E. Salpeter, Astrophys. J. , 393 (1963).[4] D. Hollenbach and E. E. Salpeter, J. Chem. Phys. , 79 (1970). [5] D. Hollenbach and E. E. Salpeter, Astrophys. J. , 155(1971).[6] D. J. Hollenbach, M. W. Werner, and E. E. Salpeter, Astrophys.J. , 165 (1971).[7] R. Smoluchowski, J. Phys. Chem. , 4229 (1983).[8] W. W. Duley and D. A. Williams, Mon. Not. R. Astron. Soc. , 177 (1986). [9] N. J. B. Green, T. Toniazzo, M. J. Pilling, D. P. Ruffle, N. Bell,and T. W. Hartquist, Astron. Astrophys. , 1111 (2001).[10] O. Biham and A. Lipshtat, Phys. Rev. E , 056103 (2002).[11] I. Lohmar and J. Krug, Mon. Not. R. Astron. Soc. , 1025(2006).[12] I. Lohmar and J. Krug, J. Stat. Phys. , 307 (2009).[13] I. Lohmar, J. Krug, and O. Biham, Astron. Astrophys. , L5(2009).[14] Q. Chang, H. M. Cuppen, and E. Herbst, Astron. Astrophys. , 599 (2005).[15] H. M. Cuppen and E. Herbst, Mon. Not. R. Astron. Soc. ,565 (2005).[16] A. G. Tsvetkov and V. I. Shematovich, Sol. Sys. Res. , 301(2009).[17] V. Mennella, Astrophys. J. Lett. , L25 (2008).[18] M. T. M. Koper, J. J. Lukkien, A. P. J. Jansen, and R. A. vanSanten, J. Phys. Chem. B , 5522 (1999).[19] E. W. Montroll, J. Math. Phys. , 753 (1969).[20] A. V. Barzykin and M. Tachiya, J. Phys. CM , 065105 (2007).[21] R. Straube, M. J. Ward, and M. Falcke, J. Stat. Phys. , 377(2007).[22] I. Langmuir, J. Am. Chem. Soc. , 1361 (1918).[23] J. W. Evans and R. S. Nord, Phys. Rev. A , 2926 (1985).[24] A. F. Voter, in Radiation Effects in Solids , edited by K. E.Sickafus, E. A. Kotomin, and B. P. Uberuaga (Springer, 2007),vol. 235 of
Nato Science Series II: Mathematics, PhysicsAnd Chemistry , URL .[25] E. W. Montroll and G. H. Weiss, J. Math. Phys. , 167 (1965).[26] O. Biham, I. Furman, V. Pirronello, and G. Vidali, Astrophys. J. , 595 (2001).[27] N. Katz, I. Furman, O. Biham, V. Pirronello, and G. Vidali, As-trophys. J. , 305 (1999).[28] H. B. Perets, Master’s thesis, The Hebrew University, Jerusalem(2004).[29] O. Biham, I. Furman, N. Katz, V. Pirronello, and G. Vidali,Mon. Not. R. Astron. Soc. , 869 (1998).[30] A. Lipshtat and O. Biham, Astron. Astrophys. , 585 (2003).[31] B. Barzel and O. Biham, J. Chem. Phys. , 144703 (2007).[32] T. Stantcheva, P. Caselli, and E. Herbst, Astron. Astrophys. ,673 (2001).[33] O. Biham, J. Krug, A. Lipshtat, and T. Michely, Small , 502(2005).[34] A. G. G. M. Tielens (1995), talk at a conference on interstellarchemistry in Leiden, The Netherlands.[35] J. Krug, Phys. Rev. E , 065102(R) (2003).[36] A. Lederhendler and O. Biham, Phys. Rev. E , 041105(2008).[37] E. Habart, F. Boulanger, L. Verstraete, C. M. Walmsley, andG. P. des Forˆets, Astron. Astrophys. , 531 (2004).[38] S. Cazaux and A. G. G. M. Tielens, Astrophys. J. , 222(2004).[39] S. Cazaux and A. G. G. M. Tielens, Astrophys. J. Lett.575