# Boundary effects in reaction-diffusion processes

\[

## Abstract

The effects of a boundary on reaction systems are examined in the framework of the general single-species reaction/coalescence process. The boundary naturally represents the reactants’ container, but is applicable to exciton dynamics in a doped TMMC crystal. We show that a density excess, which extends into the system diffusively from the boundary, is formed in two dimensions and below. This implies a surprising result for the magnetisation near a fixed spin in the coarsening of the one-dimensional critical Ising model. The universal, dimensionally-dependent functional forms of this density excess are given by an exact solution and the field-theoretic renormalisation group.

Date: .

\]

Reaction systems have been given a good deal of attention in the literature over the last few decades. They are widespread in nature and methods developed in their analysis have an applicability that extends well beyond conventional chemical systems. Moreover, they provide excellent examples of dynamical, many-body statistical processes and can exhibit a variety of interesting effects such as spontaneous symmetry breaking  and pattern formation .

In this work, we show that the imposition of impenetrable boundaries on a reaction system leads to a non-trivial spatial variation in the reactant density. Even though such a boundary would naturally represent the vessel holding the reactants, until now no such studies have appeared in the literature. In systems with many interacting degrees of freedom, boundaries often give rise to surface effects that penetrate far into the bulk. Our results demonstrate that in reaction systems such long-range effects are indeed present and have a high degree of universality.

We choose the class  of the general single-species reaction processes as the starting point for the study of these boundary effects in reaction systems, as it will provide a basis for the analysis of more complex systems. This universality class, comprising the annihilating random walk $`A+AO`$, coalescing random walk $`A+AA`$ and any combination thereof, is fundamental in the theoretical study of reaction systems and covers a broad range of physical phenomena. For example, the coalescence process is seen in the dynamics of excitonic annihilation reactions in the TMMC crystal . The predicted decay exponent from theory is in agreement with experiment for over five orders of magnitude. It should be noted that the boundary effects we introduce here are exhibited in the TMMC crystal in the presence of Mg<sup>2+</sup> or Cd<sup>2+</sup> doping. These defect ions act as perfect reflectors for the annihilating excitons .

A mapping also exists from the annihilating random walk to the domain coarsening dynamics in the critical one-dimensional Ising magnet . In this mapping the $`A`$ particles represent domain walls with an annihilation of two $`A`$ particles corresponding to a domain shrinking to zero size in a background of the opposite phase. An impenetrable boundary in the reaction system corresponds to a fixed boundary spin in the Ising magnet. Our analysis will show that domain walls are more likely to be found near the fixed spin than far away. This gives the counter-intuitive result that the absolute value of the dynamic, coarse-grained magnetisation is actually lower near a fixed spin.

In the interests of notational simplicity, we provide here the analysis specifically for the process $`A+AO`$. However, all densities given can be converted to the result for the coalescence process $`A+AA`$ by a simple factor of two. The annihilating random walk has been studied extensively in homogeneous unbounded systems, either of infinite extent or with periodic boundary conditions. Throughout the work, this will be referred to as the bulk case. It is well known  that the variation of the density $`\varrho `$ with time $`t`$ differs from the mean-field prediction $`\varrho 1/t`$ for dimensions $`d2`$. In fact the actual density decays are $`\varrho t^1\mathrm{log}t`$ ($`d=2`$), and $`\varrho t^{d/2}`$ ($`d<2`$) with a universal amplitude. In all cases, it must be stressed that the density remains uniform throughout the system.

In the following, we will first introduce the boundary into the reaction process and specify the model to be studied. The mean-field approximation will be shown to predict a homogeneous density unchanged from the bulk case. However, an argument will be presented to show that the mean-field prediction breaks down in low dimensions. The central result of this work is that a fluctuation-induced excess density develops at the boundary and extends into the system diffusively. We outline the field-theoretic renormalisation group (RG) description which we use to identify the universal quantities of this excess. These calculations were performed in real space to one-loop order and show that in two dimension and below, the density excess $`\varrho _E`$ has the following form

$`\varrho _E`$ $`=`$ $`{\displaystyle \frac{1}{(8\pi Dt)^{d/2}}}f_d\left({\displaystyle \frac{z^2}{2Dt}}\right).`$ (1)   
Here, $`D`$ is the diffusion constant of the reactants, $`z`$ is the normal distance from the boundary and $`f_d`$ is the dimensionally dependent scaling function. We were able to find the late-time scaling functions $`f_2`$ and $`f_1`$ exactly in both two and one dimensions (given in Eqs. (9) and (14) respectively). The former we derive from the RG improved field-theoretic calculation and the latter from an exact solution. Finally, in the context of surface critical phenomena, the behaviour in the related ballistic annihilation reaction system will be briefly examined and compared with the diffusive case.

We now introduce the model. The system is defined on a hypercubic, $`d`$-dimensional lattice with a lattice spacing of unity. The lattice, is infinite in $`d1`$ transverse dimensions and semi-infinite (sites $`1,2,\mathrm{},\mathrm{}`$) in what will be called the $`z`$ direction. At time $`t=0`$ the lattice is filled with an initial density $`\varrho _0`$ of identical particles that perform two types of dynamics: diffusion and mutual annihilation. The diffusion is represented by each particle hopping at a rate $`D`$ to any neighbouring site at random. A hop from site $`z=1`$ towards the boundary is disallowed. The diffusion is independent for each particle, and hence multiple occupancy of each site is possible, leading to bosonic particle statistics. However, if there are $`n2`$ particles on a particular site, a reaction can occur there with a rate $`\lambda n(n1)`$ to reduce $`n`$ by 2.

The above dynamics can be approximated by a mean-field description. This involves ignoring all possible correlations by considering a self-consistent equation involving just the average density $`\overline{\varrho }`$. In the continuum limit the boundary is enforced by a zero-current restriction, thus

$`_t\overline{\varrho }=D^2\overline{\varrho }2\lambda \overline{\varrho }^2`$ with $`_z\overline{\varrho }|_{z=0}=0.`$ (2)   
The boundary restriction is compatible with the bulk solution $`\overline{\varrho }=\varrho _0/(1+\varrho _02\lambda t)`$. Hence, in the absence of strong fluctuations, the density is uniform throughout the system.

However, correlations must be properly accounted for in low dimensions. Consider the dynamics of the model, in one dimension, up to a time $`t`$ and far from the wall. Because random walks in one-dimension are recurrent, most particles within a diffusion length $`l_b\sqrt{2Dt}`$ in the bulk, will have interacted and annihilated. This leads to a density in the bulk of the system of $`\varrho _bl_b^1c_b/\sqrt{t}`$. However, close to the wall the diffusion length is smaller, $`\varrho _wl_w^1c_w/\sqrt{t}`$ since $`c_w>c_b`$, leading to a density excess near the boundary.

Nevertheless, this argument is rather crude and a method for systematically including fluctuations is required. The RG has provided such a method for calculating bulk quantities , with the advantage of clearly identifying universal properties. We now present an overview of the generalisation to a system with a boundary: a case technically more complicated due to the lack of translational invariance. Details of the calculation will be provided elsewhere .

The field-theoretic description is obtained by first writing a master equation. This describes the flow of probability between microstates of the system and is conveniently written in second-quantised form $`_t|P>=|P>`$. The vector $`|P>`$ is the probability-state vector written in a Fock space and acted upon by the evolution operator $``$

$``$ $`=`$ $`{\displaystyle \underset{i}{}}\left[D{\displaystyle \underset{j}{}}a_i^{}(a_ia_j)\lambda (1(a_i^{})^2)a_i^2\right]`$ (3)   
where $`a^{},a`$ are the usual bosonic operators. The sum $`i`$ is over all lattice sites, and the sum $`j`$ is over all of site $`i`$’s neighbours, with the condition that both sums are restricted to the half-space.

The algebraic description $`(a,a^{})`$ is mapped onto a continuum path-integral for the action $`𝒮(\varphi ,\overline{\varphi })`$ where the complex fields $`\varphi `$ and $`\overline{\varphi }`$ are analogous to $`a`$ and $`a^{}1`$. The action $`𝒮=𝒮_D+𝒮_R+𝒮_{\varrho _0}`$ comprises diffusion $`𝒮_D`$, reaction $`𝒮_R`$ and initial-condition $`𝒮_{\varrho _0}`$ components. The diffusive part provides the propagator for the theory which is Gaussian for the transverse dimensions and has the the following “mirrored” form for the $`z`$-dimension

$`𝒢_z(z_f,z_i,t)`$ $`=`$ $`G(z_fz_i,t)+G(z_f+z_i,t)`$ (4)   
where $`G(z,t)`$ is a Gaussian with a standard deviation of $`2Dt`$. The other components in the action are

$`𝒮_R`$ $`=`$ $`2\lambda {\displaystyle _0^t}𝑑t{\displaystyle _{z>0}}d^dr\overline{\varphi }\varphi ^2+\lambda {\displaystyle _0^t}𝑑t{\displaystyle _{z>0}}d^dr\overline{\varphi }^2\varphi ^2,`$ (5)   
$`𝒮_{\varrho _0}`$ $`=`$ $`\varrho _0{\displaystyle _0^t}𝑑t{\displaystyle _{z>0}}d^dr\overline{\varphi }\delta (t).`$ (6)   
The upper critical dimension of the theory is $`d_c=2`$ and observables were rendered finite by dimensional reguralisation in $`d=2ϵ`$. The propagator is not dressed by the interactions (5), implying the boundary remains effectively reflecting on all scales. In the language of surface critical phenomena, this corresponds to the special transition  persisting at all orders. This is different from the behaviour frequently seen in equilibrium surface critical phenomena  and in related non-equilibrium systems . In fact only the reaction rate $`\lambda `$ is renormalised, with a fixed point structure identical to the bulk case . This is understandable as physically the renormalisation of $`\lambda `$ is connected to the fact that random walks in two dimensions and below are recurrent: a feature unaffected by the presence of a boundary.

To get non-trivial, $`z`$-dependent results it is clear (from the lack of an excess in the mean-field equation) that the RG improved perturbation expansion must be taken to at least one-loop order. Writing the density as an expansion in $`ϵ=2d`$ and splitting the contribution into an excess $`\varrho _E(z,t)`$ and a homogeneous, background bulk density $`\varrho _B(t)`$

$`\varrho (z,t)`$ $`=`$ $`\varrho _B(t)+\varrho _E(z,t),`$ (7)   
the homogeneous bulk density for $`d<2`$ is found to be

$`\varrho _B(t)`$ $`=`$ $`{\displaystyle \frac{1}{4\pi ϵ(Dt)^{d/2}}}\left[1+{\displaystyle \frac{ϵ}{4}}\left(2\mathrm{log}(8\pi )5\right)\right]+O(ϵ).`$ (8)   
This is exactly the result found in  as expected.

However, a fluctuation-induced density excess is also found, representing the new result from this calculation

$`\varrho _E(z,t)`$ $`=`$ $`{\displaystyle \frac{1}{8\pi (Dt)^{d/2}}}f_2\left({\displaystyle \frac{z^2}{2Dt}}\right)+O(ϵ),`$ (9)   
$`f_2(\xi ^2)`$ $`=`$ $`{\displaystyle _0^1}𝑑s{\displaystyle _0^s}𝑑r\left({\displaystyle \frac{r}{s}}\right)^2{\displaystyle \frac{\mathrm{exp}\left(\frac{\xi ^2}{(2sr)}\right)}{\left[(sr)(2sr)\right]^{1/2}}}.`$ (10)   
The function  has the asymptotics $`f_2\mathrm{exp}^{\xi ^2/2}/\xi ^3`$ and is plotted in Fig. 1. A few things should be noted about the form of $`\varrho _E`$. First, this excess is not localised at the boundary but extends into the system diffusively, by virtue of its functional dependence of $`z^2/Dt`$. Also, the excess shares the same universality as the bulk density in that it is independent of the reaction rate $`\lambda `$ and the initial density $`\varrho _0`$. Finally, for $`d<2`$ the amplitude of the excess decays with the same exponent as the bulk. Hence, the ratio of the boundary to the bulk density

$`{\displaystyle \frac{\varrho (0,t)}{\varrho _B(t)}}`$ $`=`$ $`1+ϵ\left({\displaystyle \frac{3}{4}}+{\displaystyle \frac{\pi }{2}}{\displaystyle \frac{3\pi ^2}{16}}\right)+O(ϵ^2)`$ (11)   
is a constant, universal quantity independent of all system parameters except the dimension.

The behaviour in two dimensions provides for an interesting result: the one-loop calculation Eq. (9) is the exact late-time density excess. It is independent from the renormalised reaction rate and therefore represents the universal leading order, with higher-loop corrections decaying as $`(t\mathrm{log}t)^1`$. The excess given in Eq. (9) gives surprisingly accurate results even for short times, as can be seen in Fig. 1 where the result for the density excess is compared with simulations.

The RG has provided information about the general behaviour of the density excess as a function of dimension. The universal quantities have been identified and the excess density correctly predicted for late times in $`d=2`$. Unfortunately, the $`ϵ`$ expansion gives disappointing results for the amplitude in $`d=1`$. Motivated by the mapping of the reaction dynamics onto the Ising magnet and the related excitonic coalescence process, we now provide the exact solution in one dimension. We generalise the bulk model described in  to include an impenetrable boundary. The on-site reaction rate is chosen to be infinite and hence a site in the system can be occupied by at most one particle. The dynamical rules can be written in terms of the possible evolutions of a pair of neighbouring sites. Denoting a particle on site $`k`$ by $`A_k`$ and an empty site by $`O_k`$, the allowed changes are

$`O_zA_{z+1}`$ $``$ $`A_zO_{z+1}\text{ with a rate D}`$ (12)   
$`A_zA_{z+1}`$ $``$ $`O_zO_{z+1}\text{ with a rate 2D}`$ (13)   
where the site label $`z=1,2\mathrm{}`$ is restricted to positive integers. The master equation can be written in the language of spin-half operators, and we were able to solve for the density using a similarity transformation found in  (details of the calculation will be presented elsewhere ). The density can be written as a sum of Bessel functions, of which the continuum limit gives the following density excess

$`\varrho _E(z,t)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{8\pi Dt}}}f_1\left({\displaystyle \frac{z^2}{2Dt}}\right),`$ (14)   
$`f_1(\xi ^2)`$ $`=`$ $`\sqrt{2}\text{erfc}\left({\displaystyle \frac{\xi }{\sqrt{2}}}\right)\mathrm{exp}\left({\displaystyle \frac{\xi ^2}{2}}\right)\text{erfc}\left(\xi \right).`$ (15)   
This function has the asymptotics $`f_1\mathrm{exp}^{\xi ^2/2}/\xi `$ and is plotted in Fig. 1. Again, the bulk component of the density $`\varrho _B=1/\sqrt{8\pi Dt}`$ was found to be identical to the infinite case . As expected, the form of the solution is in agreement with that predicted from the RG treatment. The universal ratio of the density near the boundary to the bulk density given in Eq. (11) is found to be exactly $`\sqrt{2}`$. It should be noted that the one-loop RG prediction for the universal ratio (given in Eq. (11) with $`ϵ=1`$) is numerically $`1.47`$ for one dimension, and therefore gives quite a fair indication of the exact value. This should be compared with the poor one-loop predictions for the amplitudes of the density itself in one dimension.

The universality (independence from the reaction rate $`\lambda `$) predicted by the RG treatment, can be seen by comparing this exact result for $`\lambda =\mathrm{}`$ with data from a simulation of a system with finite reaction rate (Fig. 1).

The result, Eq. (14), gives the time-dependent probability density of domain walls in the coarsening, one-dimensional critical Ising model near a very strong magnetic field, i.e. a fixed spin. The magnitude of the coarse-grained magnetisation is a function of the local density of domain walls - the fewer the domain walls the higher the magnetisation. Hence, the result implies that the absolute value of the magnetisation measured near a fixed spin is lower than in the bulk of the system. It would be very interesting to see if this dynamic effect is seen in other magnetic systems or in higher dimensions.

In summary, results have been presented from an analysis of a reaction-diffusion process near an impenetrable boundary. The mean-field equation was shown to predict a flat density profile. However, it was demonstrated that in two dimensions and below, a density excess develops at the boundary. In one dimension it was found that the density excess is as significant a contribution as the bulk density with both decaying as $`t^{1/2}`$. In two dimensions the excess was found to be marginally subdominant, decaying as $`t^1`$. Both these density excesses shares the same universality as the bulk density, in that they are independent of the reaction rate and the initial density. Moreover, a higher degree of universality was found in the ratio between the boundary and bulk densities: a quantity depending only on the dimension of space. The functional forms for the density excess were obtained for both one and two dimensions and a link was made to the coarsening dynamics of the one-dimensional Ising model. Most importantly, it was shown that the excess is not localised at the boundary but extends diffusively throughout the bulk of the system. It should be noted that above the critical dimensions there is a sub-dominant density excess. However, these effects are transient and quickly decay to yield the mean-field result, as predicted.

A quantity much studied recently and related to domain coarsening dynamics is the persistence exponent. This exponent describes the time dependence of the distribution of sites not yet visited by a domain wall. In the homogeneous bulk case, persistence has been examined by the RG  and an exact result found in one dimension . In light of result (14) and the fact that the excess penetrates into the system, it would be worth examining the behaviour of the persistence in the presence of a fixed spin.

Finally, it is interesting to consider the difference between the behaviour of diffusing and gas-phase ballistic reaction kinetics near a boundary. A simple realisation of the ballistic $`A+AO`$ system  shows the same bulk density-decay exponent in one dimension as the diffusive case. However, contrary to the diffusive case, it can be shown that for ballistic reactions with an (elastic) impenetrable boundary there is a lower density of reactants near the wall . This is understandable, because in late times the remaining particles tend to congregate in groups moving in the same direction. When such a group hits the elastic wall it mostly annihilates within itself, leaving few particles to return. In the context of surface critical phenomena this corresponds to the impenetrable boundary behaving as an effectively absorbing boundary on long time scales, in contrast to the diffusive case.   

Acknowledgments. We would like to thank David Mukamel and Gunter Schütz for useful discussions. Fabian Essler, Martin Evans and Klaus Oerding are also thanked for their careful reading of the manuscript. The authors acknowledge support from the Israeli Science Foundation.