format: text

text:
# DIFFUSION-LIMITED REACTION IN ONE DIMENSION: PAIRED AND UNPAIRED NUCLEATION

## 1 Introduction

Reaction rates controlled by collisions between diffusing particles depend on the distribution of distances between particles as well as on the density of particles. In particular, as Noyes stated in 1961 “Any rigorous treatment of chemical kinetics in solution must consider concentration gradients that are established by the existence of the reaction itself”.  Here, we study the dynamics of point particles in one dimension, nucleated at random positions and times, then diffusing until colliding with and annihilating another particle. Competition between nucleation and annihilation produces a statistically steady state with a well-defined mean density of particles and distribution of distances between particles. We shall contrast two types of nucleation: unpaired, in which particles are deposited at random locations at random times, and paired, in which pairs of particles are deposited at random locations. The dynamics is as follows:

1. Particles are nucleated in pairs with initial separation $`b`$;
2. Nucleation occurs at random times and positions with rate $`\mathrm{\Gamma }`$;
3. Once born, all particles diffuse independently with diffusivity $`D`$; and
4. Particles annihilate on collision.

A portion of a typical realization of these dynamics is shown in Fig. 1. For unpaired nucleation (i) and (ii) are replaced by

1. Particles are nucleated at random times and positions with rate $`Q`$.

An existing method of analysis, based on a truncated hierarchy of correlation functions, is developed and extended in this article to the case of paired nucleation, yielding expressions for the correlation functions in the steady state, and for the time scales for relaxation towards the steady state. We also introduce a different method of analysis that yields an exact explicit expression for the steady-state density and for the time dependence of the density starting from arbitrary initial conditions. Our analytical predictions are compared with the results of direct numerical simulations. In the simulations, large numbers of diffusing particles are simultaneously evolved in continuous space, with annihilation whenever two paths cross and nucleation (paired or unpaired) at random times and positions.

A striking difference between paired and unpaired nucleation is the scaling of the steady state density of particles, $`\rho _0`$, with the nucleation rate: $`\rho _0\mathrm{\Gamma }^{\frac{1}{2}}`$ (paired) versus $`\rho _0Q^{\frac{1}{3}}`$ (unpaired). Here, we shall exhibit the crossover between these two cases in terms of the following dimensionless quantity:

$$\epsilon =\left(\frac{2\mathrm{\Gamma }}{D}\right)^{\frac{1}{3}}b.$$
(1)
For $`\epsilon \mathrm{}`$, the dynamics described by (i)–(iv) is equivalent to that described by (i’), (iii)–(iv) with the replacement

$$Q2\mathrm{\Gamma }.$$
(2)
The paper is arranged as follows. In the remainder of this section we summarize published results for reaction-diffusion systems. In Sec. 2 we analyze the dynamics using a hierarchy of equations for particle density functions, called “reduced distribution functions” by Van Kampen . Derivation of the reaction kernel leads to an exact relation between the density of particles and the derivative of the correlation function. We also explore the linear response to a perturbation away from the steady state to establish the time scales for relaxation. In Sec. 3, by introducing a function that satisfies a closed linear partial differential equation, we present exact expressions for the steady-state density and for the time evolution of the density with arbitrary initial conditions. In particular, analytical results are presented describing the rapid initial annihilation that transforms an initially random distribution into one characterized by an effective repulsion between particles.

### 1.1 Unpaired nucleation

Analysis of diffusion-limited reaction dates back to M. von Smoluchowski. His Mathematische Theorie der raschen Koagulation  considered reaction between diffusing particles resulting in merger, with the reaction taken to occur immediately whenever two particles are a distance $`R`$ apart. He introduced a diffusion equation for the density of particles relative to the position of a test particle and noted that the density is zero at all times at radius $`R`$ . For many years it was assumed that the final result of a complete calculation following the procedure outlined by Smoluchowski would be an equation for the mean density of particles, $`\rho `$, of the form

$$\dot{\rho }=Qk_s\rho ^2,$$
(3)
where $`Q`$ is the rate (per unit length and time) of appearance of new particles and $`k_s`$ is constant. This would imply, for the case without nucleation ($`Q=0`$), that the density is proportional to $`t^1`$ for $`t\mathrm{}`$. However, arguments based on dimensional analysis and scaling show that this is not true in one dimension . In 1983, Torney and McConnell studied this case and published an exact solution for the mean density as a function of time . Starting from an initial random distribution of particles, they found

$$\rho (t)=\rho (0)\mathrm{exp}\left(8Dt\rho ^2(0)\right)\mathrm{erfc}\left(\rho (0)(8Dt)^{\frac{1}{2}}\right).$$
(4)
In particular, $`\rho (8\pi Dt)^{\frac{1}{2}}`$ for $`t\mathrm{}`$. A rederivation of the result of Torney and McConnell was provided by Spouge , whose insight was that an annihilation process is equivalent to a coagulation process if coagulants made up of an even number of particles are considered as diffusing “ghosts.” Derivations based on reflection principle  and field theory  methods have also been published.

In discrete models of diffusion-limited reaction, diffusion is approximated by hopping between neighboring sites on a lattice. Here, too, the density of particles without nucleation is proportional to $`t^{\frac{1}{2}}`$ for $`t\mathrm{}`$ . Moreover, with unpaired nucleation, the steady state density is proportional to the third power of the nucleation rate . This can be interpreted as evidence for a time-dependent rate constant $`k_s`$ in (3), or as requiring (3) to be replaced by an equation of the form

$$\dot{\rho }=Qk_c\rho ^3.$$
(5)
However, no polynomial equation for the density can describe both the steady state with nucleation and the long-time decay of the density without nucleation .

An exact solution has been found in one dimension for a discrete coagulation model with one fixed source. The latter solution is related to the probability that a given spin in an Ising chain with random initial conditions does not change its value before time $`t`$ . For discrete and continuous coagulation models, exact results are available not only for the density but also for the spectrum of relaxation rates,  the distribution of interparticle distances,  and correlation functions.  They are obtained by considering the function $`E(n\mathrm{\Delta }x,t)`$, defined as the probability that an arbitrarily chosen segment of $`n`$ consecutive sites contains no particles, satisfying a closed kinetic equation. It has, however, not proven possible to extend this method to the case of annihilation on contact, because the function $`E(n\mathrm{\Delta }x,t)`$ does not satisfy a closed equation .

### 1.2 Paired nucleation

The “coefficient of recombination” of two particles initially close together was introduced in the study of subatomic particles . The relative motion of two diffusing particles is equivalent to a problem of Brownian motion of one particle.

A discrete model that corresponds to paired nucleation is the Ising model, with nucleation at neighboring sites. Its dynamics was studied analytically by Glauber in 1963 ; the nucleation rate is proportional to the square of the steady state density for nucleation rates sufficiently small that excluded volume effects can be neglected . Computer simulations of a discretized reaction-diffusion model $`A+B0`$, published in 1987 , contrasted the scalings of the steady state density according to whether nucleation occurred at random sites or in pairs at neighboring sites. In the latter case, the scaling $`\mathrm{\Gamma }\rho ^2`$ was found.

A different approach to diffusion-limited reaction was recently introduced in the context of kink dynamics in a stochastic partial differential equation (PDE) . There, the dynamics was termed “mesoscopic” because it was an approximate model that ignored the internal structure of kinks and antikinks, treating them simply as particles that happen to be nucleated in pairs. The treatment was based on classifying particles according to whether they are annihilated in a collision with their nucleation partner (recombination) or with a different particle (nonrecombinant annihilation). The steady-state density $`\rho _0`$ is related to the mean lifetime of a particle, $`\tau `$, by

$$\rho _0=2\mathrm{\Gamma }\tau .$$
(6)
The mean lifetime $`\tau `$ was estimated directly by averaging over the possible histories of a pair of particles born together. This approximate analysis yielded the estimate $`\rho _0=(3b\mathrm{\Gamma }/8D)^{\frac{1}{2}}`$.

## 2 Hierarchy of distribution functions

Let $`f_n(x_1,\mathrm{}x_n;t)\mathrm{d}x_1\mathrm{}\mathrm{d}x_n`$ be the probability that there is one particle in $`(x_1,x_1+\mathrm{d}x_1)`$, one in $`(x_2,x_2+\mathrm{d}x_2)`$, $`\mathrm{}`$, and one in $`(x_n,x_n+\mathrm{d}x_n)`$ at time $`t`$, regardless of the positions of the other particles . The function $`f_1(x_1;t)`$ is the particle density at $`x_1`$ at time $`t`$. On deriving the differential equation for its time derivative, one finds that it involves $`f_2(x_1,x_2;t)`$ . Similarly, the time derivative of $`f_2(x_1,x_2;t)`$ involves $`f_3(x_1,x_2,x_3;t)`$. One is thus led to a hierarchy of differential equations for the evolution of the distribution functions.

In this section we derive the source terms appropriate for paired nucleation in the hierarchy of differential equations. We also derive the reaction terms corresponding to diffusion with annihilation on collision, without needing to introduce a reaction radius. Three parameters remain in the theory for paired nucleation: the nucleation rate of pairs $`\mathrm{\Gamma }`$, their separation at nucleation $`b`$ and the diffusivity of a particle $`D`$. For unpaired nucleation there are two parameters: the nucleation rate $`Q`$ and the diffusivity of a particle $`D`$. The annihilation process is immediate on collision and therefore does not require extra parameters. It manifests itself instead in boundary conditions on the distribution functions. We shall truncate the hierarchy of distribution functions using an ansatz for the three-point correlation function, introduced in the literature for unpaired nucleation , thus obtaining a closed pair of differential equations for the density and two-point correlation function. Their solution yields analytical approximations for the steady-state density and two-point correlation function. By examining perturbations away from the steady state, we derive the time scales for relaxation towards the steady state.

The evolution of the reduced distribution functions has a number of contributions

$`{\displaystyle \frac{}{t}}f_n(x_1,\mathrm{}x_n;t)`$ $`=`$ $`D^2f_n(x_1,\mathrm{}x_n;t)`$ (7)
$`{\displaystyle \underset{(ij)}{\overset{n}{}}}k(x_i,x_j)f_n(x_1,\mathrm{}x_n;t)`$
$`{\displaystyle \underset{i=1}{\overset{n}{}}}{\displaystyle _L^L}dx_{n+1}k(x_i,x_{n+1})f_{n+1}(x_1,\mathrm{}x_n,x_{n+1};t)`$
$`+\text{sources}.`$
First on the right is the diffusion term, due to the motion of each particle with diffusion coefficient $`D`$. The second term represents the reaction between two of the $`n`$ particles: $`k(x,x^{})`$ is the probability per unit time that a particle at $`x`$ and one at $`x^{}`$ react, and the summation is over all pairs that can be selected from the $`n`$ particles. The third term accounts for the fact that each of the $`n`$ particles may react with another that is not part of the set of $`n`$ particles, and $`2L`$ is the size of the system. The last term is a source contribution whose form is given in detail below. Equation (7) is one in an infinite hierarchy. Explicitly, the first two equations in the hierarchy are:

$`{\displaystyle \frac{}{t}}f_1(x_1;t)`$ $`=`$ $`D{\displaystyle \frac{^2}{x_1^2}}f_1(x_1;t){\displaystyle _L^L}dx^{}k(x_1,x^{})f_2(x_1,x^{};t)+q_1(x_1),`$
$`{\displaystyle \frac{}{t}}f_2(x_1,x_2;t)`$ $`=`$ $`D^2f_2(x_1,x_2;t)k(x_1,x_2)f_2(x_1,x_2;t)`$ (9)
$`{\displaystyle _L^L}dx^{}[k(x_1,x^{})+k(x_2,x^{})]f_3(x_1,x_2,x^{};t)`$
$`+q_2(x_1,x_2)+f_1(x_1;t)q_1(x_2)+f_1(x_2;t)q_1(x_1).`$
### 2.1 Source terms for paired nucleation

The term $`q_1(x_1)`$ in (2) and (9) is the probability density per unit time for the creation of a particle at $`x_1`$; the term $`q_2(x_1,x_2)`$ is the probability density per unit time for the simultaneous creation of a particle at $`x_1`$ and another at $`x_2`$. When creation of particles always occurs in pairs, these two source functions are related

$$q_1(x_1)=_L^Ldx_2q_2(x_1,x_2).$$
(10)
When the particle creation rates are independent of position and time, $`q_1(x)`$ is constant

$$q_1(x)=2\mathrm{\Gamma },$$
(11)
and $`q_2(x,x+y)`$ is independent of $`x`$. The constant $`\mathrm{\Gamma }`$ is the rate of creation of pairs per unit length.

Here, because particles are indistinguishable, $`q_2(x_1,x_2)`$ depends only on $`y=|x_1x_2|`$ and the functions $`f_n(x_1,\mathrm{},x_i,\mathrm{},x_n;t)`$ are independent of the order of the $`x_i`$. The probability that two particles initially at $`x_1`$ and $`x_2`$ react is the probability that they diffuse and collide. Since particle diffusion is isotropic and independent of position, $`k(x_1,x_2)`$ also depends only on $`y=|x_1x_2|`$. We therefore define

$$q(y)q_2(x_1,x_2)\text{ and }K(y)k(x_1,x_2).$$
(12)
The function $`q(y)`$ describes the probability density of distances between particles nucleated simultaneously. We shall use the following forms for this function, corresponding to unpaired nucleation and to paired nucleation with initial separation $`b`$:

$`q(y)=\{\begin{array}{ccc}\mathrm{\Gamma }/L\hfill & \text{unpaired};\hfill & \\ \mathrm{\Gamma }\delta (yb)\hfill & \text{paired}.\hfill & \end{array}`$ (15)
We can now rewrite (2) and (9) as follows:

$`{\displaystyle \frac{}{t}}f_1(x;t)`$ $`=`$ $`D{\displaystyle \frac{^2}{x^2}}f_1(x;t)2{\displaystyle _0^L}dyK(y)f_2(x,x+y;t)+2\mathrm{\Gamma },`$
$`{\displaystyle \frac{}{t}}f_2(x,x+y;t)`$ $`=`$ $`D^2f_2(x,x+y;t)K(y)f_2(x,x+y;t)`$ (17)
$`{\displaystyle _L^L}dz\left(K(|z|)+K(|zy|)\right)f_3(x,x+y,x+z;t)`$
$`+q(y)+2\mathrm{\Gamma }(f_1(x;t)+f_1(x+y;t)).`$
If the initial conditions are homogeneous, then the functions $`f_1`$, $`f_2`$, $`\mathrm{}`$ will be homogeneous at all times. In particular, $`f_1(x;t)`$ will be independent of $`x`$ at every $`t`$. Let

$`\rho (t)`$ $``$ $`f_1(x;t),`$
$`F_n(x_2,\mathrm{},x_n;t)`$ $``$ $`\rho ^n(t)f_n(x,x_2x,\mathrm{},x_nx;t).`$ (18)
We shall in particular be interested in the dimensionless correlation function defined by

$$g(y,t)F_2(y;t).$$
(19)
The function $`g(y,t)`$ is the probability density at time $`t`$ of particles at a distance $`y`$ from a reference particle, divided by the overall density of particles. It is constructed numerically as follows. Choose a sample of $`N`$ reference particles, located at $`\{x_i,i=1,\mathrm{},N\}`$ at time $`t`$. For each $`x_i`$, construct

$$G_i(y,t)=\{\text{number of particles between}x_i\text{and}x_i+y\}$$

for $`y>0`$. Then $`G(y,t)`$ is the average over the $`N`$ particles of the $`G_i(y,t)`$ and

$$g(y,t)=(\rho (t))^1\frac{}{y}G(y,t).$$
(20)
If there is no correlation between particles at time $`t`$, then $`g(y,t)`$=1 for all $`y0`$. In all the situations considered here, the total length $`2L`$ of the system is sufficently large compared to the correlation length so that

$$\underset{y\mathrm{}}{lim}g_0(y)=1,$$
(21)
where $`g_0(y)`$ denotes the steady state correlation function.

In terms of $`\rho (t)`$ and $`g(y,t)`$, Eqs. (2.1) and (17) now simplify to the pair of equations

$`\dot{\rho }(t)`$ $`=`$ $`2\rho ^2(t){\displaystyle _0^L}dyK(y)g(y,t)+2\mathrm{\Gamma },`$ (22)
$`{\displaystyle \frac{}{t}}g(y,t)`$ $`=`$ $`2D{\displaystyle \frac{^2}{y^2}}g(y,t)K(y)g(y,t)`$ (23)
$`\rho (t){\displaystyle _L^L}dz\left(K(|z|)+K(|zy|)\right)F_3(y,z;t)`$
$`+4\mathrm{\Gamma }\rho ^1(t)+\rho ^2(t)q(y)2\rho ^1(t)g(y,t)\dot{\rho }(t).`$
### 2.2 The reaction kernel

To complete the description of the dynamics of the system, we consider the reaction terms for the case where particles diffuse with diffusivity $`D`$ and annihilate on collision. We shall see that a consequence of annihilation on collision is that $`g(0)=0`$ for all $`t>0`$, where $`g(y)`$ is the correlation function defined in (19). We derive an exact relation between $`g^{}(0)`$ and the rate of collisions between particles.

Let $`s(y,\mathrm{\Delta }t)`$ be the probability that two particles, with initial separation $`y`$, collide before $`\mathrm{\Delta }t`$. Then the reaction kernel $`K(y)`$, defined in (12), is given by

$$K(y)=\underset{\mathrm{\Delta }t0}{lim}\frac{1}{\mathrm{\Delta }t}s(y,\mathrm{\Delta }t).$$
(24)
If both particles diffuse with diffusivity $`D`$ then

$$s(y,\mathrm{\Delta }t)=\mathrm{erfc}\left(\frac{y}{(8D\mathrm{\Delta }t)^{\frac{1}{2}}}\right),$$
(25)
where we assume $`L(D\mathrm{\Delta }t)^{\frac{1}{2}}`$. To calculate the frequency of collisions between particles, we consider a time interval $`t,t+\mathrm{\Delta }t`$. Given the density $`\rho (t)`$ and the correlation function $`g(y,t)`$ defined in (19), we can imagine following the paths of all the particles from time $`t`$ to time $`t+\mathrm{\Delta }t`$ without removing those that collide. Then, the probability that a particle chosen at random undergoes a collision between time $`t`$ and time $`t+\mathrm{\Delta }t`$ is $`P(t,\mathrm{\Delta }t)`$, where

$$P(t,\mathrm{\Delta }t)=2\rho (t)_0^Ldys(y,\mathrm{\Delta }t)g(y,t).$$
(26)
The expression (26) overestimates the number of collisions in the system with annihilation on collision due to the possibility that the same particle undergoes two (or more) collisions in the interval $`t,t+\mathrm{\Delta }t`$. However, this latter probability is proportional to $`(\mathrm{\Delta }t)^2`$ as $`\mathrm{\Delta }t0`$, and so (26) is valid for our system in the limit $`\mathrm{\Delta }t0`$.

Next, consider the dynamics of the system as a whole. The mean number of distinct collisions between time $`t`$ and time $`t+\mathrm{\Delta }t`$ is given by $`L\rho (t)P(t,\mathrm{\Delta }t)`$. We can, therefore, write

$`{\displaystyle _0^L}dys(y,\mathrm{\Delta }t)g(y,t)`$
$`=`$ $`{\displaystyle _0^L}dys(y,\mathrm{\Delta }t)\left(g(0,t)+yg^{}(0^+,t)+{\displaystyle \frac{1}{2}}y^2g^{\prime \prime }(0^+,t)+\mathrm{}\right)`$
$`=`$ $`\left({\displaystyle \frac{1}{\sqrt{\pi }}}(8D\mathrm{\Delta }t)^{\frac{1}{2}}g(0,t)+{\displaystyle \frac{1}{4}}(8D\mathrm{\Delta }t)g^{}(0^+,t)+𝒪(D\mathrm{\Delta }t)^{\frac{3}{2}}\right).`$
The number of collisions between time $`t`$ and $`t+\mathrm{\Delta }t`$ is proportional to $`\mathrm{\Delta }t`$ if

$$g(0,t)=0.$$
(28)
Because the number of nucleation events between time $`t`$ and time $`t+\mathrm{\Delta }t`$ is proportional to $`2L\mathrm{\Delta }t`$, the condition (28) is necessary if there is to be a steady-state balance between nucleation and annihilation. More generally, it is necessary if $`\rho (t)`$ is to obey a differential equation. It is, of course, possible to construct initial conditions that do not satisfy (28): a random distribution of particles, for example. Then the number of annihilation events will initially be proportional to $`t^{\frac{1}{2}}`$; this period of rapid annihilation creates a “depletion zone”  in $`g(y,t)`$, which thereafter satisfies (28). That $`g(y,t)<1`$ for $`y0`$ implies an effective repulsion: particles are less likely to be found close to a reference particle than a large distance from it.

Using (LABEL:kercalc) and (28) gives exact expressions for the evolution of the density:

$$\{\begin{array}{ccc}\dot{\rho }(t)\hfill & =\hfill & 4D\rho ^2(t)g^{}(0^+,t)+2\mathrm{\Gamma },\text{paired};\hfill \\ \dot{\rho }(t)\hfill & =\hfill & 4D\rho ^2(t)g^{}(0^+,t)+Q,\text{unpaired}.\hfill \end{array}$$
(29)
In particular, we have the following relationship between the steady state density and the derivative of the correlation function. Let $`\rho _0`$ and $`g_0(y)`$ denote the steady state density and correlation function. Then

$$\{\begin{array}{ccc}\mathrm{\Gamma }=2D\rho _0^2g_0^{}(0^+)\text{paired};\hfill & & \\ Q=4D\rho _0^2g_0^{}(0^+)\text{unpaired}.\hfill & & \end{array}$$
(30)
The reaction kernel $`K(y)`$ is a singular function

$$_0^LdyK(y)g(y,t)=\underset{\mathrm{\Delta }t0}{lim}\frac{1}{\mathrm{\Delta }t}_0^Ldys(y,\mathrm{\Delta }t)g(y,t)=2Dg^{}(0^+,t).$$
(31)
We have assumed that (28) holds. The derivative of $`g(y)`$ is one-sided

$$g^{}(0^+,t)\underset{a0^+}{lim}\frac{g(a,t)}{a},$$
(32)
because $`g(y)`$ is only defined for $`y>0`$. In other words

$$K(y)=2D\underset{a0^+}{lim}\frac{\delta (ya)}{a}.$$
(33)
It is interesting to compare (33) with the form of $`K(y)`$ used, for example, by Lindenberg et al.

$$K(y)=k\delta (ya),$$
(34)
which introduced a reaction radius $`a`$ and a rate coefficient $`k`$. There it was assumed that these constants are connected to the diffusivity via the Smoluchowski relation $`k=2D/a`$ and that in the limit $`a0`$, $`k\mathrm{}`$ their product remains finite. In the form used here, by contrast, we are able to explicitly take the limit of zero reaction radius: $`a0`$.

Similarly to (31)

$$_L^LdzK(|z|)F_3(y,z;t)=4DF_3^{}(y,0^+;t),$$
(35)
where

$$F_3^{}(y,0^+;t)=\underset{a0^+}{lim}a^1F_3(y,a;t),$$
(36)
and

$$_L^LdzK(|zy|)F_3(y,z;t)=4DF_3^{}(y,y^+;t).$$
(37)
A similar expression was derived for a discrete coagulation model without nucleation by Lin, Doering, and ben-Avraham . Since $`F_3(y,0;t)=F_3(y,y;t)`$, Eqs. (22) and (23) simplify to the pair of equations

$`\dot{\rho }(t)`$ $`=`$ $`4D\rho ^2(t)g^{}(0^+,t)+2\mathrm{\Gamma },`$ (38)
$`{\displaystyle \frac{}{t}}g(y,t)`$ $`=`$ $`2D{\displaystyle \frac{^2}{y^2}}g(y,t)K(y)g(y,t)8D\rho (t)F_3^{}(y,0^+;t)`$ (39)
$`+4\mathrm{\Gamma }\rho ^1(t)+\rho ^2(t)q(y)2\rho ^1(t)g(y,t)\dot{\rho }(t).`$
Annihilation on collision is described by the terms involving the reaction kernel $`K(y)`$.

### 2.3 Truncation of the hierarchy

We have obtained exact expressions for the evolution of the density. However, to obtain a closed set of equations, we truncate the hierarchy of distribution functions via an approximation. Various methods have been used to break hierarchies resulting from reaction-diffusion systems . We shall restrict ourselves to the simplest. In the hierarchy that begins with (38) and (39) we make the ansatz

$$F_3^{}(y,0^+;t)=g(y,t)g^{}(0^+,t).$$
(40)
This choice, which would be exact if successive interparticle spacings were independent , is not per se the most compelling, but it has been shown to produce excellent results (when compared with simulations) for batch reactions and in the steady state with unpaired nucleation . In Sec. 3 we shall compare the steady-state density obtained with this closure to the exact result.

With the approximation (40) we find $`4D\rho ^2(t)F_3^{}(y,0^+)=g(y,t)(\dot{\rho }(t)2\mathrm{\Gamma })`$, and so (38) and (39) reduce to the following closed set of equations, linear in $`g(y,t)`$:

$`\dot{\rho }(t)`$ $`=`$ $`4D\rho ^2(t)g^{}(0^+,t)+2\mathrm{\Gamma },`$ (41)
$`{\displaystyle \frac{}{t}}g(y,t)`$ $`=`$ $`2D{\displaystyle \frac{^2}{y^2}}g(y,t)2K(y)g(y,t)+{\displaystyle \frac{4\mathrm{\Gamma }}{\rho (t)}}(1g(y,t))+{\displaystyle \frac{q(y)}{\rho ^2(t)}},`$
plus the condition (28). In the case of paired nucleation, $`q(y)=\mathrm{\Gamma }\delta (yb)`$. In the case of unpaired nucleation in the (thermodynamic) limit, $`L\rho (t)^1`$, (LABEL:gdott) reduces to

$`{\displaystyle \frac{}{t}}g(y,t)`$ $`=`$ $`2D{\displaystyle \frac{^2}{y^2}}g(y,t)2K(y)g(y,t)+{\displaystyle \frac{2Q}{\rho (t)}}\left(1g(y,t)\right),`$
with, as before, $`Q2\mathrm{\Gamma }`$. Note that no further simplifications can be obtained by assuming low density. In particular, a low density expansion cannot be used to justify the truncation (40).

### 2.4 Steady states

The density and correlation function in the steady state, $`\rho _0`$ and $`g_0(y)`$, are found by setting to zero the time derivatives on the left-hand side of (41) and of (LABEL:gdott) or (LABEL:pairnoaunc). There thus results a second-order equation for $`g_0(y)`$, with the two relations (28) and (30).

For unpaired nucleation one finds

$`\rho _0^{\mathrm{tu}}`$ $`=`$ $`\left({\displaystyle \frac{Q}{16D}}\right)^{\frac{1}{3}},`$ (44)
$`g_0^{\mathrm{tu}}(y)`$ $`=`$ $`1e^{2(Q/2D)^{1/3}y},`$ (45)
where we have introduced the superscript “t” to indicate that the result is obtained from the truncation (40) and “u” to denote “unpaired nucleation.” In Fig. 2 we compare numerical results for the correlation function with (45).

For paired nucleation the steady-state equation for the correlation function is

$`0`$ $`=`$ $`2D{\displaystyle \frac{^2}{y^2}}g_0(y)2K(y)g_0(y)+4{\displaystyle \frac{\mathrm{\Gamma }}{\rho _0}}[1g_0(y)]+{\displaystyle \frac{\mathrm{\Gamma }}{\rho _0^2}}\delta (yb).`$ (46)
The solution of (46) is derived in Appendix A. When $`b`$ is sufficiently large the results are equivalent to (45) with (44). Of interest here is the opposite situation: $`\epsilon 0`$, with $`ϵ`$ defined in Eq. (1). The separation $`b`$ in the latter case is much smaller than the length scale defined by the inverse density, and the steady-state density is given by

$$\rho _0^{\mathrm{tp}}=\left(\frac{\mathrm{\Gamma }b}{2D}\right)^{1/2}.$$
(47)
The correlation function in the same limit is

$$g_0^{\mathrm{tp}}(y)=\{\begin{array}{ccc}\frac{y}{b}\hfill & 0y<b;\hfill & \\ 1\hfill & yb.\hfill & \end{array}$$
(48)
Corrections to (48) are proportional to $`\epsilon ^{\frac{3}{4}}`$. In Fig. 3 this correlation function is compared with numerical results.

### 2.5 Relaxation to the steady state

In order to study the relaxation to the steady state, we decompose the functions $`\rho (t)`$ and $`g(y,t)`$ as follows:

$`\rho (t)`$ $`=`$ $`\rho _0+\delta \rho (t),`$ (49)
$`g(y,t)`$ $`=`$ $`g_0(y)+\delta g(y,t),`$ (50)
with $`\rho _0`$ and $`g_0(y)`$ the steady-state density and the steady-state correlation function, respectively. This decomposition is valid for both unpaired and paired nucleation. Assuming that we are close to the steady state, we can obtain linearized equations for the deviations $`\delta \rho `$ and $`\delta g`$ from their steady-state values. For paired nucleation

$`{\displaystyle \frac{}{t}}\delta \rho (t)`$ $`=`$ $`4D\rho _0^2\delta g^{}(0^+,t)8D\rho _0g_0^{}(0^+)\delta \rho (t),`$ (51)
$`{\displaystyle \frac{}{t}}\delta g(y,t)`$ $`=`$ $`2D{\displaystyle \frac{^2}{y^2}}\delta g(y,t)2K(y)\delta g(y,t)`$
$`{\displaystyle \frac{4\mathrm{\Gamma }}{\rho _0}}\left[\delta g(y,t)+{\displaystyle \frac{1g_0(y)}{\rho _0}}\delta \rho (t)\right]{\displaystyle \frac{2\mathrm{\Gamma }}{\rho _0^3}}\delta (yb)\delta \rho (t).`$
For unpaired nucleation the last term in (LABEL:lgho1) is absent.

Formal solution of these coupled linear equations is presented in Appendix B. Explicit solution for all times is in fact possible for the unpaired nucleation case (and presented in the appendix). Ultimately we are interested in their asymptotic relaxation behavior. If the relaxation processes each involve a single exponential decay

$`\delta \rho (t)`$ $`\stackrel{t\mathrm{}}{}`$ $`Ae^{\alpha t},`$ (53)
$`\delta g^{}(0^+,t)`$ $`\stackrel{t\mathrm{}}{}`$ $`Be^{\beta t},`$ (54)
then the density and correlation function decay on the same time scale, i.e., $`\beta =\alpha `$. For unpaired nucleation we find from the exact result (B.22) that the asymptotic decay is indeed exponential, with

$`\alpha ^\mathrm{u}`$ $`=`$ $`(5+\sqrt{5})(D\mathrm{\Gamma }^2)^{1/3}=7.236\mathrm{}(D\mathrm{\Gamma }^2)^{1/3}`$ (55)
$`=`$ $`{\displaystyle \frac{5+\sqrt{5}}{4^{1/3}}}(DQ^2)^{1/3}=4.558\mathrm{}(DQ^2)^{1/3}.`$
In the case of paired nucleation, if we assume exponential decay we find for the inverse time scale

$$\alpha ^\mathrm{p}=\left(\frac{32D\mathrm{\Gamma }}{b}\right)^{1/2}.$$
(56)
However, on the basis of the exact results reported below (and also in parallel work ) there is reason to suspect that the decay may not be purely exponential in the paired nucleation case.

## 3 Exact results

In this section we derive exact expressions for the density of particles, using a function that obeys a linear partial differential equation. The function is similar in interpretation to the pair-pair correlation function in the Ising model . The methodology is also similar to that used to obtain exact results for models of diffusion-limited coagulation . Here, we obtain explicit exact expressions for the density of particles, in steady state and nonsteady state, for paired and unpaired nucleation. Previous exact results for diffusion-limited reaction with annihilation have been limited to the case of no nucleation .

Let the function $`r(x,t)`$ be defined as follows:

$`r(x,t)=`$ $`\{\text{probability that the number of particles}`$
$`\text{between }0\text{ and }x\text{ at time }t\text{ is even}\}.`$
Note that, by translational invariance, we can replace the interval $`(0,x)`$ by $`(X,X+x)`$ for any $`X`$. The value of $`r(x,t)`$ changes due to diffusion of particles in or out of of the region $`(0,x)`$, and due to nucleation of a single particle in the region.

In Appendix C we derive equations for the space and time derivatives of $`r(x,t)`$. At any time $`t`$ the density $`\rho (t)`$ is given by

$$\rho (t)=\frac{}{x}r(x,t)|_{x=0^+}.$$
(58)
To describe the time evolution of $`r(x,t)`$, we distinguish between the cases of unpaired and paired nucleation using the superscripts $`\mathrm{u}`$ and $`\mathrm{p}`$.

In the case of unpaired nucleation, $`r^\mathrm{u}(x,t)`$ satisfies

$`{\displaystyle \frac{}{t}}r^\mathrm{u}(x,t)`$ $`=`$ $`2D{\displaystyle \frac{^2}{x^2}}r^\mathrm{u}(x,t)xQr^\mathrm{u}(x,t)+xQ(1r^\mathrm{u}(x,t))`$ (59)
$`=`$ $`2D{\displaystyle \frac{^2}{x^2}}r^\mathrm{u}(x,t)+xQ(12r^\mathrm{u}(x,t)),`$
with the boundary conditions

$$r^\mathrm{u}(0,t)=1\text{and}\underset{x\mathrm{}}{lim}r^\mathrm{u}(x,t)=\frac{1}{2},t>0.$$
(60)
In the case of paired nucleation $`r^\mathrm{p}(x,t)`$ satisfies

$`{\displaystyle \frac{}{t}}r^\mathrm{p}(x,t)=\{\begin{array}{ccc}2D{\displaystyle \frac{^2}{x^2}}r^\mathrm{p}(x,t)+2x\mathrm{\Gamma }(12r^\mathrm{p}(x,t))\hfill & xb;\hfill & \\ 2D{\displaystyle \frac{^2}{x^2}}r^\mathrm{p}(x,t)+2b\mathrm{\Gamma }(12r^\mathrm{p}(x,t))\hfill & x>b,\hfill & \end{array}`$ (63)
with the boundary conditions

$$r^\mathrm{p}(0,t)=1\text{and}\underset{x\mathrm{}}{lim}r^\mathrm{p}(x,t)=\frac{1}{2},t>0.$$
(64)
### 3.1 Steady state: Unpaired nucleation

The steady-state solution of (59) will be denoted by $`r_0^\mathrm{u}(x)`$. It satisfies

$$2D\frac{^2}{x^2}r_0^\mathrm{u}(x)+xQ(12r_0^\mathrm{u}(x))=0.$$
(65)
The solution is

$$r_0^\mathrm{u}(x)=\frac{1}{2}\left(1+\frac{\mathrm{Ai}((Q/D)^{\frac{1}{3}}x)}{\mathrm{Ai}(0)}\right)$$
(66)
and is shown in Fig. 4. Thus the exact steady-state density for unpaired nucleation is

$$\rho _0^\mathrm{u}=\frac{}{x}r_0^\mathrm{u}(x)|_{x=0^+}=\frac{1}{2}\left(\frac{Q}{D}\right)^{\frac{1}{3}}\frac{|\mathrm{Ai}^{}(0)|}{\mathrm{Ai}(0)}=\left(\frac{Q}{16D}\right)^{\frac{1}{3}}0.9186\mathrm{}.$$
(67)
Note that the exact result for the steady-state density is $`0.9186`$ of the density (44) predicted from the truncated hierarchy. For comparison, in a discrete model with nucleation rate $`R`$, where collision of particles produces coagulation rather than annihilation, the steady-state density is given by

$$\rho _0^\mathrm{u}=\frac{1}{2}\left(\frac{R}{2D}\right)^{\frac{1}{3}}\frac{|\mathrm{Ai}^{}(0)|}{\mathrm{Ai}(0)}\text{(coagulation)}.$$
(68)
### 3.2 Steady state: Paired nucleation

The steady-state solution of (63) is

$$r_0^\mathrm{p}(x)=\{\begin{array}{ccc}\frac{1}{2}\left[c_1\mathrm{Ai}\left(\left(\frac{2\mathrm{\Gamma }}{D}\right)^{\frac{1}{3}}x\right)+c_2\mathrm{Bi}\left(\left(\frac{2\mathrm{\Gamma }}{D}\right)^{\frac{1}{3}}x\right)+1\right]\hfill & xb;\hfill & \\ \frac{1}{2}\left[c_3\mathrm{exp}((\frac{2\mathrm{\Gamma }b}{D})^{\frac{1}{2}}x)+1\right]\hfill & x>b\hfill & \end{array}$$
(69)
and is shown in Fig. 5. We have used the second of the boundary conditions (64) to rule out increasing exponential solutions for $`x>b`$. The constants $`c_1`$, $`c_2`$, and $`c_3`$ are fixed by requiring $`r_0^\mathrm{p}(0)=0`$ and imposing continuity of $`r_0^\mathrm{p}(x)`$ and $`\frac{\mathrm{d}}{\mathrm{d}x}r_0^\mathrm{p}(x)`$ at $`x=b`$.

The density $`\rho (t)`$ is given by (58). In the steady state

$`\rho _0^\mathrm{p}`$ $`=`$ $`{\displaystyle \frac{\mathrm{d}}{\mathrm{d}x}}r_0^\mathrm{p}(x)|_{x=0^+}`$ (70)
$`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{2\mathrm{\Gamma }}{D}}\right)^{\frac{1}{3}}(c_1\mathrm{Ai}^{}(0)+c_2\mathrm{Bi}^{}(0))`$
$`=`$ $`\left({\displaystyle \frac{\mathrm{\Gamma }}{4D}}\right)^{\frac{1}{3}}{\displaystyle \frac{|\mathrm{Ai}^{}(0)|}{\mathrm{Ai}(0)}}\left({\displaystyle \frac{\mathrm{Bi}^{}(\epsilon )+\sqrt{3\mathrm{A}\mathrm{i}^{}}(\epsilon )+\epsilon (\mathrm{Bi}(\epsilon )+\sqrt{3\mathrm{A}\mathrm{i}}(\epsilon ))}{\mathrm{Bi}^{}(\epsilon )\sqrt{3\mathrm{A}\mathrm{i}^{}}(\epsilon )+\epsilon (\mathrm{Bi}(\epsilon )\sqrt{3\mathrm{A}\mathrm{i}}(\epsilon ))}}\right),`$
in terms of the dimensionless quantity $`\epsilon `$ defined in Eq. (1). The function (70) is plotted in Fig. 6.

In the limit $`\epsilon 0`$, $`\mathrm{Bi}(\epsilon )\sqrt{3\mathrm{A}\mathrm{i}}(\epsilon )`$ and $`\mathrm{Bi}^{}(\epsilon )\sqrt{3\mathrm{A}\mathrm{i}^{}}(\epsilon )`$, so

$`\rho _0^\mathrm{p}`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{\Gamma }}{4D}}\right)^{\frac{1}{3}}\left(\epsilon ^{\frac{1}{2}}{\displaystyle \frac{1}{2}}\epsilon ^2+\mathrm{}\right)`$ (71)
$`=`$ $`\left({\displaystyle \frac{b\mathrm{\Gamma }}{2D}}\right)^{\frac{1}{2}}(1+𝒪(\epsilon ^{\frac{3}{2}})).`$
For $`\epsilon \mathrm{}`$

$$\rho _0^\mathrm{p}\frac{1}{2}\left(\frac{2\mathrm{\Gamma }}{D}\right)^{\frac{1}{3}}\frac{|\mathrm{Ai}^{}(0)|}{\mathrm{Ai}(0)},$$
(72)
and we regain the result (67).

### 3.3 Time-dependent statistics: Unpaired nucleation

Let us introduce

$$h^\mathrm{u}(x,t)=r^\mathrm{u}(x,t)r_0^\mathrm{u}(x).$$
(73)
Then, $`h^\mathrm{u}(x,t)`$ satisfies

$$\frac{}{t}h^\mathrm{u}(x,t)=2D\frac{^2}{x^2}h^\mathrm{u}(x,t)2Qxh^\mathrm{u}(x,t),$$
(74)
with the boundary conditions

$$h^\mathrm{u}(0,t)=0\text{and}\underset{x\mathrm{}}{lim}h^\mathrm{u}(x,t)=0,$$
(75)
for all $`t`$. We can expand the general solution as follows:

$$h^\mathrm{u}(x,t)=\underset{i=1}{\overset{\mathrm{}}{}}c_ih_i^\mathrm{u}(x)\mathrm{e}^{\lambda _it},$$
(76)
where the eigenfunctions $`h_i^\mathrm{u}(x)`$ satisfy

$$\frac{\mathrm{d}^2}{\mathrm{d}x^2}h_i^\mathrm{u}(x)\frac{Q}{D}xh_i^\mathrm{u}(x)=\lambda _ih_i^\mathrm{u}(x),$$
(77)
with the boundary conditions

$$h_i^\mathrm{u}(0)=0\text{and}\underset{x\mathrm{}}{lim}h_i^\mathrm{u}(x)=0.$$
(78)
The eigenfunctions $`h_i^\mathrm{u}(x)`$ are thus given by

$$h_i^\mathrm{u}(x)=N_i\mathrm{Ai}\left(\left(\frac{Q}{D}\right)^{\frac{1}{3}}(x\frac{\lambda _i}{2Q})\right).$$
(79)
The eigenfunctions are normalized by choosing

$$N_i^2=_0^{\mathrm{}}dx\mathrm{Ai}^2\left(\left(\frac{Q}{D}\right)^{\frac{1}{3}}(x\frac{\lambda _i}{2Q})\right)=\left(\frac{D}{Q}\right)^{\frac{1}{3}}_{a_i}^{\mathrm{}}dz\mathrm{Ai}^2(z).$$
(80)
The eigenvalues $`\lambda _i`$ are related to the zeros of the Airy function (all on the negative real axis)

$$\lambda _i=(8DQ^2)^{\frac{1}{3}}a_i,$$
(81)
where $`a_i`$ is the $`i`$th zero counting away from $`0`$. Relaxation towards the steady state is determined for late times by the smallest eigenvalue

$$\lambda _1=(8DQ^2)^{\frac{1}{3}}a_1.$$
(82)
An explicit analytical solution for the density as a function of time is obtained once the constants $`c_i`$ are determined from the initial condition $`r^\mathrm{u}(x,0)`$

$$c_i=_0^{\mathrm{}}dxh_i^\mathrm{u}(x)(r^\mathrm{u}(x,0)r_0^\mathrm{u}(x)).$$
(83)
Thus,

$`h^\mathrm{u}(x,t)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{_0^{\mathrm{}}dz\mathrm{Ai}(z+a_i)\left[r(z(\frac{D}{Q})^{\frac{1}{3}},0)\frac{1}{2}\left(1+\frac{\mathrm{Ai}\left(z\right)}{\mathrm{Ai}\left(0\right)}\right)\right]}{_{a_i}^{\mathrm{}}dz\mathrm{Ai}^2(z)}}`$ (84)
$`\times \mathrm{Ai}\left(({\displaystyle \frac{Q}{D}})^{\frac{1}{3}}(x{\displaystyle \frac{\lambda _i}{2Q}})\right)\mathrm{e}^{\lambda _it},`$
and

$`\rho (t)`$ $`=`$ $`\rho _0^\mathrm{u}{\displaystyle \frac{}{x}}h^\mathrm{u}(x,t)|_{x=0^+}`$ (86)
$`=`$ $`\rho _0^\mathrm{u}{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{Q}{D}}\right)^{\frac{1}{3}}{\displaystyle \underset{i=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{_0^{\mathrm{}}dz\mathrm{Ai}(z+a_i)\left[r^\mathrm{u}(z(\frac{D}{Q})^{\frac{1}{3}},0)\frac{1}{2}\left(1+\frac{\mathrm{Ai}\left(z\right)}{\mathrm{Ai}\left(0\right)}\right)\right]}{_{a_i}^{\mathrm{}}dz\mathrm{Ai}^2(z)}}`$
$`\times \mathrm{Ai}^{}(a_i)\mathrm{e}^{\lambda _it}.`$
#### 3.3.1 Zero initial density

If $`\rho (0)=0`$, then $`r^\mathrm{u}(x,0)=1`$ for all $`x>0`$ and

$$c_i=\frac{1}{2}N_i\left(\frac{D}{Q}\right)^{\frac{1}{3}}_0^{\mathrm{}}dz\mathrm{Ai}(z+a_i)\left(1\frac{\mathrm{Ai}(z)}{\mathrm{Ai}(0)}\right).$$
(87)
Thus,

$`\rho (t)`$ $`=`$ $`\rho _0^\mathrm{u}{\displaystyle \frac{}{x}}h^\mathrm{u}(x,t)|_{x=0^+}`$
$`=`$ $`\rho _0^\mathrm{u}`$
$``$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{Q}{D}}\right)^{\frac{1}{3}}{\displaystyle \underset{i=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{_0^{\mathrm{}}dz\mathrm{Ai}(za_i)\left(1\frac{\mathrm{Ai}\left(z\right)}{\mathrm{Ai}\left(0\right)}\right)}{_{a_i}^{\mathrm{}}dz\mathrm{Ai}^2(z)}}\mathrm{Ai}^{}(a_i)\mathrm{e}^{\lambda _it}.`$ (89)
In Fig. 7, the exact time evolution is compared with numerical results, obtained with $`L=3\times 10^6`$. The lower dotted line is obtained by plotting only the first term of the sum in (89), using the values from Table 1. Explicitly, the first eigenvalue that determines the long-time approach to the steady state is

$$\lambda _1=4.676\mathrm{}(DQ^2)^{1/3}=7.4227\mathrm{}(D\mathrm{\Gamma }^2)^{1/3}.$$
(90)
#### 3.3.2 Random initial density

An interesting case is provided by starting the system with the exact steady-state density $`\rho (0)=\rho _0^\mathrm{u}`$, but with a random initial distribution of particles. There is an initial period of rapid annihilation that reduces the density, followed by a slower relaxation back to the steady-state value.

For a random initial distribution of particles with density $`\rho `$, the number of particles in $`(0,x)`$ is a Poisson random variable with mean $`\rho x`$. The function $`r^\mathrm{u}(x,0)`$ can be calculated as follows:

$`r^\mathrm{u}(x,0)`$ $`=`$ $`𝒫[0\text{ particles between }0\text{ and }x]`$ (91)
$`+𝒫[2\text{ particles between }0\text{ and }x]+\mathrm{}`$
$`=`$ $`\mathrm{e}^{\rho x}+\mathrm{e}^{\rho x}{\displaystyle \frac{(\rho x)^2}{2!}}+\mathrm{}`$
$`=`$ $`{\displaystyle \frac{1}{2}}\left(1+\mathrm{e}^{2\rho x}\right).`$
Figure 8 shows data from a numerical simulation, performed with $`L=2\times 10^5`$, along with the results of the calculation of the coefficients in (86), using (91).

### 3.4 Time-dependent statistics: Paired nucleation

Let us introduce

$$h^\mathrm{p}(x,t)=r^\mathrm{p}(x,t)r_0^\mathrm{p}(x).$$
(92)
Then $`h^\mathrm{p}(x,t)`$ satisfies

$$\frac{}{t}h^\mathrm{p}(x,t)=h^\mathrm{p}(x,t),$$
(93)
where the operator $``$ is defined by

$`f(x)=\{\begin{array}{ccc}2D{\displaystyle \frac{\mathrm{d}^2}{\mathrm{d}x^2}}f(x)4x\mathrm{\Gamma }f(x)\hfill & xb;\hfill & \\ 2D{\displaystyle \frac{\mathrm{d}^2}{\mathrm{d}x^2}}f(x)4b\mathrm{\Gamma }f(x)\hfill & x>b.\hfill & \end{array}`$ (96)
The boundary conditions on $`h^\mathrm{p}(x,t)`$ are

$$h^\mathrm{p}(0,t)=0\text{and}h^\mathrm{p}(x,t)\text{bounded as}x\mathrm{},$$
(97)
for $`t>0`$.

Let us introduce

$$\alpha _i=\frac{\lambda _i}{4b\mathrm{\Gamma }}.$$
(98)
For $`\alpha _i>1`$, the eigenvalue equation

$$h_i^\mathrm{p}(x)=\lambda _ih_i^\mathrm{p}(x),$$
(99)
has a continuous spectrum of solutions

$$h_i^\mathrm{p}(x)=\{\begin{array}{ccc}c_1(\alpha _i)\left[\sqrt{3\mathrm{Ai}}\left(\epsilon (\frac{x}{b}\alpha _i)\right)\mathrm{Bi}\left(\epsilon (\frac{x}{b}\alpha _i)\right)\right]\hfill & xb;\hfill & \\ c_2(\alpha _i)\mathrm{sin}\left[(\alpha _i1)^{\frac{1}{2}}(\frac{2\mathrm{\Gamma }b}{D})^{\frac{1}{2}}(xb)\right]\hfill & & \\ +c_3(\alpha _i)\mathrm{cos}\left[(\alpha _i1)^{\frac{1}{2}}(\frac{2\mathrm{\Gamma }b}{D})^{\frac{1}{2}}(xb)\right]\hfill & x>b.\hfill & \end{array}$$
(100)
When $`\epsilon `$ is sufficiently large, there are also discrete eigenvalues at values of $`\alpha _i<1`$ satisfying

$`\mathrm{Ai}(\epsilon \alpha _i)\mathrm{Bi}^{}(\epsilon (1\alpha _i))\mathrm{Bi}(\epsilon \alpha _i)\mathrm{Ai}^{}(\epsilon (1\alpha _i))`$
$`+\left(\epsilon (1\alpha _i)\right)^{\frac{1}{2}}\left(\mathrm{Ai}(\epsilon \alpha _i)\mathrm{Bi}(\epsilon (1\alpha _i))\mathrm{Bi}(\epsilon \alpha _i)\mathrm{Ai}(\epsilon (1\alpha _i))\right)`$
$`=0.`$ (101)
The eigenfunctions in this case are

$`h_i^\mathrm{p}(x)=\{\begin{array}{ccc}c_1(\alpha _i)\left[\sqrt{3\mathrm{Ai}}\left(\epsilon (\frac{x}{b}\alpha _i)\right)\mathrm{Bi}\left(\epsilon (\frac{x}{b}\alpha _i)\right)\right]\hfill & xb;\hfill & \\ c_4(\alpha _i)\mathrm{exp}[(\frac{2\mathrm{\Gamma }b}{D})^{\frac{1}{2}}(1\alpha _i)^{\frac{1}{2}}x]\hfill & x>b.\hfill & \end{array}`$ (104)
Thus, for all finite $`\epsilon `$, there is a continuous spectrum of eigenvalues with $`\lambda _i4b\mathrm{\Gamma }`$. For $`\epsilon `$ smaller than $`\epsilon _c`$ these are the only eigenvalues. The critical value $`\epsilon _c`$ satisfies

$$\mathrm{Ai}(\epsilon _c)\mathrm{Bi}^{}(0)\mathrm{Bi}(\epsilon _c)\mathrm{Ai}^{}(0)=0,$$
(105)
so

$$\epsilon _c=1.98\mathrm{}.$$
(106)
Discrete eigenvalues appear for larger values of $`\epsilon `$ (Fig. 9). The unpaired limit is regained as $`\epsilon \mathrm{}`$, when $`\epsilon \alpha _ia_i`$, so that $`\lambda _i2.338(32D\mathrm{\Gamma }^2)^{\frac{1}{3}}`$.

In Figs. 10 and 11 we compare the exact expressions for the steady-state density and exponents characterizing relaxation toward the steady-state density with numerical results. In each case the curved dotted line is obtained using the lowest exponent available, but the coefficient is obtained from a best fit. In the case depicted in Fig. 10, there is a discrete eigenvalue $`\lambda _i<4b\mathrm{\Gamma }`$; in the case depicted in Fig. 11, there is only the continuum of eigenvalues $`\lambda _i4b\mathrm{\Gamma }`$, and the relaxation process may not be purely exponential.

### 3.5 Time-dependent statistics: No nucleation

In the absence of nucleation, $`r^\mathrm{n}(x,t)`$ satisfies the heat equation

$`{\displaystyle \frac{}{t}}r^\mathrm{n}(x,t)`$ $`=`$ $`2D{\displaystyle \frac{^2}{x^2}}r^\mathrm{n}(x,t),`$ (107)
with the boundary conditions

$$r^\mathrm{n}(0,t)=1\text{and}\underset{x\mathrm{}}{lim}r^\mathrm{n}(x,t)=\frac{1}{2},$$
(108)
for all $`t>0`$. The solution of (107) is given by

$`r^\mathrm{n}(x,t)`$ $`=`$ $`1+\left(8\pi Dt\right)^{\frac{1}{2}}({\displaystyle _0^{\mathrm{}}}\mathrm{d}y(1r^\mathrm{n}(y,0))\mathrm{e}^{(xy)^2/8Dt}`$ (109)
$`+{\displaystyle _{\mathrm{}}^0}\mathrm{d}y(r^\mathrm{n}(y,0)1)\mathrm{e}^{(xy)^2/8Dt}).`$
If the initial distribution of particles is random with density $`\rho (0)`$, then

$$r^\mathrm{n}(x,0)=\frac{1}{2}\left(1\mathrm{e}^{2\rho (0)x}\right).$$
(110)
Now, using (58), we derive the density of particles as a function of time for random initial conditions

$`\rho (t)`$ $`=`$ $`\left(8\pi Dt\right)^{\frac{1}{2}}({\displaystyle _0^{\mathrm{}}}\mathrm{d}y\mathrm{e}^{2\rho (0)t}{\displaystyle \frac{y}{8Dt}}\mathrm{e}^{y^2/8Dt}`$ (111)
$`+{\displaystyle _{\mathrm{}}^0}\mathrm{d}y\mathrm{e}^{2\rho (0)t}{\displaystyle \frac{y}{8Dt}}\mathrm{e}^{y^2/8Dt})`$
$`=`$ $`2\rho (0)\left(8\pi Dt\right)^{\frac{1}{2}}{\displaystyle _0^{\mathrm{}}}dy\mathrm{e}^{2\rho (0)y}\mathrm{e}^{y^2/8Dt}`$
$`=`$ $`\rho (0)\mathrm{exp}\left(8Dt\rho ^2(0)\right)\mathrm{erfc}\left(\rho (0)(8Dt)^{\frac{1}{2}}\right).`$
We thus reproduce the result (4) of Torney and McConnell .

## 4 Discussion

In the case of unpaired nucleation, there is only one length scale, proportional to $`(D/Q)^{\frac{1}{3}}`$, and only one time scale, proportional to $`(DQ^2)^{\frac{1}{3}}`$. The relaxation time to equilibrium and the mean lifetime of a particle are proportional to one another. This is made clear in the mesoscopic approach , although in the approaches detailed here we do not make this explicit distinction. It is noteworthy that the reaction-diffusion approach yields a steady state that is within 9% of the correct one and a relaxation rate that differs from the exact result by only 2.3%.

Although the case of paired nucleation is more complicated in terms of the time scales associated with its dynamics, its steady-state distribution of particles is closer to a classical equilibrium random distribution than the corresponding distribution produced by unpaired nucleation. Note that the truncated hierarchy approach in this case leads to the exact steady state density. The underlying reason is that the dynamics produced by paired nucleation is close to time-reversal invariant. For comparison, an ensemble of noninteracting diffusing particles has a two-point function $`g(y)`$ identically equal to 1 . We can imagine producing a spacetime diagram such as shown in Fig. 1 from a diagram associated with noninteracting particles in two steps. First, when two particles collide, move them to a different, randomly chosen part of the system. Second, separate them by a distance $`b`$. The first step does not affect the correlation function or time-reversal invariance. The second step directly changes the correlation function for separations smaller than $`b`$. Thus, for diffusion-limited annihilation with paired nucleation, as the parameter for $`\epsilon `$ that measures the distance between newly nucleated pairs tends to $`0`$, the two-point function in the steady state is appreciably different from $`1`$ only in a region whose width is proportional to $`b`$.

Results obtained for the steady-state density and relaxation rate are summarized in Table 2.

For the evaluation of time scales, the case of paired nucleation is more complicated. The different approaches indicate the occurrence of multiple length and time scales. The mesoscopic approach  leads to a characteristic time for approach to equilibrium proportional to $`(b\mathrm{\Gamma })^1`$, and a distinct mean lifetime of a particle proportional to $`(b/D\mathrm{\Gamma })^{\frac{1}{2}}`$. These two time scales were identified as corresponding, respectively, to recombination (two particles created a distance $`b`$ apart collide and annihilate) and to nonrecombinant annihilation (collision between two particles nucleated at different times). Unpaired annihilation is less frequent than paired annihilation, but both time scales are important in the dynamics of the system . The exact approach leads to the former time scale as an upper bound of a continuum of scales. The truncated reaction-diffusion hierarchy leads to the latter time scale under the assumption of exponential decay, which may not be valid. Understanding the time scales in the case of paired nucleation requires further research.

The theoretical approach based on truncation of a hierarchy of distribution functions thus permits the calculation of steady-state densities and correlation functions that are in fair agreement with simulations. It gives the exact result for the steady-state density in the limit $`\epsilon 0`$, where the statistical distribution of particles is close to random. As pointed out by van Kampen , an approach using a truncation can be made systematic if it is based on an expansion in a small parameter. However, his suggestion that the small parameter be the density of particles is not applicable to the case of unpaired nucleation of point particles because there is no other quantity with the dimension of length, i.e., nothing for the density to be small compared to. The exact approach based on the function $`r(x,t)`$ sidesteps these difficulties by providing a direct method to calculate the density of particles. For any value of $`\epsilon `$, the density in the steady state and its time-dependent statistics can be exactly calculated. However, the method has not yet been extended to exact calculation of the full distribution of interparticle distances.

## 5 Acknowledgments

The authors acknowledge support from IGPP under Project Los Alamos/DOE 822AR. K. L. gratefully acknowledges support from the Engineering Research Program of the Office of Basic Energy Sciences at the U. S. Department of Energy under Grant No. DE-FG03-86ER13606. C.M.-P. wishes to acknowledge support from Centro de Astrobiología. S. H., G. L., and C. M.-P. wish to thank K.L. for her hospitality while visiting USCD, where part of this work was carried out. The authors are grateful to Ubbo Felderhof for his insightful comments and corrections, and for the exact inversion of Eq. (B.17).

APPENDIX A: Steady state solution of truncated hierarchy

In this Appendix we find the steady state solutions of (LABEL:gdott) and (LABEL:pairnoaunc).

##### Unpaired nucleation

We solve (LABEL:pairnoaunc) with the left-hand side set to zero. Fourier transforming the quasilinear equation according to

$$\widehat{g}_n=\frac{2}{L^{1/2}}_0^Ldyg_0(y)\mathrm{cos}\frac{2\pi ny}{L},$$
(A.1)
leads for $`n0`$ to

$$\widehat{g}_n=\frac{\frac{2\mathrm{\Gamma }}{L^{1/2}\rho _0^2}}{\frac{8D\pi ^2n^2}{L^2}+\frac{4\mathrm{\Gamma }}{\rho _0}},$$
(A.2)
where Eq. (30) has been used. Normalization sets $`\widehat{g}_0=L^{1/2}`$. Fourier inversion according to

$$g(y)=\frac{1}{L^{1/2}}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}\widehat{g}_n\mathrm{cos}\frac{2\pi ny}{L},$$
(A.3)
can be done by separating out the $`n=0`$ contribution explicitly and changing the resulting sum to an integral (valid as $`L\mathrm{}`$)

$`g_0(y)`$ $`=`$ $`1{\displaystyle \frac{\mathrm{\Gamma }}{\pi \rho _0^2}}{\displaystyle _0^{\mathrm{}}}dq{\displaystyle \frac{\mathrm{cos}qy}{q^2+{\displaystyle \frac{2\mathrm{\Gamma }}{D\rho _0}}}}`$ (A.4)
$`=`$ $`1{\displaystyle \frac{1}{4}}\left({\displaystyle \frac{2\mathrm{\Gamma }}{D\rho _0^3}}\right)^{1/2}\mathrm{e}^{y(2\mathrm{\Gamma }/D\rho _0)^{1/2}}.`$
The requirement (28) must also hold in the steady state. Together with (A.4) this leads to Eqs. (44) and (45).

##### Paired nucleation

We now solve Eq. (46), where the steady-state density $`\rho _0`$ is related to the derivative of $`g_0(y)`$ as $`y0`$

$$\mathrm{\Gamma }=2Dg_0^{}(0^+)\rho _0^2.$$
(A.5)
The $`\delta `$-function contribution in the last term of (46) leads to a discontinuity in the derivative $`g_0^{}(y)`$ at $`y=b`$. This leads to the search for a solution of the form

$$g_0(y)=\{\begin{array}{ccc}1+(S1)e^{(2\mathrm{\Gamma }/D\rho _0)^{1/2}y}\hfill & 0y<b;\hfill & \\ 1+Pe^{(2\mathrm{\Gamma }/D\rho _0)^{1/2}y}+Se^{(2\mathrm{\Gamma }/D\rho _0)^{1/2}y}\hfill & yb.\hfill & \end{array}$$
(A.6)
The constant $`S`$ is determined from the condition (A.5)

$$12S=\left(\frac{\mathrm{\Gamma }}{8D\rho _0^3}\right)^{\frac{1}{2}}.$$
(A.7)
Now, the constant $`P`$ is determined by enforcing continuity of the solution $`g_0(y)`$ at $`y=b`$

$$P=S(1e^{2(2\mathrm{\Gamma }/D\rho _0)^{1/2}b})1.$$
(A.8)
The discontinuity in the derivative $`g_0^{}(y)`$ at $`y=b`$ is

$$g_0^{}(b^+)g_0^{}(b^{})=\frac{\mathrm{\Gamma }}{2D\rho _0^2}.$$
(A.9)
Using (A.6) to evaluate the left-hand side of (A.9) and rearranging gives an implicit expression for the steady-state density

$$\left(\frac{\mathrm{\Gamma }}{8D}\right)^{1/2}=\left(\left(\frac{\mathrm{\Gamma }}{8D}\right)^{1/2}\rho _0^{3/2}\right)e^{(2\mathrm{\Gamma }/D\rho _0)^{1/2}b},$$
(A.10)
or, rearranging again

$$(b\rho _0)^{3/2}=\frac{\sigma ^2}{8}\left(1\mathrm{e}^{\sigma ^2/2(b\rho _0)^{1/2}}\right),$$
(A.11)
where $`\sigma `$ is defined by

$$\sigma =\left(\frac{8b^3\mathrm{\Gamma }}{D}\right)^{\frac{1}{4}}\sqrt{2}\epsilon ^{3/4}.$$
(A.12)
While we cannot invert (A.11) explicitly for the steady-state density, we can examine the limits in the dimensionless parameter $`\sigma `$. In the limit $`\sigma 1`$, corresponding to large initial separation, we find

$$(b\rho _0)^{3/2}\frac{\sigma ^2}{8},$$
(A.13)
so $`\rho _0^3\mathrm{\Gamma }/8D`$, as obtained in (44) for the unpaired nucleation case. Morevover, $`S0`$ and $`g(y)1\mathrm{e}^{4\rho _0y}`$ as in (45).

The limit $`\sigma 0`$ corresponds to small initial separation. Then (A.11) reduces to

$$b\rho _0=\frac{\sigma ^2}{4}\left(1\frac{1}{4}\sigma +𝒪(\sigma ^2)\right).$$
(A.14)
Note that $`\sigma 0`$ corresponds to $`b\rho _00`$. Expanding in powers of $`\sigma `$, we find

$$\rho _0=\left(\frac{b\mathrm{\Gamma }}{2D}\right)^{\frac{1}{2}}(1\frac{1}{4}\sigma +𝒪(\sigma ^2)),S=\frac{1}{2}\sigma ^1+\frac{5}{16}+𝒪(\sigma )\text{and}P=\frac{1}{2}\sigma +𝒪(\sigma ^2).$$
(A.15)
The correlation function $`g_0`$ can be expanded as

$$g_0(y)=\{\begin{array}{ccc}\frac{y}{b}(1+\frac{1}{2}\sigma )+𝒪(\sigma ^2)\hfill & 0y<b;\hfill & \\ 1+\frac{1}{2}\sigma e^{\sigma y/b}+𝒪(\sigma ^2)\hfill & yb.\hfill & \end{array}$$
(A.16)
These results as $`\sigma 0`$ are reported in Eqs. (47) and (48).

APPENDIX B: Relaxation to the Steady State

Here we detail the calculation of approach to the steady state in the truncated hierarchy approach.

##### Unpaired nucleation

It is convenient to introduce the symbols

$$\chi =(D\mathrm{\Gamma }^2)^{1/3},\gamma =(D^2\mathrm{\Gamma })^{1/3}.$$
(B.1)
The linearized perturbation equations (51) and (LABEL:lgho1) in the unpaired case are

$`{\displaystyle \frac{}{t}}\delta \rho (t)`$ $`=`$ $`\chi \left[\delta g^{}(0^+,t)+8\delta \rho (t)\right],`$ (B.2)
$`{\displaystyle \frac{}{t}}\delta g(y,t)`$ $`=`$ $`2D{\displaystyle \frac{^2}{y^2}}\delta g(y,t)2K(y)\delta g(y,t)8\chi \delta g(y,t)`$ (B.3)
$`16\gamma \left[1g_0^\mathrm{u}(y)\right]\delta \rho (t),`$
where the specific form (44) has been implemented. As initial conditions we choose

$$\delta g(y,t=0)=0,$$
(B.4)
and an arbitrary $`\delta \rho (0)`$. Note that having implemented the condition (21) on the steady-state solution implies that $`\delta \widehat{g}_{n=0}(t)=0`$.

The solution of Eq. (B.2) is

$$\delta \rho (t)=e^{8\chi t}\delta \rho (0)\chi _0^td\tau e^{8\chi (t\tau )}\delta g^{}(0^+,\tau ).$$
(B.5)
Transforming (B.3) according to (A.1) gives, for $`n0`$

$$\frac{}{t}\delta \widehat{g}_n(t)=\frac{8D\pi ^2n^2}{L^2}\delta \widehat{g}_n(t)8\chi \delta \widehat{g}_n(t)\frac{8D}{L^{1/2}}\delta g^{}(0^+,t)+16\gamma \widehat{g}_n\delta \rho (t).$$
(B.6)
We can formally solve this equation as well, to obtain

$$\delta \widehat{g}_n(t)=_0^td\tau e^{8\left(\frac{D\pi ^2n^2}{L^2}+\chi \right)(t\tau )}\left[\frac{8D}{L^{1/2}}\delta g^{}(0^+,\tau )+16\gamma \widehat{g}_n\delta \rho (\tau )\right],$$
(B.7)
and hence its inverse Fourier transform

$$\delta g(y,t)=_0^td\tau 𝒦(y,t\tau )\delta g^{}(0^+,\tau )+8_0^td\tau 𝒢(y,t\tau )\delta \rho (\tau ).$$
(B.8)
We have introduced the following functions (and taken the limit $`L\mathrm{}`$):

$`𝒦(y,t)`$ $`=`$ $`{\displaystyle \frac{8D}{L}}e^{8\chi t}{\displaystyle \stackrel{}{}_{n=\mathrm{}}^+\mathrm{}}\mathrm{cos}{\displaystyle \frac{2\pi ny}{L}}e^{\frac{8D\pi ^2n^2}{L^2}t}`$ (B.9)
$`=`$ $`{\displaystyle \frac{8D}{(2\pi Dt)^{1/2}}}e^{8\chi t}e^{\frac{y^2}{2Dt}},`$
$`𝒢(y,t)`$ $`=`$ $`{\displaystyle \frac{2\gamma }{L^{1/2}}}e^{8\chi t}{\displaystyle \stackrel{}{}_{n=\mathrm{}}^+\mathrm{}}\mathrm{cos}{\displaystyle \frac{2\pi ny}{L}}e^{\frac{8D\pi ^2n^2}{L^2}t}\widehat{g}_n`$ (B.10)
$`=`$ $`{\displaystyle \frac{8\gamma ^2}{\pi }}e^{8\chi t}{\displaystyle _{\mathrm{}}^{\mathrm{}}}dqe^{2Dq^2t}{\displaystyle \frac{1}{8\chi +2Dq^2}}\mathrm{cos}qy,`$
where the prime on the sums indicates omission of the $`n=0`$ term and we have used Eq. (A.2).

Since the unknowns appear on both sides of (B.5) and (B.8), these are only formal solutions. To proceed, we Laplace transform them (indicated by a tilde, $`f(t)=\stackrel{~}{f}(s)`$, and solve the resulting set self-consistently. The limit (32) must be handled carefully and not implemented prematurely. We find

$`\delta \stackrel{~}{\rho }(s)`$ $`=`$ $`{\displaystyle \frac{1}{8\chi (u+1)}}\left(\delta \rho (0)\chi \delta \stackrel{~}{g}^{}(0^+,s)\right),`$ (B.11)
$`\delta \stackrel{~}{g}(y,s)`$ $`=`$ $`\stackrel{~}{𝒦}(y,s)\delta \stackrel{~}{g}^{}(0^+,s)+8\stackrel{~}{𝒢}(y,s)\delta \stackrel{~}{\rho }(s).`$ (B.12)
From these two equations we can obtain an expression for $`\delta \stackrel{~}{g}^{}(0^+,s)`$ as follows. First, set $`y=0`$ in (B.12). Since $`g(0,t)=0`$ and since $`g_0(0)=0`$, it follows from Eq. (50) that $`\delta g(0,t)=0`$ for all $`t`$ and therefore $`\delta \stackrel{~}{g}(0,s)=0`$ for all $`s`$. Thus, we find

$$0=\stackrel{~}{𝒦}(0,s)\delta \stackrel{~}{g}^{}(0^+,s)+8\stackrel{~}{𝒢}(0,s)\delta \stackrel{~}{\rho }(s).$$
(B.13)
Subsitution of (B.11) into this result immediately leads to

$$\chi \delta \stackrel{~}{g}^{}(0^+,s)=\frac{\delta \rho (0)\stackrel{~}{𝒢}(0,s)}{(u+1)\stackrel{~}{𝒦}(0,s)+\stackrel{~}{𝒢}(0,s)},$$
(B.14)
where we have set $`us/8\chi `$. This, together with Eqs. (B.11) and (B.12), constitutes a complete solution of the Laplace transform of the problem. The indicated transforms can be calculated explicitly

$`\stackrel{~}{𝒦}(0,s)`$ $`=`$ $`{\displaystyle \frac{2(D/\mathrm{\Gamma })^{1/3}}{\sqrt{u+1}}},`$ (B.15)
$`\stackrel{~}{𝒢}(0,s)`$ $`=`$ $`{\displaystyle \frac{(D/\mathrm{\Gamma })^{1/3}}{4u}}\left(1{\displaystyle \frac{1}{\sqrt{u+1}}}\right),`$ (B.16)
and one readily obtains

$$\delta \stackrel{~}{\rho }(s)=\frac{\delta \rho (0)}{\chi }\left(\frac{u}{1\sqrt{1+u}+8u(1+u)}\right).$$
(B.17)
It is possible to Laplace invert these expressions exactly . In particular, Eq. (B.17) can be rewritten as

$$\delta \stackrel{~}{\rho }(s)=\frac{\delta \rho (0)}{\chi }\left(\frac{A_1}{\sqrt{u+1}y_1}+\frac{A_2}{\sqrt{u+1}y_2}+\frac{A_3}{\sqrt{u+1}y_3}\right).$$
(B.18)
where

$$A_1=\frac{5+\sqrt{5}}{40}A_2=\frac{5\sqrt{5}}{40}A_3=\frac{1}{4}$$
(B.19)
and

$$y_1=\frac{1}{4}\left(\sqrt{5}1\right)y_2=\frac{1}{4}\left(\sqrt{5}+1\right)y_3=\frac{1}{2}.$$
(B.20)
The inversion

$$^1\left(\frac{1}{\sqrt{s+A}B}\right)=e^{At}\left(\frac{1}{\sqrt{\pi t}}+Be^{B^2t}\mathrm{erfc}(B\sqrt{t})\right)$$
(B.21)
can then be applied to obtain

$$\delta \rho (t)=8\delta \rho (0)e^{8\chi t}\underset{i}{}A_iy_ie^{8\chi y_i^2t}\mathrm{erfc}(y_i\sqrt{8\chi t}).$$
(B.22)
It is noteworthy that this solution is exact within the truncation approximations of the model, that is, it represents the full time-dependent solution for the model.

Asymptotic analysis of the exact result yields pure exponential decay as indicated in Eqs. (53) and (54), with

$$\alpha =(5+\sqrt{5})\chi =7.236\mathrm{}\chi .$$
(B.23)
The proportionality of $`\delta \stackrel{~}{g}`$ and $`\delta \stackrel{~}{\rho }`$ clearly leads to the same decay rate for $`\delta g(y,t)`$ as for $`\delta \rho (t)`$.

##### Paired nucleation

Here it is convenient to introduce the symbols

$$\mathrm{\Omega }=\frac{4\mathrm{\Gamma }}{\rho _0},\theta =\mathrm{\Gamma }b.$$
(B.24)
The linearized perturbation equations (51) and (LABEL:lgho1) in the paired case are

$`{\displaystyle \frac{}{t}}\delta \rho (t)`$ $`=`$ $`4D\rho _0^2\delta g^{}(0^+,t)\mathrm{\Omega }\delta \rho (t),`$ (B.25)
$`{\displaystyle \frac{}{t}}\delta g(y,t)`$ $`=`$ $`2D{\displaystyle \frac{^2}{y^2}}\delta g(y,t)2K(y)\delta g(y,t)\mathrm{\Omega }\delta g(y,t)`$ (B.26)
$`{\displaystyle \frac{4\mathrm{\Gamma }}{\rho _0}}\left[1g_0^\mathrm{p}(y)\right]\delta \rho (t){\displaystyle \frac{2\mathrm{\Gamma }}{\rho _0^3}}\delta (yb)\delta \rho (t).`$
With the same initial condition as in the unpaired case, the solution of Eq. (B.25) is formally given by

$$\delta \rho (t)=e^{\mathrm{\Omega }t}\delta \rho (0)4D\rho _0^2_0^td\tau e^{\mathrm{\Omega }(t\tau )}\delta g^{}(0^+,\tau ).$$
(B.27)
Transforming (B.26) according to (A.1) gives for $`n0`$

$`{\displaystyle \frac{}{t}}\delta \widehat{g}_n(t)`$ $`=`$ $`{\displaystyle \frac{8D\pi ^2n^2}{L^2}}\delta \widehat{g}_n(t)\mathrm{\Omega }\delta \widehat{g}_n(t){\displaystyle \frac{8D}{L^{1/2}}}\delta g^{}(0^+,t)`$ (B.28)
$`+{\displaystyle \frac{\mathrm{\Omega }}{\rho _0}}\widehat{g}_n\delta \rho (t){\displaystyle \frac{4\mathrm{\Gamma }}{\rho _0^3L^{1/2}}}\mathrm{cos}\left({\displaystyle \frac{2\pi nb}{L}}\right)\delta \rho (t).`$
We can formally solve this equation as well, to obtain

$`\delta \widehat{g}_n(t)`$ $`=`$ $`{\displaystyle _0^t}d\tau e^{\left(\frac{8D\pi ^2n^2}{L^2}+\mathrm{\Omega }\right)(t\tau )}`$
$`\times \left\{{\displaystyle \frac{8D}{L^{1/2}}}\delta g^{}(0^+,\tau )+\left[{\displaystyle \frac{\mathrm{\Omega }}{\rho _0}}\widehat{g}_n{\displaystyle \frac{4\mathrm{\Gamma }}{\rho _0^3L^{1/2}}}\mathrm{cos}\left({\displaystyle \frac{2\pi nb}{L}}\right)\right]\delta \rho (\tau )\right\},`$
and hence its inverse Fourier transform (notice that $`\delta g_{n=0}(t)=0`$ for all $`t`$)

$`\delta g(y,t)`$ $`=`$ $`{\displaystyle _0^t}d\tau 𝒦(y,t\tau )\delta g^{}(0^+,\tau )+{\displaystyle _0^t}d\tau 𝒢(y,t\tau )\delta \rho (\tau )`$ (B.29)
$`{\displaystyle \frac{4\mathrm{\Gamma }}{\rho _0^3L^{1/2}}}{\displaystyle _0^t}d\tau (y,t\tau )\delta \rho (\tau ).`$
We have introduced the following functions:

$`𝒦(y,t)`$ $`=`$ $`{\displaystyle \frac{8D}{L}}e^{\mathrm{\Omega }t}{\displaystyle \stackrel{}{}_{n=\mathrm{}}^+\mathrm{}}\mathrm{cos}{\displaystyle \frac{2\pi ny}{L}}e^{\frac{8D\pi ^2n^2}{L^2}t}`$ (B.30)
$`=`$ $`{\displaystyle \frac{8D}{(2\pi Dt)^{1/2}}}e^{\mathrm{\Omega }t}e^{\frac{y^2}{2Dt}},`$
$`𝒢(y,t)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }}{\rho _0L^{1/2}}}e^{\mathrm{\Omega }t}{\displaystyle \stackrel{}{}_{n=\mathrm{}}^+\mathrm{}}\mathrm{cos}{\displaystyle \frac{2\pi ny}{L}}e^{\frac{8D\pi ^2n^2}{L^2}t}\widehat{g}_n,`$ (B.31)
$`(y,t)`$ $`=`$ $`{\displaystyle \frac{1}{L^{1/2}}}e^{\mathrm{\Omega }t}{\displaystyle \stackrel{}{}_{n=\mathrm{}}^+\mathrm{}}\mathrm{cos}{\displaystyle \frac{2\pi ny}{L}}\mathrm{cos}{\displaystyle \frac{2\pi nb}{L}}e^{\frac{8D\pi ^2n^2}{L^2}t}.`$ (B.32)
We Laplace transform (B.27) and (B.29) and solve the resulting set self-consistently, to obtain

$`\delta \stackrel{~}{\rho }(s)`$ $`=`$ $`{\displaystyle \frac{1}{s+\mathrm{\Omega }}}\delta \rho (0)4D\rho _0^2{\displaystyle \frac{\delta \stackrel{~}{g}^{}(0^+,s)}{s+\mathrm{\Omega }}},`$
$`\delta \stackrel{~}{g}(y,s)`$ $`=`$ $`\stackrel{~}{𝒦}(y,s)\delta \stackrel{~}{g}^{}(0^+,s)+\stackrel{~}{𝒢}(y,s)\delta \stackrel{~}{\rho }(s){\displaystyle \frac{4\mathrm{\Gamma }}{\rho _0^3L^{1/2}}}\stackrel{~}{}(y,s)\delta \stackrel{~}{\rho }(s),`$
from which in turn we find

$`\delta \stackrel{~}{g}(y,s)`$ $`=`$ $`\stackrel{~}{𝒦}(y,s)\delta \stackrel{~}{g}^{}(0^+,s)+\left[\stackrel{~}{𝒢}(y,s){\displaystyle \frac{4\mathrm{\Gamma }}{\rho _0^3L^{1/2}}}\stackrel{~}{}(y,s)\right]`$ (B.33)
$`\times \left[{\displaystyle \frac{1}{s+\mathrm{\Omega }}}\delta \rho (0)4D\rho _0^2{\displaystyle \frac{\delta \stackrel{~}{g}^{}(0^+,s)}{s+\mathrm{\Omega }}}\right].`$
We are interested in obtaining $`\delta \stackrel{~}{g}^{}(0^+,s)`$. This can readily be done by evaluating the previous equation at $`y=0`$, to obtain (notice that we have chosen $`\delta g(y=0,t)=0`$ for all $`t`$)

$`0`$ $`=`$ $`\stackrel{~}{𝒦}(0,s)\delta \stackrel{~}{g}^{}(0^+,s)+\left[\stackrel{~}{𝒢}(0,s){\displaystyle \frac{4\mathrm{\Gamma }}{\rho _0^3L^{1/2}}}\stackrel{~}{}(0,s)\right]`$ (B.34)
$`\times \left[{\displaystyle \frac{1}{s+\mathrm{\Omega }}}\delta \rho (0)4D\rho _0^2{\displaystyle \frac{\delta \stackrel{~}{g}^{}(0^+,s)}{s+\mathrm{\Omega }}}\right].`$
The previously introduced functions evaluated at the origin become

$`𝒦(0,t)`$ $`=`$ $`{\displaystyle \frac{8D}{(2\pi Dt)^{1/2}}}e^{\mathrm{\Omega }t},`$ (B.35)
$`𝒢(0,t)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }}{\rho _0L^{1/2}}}e^{\mathrm{\Omega }t}{\displaystyle \stackrel{}{}_{n=\mathrm{}}^+\mathrm{}}e^{\frac{8D\pi ^2n^2}{L^2}t}\widehat{g}_n,`$ (B.36)
$`(0,t)`$ $`=`$ $`{\displaystyle \frac{1}{L^{1/2}}}e^{\mathrm{\Omega }t}{\displaystyle \stackrel{}{}_{n=\mathrm{}}^+\mathrm{}}\mathrm{cos}{\displaystyle \frac{2\pi nb}{L}}e^{\frac{8D\pi ^2n^2}{L^2}t}.`$ (B.37)
We are interested in obtaining the Laplace transform of these functions. It is easy to obtain

$`\stackrel{~}{𝒦}(0,s)`$ $`=`$ $`{\displaystyle \frac{2(D/\mathrm{\Gamma })^{1/3}}{\sqrt{s/\mathrm{\Omega }+1}}},`$ (B.38)
$`(0,t)`$ $`=`$ $`{\displaystyle \frac{L^{1/2}}{\pi }}e^{\mathrm{\Omega }t}{\displaystyle _0^+\mathrm{}}dqe^{2Dq^2t}\mathrm{cos}qb,`$ (B.39)
so that

$`{\displaystyle \frac{4\mathrm{\Gamma }}{\rho _0^3L^{1/2}}}(0,t)`$ $`=`$ $`{\displaystyle \frac{2\mathrm{\Gamma }}{\rho _0^3\sqrt{2\pi Dt}}}e^{\mathrm{\Omega }t}e^{b^2/(8Dt)},`$ (B.40)
$`{\displaystyle \frac{4\mathrm{\Gamma }}{\rho _0^3L^{1/2}}}\stackrel{~}{}(0,s)`$ $`=`$ $`{\displaystyle \frac{2\mathrm{\Gamma }}{\rho _0^3\sqrt{2D}}}{\displaystyle \frac{1}{\sqrt{s+\mathrm{\Omega }}}}e^{b\sqrt{\frac{s+\mathrm{\Omega }}{2D}}}.`$ (B.41)
The Fourier components for the two-point correlation function associated with Eq. (A.6) are

$`\widehat{g}_n`$ $`=`$ $`{\displaystyle \frac{2}{L^{1/2}}}{\displaystyle \frac{\sqrt{2\mathrm{\Gamma }/(D\rho _0)}}{\frac{2\mathrm{\Gamma }}{D\rho _0}+\frac{4\pi ^2n^2}{L^2}}}\left[(2S1)2Se^{b\sqrt{2\mathrm{\Gamma }/(D\rho _0)}}\mathrm{cos}{\displaystyle \frac{2\pi nb}{L}}\right].`$ (B.42)
This form in turn leads to the Laplace transform of $`𝒢(0,t)`$:

$`\stackrel{~}{𝒢}(0,s)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }}{\pi }}\sqrt{{\displaystyle \frac{2\mathrm{\Gamma }}{D\rho _0^3}}}2(2S1){\displaystyle _0^+\mathrm{}}dq{\displaystyle \frac{1}{(q^2+2\mathrm{\Gamma }/(D\rho _0))(s+\mathrm{\Omega }+2Dq^2)}}`$ (B.43)
$`{\displaystyle \frac{\mathrm{\Omega }}{\pi }}\sqrt{{\displaystyle \frac{2\mathrm{\Gamma }}{D\rho _0^3}}}2Se^{b\sqrt{2\mathrm{\Gamma }/(D\rho _0)}}`$
$`\times 2{\displaystyle _0^+\mathrm{}}dq{\displaystyle \frac{\mathrm{cos}qb}{(q^2+2\mathrm{\Gamma }/(D\rho _0))(s+\mathrm{\Omega }+2Dq^2)}}.`$
It is now a straightforward matter to 1) collect the various Laplace transform expressions to solve for $`\delta \stackrel{~}{g}^{}(0^+,s)`$ using Eq. (B.34), 2) substitute this result into Eq. (B.11), and explore the poles of the denominator of the resulting $`\delta \stackrel{~}{\rho }(s)`$ in the limit $`ϵ0`$. The procedure is tedious but leads to the inverse time scale (56). The proportionality of $`\delta \stackrel{~}{g}`$ and $`\delta \stackrel{~}{\rho }`$ clearly leads to the same decay rate for $`\delta g(y,t)`$ as for $`\delta \rho (t)`$. We do note that it is not clear from this procedure that the relaxation process is actually exponential in time. If it is exponential (and there is reason to question this from the results of the exact and mesoscopic procedures), then it is necessary to perform the inverse Laplace transform more carefully. This is possible, but beyond the scope of this paper.

APPENDIX C: The function $`r(x,t)`$ and its derivatives

We first calculate the derivative of $`r(x,t)`$ with respect to $`x`$ by considering the intervals shown in Fig. 12.

* Let $`𝒫_{\mathrm{e0}}(x,\mathrm{\Delta }x,t)`$ be the probability that there is an even number of particles in $`(0,x)`$ and no particle in $`(x,x+\mathrm{\Delta }x)`$ at time $`t`$.
* Let $`𝒫_{\mathrm{e1}}(x,\mathrm{\Delta }x,t)`$ be the probability that there is an even number of particles in $`(0,x)`$ and one particle in $`(x,x+\mathrm{\Delta }x)`$ at time $`t`$.
* Let $`𝒫_{\mathrm{o1}}(x,\mathrm{\Delta }x,t)`$ be the probability that there is an odd number of particles in $`(0,x)`$ and one particle in $`(x,x+\mathrm{\Delta }x)`$ at time $`t`$.

The function $`r(x,t)`$ defined in Eq. (3) can be expressed in terms of these quantities as follows:

$`r(x+\mathrm{\Delta }x,t)`$ $`=`$ $`𝒫_{\mathrm{e0}}(x,\mathrm{\Delta }x,t)+𝒫_{\mathrm{o1}}(x,\mathrm{\Delta }x,t)+𝒪(\mathrm{\Delta }x^2)`$
$`r(x,t)`$ $`=`$ $`𝒫_{\mathrm{e0}}(x,\mathrm{\Delta }x,t)+𝒫_{\mathrm{e1}}(x,\mathrm{\Delta }x,t)+𝒪(\mathrm{\Delta }x^2).`$ (C.1)
Thus

$$r(x+\mathrm{\Delta }x,t)r(x,t)=𝒫_{\mathrm{o1}}(x,\mathrm{\Delta }x,t)𝒫_{\mathrm{e1}}(x,\mathrm{\Delta }x,t)+𝒪(\mathrm{\Delta }x^2).$$
(C.2)
In particular, by choosing $`x=0`$,

$`r(\mathrm{\Delta }x,t)r(0,t)`$ $`=`$ $`0𝒫_{\mathrm{e1}}(0,\mathrm{\Delta }x,t)+𝒪(\mathrm{\Delta }x^2)`$ (C.3)
$`=`$ $`\mathrm{\Delta }x\rho (t)+𝒪(\mathrm{\Delta }x^2).`$
Thus, the density $`\rho (t)`$ is given by

$$\rho (t)=\frac{}{x}r(x,t)|_{x=0^+},$$
(C.4)
which proves Eq. (58).

Next, we consider the intervals $`(0,x\mathrm{\Delta }x)`$, $`(x\mathrm{\Delta }x,x)`$, and $`(x,x+\mathrm{\Delta }x)`$. Let $`𝒫_{\mathrm{e01}}(x,\mathrm{\Delta }x,t)`$ be the probability that there is an even number of particles in $`(0,x\mathrm{\Delta }x)`$, no particle in $`(x\mathrm{\Delta }x,x)`$, and one particle in $`(x,x+\mathrm{\Delta }x)`$ at time $`t`$. Let $`𝒫_{\mathrm{e00}}(x,\mathrm{\Delta }x,t)`$, $`𝒫_{\mathrm{e10}}(x,\mathrm{\Delta }x,t)`$, $`𝒫_{\mathrm{e11}}(x,\mathrm{\Delta }x,t)`$, $`𝒫_{\mathrm{o00}}(x,\mathrm{\Delta }x,t)`$, $`𝒫_{\mathrm{o01}}(x,\mathrm{\Delta }x,t)`$, $`𝒫_{\mathrm{o10}}(x,\mathrm{\Delta }x,t)`$, and $`𝒫_{\mathrm{o11}}(x,\mathrm{\Delta }x,t)`$ be defined in the obvious way. The appropriate intervals are shown in Fig. 13.

Because $`g(0,t)=0`$ for all $`t>0`$, the probability that there are two particles in $`(x\mathrm{\Delta }x,x+\mathrm{\Delta }x)`$ is proportional to $`\mathrm{\Delta }x^3`$ as $`\mathrm{\Delta }x0`$. Thus,

$`r(x+\mathrm{\Delta }x,t)+r(x\mathrm{\Delta }x,t)2r(x,t)`$
$`=𝒫_{\mathrm{e10}}(x,\mathrm{\Delta }x,t)𝒫_{\mathrm{o10}}(x,\mathrm{\Delta }x,t)`$
$`𝒫_{\mathrm{e01}}(x,\mathrm{\Delta }x,t)+𝒫_{\mathrm{o01}}(x,\mathrm{\Delta }x,t)+𝒪(\mathrm{\Delta }x^3).`$ (C.5)
We derive the contribution due to diffusion of particles to the partial differential equations (59) and (63) for the evolution of $`r(x,t)`$ by considering the probability that a particle at $`x\mathrm{\Delta }x`$ at time $`t`$ diffuses out of the region $`(0,x)`$ before time $`t+\mathrm{\Delta }t`$, and the probability that a particle at $`x+\mathrm{\Delta }x`$ at time $`t`$ diffuses out of the region $`(0,x)`$ before time $`t+\mathrm{\Delta }t`$.

The probability that a particle, at $`x+\mathrm{\Delta }x`$ at time $`t`$, is in $`(0,x)`$ at time $`t+\mathrm{\Delta }t`$ is given by

$`Q(\mathrm{\Delta }x,\mathrm{\Delta }t)=\left(4\pi D\right)^{\frac{1}{2}}{\displaystyle _{\mathrm{\Delta }x}^{\mathrm{}}}dx\mathrm{exp}\left(x^2/4D\mathrm{\Delta }t\right)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{erfc}\left({\displaystyle \frac{\mathrm{\Delta }x}{(4D\mathrm{\Delta }t)^{\frac{1}{2}}}}\right)`$
Let

$$R_{\mathrm{e01}}(x,\mathrm{\Delta }x,t)=\frac{}{\mathrm{\Delta }x}𝒫_{\mathrm{e01}}(x,\mathrm{\Delta }x,t).$$
(C.7)
As $`\mathrm{\Delta }x0`$, this quantity is the probability density for finding a particle at $`x+\mathrm{\Delta }x`$ at time $`t`$, given that the number of particles in $`(0,x)`$ is even. Similarly, let

$$R_{\mathrm{e10}}(x,\mathrm{\Delta }x,t)=\frac{}{\mathrm{\Delta }x}𝒫_{\mathrm{e10}}(x,\mathrm{\Delta }x,t),R_{\mathrm{o01}}(x,\mathrm{\Delta }x,t)=\frac{}{\mathrm{\Delta }x}𝒫_{\mathrm{o01}}(x,\mathrm{\Delta }x,t),$$
(C.8)
and so on.

The time derivative of $`r(x,t)`$

$$\frac{}{t}r(x,t)=\underset{\mathrm{\Delta }t0}{lim}\frac{1}{\mathrm{\Delta }t}(r(x,t+\mathrm{\Delta }t)r(x,t)),$$
(C.9)
is found using Eq. (C.5), integrating over $`\mathrm{\Delta }x`$, and taking the limit $`\mathrm{\Delta }t0`$

$`r(x,t+\mathrm{\Delta }t)r(x,t)`$ $`=`$ $`2{\displaystyle _0^{\mathrm{}}}d\mathrm{\Delta }xQ(\mathrm{\Delta }x,\mathrm{\Delta }t)\left(R_{\mathrm{e10}}(x,\mathrm{\Delta }x,t)R_{\mathrm{o10}}(x,\mathrm{\Delta }x,t)\right)`$ (C.10)
$`+2{\displaystyle _0^{\mathrm{}}}d\mathrm{\Delta }xQ(\mathrm{\Delta }x,\mathrm{\Delta }t)\left(R_{\mathrm{e01}}(x,\mathrm{\Delta }x,t)+R_{\mathrm{o01}}(x,\mathrm{\Delta }x,t)\right)`$
$`=`$ $`2{\displaystyle _0^{\mathrm{}}}\mathrm{d}\mathrm{\Delta }xQ(\mathrm{\Delta }x,\mathrm{\Delta }t){\displaystyle \frac{}{\mathrm{\Delta }x}}(𝒫_{\mathrm{e10}}(x,\mathrm{\Delta }x,t)𝒫_{\mathrm{e10}}(x,\mathrm{\Delta }x,t)`$
$`𝒫_{\mathrm{e01}}(x,\mathrm{\Delta }x,t)+𝒫_{\mathrm{o01}}(x,\mathrm{\Delta }x,t))`$
$`=`$ $`2{\displaystyle _0^{\mathrm{}}}d\mathrm{\Delta }xQ(\mathrm{\Delta }x,\mathrm{\Delta }t){\displaystyle \frac{}{\mathrm{\Delta }x}}\left({\displaystyle \frac{^2}{x^2}}r(x,t)\mathrm{\Delta }x^2+𝒪(\mathrm{\Delta }x^3)\right)`$
$`=`$ $`{\displaystyle _0^{\mathrm{}}}d\mathrm{\Delta }x\mathrm{erfc}\left({\displaystyle \frac{\mathrm{\Delta }x}{(4D\mathrm{\Delta }t)^{\frac{1}{2}}}}\right)\left({\displaystyle \frac{^2}{x^2}}r(x,t)\mathrm{\Delta }x+𝒪(\mathrm{\Delta }x^2)\right)`$
$`=`$ $`2D{\displaystyle \frac{^2}{x^2}}r(x,t)\mathrm{\Delta }t+𝒪(\mathrm{\Delta }t^2).`$