{'fraction_non_alphanumeric': 0.028276720611773075, 'fraction_numerical': 0.012235461497421305, 'mean_word_length': 4.581893984514592, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fluctuation-driven heterogeneous chemical processes

We explore a new framework for describing the kinetics of a heterogeneous chemical reaction where two particles of the same chemical species form a reaction product of another chemical species on the surface of a seed particle. Traditional treatments neglect the effect of statistical fluctuations in populations. We employ techniques in a manner analogous to the treatment of quantum systems to develop a stochastic description of processes beyond the mean field approximation.


Introduction
Traditional models employing the evolution of the mean population of species in a system provide a good enough description of various physical processes such as heterogeneous chemical reactions and heterogeneous nucleation of aerosols as long as the average particle number is large. However, when the mean populations involved are small we expect to see a significant deviation from the solution to the classical equations.
As an example consider the chemical reaction of two hydrogen atoms on the surface of a dust particle in interstellar space [13]. Atoms can be adsorbed onto and evaporated from the grain particle. The adsorbed molecules will diffuse on the surface of the seed and eventually collide with another reactant to form diatomic hydrogen. Under interstellar conditions the rate of the adsorption of atoms will be small compared to the reaction rate. If the mean population of hydrogen atoms is small, statistical fluctuations brought about by gain and loss, and by the random diffusion of atoms on the surface of the dust particle, are important. Yet the traditional model does not include the treatment of such fluctuations. Several attempts have been introduced to resolve that problem; among them the Monte Carlo approach, the modified rate approach, the direct master equation approach and the Gauge Poisson representation approach -see [9,10] and [7]. Related studies have been carried out in [4,11]. Analytical solutions to a master equation approach have been found for the steady state case [1,8,16]. Our work is a first step to an alternative approach for describing the kinetics of a heterogeneous chemical reaction which can be extended to other areas where a similar problem occurs, for example when computing the rate of nucleation processes taking place on small particles. Again, the number of adsorbed molecules is small and statistical fluctuations need to be taken into account. In [16] the steady state solution was studied employing a description using master equations instead of mean population dynamics. A nucleation rate lower than the classical equations predicted was obtained. Therefore, it is important to develop a model by which we replace the classical equations with a set of stochastic equations.
For the correct treatment of population fluctuations we employ methods of Quantum Field Theory [5,6]. Starting from a master equation we introduce a spatial lattice where the microstates of the system correspond to a set of occupation numbers at each lattice site. A Fock space is constructed using annihilation and creation operators at each lattice site. By means of this set-up it is easy to show that the master equation is equivalent to a Schrödinger equation with imaginary time. This enables us to employ techniques originally developed in order to describe a quantum mechanical system where the fluctuations are due to a quantum uncertainty. We obtain the average particle population of the classical many-body system by developing a mechanism for computing expectation values of observables analogous to Feynman's path integral formulation [2,18]. In this paper we will not concern ourselves with renormalisation group analysis although studies in that direction have been undertaken by various people, for example [12,15,21]. Introducing a stochastic variable [20], a mathematical trick -for a nice review paper we refer to [23]helps with the evaluation of the expression for the expectation values. The complex fluctuating solutions to a set of constraint equations, which are stochastic partial differential equations, are then averaged over all realisations of the stochastic noise. For numerical investigations, the solutions to the constraint equations can be generated by various numerical schemes [14]. The path integral average is computed using Monte Carlo methods [17]. The Code is written in C and computations do not take longer than a few seconds up to a few minutes running on a standard laptop.
In the sequel we will concentrate on a heterogeneous chemical reaction where two reaction partners of the same sort of particle type react on the surface of a seed particle to form a reaction product.
For readers more interested in the physical content than in the mathematical details of the formalism we recommend to skip the first few pages and instead have a look at the most important equations on page 8 and continue from there. Equation (42) which is a path integral average (PIA) gives the average particle density of the various chemical species once the complex, fluctuating solutions to the constraint equations (40) and (41) have been found and inserted into the PIA. Together with the correlations (43) for the stochastic noise occuring in the PIA this forms a complete set of equations.

Master Equation And Schrödinger-like Equation
We concentrate on chemical reactions of type A + A −→ C, that is situations in which the atoms adsorb onto grain particles where they can associate with themselves to produce diatomic molecules. If the number of reactive species on an individual grain is small, traditional rate equations will fail to accurately describe the diffusive chemistry occuring on the surface of the grain particle. We start our investigations by determining a master equation that describes such a heterogeneous chemical process.
A general master equation can be written as where T n→m represents the transition amplitude or propagator from a microstate n to a microstate m and P (m) is the probability to find the system in state m. Considering a d-dimensional lattice L with lattice constant l, the microstates correspond to the occupation numbers {N i } = {N 1 , N 2 , ...} at each lattice site i.
The chemical reaction of pairs of species A to form a product of species C on a particle or droplet , A + A −→ C, is modelled by the following master equation This equation describes the evolution of the probability distribution P ({N A }, {N C }; t) for the total number of adsorbed molecules {N A } of species A and for the number of reaction products {N C } of species C. The symbols N Ai,Ci denote the numbers of A or C molecules at lattice site i, respectively. The constants j A,C , κ and λ A,C are rate coefficients, D A,C is the diffusion constant and V stands for the volume of the droplet. The rate coefficient j A,C is called source rate, the rate coefficient λ A,C is called evaporation rate and κ is known as the reaction rate.
In our model, the chemical reaction is taking place on a d-dimensional lattice, allowing for multiple occupancy on each site. This configuration is also called bosonic representation.
The changes in population which we consider are caused by: The summation in the fourth and fifth line of the above equation (2) is taken over nearest neighbour sites only. The factors (N A,C + 2), (N A,C + 1), N A,C , (N A,C − 1) describe the number of ways of choosing the particles involved in the considered process. In the continuum limit the particle hopping from one site to another corresponds to the diffusion of the particles. The initial condition is chosen corresponding to a Poissonian distribution on each site where n A,C (0) is the average occupation number per lattice site for the A or C particles respectively.
In the next step we will apply the methods of second quantisation [5,6]. We will rewrite the master equation as a Schrödinger-like equation for a many-body wave function. This approach can be justified by noting that, first of all, the master equation is a differential equation of first order with respect to time. The second reason to suggest the treatment of the master equation according to the second quantisation is that the master equation is linear in the probability.
In order to simplify the notation we will suppress the dependence on space coordinates x = (x 1 , x 2 , ..., x d ). We will be working in an appropriate space, the Fock space. A Fock space F ν (H) is a Hilbert space made from the direct sum of tensor products of single-particle Hilbert spaces H with S ν a symmetrising (for the case of bosons) or antisymmetrising (for the case of fermions) operator. The Fock space is constructed by introducing the following operators at each lattice site i + a i , i ∈ L : creation operator, − a i , i ∈ L : annihilation operator, which satisfy the commutation relationships The vacuum state |{0} is defined by with where |0 j denotes the vacuum state in a single-particle Hilbert space.
The master equation (2) is equivalent to the Schrödinger-like equation -a Schrödinger equation with imaginary time- with the many-body wave function and the Hamiltonian operator For verification of the above statement one has to insert the states i ( For the timebeing let us suppress the subindices that identify the particle type. The form of the many-body wave function (9) can be made plausible when considering the state vector |N i at site i ∈ L, namely It holds that

Expectation Values Of Observables
We are interested in obtaining the expectation values for various observables, especially the average number density. The expectation values of observables O are given by We want the expectation values of the observables -diagonal in the occupation number basisto be linear in the probabilities. After some straightforward computation -see for example [23]one can see that the above expression is equivalent to where The projection state obeys the relation By definition, the projection state is a left eigenstate of all creation operators with unit eigenvalue Furthermore, {P}|Ψ(t) = 1. Conservation of probability requires that {P}|H = {0}|.
We break the time interval [t 0 , t T ] into T short slices of duration ∆t = tT −t0 T . At each time slice we insert a complete set of coherent states -see, for example, [2]. Coherent states, |C i (t) , are right eigenstates of the annihilation operator where the eigenvalue Φ i is a complex function. The duals C i (t)| are left eigenstates of the creation operator We have The coherent states are over-complete. Still, we can use them to create the identity for a single lattice site i ∈ L, and for multiple lattice sites accordingly with |{C} = j |C j . Let us recall the formula for the expectation values of observables where the initial many-body wave function takes the form -see equations (3) and (9)- We observe the following proportionalities for all admissible values of j. Therefore, one can recast the equation for the expectation values where where we have labeled each time slice by a time index τ ∈ [0, ∆t, 2∆t, ..., T ]. The normalisation constant has to be determined lateron. The consideration of the lattice expectation value is not sufficient if one is interested in long wavelength properties. In the formal continuum limit, we obtain Next, we expand the exponential function for small ∆t, neglect higher order terms in ∆t and recast the continuous average. When interested in inclusive probabilities -e.g. the average number of particles at a given lattice site irrespective of the number of particles elsewhere-it is convenient to commute the factor of e The operators are now normal ordered. It is valid that the operator O and its normal ordered counterpart have the same expectation value if all creation operators occuring in the normal ordered operator are replaced by unity -see for example [23].
In particular, the density operator + a − a reduces to the annihilation operator − a. Therefore, in the continuum limit the average particle density of the A particles is given by the stochastic average of the complex eigenvalues of the coherent state vectors under the annihilation operator, that is one chooses the operator O to be Φ A (x, t).
We take the continuum limit -the dimensions of the constants are chosen by looking at the discrete Hamiltonian operator (10) In the continuum limit, the average particle density in the stochastic model is then given by where D denotes the measure of the functional integral andΦ A,C is the shifted eigenvalue of the dual of the coherent state under the creation operator defined byΦ A,C := Φ * A,C − 1. Accordingly, all fields that incorporate shifted eigenvalues instead of the original eigenvalues will be denoted by a twiddle in the sequel. Note that the average (31) is performed taking into account the dynamics and the initial conditions -for a more detailed discussion see [23]. We already have incorporated the shifted initial state in the shifted actionS. The symbolS A+A→C represents the shifted action which is defined as follows with the shifted HamiltonianH A+A→C . The shifted action for the chemical reaction A + A −→ C takes the form We want to untangle the quadratic termΦ 2 A . A linear expression inΦ A can be obtained by means of a Gaussian transformation where P[η] A+A→C is the probability distribution for a white noise η(x, t). The Gaussian distribution reads Now the shifted actionS is linear inΦ A and one can easily integrate out overΦ A (x, t) and Φ C (x, t) in (31). One obtains where δ[F A,C ] is a functional Dirac delta distribution. In its generalised Fourier representation it is defined by with z(y) being a multicomponent field satisfying the constraint It follows from equation (37) that the average particle density for the A respectively C molecules is now given by The stochastic noise η has zero mean value and a correlation given by This is obvious when considering the Gaussian distribution (36). Note that the above average is no longer taken over the initial conditions. The constraint equation (40) is an inhomogeneous partial stochastic differential equation with additive noise for a complex fluctuating unknown field in the Itô calculus. It resembles the deterministic partial differential equation that describes the evolution of the mean particle density in the classical theory. But despite the temptation for an intuitive interpretation it is very important to keep in mind that in equation (40) we are confronted with a complex fluctuating quantity that has as such no physical interpretation. Only if the path integral average (PIA) (42) of a solution to (40) or (41) is taken over all possible realisations of the stochastic noise that appears in the constraint equation can one interpret the outcome of this computation as a mean particle density.

Case A: Vanishing Source Rate
In the remainder of this paper, let us concentrate on a single spatial site model. We will now compare the results in the stochastic model to the observations made in the traditional approach.
The classical equation for the evolution of the mean particle density in the single spatial site model reads where n A (t) denotes the mean particle density of the A molecules in the mean field approximation. The classical evolution equation (44) is solved by The stochastic constraint equation for the complex fluctuating fieldΦ A (t) in zero spatial dimensions with vanishing source rate takes the form The traditional equation for the average particle density of the A molecules in the mean field approach (44) and the stochastic constraint equation associated with the A molecules (46) resemble each other at first sight. But as mentioned before the solution of the stochastic differential equation (46) is a complex, fluctuating field that can only be interpreted as an average particle density after it has been averaged in the sense of equation (42). For vanishing source rate it is fairly easy to find an analytic solution of equation (46). The stochastic constraint equation (46) -because of the continuous but not smooth nature of a stochastic process-has to be understood in terms of a stochastic integral equation where The stochastic noise η(t) is rewritten in terms of the Wiener process W (t) For vanishing source ratej A the stochastic constraint equation for the A particle density (40) reduces to the following equation The above equation is a nonlinear reducible stochastic differential equation with polynomial drift of degree two in the Itô picture. In contrary to a Stratonovich stochastic differential equation an Itô stochastic differential equation can not be solved directly by methods of classical calculus 1 .
For an analytical solution of an Itô stochastic differential equation one has to use a modified version of the drift coefficient where the derivative in the last term is a functional derivative. Equation (50) is a stochastic version of a Verhulst-like equation -see [14]. It can be reduced to a linear stochastic differential equation with multiplicative noise. We obtain the solution to the first stochastic constraint equation (40) for vanishing source rate, namelȳ Inserting the above solution into the path integral average (42) one obtains the average particle density for the A molecules in the stochastic picture.
The stochastic constraint equation for the reaction product, the C particles, in zero dimensions looks -in its form-identical to the classical evolution equation for the mean density of C particles In the single spatial site model, the full solution of the second constraint equation -simply takē Φ C (t) instead of n C (t) in the above equation (53)-can be obtained even for non-vanishing source rate and readsΦ The stochasticity of the above equation is hidden in the first term employing the fluctuating solutionΦ A (t) of the first constraint equation (46).
Once the solutions to the stochastic constraint equations (52) and (54) are known, one has to insert either of the solutions into the path integral average (42) and compute the path integral by means of a Monte Carlo calculation in order to obtain the average particle density for the A or C particle population, respectively. Random samples are generated according to the Gaussian probability distribution (36); that is we generate Wiener processes. We estimate the path integral (42) by summing a large number of solutions of the constraint equations associated to the set of random samples generated in the above sense and divide the sum by the number of random samples. The Monte Carlo method displays a convergence of 1 √ N where N is the number of random samples -see [19].
Instead of using the expressions of the analytical solutions, equation (52) and equation (54), to generate solutions to the stochastic constraint equations one can alternatively compute the paths directly from the stochastic differential equations (40) and (41). The latter method turns out to be less time consuming. The stochastic differential equation (40) in zero dimensions can be converted into where X A,n := Φ A (t n ) in discretised time t n for n = 0, .., N , ∆ n := t n+1 −t n and ∆W n := W tn+1 − W tn . The ∆W n is generated by two uniformly distributed independently random variables via the Box-Muller transformation -see, for example, [14]. The numerical scheme (55) is called the Euler scheme and is the most straightforward approach to undertake some numerical investigations. Accordingly, the second constraint equation in the single spatial site model and for vanishing source rate takes the following form As stochastic differential equations are extremely sensitive one has to convince oneself that the code is stable and converging as it should be. Other schemes we used that might be more accurate or stable than the Euler method are the Milstein scheme, the simplified order 2.0 weak Taylor scheme, the implicit order 1.0 strong Runge-Kutta scheme, the predictor-corrector method of order 1 with modified trapezoidal method weak order 1.0 -see [14].
For the numerical evaluation we employ values for the rate coefficients that can be found in realistic physical set-ups. Instead of using the coefficients introduced in the continuum limit, namelyλ A,C andκ, we employ the traditional rate coefficients which we denote by L A,C and K A,C and which have dimensions per unit time. Dimensional analysis shows that using these new constants we are now calculating an average particle population instead of an average particle density -this can be easily verified in equations (46) and (53).
The following plots were generated for the situation where two hydrogen atoms react on the surface of an interstellar dust particle. According to [3,22] the reaction rate takes the value of K = 1.45 × 10 5 s −1 , the evaporation rate for the hydrogen atoms L H = 1.88 × 10 −3 s −1 and the evaporation rate for the reaction product L H2 = 6.9 × 10 −8 s −1 . As initial values we used Φ A (0) =Φ C (0) = 6.
In Figure 1 we generate the real part of one solution to the stochastic constraint equation for the hydrogen atoms under the above conditions. Figure 2 shows the imaginary part of the same solution of the stochastic constraint equation for the reaction partners. If one compares these plots to Figure 3 and Figure 4 where the real and imaginary part of the path integral average over 1000 realisations of the white Gaussian noise for the H atoms are given one observes that the real part of the path integral average smoothes out and the fluctuations in the imaginary part decrease in intensity. For increasing number of paths employed in the path integral average the imaginary part of the PIA tends to zero. Therefore, it is safe to interpret the real part of the path integral average as the average particle population. Together with Figure 7 one can interpret the path integral average in Figure 3 in the following way: according to Figure 7 the chemical reaction stops after certain transient processes. The average population of the diatomic hydrogen is constant. The intuitive physical reason for this is because all the potential reaction partners H have already been used to form the reaction product H 2 . On the other hand, from Figure 3 one sees how the plot for the hydrogen atoms approaches asymptotically the value one half which could be interpreted as a state consisting of either zero or one particle. The value of Φ A (t → ∞) = 1/2 is the lowest possible eigenvalue of the coherent states under the annihilation operator once the system has reaches its equilibrium. The rate coefficients K and L A,C have to be understood in a probabilistic sense, in a similar fashion as it is done with, say, the mean life expectancy of a radioactive isotope. For the specific values we used in our calculations the reaction rate dominates over the evaporation rate.
The discussion above can be compared with the results obtained from the solution to the classical evolution equation (45) which predicts an asymptotic value n A (t) −→ 0 as t −→ ∞ for L A = 0. That is, the classical model predicts the extinction of all the reactants.

Case B: non-vanishing source rate
As in the previous section, we generate solutions to the constraint equations (40) and (41) in zero dimensions by the numerical schemes discussed above but now with a non-zero source rate. We compare the results to the solutions of the classical evolution equations which read where α := 8κj A +λ 2 A and β := 2/α arctanh (4κn A (0) +λ A )/α and n C (t) = e −λC t t 0 eλ C s (j C +κn 2 A (s))ds + n C (0) .
We now analyse the influence of the source rate coefficient associated with the adsorption of reactants onto the surface of the seed particle. For convenience, we will leave the source rate for the reaction products at zero as it will not have significant influence on the outcome of our discussion. The other rate coefficients were chosen as before, J A,C , L A,C and K A,C denoting the rate coefficients per unit time. The plots in Figures 9 to 12 were obtained for a source rate that              (Figure 9) and in the real part of one path associated with the second constraint equation (Figure 11), Figures 9 to 12 reproduce deterministic behaviour. The real part of the path integral average of the reactants (Figure 10), as well as the real part of the path integral average of the reaction products ( Figure 12) over 1000 possible realisations of the stochastic noise, coincide with the results of the associated classical equations. This accordance between classical and stochastic model is no longer valid for Figures  13 to 16. As an example, we generated one single path for each particle population, the reaction partners and the reaction products, and plotted their real part in Figure 13 and Figure 15, respectively.
Let us now compare the average particle populations of the hydrogen atoms and the diatomic hydrogen for large source rate ( Figure 10 and Figure 12) to the average particle populations of the reactants and the reaction product for small source rate ( Figure 14 and Figure 16). In the deterministic case, that is for large source rate, the chemical reaction does not die out after a certain period of time -see Figure 12-in contrary to the observations made from Figure 16 where the real part of the path integral average of the diatomic hydrogen over 1000 paths is plotted for a small source rate. As can be seen from Figure 10 and Figure 14 respectively, for a source rate of J H = 10 8 s −1 the average particle population of the H atoms reaches an asymptotic value of 22.28 whereas for a source rate of J H = 10 −7 s −1 the average particle population of the H atoms is 1/2 as was the case for vanishing source rate in the latter section. Figure 12 and Figure 16 give the average particle population for the H 2 atoms at time t = 8 × 10 −5 s −1 , namely 3996.89 for large source rate J H = 10 8 s −1 and 8.73 for small source rate J H = 10 −7 s −1 . To determine the transition from a deterministic to a stochastic behaviour we computed the average particle population of the reaction partners and products for a source rate of the hydrogen atoms in the range of J H ∈ [10 8 s −1 , 10 −7 s −1 ] for each order of magnitude and generated Figure 17 to Figure  19. Already when the reaction rate and the source rate are of the same order of magnitude one can observe deviations from the classical behaviour, in the sense that the equilibrium value of the average particle population obtained by the stochastic methods is higher than that predicted from the classical equation (57). One also observes from Figure 17 that for a source rate between J H = 10 3 s −1 and J H = 10 2 s −1 the chemical reaction dies out because there is not one pair of reaction partners left in order to initiate a chemical reaction.

Conclusions
We have argued that the classical evolution equations for the mean particle population of a chemical species involved in a heterogeneous chemical reaction do not give the right results for small systems. Instead, we developed a stochastic model that includes statistical fluctuations and showed in our numerical investigations that those fluctuations can not be ignored for low rates of particle adsorption onto the surface of a grain particle. Although one starts from an apparently classical system the introduction of a quantum field theoretical formalism for its description seems to force us to adopt a "quantum mechanical-like" interpretation of the results.
It is possible to extend this work to other chemical reactions, for example of the type A+B −→ C. In a next step, we shall consider a network of chemical reactions in which several reactions compete against each other. One would expect that it will take considerably longer to reach an asymptotic value for the average particle population of a certain species.

Acknowledgements
This work was supported by the Leverhulme Trust under grant F/07134/BV and partly supported by the C N Davies Award of the Aerosol Society.                     Figure 19: Same curve as shown in Figure 18 but zoomed to resolve the plot close to the horizontal axis.