text:
# Perturbation theory for the one-dimensional trapping reaction

## 1 Introduction

Our understanding of the nonequilibrium dynamics of many-body systems has been greatly advanced through the study of particle reaction models. Typically these models involve a number of particles of different species that undergo reactions which lead to a lack of conservation of particle number. Whilst the most obvious physical application is to the kinetics of chemical reactions , these models also enjoy mappings to a range of phenomena. These include interfacial growth , domain coarsening in magnetic media and fluids  and aggregation . Reaction-diffusion systems have also held their own as prototypes for the development of theoretical tools, such as the field theoretic renormalisation group  and exact methods in low dimensions .

In this work we consider the dynamics of the $`A+B\mathrm{}`$ reaction, i.e. a system comprising two particle species that mutually annihilate (or form an inert product) on contact. This problem was introduced in the seminal paper of Toussaint and Wilczek  as a model of monopole-antimonopole annihilation in the early universe. Since then, there have been applications to chemical kinetics and condensed matter physics—see, for example,  for a review.

In a reaction system, such as $`A+B\mathrm{}`$, one is chiefly interested in the concentration of $`A`$ and $`B`$ particles after time $`t`$ given some initial condition. The ‘traditional’ approach to the problem is to write down differential rate equations for the density of each particle species, denoted $`\rho _A(\stackrel{}{x},t)`$ and $`\rho _B(\stackrel{}{x},t)`$ respectively . One has for $`\rho _A`$

$$\frac{\rho _A(\stackrel{}{x},t)}{t}=D_A^2\rho _Ar\rho _A\rho _B$$
(1)
and the same equation for $`\rho _B`$ if one exchanges the labels $`A`$ and $`B`$. In words this equation says that both particle species diffuse with diffusion constants $`D_A`$ and $`D_B`$ and particles are removed at a rate proportional to the reaction constant $`r`$ and the joint probability to find an $`A`$-$`B`$ pair at coordinate $`\stackrel{}{x}`$ at time $`t`$. Under the assumption that the particles are well-mixed at all times, $`\rho _{A,B}(\stackrel{}{x},t)=\rho _{A,B}(t)`$ and one can neglect the Laplacian in the rate equations. In physics parlance, this is equivalent to making a mean-field approximation which becomes exact in the “reaction-limited” regime where the reaction rate is slow compared to the diffusion rate. Here we are interested in the opposite “diffusion-limited” regime.

The long-time form of $`\rho _A`$ and $`\rho _B`$ depends on whether the initial numbers of $`A`$ and $`B`$ particles are equal or not. Consider first the case $`\rho _A(0)=\rho _B(0)`$. Since each reaction event removes both an $`A`$ and a $`B`$ particle, their densities remain equal at all times. Then, we have for $`\rho (t)=\rho _A(t)=\rho _B(t)`$

$$\frac{\rho (t)}{t}=r\rho ^2$$
(2)
whose solution at large times is $`\rho (t)1/(rt)`$. It turns out that this behaviour holds only in more than four dimensions. Through rigorous bounding arguments  and renormalisation group analysis  it is now well understood that in dimensions $`d<4`$ the density decays as $`\rho (t)1/t^{d/4}`$. These results serve as a reminder of the importance of density fluctuations and correlations in low dimensions, which were explicitly ignored in the mean-field approach. Also we note that, in this problem, there exists a dynamical scaling regime in which the characteristic length-scale in the problem (here the interparticle separation $`\mathrm{}=\rho ^{1/d}`$) is related to a characteristic time-scale through a power-law.

When the initial densities of $`A`$ and $`B`$ particles are not equal, the situation is somewhat different. Consider an initial condition that has $`\rho _A(0)<\rho _B(0)`$. Then, as time progresses the $`A+B\mathrm{}`$ recation implies that the ratio $`\rho _A(t)/\rho _B(t)`$ decreases. Eventually one has a small number of $`A`$ particles in a sea of $`B`$ particles whose density remains effectively constant. In the rate-equation formalism (1) we have, under the assumption that particles are well-mixed and $`\rho _B(t)=\text{const}=\rho _B`$, that the density of $`A`$ particles decays as $`\rho _A(t)\mathrm{exp}(r\rho _Bt)`$. Again it is known that this result holds only in a space of suitably large dimensionality. In particular Bramson and Lebowitz  used bounding arguments to show that, at long times, the $`A`$ particles experience a stretched-exponential decay $`\rho _A(t)\mathrm{exp}(\lambda _dt^{d/2})`$ in fewer than two dimensions, and $`\rho _A(t)\mathrm{exp}(\lambda _2t/\mathrm{ln}t)`$ in the critical dimension $`d=2`$. Values of the constants $`\lambda _d`$ for $`d2`$ were recently reported by us  with the details to be presented elsewhere .

Clearly the special case of a single $`A`$ particle in the presence of a sea of $`B`$ particles will be governed by these asymptotics. Then, one can interpret the $`\rho _A(t)`$ as the survival probability of the $`A`$ particle. If the distribution of $`B`$ particles is homogeneous over all space, and the diffusion constants of the $`A`$ and $`B`$ particles are the same, one can also view $`\rho _A(t)`$ as the fraction of particles that have not met any other particles. Thus the reaction $`A+B\mathrm{}`$ in the limit of a low density of $`A`$ particles has been discussed under the guises of uninfected walkers  in which random walkers infect each other on contact, diffusion in the presence of traps  in which the $`B`$ particles are considered as traps for the $`A`$ particles, and predator-prey models  in which one asks for the survival of a prey (the $`A`$ particle) being ‘chased’ by diffusing predators (the $`B`$ particles). To avoid confusion, we shall use only the trapping terminology, and in this work we restrict ourselves to the case of a single $`A`$ particle and a finite number of diffusing traps ($`B`$ particles) in one spatial dimension. When we make no distinction between the particle and the traps, we shall use the generic term walkers.

Before we formally introduce our model and review some exactly known results, let us make a few comments about the relationship between the trapping reaction and persistence in diffusive systems. Persistence is a property that can be defined for any stochastic dynamical system as it is simply the probability that a random variable does not change sign between time $`0`$ and a time $`t`$. In many cases, a persistence probability $`Q(t)`$ decreases as a power-law in time, $`Q(t)t^\theta `$ in which $`\theta `$ is a persistence exponent (see, e.g.,  for a review). In the trapping reaction we can define two distinct persistence properties: (i) a site persistence, which measures the probability that a particular point in space (or site on a lattice) has not been crossed by any walker; and (ii) a walker persistence, which is the probability that a particular walker (here the $`A`$ particle) has not yet met another walker. In this terminology, one can make the transition from walker to site persistence by decreasing the particle’s diffusion constant to zero. Then the particle is static, and the probability that it survives is equal to the probability that a particular point has not been crossed.

In the following section we will show that the site persistence probability in the presence of $`N`$ traps in one dimension is known exactly to decay with an exponent $`\theta =N/2`$. On the other hand the more general walker persistence probability, where the particle has a nonzero diffusion constant, is not known exactly. The principal purpose of this paper is present the full details of a perturbation theory for the walker persistent exponent that was outlined in a previous work . As well as describing the approach in greater detail, we also extend the calculation to second order in (a quantity closely related to) the particle’s diffusion constant. We find that the persistence exponent depends on a quantity that characterises the level of asymmetry between the number of traps initially positioned to each side of the particle.

To close this introductory section, we make a few remarks as to why the case where both the particle and the traps diffuse is considerably more complicated than that in which either species is static (the case of static traps is sometimes referred to as the Donsker-Varadhan problem  for which the walker persistence decays as a stretched exponential). When the walker is static, the probability that none of the $`N`$ independently diffusing traps has crossed a particular point is simply the product of probabilities for each of the traps separately not to have crossed that point. However, any movement of the particle (in the subset of trajectories for which the particle survives) correlates the motion of the traps, and so the probabilities for each of the traps not to have met the particle are no longer independent. To understand this point more clearly, consider a set of walkers on the one-dimensional lattice. When the particle takes a step to the left (say), all the traps to the left hop one site closer to the particle in concert, whereas those to the right all hop one site away. It is this collective motion of the traps in the particle’s frame of reference that makes the problem very difficult mathematically, and as far as we are aware, exact expressions for the walker persistence exponent are known only when $`N2`$.

In the following section, we discuss these known results more quantitatively, and make the connection to isotropic diffusion in an $`N`$-dimensional ‘hyperwedge’ geometry. Then, in section 3, we describe how one sets up a perturbation theory using a path-integral approach. A judicious choice of timescale reveals that the natural parameter in the perturbation expansion has a geometric significance in the hyperwedge problem. In section 4 we go on to calculate the coefficients of the first and second order terms in the perturbation expansion, before concluding with some discussion and open questions.

## 2 Model definition and review of known results

The model we consider is defined as follows. At time $`t=0`$, a particle is placed at the origin of a one-dimensional space $`x_0(0)=0`$. Additionally, a set of mobile traps are placed at positions $`x_i(0)=y_i`$ with $`i=1,2,\mathrm{},N`$ such that $`N`$ is the number of traps. We will at times need to distinguish between those traps initially placed to the left and right of the particle. To this end we introduce the shorthand $`\sigma _i=\mathrm{sgn}(y_i)`$ and the quantities $`N_L`$ and $`N_R`$ which are the number of particles that have $`\sigma _i<0`$ and $`\sigma _i>0`$ respectively.

The system evolves through the independent diffusion of the particle and traps. We shall take the diffusion constant of the particle to be $`D^{}`$ and that of the traps to be $`D`$. Then, we can write a Langevin equation

$$\dot{x}_i=\eta _i(t)$$
(3)
in which $`\eta _i(t)`$ is a Gaussian white noise with zero mean and correlator

$$\eta _i(t)\eta _j(t^{})=2D_i\delta _{i,j}\delta (t^{}t),$$
(4)
where $`D_0=D^{}`$ and $`D_i=D`$ for $`i0`$.

If, at any time, a trap coordinate $`x_i(t)`$ coincides with the particle coordinate $`x_0(t)`$, the particle becomes ‘trapped’. We seek an expression for the survival probability $`Q(t)`$, i.e. the probability that the particle has not yet fallen into a trap.

It is convenient to work in the particle’s frame of reference and so we introduce the relative coordinates $`X_i=x_ix_0`$ and $`Y_i=y_i`$. Then the connection to persistence becomes clearer, since $`Q(t)`$ is the probability that none of the variables $`X_i`$ has changed sign until a time $`t`$. From (3) and (4) we obtain the underlying Langevin equations in the variables $`X_i`$. These read

$$\dot{X}_i=\xi _i(t)$$
(5)
in which the noise $`\xi _i(t)`$ has mean zero and correlator

$$\xi _i(t)\xi _j(t^{})=2(D\delta _{i,j}+D^{})\delta (t^{}t)2D_{ij}\delta (t^{}t).$$
(6)
Here we have introduced the matrix $`D_{ij}=D\delta _{i,j}+D^{}`$ which expresses the way in which the trap’s trajectories are correlated in the frame of reference of a diffusing particle (i.e. when $`D^{}0`$).

In the introduction, we remarked that the survival probability $`Q(t)`$ can be calculated exactly when $`D^{}=0`$. Although this is a classic result (see, e.g., ), we present the details of the calculation here as it is will be needed when calculating the coefficients in the perturbation expansion that we construct for general $`D^{}`$ in the next section.

We begin by writing down the Fokker-Planck equation for the probability distribution function $`G_N(\{X\},t|\{Y\})`$ that holds when $`D^{}=0`$. It reads

$$\frac{}{t}G_N(\{X\},t|\{Y\})=D\underset{i=1}{\overset{N}{}}\frac{^2}{X_i^2}G_N(\{X\},t|\{Y\})$$
(7)
and is subject to the absorbing boundary condition $`G_N=0`$ if any of the $`X_i`$ are zero (i.e., when a trap meets the stationary particle). Clearly, since the traps diffuse independently, the solution of this equation is just

$$G_N(\{X\},t|\{Y\})=\underset{i=1}{\overset{N}{}}G_1(X_i,t|Y_i)$$
(8)
where $`G_1(X,t;Y)`$ is the probability for a walker starting at $`Y`$ to be at $`X`$ and not to have crossed the origin up to time $`t`$. For a particle starting at $`Y>0`$, this quantity is given by

$$G_1(X,t|Y)=\frac{1}{\sqrt{4\pi Dt}}\left[\mathrm{exp}\left(\frac{(XY)^2}{4Dt}\right)\mathrm{exp}\left(\frac{(X+Y)^2}{4Dt}\right)\right]$$
(9)
for $`X>0`$. We note that this is the correct solution since the two terms on the right-hand side separately obey (7) with $`N=1`$ and this particular combination of solutions satisfies the boundary condition $`G_1(0,t|Y)=0`$.

To obtain the probability for a walker starting at $`Y`$ not to have crossed the origin until time $`t`$ we simply integrate (9) over all $`X>0`$ to find

$$Q_1(t|Y)=\mathrm{erf}\left(\frac{Y}{2\sqrt{Dt}}\right).$$
(10)
Note that, by symmetry, we must have $`Q_1(t|Y)=Q_1(t|Y)`$ so the general result is

$$Q_1(t|Y)=\mathrm{erf}\left(\frac{|Y|}{2\sqrt{Dt}}\right).$$
(11)
Finally, using (8) we find the probability that none of the traps has met the particle is

$$Q_N(t)=\underset{i=1}{\overset{N}{}}\mathrm{erf}\left(\frac{|Y_i|}{2\sqrt{Dt}}\right)\left(\frac{1}{\pi Dt}\right)^{N/2}\underset{i=1}{\overset{N}{}}|Y_i|$$
(12)
in which we have used the fact that

$$\mathrm{erf}x=\frac{2}{\sqrt{\pi }}\left(z\frac{z^3}{3}+\mathrm{}\right)$$
(13)
to obtain a form for $`Q_N(t)`$ valid at large times. The persistence exponent $`\theta `$ associated with the survival probability $`Q(t)`$ can be defined formally as

$$\theta =\underset{t\mathrm{}}{lim}t\frac{}{t}\mathrm{ln}Q(t)$$
(14)
and so for the case $`D^{}=0`$ we have from (12) that $`\theta =N/2`$.

It is instructive to consider an approach that has been used (see  for a more in-depth overview) to find the persistence exponent $`\theta `$ exactly for the case $`N=2`$ and $`D^{}0`$ since it will reveal a connection with the perturbation series that we derive in the next section.

The essence of this approach is to perform a second coordinate transformation that renders the correlated diffusion described by equations (5) and (6) isotropic in the $`N`$-dimensional space spanned by the coordinates $`X_i`$. To perform this transformation one considers the diffusion matrix $`D_{ij}=D\delta _{i,j}+D^{}`$. The eigenvectors of this matrix indicate the directions in which the $`N`$-dimensional diffusion is independent whilst the corresponding eigenvalues give the diffusion constants in those directions. One finds that one of the eigenvectors is $`(1,1,1,\mathrm{},1)`$ and has an eigenvalue $`D+ND^{}`$. The remaining eigenvectors are degenerate with eigenvalue $`D`$ and hence form a basis set in the $`(N1)`$-dimensional space perpendicular to $`(1,1,1,\mathrm{},1)`$. In order to make the diffusion isotropic, one must scale in the $`(1,1,1,\mathrm{},1)`$ direction by a factor $`(1+ND^{}/D)^{1/2}`$.

Whilst this coordinate transformation gives rise to simple, isotropic diffusion in the $`N`$-dimensional space, the boundary conditions become more complicated. In the space spanned by the variables $`X_i`$, the absorbing boundary is constructed from the set of orthogonal planes that have $`X_i=0`$ with $`i=1,2,\mathrm{},N`$. The coordinate transformation just described rotates these planes so that they are no longer orthogonal. It is a straightforward exercise in linear algebra to determine that the angle between two planes $`i`$ and $`j`$ in the transformed coordinate system is $`\mathrm{\Theta }_{ij}=\mathrm{arccos}(\sigma _i\sigma _j\alpha _N)`$ in which the parameter $`\alpha _N`$ takes the form

$$\alpha _N=\frac{D^{}}{D+(N1)D^{}}.$$
(15)
Hence two planes that correspond to traps that start on opposite sides of the particle close up whereas two planes that correspond to traps starting on the same side of the particle open out.

To determine the persistence exponent $`\theta `$ given by equation (14) it is necessary to solve the diffusion equation in an $`N`$-dimensional wedge geometry. As far as we know, results are available only for the case of a two-dimensional wedge . Then, the particle’s survival probability decays with the exponent

$$\theta =\frac{\pi }{2\mathrm{\Theta }_{12}}.$$
(16)
For the case where both traps initially surround the particle, we have $`\sigma _1\sigma _2=1`$ and so

$$\theta =\frac{\pi }{2\mathrm{arccos}\alpha _2}.$$
(17)
Similarly when both traps are initially positioned to one side of the particle, we have $`\sigma _1\sigma _2=1`$ and thus

$$\theta =\frac{\pi }{2(\pi \mathrm{arccos}\alpha _2)}$$
(18)
since $`\mathrm{arccos}(x)=\pi \mathrm{arccos}(x)`$.

Although this approach has thus far proved intractable for $`N>2`$, we have spent some time outlining it because, in the perturbation theory that we construct below, it turns out that the natural expansion parameter is the quantity $`\alpha _N`$ introduced above. Note that $`\alpha _N=0`$ corresponds to the case where the particle is static and the $`N`$-dimensional hyperwedge is in fact a corner of an $`N`$-dimensional hypercube. Thus the expansion for small $`\alpha _N`$ derived below using an alternative path-integral approach gives results for diffusion in a geometry perturbed from an $`N`$-dimensional hypercube.

## 3 Perturbation expansion via a path-integral approach

The starting point in the derivation of a perturbation expansion for the persistence exponent $`\theta `$ is to note that the statistical weight a particular set of trajectories $`\stackrel{}{X}(t)=[X_1(t),X_2(t),\mathrm{},X_N(t)]`$ governed by equations (5) and (6) can be written in the form $`\mathrm{exp}(S[\stackrel{}{X}])`$. In this expression, $`S[\stackrel{}{X}]`$ is an ‘action’ functional defined as

$$S[\stackrel{}{X}]=\frac{1}{4}\underset{i,j=1}{\overset{N}{}}[D^1]_{ij}_0^tt^{}\dot{X}_i(t^{})\dot{X}_j(t^{})$$
(19)
in which $`D^1`$ is the inverse of the matrix $`D`$ in equation (6). It is straightforward to show that

$$[D^1]_{ij}=\frac{1}{D}\left(\delta _{ij}\frac{D^{}}{D+ND^{}}\right).$$
(20)
It is helpful at this stage to introduce a rescaled time variable

$$\tau =\frac{D+ND^{}}{D+(N1)D^{}}Dt$$
(21)
such that the action becomes

$$S[\stackrel{}{X}]=S_0\frac{D^{}}{D+(N1)D^{}}S_1=S_0\alpha _NS_1$$
(22)
in which the diagonal part $`S_0`$ and off-diagonal part $`S_1`$ of the action are given by

$`S_0[\stackrel{}{X}]`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \underset{i}{}}{\displaystyle _0^\tau }\tau ^{}[\dot{X}_i(\tau ^{})]^2`$ (23)
$`S_1[\stackrel{}{X}]`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \underset{ij}{}}{\displaystyle _0^\tau }\tau ^{}\dot{X}_i(\tau ^{})\dot{X}_j(\tau ^{})`$ (24)
where now a dot denotes differentiation to the rescaled time variable $`\tau `$. Note that when $`\alpha _N=0`$ the action $`S[\stackrel{}{X}]=S_0[\stackrel{}{X}]`$ which gives the statistical weight for $`N`$ walkers executing independent diffusive motion with diffusion constant unity. We will treat the off-diagonal part (that describes the correlations between the walkers) as a perturbation with the quantity $`\alpha _N`$ given by equation (15) as the expansion parameter, as advertised.

To this end, observe that out of all possible trajectories $`\stackrel{}{X}(\tau )`$, only a subset corresponds to the particle surviving until time $`\tau `$. These trajectories are defined through the persistence condition that $`\mathrm{sgn}X(\tau ^{})=\mathrm{sgn}X(0)`$ for $`0<\tau ^{}\tau `$. We shall use the symbol $`_R𝒟\stackrel{}{X}(\tau )`$ to denote integration over this restricted set of surviving trajectories. Then, the probability that the particle survives is

$$Q(\tau )=\frac{_R𝒟\stackrel{}{X}(\tau )\mathrm{exp}(S[\stackrel{}{X}])}{𝒟\stackrel{}{X}(\tau )\mathrm{exp}(S[\stackrel{}{X}])}$$
(25)
in which the integral of the weights over the restricted trajectories is normalised by the integral over all possible trajectories.

Now, note that the survival probability $`Q_N(\tau )`$ of a stationary particle in the presence of $`N`$ traps diffusing independently with diffusion constant unity is given by the path-integral expression

$$Q_N(\tau )=\frac{_R𝒟\stackrel{}{X}(\tau )\mathrm{exp}(S_0[\stackrel{}{X}])}{𝒟\stackrel{}{X}(\tau )\mathrm{exp}(S_0[\stackrel{}{X}])}.$$
(26)
Combining this with (25) we find that

$`Q(\tau )=Q_N(\tau ){\displaystyle \frac{_R𝒟\stackrel{}{X}(\tau )\mathrm{exp}(S[\stackrel{}{X}])}{_R𝒟\stackrel{}{X}(\tau )\mathrm{exp}(S_0[\stackrel{}{X}])}}{\displaystyle \frac{𝒟\stackrel{}{X}(\tau )\mathrm{exp}(S_0[\stackrel{}{X}])}{𝒟\stackrel{}{X}(\tau )\mathrm{exp}(S[\stackrel{}{X}])}}`$
$`=Q_N(\tau ){\displaystyle \frac{\mathrm{}^{\alpha _NS_1[\stackrel{}{X}]}_R}{\mathrm{}^{\alpha _NS_1[\stackrel{}{X}]}_U}}`$ (27)
in which the notation $``$ indicates an average over paths weighted by $`\mathrm{exp}(S_0[\stackrel{}{X}])`$ and the subscripts $`R`$ and $`U`$ denote the restricted and unrestricted ensembles of paths. (Recall that paths that cross the origin are excluded from the former).

To obtain a perturbative expression for the persistence exponent $`\theta `$ we make use of (14) and (3) to find, in the rescaled time variable $`\tau `$,

$`\theta =\underset{\tau \mathrm{}}{lim}\tau {\displaystyle \frac{}{\tau }}\mathrm{ln}Q(\tau )`$
$`=\underset{\tau \mathrm{}}{lim}\tau {\displaystyle \frac{}{\tau }}(\mathrm{ln}Q_N(\tau )+\mathrm{ln}\mathrm{}^{\alpha _NS_1[\stackrel{}{X}]}_R\mathrm{ln}\mathrm{}^{\alpha _NS_1[\stackrel{}{X}]}_U).`$ (28)
The asymptotic form of $`Q_N(\tau )`$ is known from equation (12). To compute the averages of $`S_1[\stackrel{}{X}]`$ perturbatively, we make use of the fact that $`\mathrm{ln}\mathrm{exp}\lambda x`$ defines the generating function of the cumulants of $`x`$. Specifically

$$\mathrm{ln}\mathrm{}^{\lambda x}=\lambda x+\frac{\lambda ^2}{2}\left(x^2x^2\right)+\mathrm{O}\left(\lambda ^3\right).$$
(29)
Putting all this together, we find that $`\theta `$ has a series expansion

$$\theta =\frac{N}{2}\underset{\tau \mathrm{}}{lim}\tau \frac{}{\tau }\left[\alpha _N\left(S_1_RS_1_U\right)+\frac{\alpha _N^2}{2}\left(S_1^2_RS_1_R^2S_1^2_U+S_1_U^2\right)\right]$$
(30)
to second order in the parameter $`\alpha _N`$ defined in equation (15). Of course, one could go to higher order by including more terms from the cumulant expansion (29).

## 4 Walker persistence exponent to second order

As a first step in computing the mean and variance of $`S_1`$ in the two ensembles we shall determine the structure of these quantities. For either the restricted or the unrestricted average we have

$$S_1=\frac{1}{4}_{}^{}{}_{i,j}{}^{}_0^\tau \tau ^{}\dot{X}_i(\tau ^{})\dot{X}_j(\tau ^{})$$
(31)
in which the notation $`^{}`$ means the sum over combinations of the indices such that each index is different. (Note that all permutations of a set of distinct indices are included in this sum). Since we are averaging over ensembles of trajectories of independently diffusing particles, and the indices $`i`$ and $`j`$ are always different, we have

$$S_1=\frac{1}{4}_{}^{}{}_{i,j}{}^{}_0^\tau \tau ^{}\dot{X}_i(\tau ^{})\dot{X}_j(\tau ^{}).$$
(32)
Hence to calculate the mean of $`S_1`$ we must find the mean velocity of a random walk with diffusion constant unity in each of the ensembles.

For the variance of $`S_1`$ we have

$`S_1^2S_1^2={\displaystyle \frac{1}{16}}{\displaystyle _{}^{}{}_{i,j}{}^{}}{\displaystyle _{}^{}{}_{k,\mathrm{}}{}^{}}{\displaystyle _0^\tau }\tau _1{\displaystyle _0^\tau }\tau _2(\dot{X}_i(\tau _1)\dot{X}_j(\tau _1)\dot{X}_k(\tau _2)\dot{X}_{\mathrm{}}(\tau _2)`$
$`\dot{X}_i(\tau _1)\dot{X}_j(\tau _1)\dot{X}_k(\tau _2)\dot{X}_{\mathrm{}}(\tau _2)).`$ (33)
The terms in the double summation can be divided into three groups. In the first, all four indices are different, and two averages factorise to cancel. A second set of terms has two indices the same, the other two different (there are four ways this property can be satisfied). Finally there are two ways to arrange for the indices to comprise two pairs the same. Thus the variance can be written as

$`S_1^2S_1^2={\displaystyle \frac{1}{16}}{\displaystyle _0^\tau }\tau _1{\displaystyle _0^\tau }\tau _2`$
$`[4{\displaystyle _{}^{}{}_{i,j,k}{}^{}}(\dot{X}_i(\tau _1)\dot{X}_i(\tau _2)\dot{X}_i(\tau _1)\dot{X}_i(\tau _2))\dot{X}_j(\tau _1)\dot{X}_k(\tau _2)`$
$`+2{\displaystyle _{}^{}{}_{i,j}{}^{}}(\dot{X}_i(\tau _1)\dot{X}_i(\tau _2)\dot{X}_j(\tau _1)\dot{X}_j(\tau _2)\dot{X}_i(\tau _1)\dot{X}_i(\tau _2)\dot{X}_j(\tau _1)\dot{X}_j(\tau _2))]`$ (34)
which reveals that we must calculate not just the mean velocity $`\dot{X}_i(\tau ^{})`$ but also the velocity correlation function $`\dot{X}_i(\tau _1)\dot{X}_i(\tau _2)`$ over the two ensembles.

The unrestricted averages are easily found. Since each of the walkers performs diffusion with a diffusion constant of unity in the rescaled time $`\tau `$, we have $`\dot{X}_i(\tau ^{})=\eta _i(\tau ^{})`$ where $`\eta _i(\tau ^{})`$ is a Gaussian white noise with zero mean and a correlator $`\eta _i(\tau _1)\eta _i(\tau _2)=2\delta (\tau _2\tau _1)`$. Hence

$$\dot{X}_i(\tau ^{})_U=0\text{and}\dot{X}_i(\tau _1)\dot{X}_i(\tau _2)_U=2\delta (\tau _2\tau _1).$$
(35)
The restricted averages involve more work, since we must calculate the mean velocity of a walker at time $`\tau ^{}`$ given that it survives until a later time $`\tau `$; similarly we must include the fact that the walker survives until a specified time $`\tau `$ in the calculation of the velocity correlation function. Before getting into the details, therefore, let us propose their general form and the implications of such a form in context of the expression (30). On dimensional grounds, we suggest that

$`\dot{X}_i(\tau ^{})_R={\displaystyle \frac{1}{\sqrt{\tau ^{}}}}F_i(\tau _0/\tau ^{},\tau ^{}/t)`$ (36)
$`\dot{X}_i(\tau _1)\dot{X}_i(\tau _2)_R=2\left[\delta (\tau _2\tau _1)+{\displaystyle \frac{1}{\tau _+}}G_i(\tau _0/\tau _{},\tau _{}/\tau _+,\tau _+/\tau )\right]`$ (37)
in which $`\tau _0`$ is some early timescale in the problem and $`\tau _{}`$ ($`\tau _+`$) are the smaller (respectively larger) of $`\tau _1`$ and $`\tau _2`$. Guided by (35), and confirmed by explicit calculation to be presented below, we have included a delta-function contribution in the two-time correlation function.

The upshot of this is that all the integrals in (32) and (4) reduce to one of two forms. These are

$`{\displaystyle _0^\tau }{\displaystyle \frac{\tau ^{}}{\tau ^{}}}f(\tau _0/\tau ^{},\tau ^{}/\tau )`$ (38)
$`{\displaystyle _0^\tau }{\displaystyle \frac{\tau _+}{\tau _+^2}}{\displaystyle _0^{\tau _+}}\tau _{}g(\tau _0/\tau _{},\tau _{}/\tau _+,\tau _+/\tau )`$ (39)
in which the functions $`f`$ and $`g`$ are various combinations of the velocity and velocity-correlation functions. It is not the values of these integrals that we are interested in, but the limit that appears in the perturbation expansion (30). It is straightforward to show that

$`\underset{\tau \mathrm{}}{lim}\tau {\displaystyle \frac{}{\tau }}{\displaystyle _0^\tau }{\displaystyle \frac{\tau ^{}}{\tau ^{}}}f(\tau _0/\tau ^{},\tau ^{}/\tau )=f(0,0)`$ (40)
$`\underset{\tau \mathrm{}}{lim}\tau {\displaystyle \frac{}{\tau }}{\displaystyle _0^\tau }{\displaystyle \frac{\tau _+}{\tau _+^2}}{\displaystyle _0^{\tau _+}}\tau _{}g(\tau _0/\tau _{},\tau _{}/\tau _+,\tau _+/\tau )={\displaystyle _0^1}ug(0,u,0).`$ (41)
The fact that zeros appear in the arguments of the functions $`f`$ and $`g`$ implies that when calculating $`\dot{X}_i(\tau ^{})_R`$ and $`\dot{X}_i(\tau _1)\dot{X}_i(\tau _2)_R`$ we need only determine their forms in the regime $`\tau _0\tau ^{},\tau _1,\tau _2\tau `$ to obtain the expansion (30).

We now indicate how to calculate these quantities. First, note that $`\dot{X}_i(\tau ^{})_R=\frac{}{\tau ^{}}X_i(\tau ^{})_R`$. This reduces the problem to the calculation of the mean position of a walker at time $`\tau ^{}`$ taking into account that it does not cross the origin (survives) until at least a time $`\tau `$. For the case where a walker’s initial position $`Y_i>0`$, this quantity can be expressed as

$$X_i(\tau ^{})_R=_0^{\mathrm{}}X_iX_i\frac{P_1(\tau ;X_i,\tau ^{}|Y_i)}{Q_1(\tau |Y_i)}$$
(42)
in which $`P_1(\tau ;X_i,\tau ^{}|Y_i)`$ is the probability that a single walker visits $`X_i`$ at time $`\tau ^{}`$ and does not visit the origin at any time $`0\tau ^{}\tau `$. In order to average $`X_i`$ with respect to this distribution, it has been normalised through the factor $`Q_1(\tau |Y_i)`$, which is the probability the walker survives until time $`\tau `$ and is given by equation (11).

Since the diffusion process is Markovian, we can perform the factorisation $`P_1(\tau ;X_1,\tau ^{}|Y_i)=Q_1(\tau \tau ^{}|X_i)G_1(X_i,\tau ^{}|Y_i)`$, where we recall that $`G_1(X_i,\tau ^{}|Y_i)`$ is that probability that a walker is at $`X_i`$ at time $`\tau ^{}`$ given that it started at $`Y_i`$ and whose form is given by equation (9). Inserting the explicit expressions for $`Q_1(\tau |Y)`$ and $`G_1(X,\tau |Y)`$ into (42) we find

$$X_i(\tau ^{})_R=_0^{\mathrm{}}X_iX_i\frac{\mathrm{erf}\left(\frac{X_i}{2\sqrt{\tau \tau ^{}}}\right)\left[\mathrm{exp}\left(\frac{(X_iY_i)^2}{4\tau ^{}}\right)\mathrm{exp}\left(\frac{(X_i+Y_i)^2}{4\tau ^{}}\right)\right]}{\sqrt{4\pi \tau ^{}}\mathrm{erf}\left(\frac{Y_i}{2\sqrt{\tau }}\right)}.$$
(43)
Now we can use the fact that when calculating the exponent $`\theta `$ using (30) it is sufficient to know $`X_i(\tau ^{})_R`$ in the regime $`\tau _0\tau ^{}\tau `$. We take $`\tau _0=Y_i^2`$ (recalling that the rescaled time variable has units of length squared) and on expanding the integrand in the previous expression, we find

$$X_i(\tau ^{})_R=\frac{1}{2\sqrt{\pi }(\tau ^{})^{3/2}}_0^{\mathrm{}}X_iX_i^3\mathrm{exp}\left(\frac{X_i^2}{4\tau ^{}}\right)=4\sqrt{\frac{\tau ^{}}{\pi }}$$
(44)
in this regime. Finally, we differentiate to obtain

$$\dot{X}_i(\tau ^{})_R=\frac{2}{\sqrt{\pi \tau ^{}}}$$
(45)
when $`\tau _0\tau ^{}\tau `$ and $`Y_i>0`$. Since we are working at times much larger than $`\tau _0=Y_i^2`$, the dependence on the magnitude of $`Y_i`$ has dropped out of this expression. However, the sign of the average velocity clearly changes if we consider a walker with initial position $`Y_i<0`$. Hence the correct limiting form of the function $`F_i(\tau _0/\tau ^{},\tau ^{}/\tau )`$ in equation (36) is

$$F_i(0,0)=\frac{2\sigma _i}{\sqrt{\pi }}$$
(46)
in which $`\sigma _i=\mathrm{sgn}(Y_i)`$.

The two-point velocity correlation function $`\dot{X}_i(\tau _1)\dot{X}_i(\tau _2)_R`$ is obtained in a similar way. We defer the details of the calculation to the appendix, quoting here only the result over the required range $`\tau _0\tau _1,\tau _2\tau `$ which is

$$\dot{X}_i(\tau _1)\dot{X}_i(\tau _2)_R=2\left[\delta (\tau _2\tau _1)+\frac{2}{\pi \tau _+}\sqrt{\frac{\tau _+\tau _{}}{\tau _{}}}\right].$$
(47)
Recall that $`\tau _{}=\mathrm{min}\{\tau _1,\tau _2\}`$ and $`\tau _+=\mathrm{max}\{\tau _1,\tau _2\}`$. Note that this expression is in agreement with the proposal (37), and reveals the limiting form of the function $`G`$ to be

$$G_i(0,\tau _{}/\tau _+,0)=\frac{2}{\pi }\sqrt{\frac{\tau _+}{\tau _{}}1}.$$
(48)
Once again, the dependence of this quantity on the walker’s initial position $`Y_i`$ has dropped out; furthermore, the two-point correlation function is independent of the sign of $`Y_i`$, so we can drop the subscript $`i`$ on $`G`$.

We now consider the first order term in (30). We have

$`S_1_RS_1_U`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle _{}^{}{}_{i,j}{}^{}}{\displaystyle _0^\tau }\tau ^{}\dot{X}_i(\tau ^{})_R\dot{X}_j(\tau ^{})_R`$ (49)
$`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle _{}^{}{}_{i,j}{}^{}}{\displaystyle _0^\tau }{\displaystyle \frac{\tau ^{}}{\tau ^{}}}F_i(\tau _0/\tau ^{},\tau ^{}/t)F_j(\tau _0/\tau ^{},\tau ^{}/t).`$
The unrestricted averages have disappeared from the right-hand side of this expression because $`\dot{X}_i_U=0`$. Now, using equation (40), we find that the limiting form of this expression required in the expansion (30) reads

$$\underset{\tau \mathrm{}}{lim}\tau \frac{}{\tau }\left(S_1_RS_1_U\right)=\frac{1}{4}_{}^{}{}_{i,j}{}^{}F_i(0,0)F_j(0,0)=\frac{1}{\pi }_{}^{}{}_{i,j}{}^{}\sigma _i\sigma _j.$$
(50)
To perform the remaining summation, we note that $`\sigma _i\sigma _j`$ is equal to $`+1`$ if walkers $`i,j`$ both start on the same side of the origin, and is equal to $`1`$ if they start on opposite sides of the origin. If there are $`N_L`$ ($`N_R`$) walkers initially to the left (right), there are $`N_L(N_L1)+N_R(N_R1)`$ positive terms in the summation and $`2N_LN_R`$ negative terms. Thus one finds that the coefficient of $`\alpha _N`$ in (30) is

$$\underset{\tau \mathrm{}}{lim}\tau \frac{}{\tau }\left(S_1_RS_1_U\right)=\frac{1}{\pi }(\mathrm{\Delta }^2N)$$
(51)
in which $`\mathrm{\Delta }=N_LN_R`$.

Let us now turn to the second-order term in (30). From (4) we have the general form

$`{\displaystyle \frac{1}{2}}\left(S_1^2_RS_1_R^2S_1^2_U+S_1_U^2\right)=`$
$`{\displaystyle \frac{1}{32}}{\displaystyle _0^\tau }\tau _1{\displaystyle _0^\tau }\tau _2\left[4{\displaystyle _{}^{}{}_{i,j,k}{}^{}}A_{ijk}(\tau _0/\tau _1,\tau _1,\tau _2,\tau _2/\tau )+2{\displaystyle _{}^{}{}_{i,j}{}^{}}B_{ij}(\tau _0/\tau _1,\tau _1,\tau _2,\tau _2/\tau )\right].`$ (52)
In the regime $`\tau _0\tau _1,\tau _2t`$ we have from (35), (45) and (47)

$`A_{ijk}(0,\tau _{},\tau _+,0)={\displaystyle \frac{8\sigma _j\sigma _k}{\pi ^2}}\left[{\displaystyle \frac{\pi \delta (\tau _2\tau _1)}{\tau _2}}+{\displaystyle \frac{2}{\tau _+^2}}{\displaystyle \frac{\tau _+}{\tau _{}}}\left(\sqrt{1{\displaystyle \frac{\tau _{}}{\tau _+}}}1\right)\right]`$ (53)
$`B_{ij}(0,\tau _{},\tau _+,0)={\displaystyle \frac{16}{\pi ^2\tau _+^2}}.`$ (54)
As with the first order term, we calculate the limit in (30) by making use of the results (40) and (41). This procedure yields the result

$`\underset{\tau \mathrm{}}{lim}\tau {\displaystyle \frac{}{\tau }}{\displaystyle \frac{1}{2}}\left(S_1^2_RS_1_R^2S_1^2_U+S_1_U^2\right)=`$
$`{\displaystyle \frac{1}{\pi ^2}}\left[\left(\pi +4{\displaystyle _0^1}u{\displaystyle \frac{\sqrt{1u}1}{u}}\right){\displaystyle _{}^{}{}_{i,j,k}{}^{}}\sigma _j\sigma _k2{\displaystyle _{}^{}{}_{i,j}{}^{}}{\displaystyle _0^1}u\right].`$ (55)
To obtain the final expression for the expansion (30) we use the fact that

$$_0^1u\frac{\sqrt{1u}1}{u}=2(\mathrm{ln}21).$$
(56)
Also, we must consider the possible ways of choosing three particles labelled $`i,j,k`$ from a set of $`N`$ particles comprising a number $`N_L`$ with $`\sigma =1`$ and $`N_R`$ with $`\sigma =1`$. These considerations lead us to

$$_{}^{}{}_{i,j,k}{}^{}\sigma _j\sigma _k=(N2)(\mathrm{\Delta }^2N).$$
(57)
Similarly the number of terms in the sum $`_{i,j}^{}`$ is just $`N(N1)`$ which finally allows us to write down an expression for the exponent $`\theta `$ which constitutes the main result of this paper. It reads

$`\theta ={\displaystyle \frac{N}{2}}+{\displaystyle \frac{1}{\pi }}(N\mathrm{\Delta }^2)\alpha _N`$
$`+{\displaystyle \frac{1}{\pi ^2}}\left[(N2)(N\mathrm{\Delta }^2)(\pi 8(1\mathrm{ln}2))+2N(N1)\right]\alpha _N^2.`$ (58)
As a check of this formula, we consider the exactly solvable cases of $`N=1`$ and $`N=2`$. When $`N=1`$, one can simply transform to a frame in which the particle is stationary. Then, its survival probability is given by equation (11) with $`D`$ replaced by $`D+D^{}`$ which implies a persistence exponent of $`\theta =\frac{1}{2}`$. Since $`\mathrm{\Delta }^2=1`$ when $`N=1`$, we find that the perturbative expression for $`\theta `$ (4) also gives the exact result $`\theta =\frac{1}{2}`$.

When $`N=2`$, there are two possible arrangements of the walkers. Either they surround the particle, in which case $`\mathrm{\Delta }=0`$ and the exact exponent is given by equation (17), or they are both on one side of the particle, which has $`\mathrm{\Delta }=\pm 2`$ and $`\theta `$ is given by equation (18). Expanding these exact expressions as series in $`\alpha _2`$ one finds

$$\theta =1\pm \frac{2}{\pi }\alpha _2+\frac{4}{\pi ^2}\alpha _2^2+\mathrm{}$$
(59)
where the plus sign is taken for the case $`\mathrm{\Delta }=0`$ and the minus sign for $`\mathrm{\Delta }=\pm 2`$. Thus to second order in $`\alpha _2`$ we find agreement between the perturbation series (4) and the exactly known results.

## 5 Discussion and conclusion

In this work we have studied the survival probability of a diffusing particle in the presence of a finite number of mobile traps. In the absence of a complete exact solution for the problem, we have devised a method for calculating the persistence exponent $`\theta `$ \[defined by equation (14)\] valid when the particle’s diffusion constant is small. In order to calculate the coefficients in the expansion, it is necessary to find velocity correlation functions for a random walker in the presence of an absorbing boundary. In this paper, we calculated the exponent $`\theta `$ to second order, culminating in the expression (4). We stress that, in principle, one could go to higher order, although the amount of work required to obtain $`n`$-time velocity correlation functions is likely to increase rapidly with $`n`$.

We wish to make a few remarks about our result (4). First, as we have already noted, the expansion parameter $`\alpha _N`$ defined by equation (15) has an interpretation in the context of the (hyper)wedge geometry that has been previously used to obtain results for this problem. Although the connection appears at first sight intriguing, it is probably no coincidence since in the two approaches the parameter arises through a rescaling that gives rise to uncorrelated diffusion of the traps with diffusion constant unity. However we do see that a system with perturbed equations of motion (treated using the path-integral approach) can be related to one with the original equations of motion but perturbed boundaries (i.e. the hyperwedge geometry).

Secondly, we find that the expansion to second order (4) depends on two parameters: the total number of traps in the system $`N`$ and the square difference $`\mathrm{\Delta }^2`$ in the number of traps initially placed on either side of the particle. The origin of this observation lies in the fact that at sufficiently large times the initial position of a particular trap is unimportant. Of course, the number of traps initially positioned to each side of the particle plays an important role, and enters into the perturbation expansion through summations in which some terms are positive and some negative depending on the initial arrangement of the traps. We would expect, therefore, the combination $`N\mathrm{\Delta }^2`$ to appear in higher-order terms in the expansion.

We also note from the expansion (4) that when the asymmetry is small, i.e. $`\mathrm{\Delta }^2<N`$, the exponent $`\theta `$ is an increasing function of the particle’s diffusion constant $`D^{}`$ at fixed trap diffusion constant $`D`$. Hence, in this near-symmetric situation, the particle is more likely to survive longer if it is at rest than if it is slowly diffusing. Of course, we do not know if $`\theta `$ continues to increase for large $`D^{}`$ since we only have the first two terms in the perturbation expansion. However, our intuition suggests that this be the case.

The case of large asymmetry, $`\mathrm{\Delta }^2>N`$ is more subtle since one of the terms in the coefficient of $`\alpha _N^2`$ is always positive. We learn a little by defining the critical asymmetry $`\mathrm{\Delta }_c`$ at fixed $`N`$ and $`\alpha _N`$ to be that for which the persistence exponent $`\theta =\frac{N}{2}`$. Then for $`\mathrm{\Delta }>\mathrm{\Delta }_c`$ one has a regime in which the particle’s survival probability is not (locally) maximised by remaining still. One can understand this phenomenon from the extreme case of all traps being on one side of the particle ($`\mathrm{\Delta }=\pm N`$). Then, there is some probability for the particle to diffuse into the region which is devoid of traps, thereby increasing its survival probability. By rearranging the expansion (4) one finds for small $`\alpha _N`$ that

$$\mathrm{\Delta }_c^2=N+\frac{2N(N1)}{\pi }\alpha _N+\mathrm{O}(\alpha _N^2).$$
(60)
Of course, it is again difficult to make any statements for large $`\alpha _N`$ given only the first two terms of the perturbation series.

The final comment we wish to make about the form of (4) is that, although the quantity $`\alpha _N`$ seems to be the natural expansion parameter in the problem, the presence of a nontrivial coefficient in the second order term is not suggestive of there being a simple, closed expression valid for all values of the parameters $`N`$, $`\mathrm{\Delta }`$, $`D`$ and $`D^{}`$. Nevertheless, one might find that the marginal case $`N=\mathrm{\Delta }^2`$ is in some way much easier to understand.

The trapping reaction studied in this work suggests a couple of generalisations. Firstly there is the question of how to treat higher dimensions. In a forthcoming work , we extend the methods introduced in  for the case of an infinite sea of traps to dimensions greater than one. It would be interesting to see if the perturbation theory could be applied to higher dimensions to treat the case of a finite number of traps. Finally, we have studied in this work the trapping reaction on an infinite line. It appears that the problem in the periodic system is much harder to treat than the case considered here.

RAB acknowledges financial support from the EPSRC under grant GR/R53197.

## Appendix A Two-point velocity function for surviving walks

In this appendix we outline the derivation of the correlation function $`\dot{X}(\tau _1)\dot{X}(\tau _2)_R`$ as given by equation (47). Recall that this is the two-point velocity correlation function for a single random walker that has a diffusion constant $`D_0`$ and does not cross the origin at any time $`0\tau ^{}\tau `$, where $`\tau >\tau _1,\tau _2`$. We shall assume, with no loss of generality, that the walker’s initial position $`Y>0`$ and that $`\tau _1<\tau _2`$.

The starting point is to consider the probability $`P_1(\tau ;X_1,\tau _1;X_2,\tau _2|Y)`$ for the walker to be at $`X_1`$ at time $`\tau _1`$, $`X_2`$ at time $`\tau _2`$ and to survive until time $`t`$ given that it started at $`Y`$. Then,

$`X(\tau _1)X(\tau _2)_R={\displaystyle _0^{\mathrm{}}}X_1{\displaystyle _0^{\mathrm{}}}X_2X_1X_2{\displaystyle \frac{P_1(\tau ;X_1,\tau _1;X_2,\tau _2|Y)}{Q_1(\tau |Y)}}`$
$`={\displaystyle _0^{\mathrm{}}}X_1{\displaystyle _0^{\mathrm{}}}X_2X_1X_2{\displaystyle \frac{Q_1(\tau \tau _2|X_2)G_1(X_2,\tau _2\tau _1|X_1)G_1(X_1,\tau _1|Y)}{Q_1(\tau |Y)}},`$ (61)
in which the Markovian property of the diffusion process has been used to write the joint probability distribution function $`P_1(\tau ;X_1,\tau _1;X_2,\tau _2|Y)`$ in terms of the functions $`G_1(X,\tau |Y)`$ and $`Q_1(\tau |Y)`$ defined through equations (9) and (11). Note that it was also necessary to normalise the distribution according to the probability that a walker starting at $`Y`$ survives a time $`\tau `$.

Assuming that $`\tau _0=Y^2\tau _1,\tau _2\tau `$, we can expand both the numerator and denominator in (A) to first order in $`Y/\sqrt{\tau _1}`$ and $`X_2/\sqrt{\tau }`$ to find

$`X(\tau _1)X(\tau _2)_R={\displaystyle \frac{1}{4\pi \tau _1\sqrt{\tau _1(\tau _2\tau _1)}}}{\displaystyle _0^{\mathrm{}}}X_1{\displaystyle _0^{\mathrm{}}}X_2X_1^2X_2^2\mathrm{exp}\left({\displaystyle \frac{X_1^2}{4\tau _1}}\right)`$
$`\times \left[\mathrm{exp}\left({\displaystyle \frac{(X_2X_1)^2}{4(\tau _2\tau _1)}}\right)\mathrm{exp}\left({\displaystyle \frac{(X_2+X_1)^2}{4(\tau _2\tau _1)}}\right)\right].`$ (62)
This integral can be evaluated without recourse to further approximations. First one integrates over $`X_2`$, which yields

$`X(\tau _1)X(\tau _2)_R={\displaystyle \frac{1}{2\sqrt{\pi }\tau _1^{3/2}}}{\displaystyle _0^{\mathrm{}}}X_1X_1^2\mathrm{exp}({\displaystyle \frac{X_1^2}{4\tau _1}})[(X_1^2+2(\tau _2\tau _1))`$
$`\times \mathrm{erfc}\left({\displaystyle \frac{X_1}{\sqrt{4(\tau _2\tau _1)}}}\right)+\sqrt{{\displaystyle \frac{4(\tau _2\tau _1)}{\pi }}}X_1\mathrm{exp}({\displaystyle \frac{X_1^2}{4(\tau _2\tau _1)}})].`$ (63)
To perform the remaining integration, it is helpful to make a change of variable $`X_1=2\sqrt{\tau _2\tau _1}u`$ which leads to

$$X(\tau _1)X(\tau _2)_R=\frac{16(\tau _2\tau _1)^{5/2}}{\sqrt{\pi }\tau _1^{3/2}}_0^{\mathrm{}}u\left[\frac{1}{\sqrt{\pi }}u^3\mathrm{e}^{\alpha u^2}+\frac{1}{2}u^2\mathrm{e}^{\beta u^2}\mathrm{erf}u+u^4\mathrm{e}^{\beta u^2}\right]$$
(64)
in which $`\alpha =\tau _2/\tau _1`$ and $`\beta =(\tau _2\tau _1)/\tau _1`$. Then, using

$`{\displaystyle _0^{\mathrm{}}}uu^3\mathrm{e}^{\alpha u^2}={\displaystyle \frac{1}{2\alpha ^2}}`$ (65)
$`{\displaystyle _0^{\mathrm{}}}uu^{2n}\mathrm{e}^{\beta u^2}\mathrm{erf}u=\left({\displaystyle \frac{}{\beta }}\right)^n{\displaystyle _0^{\mathrm{}}}u\mathrm{e}^{\beta u^2}\mathrm{erf}u=\left({\displaystyle \frac{}{\beta }}\right)^n{\displaystyle \frac{\mathrm{arctan}(1/\sqrt{\beta })}{\sqrt{\pi \beta }}}`$ (66)
we obtain, after a lot of manipulation,

$$X(\tau _1)X(\tau _2)_R=\frac{4\tau _1}{\pi }\left[\left(2+\frac{\tau _2}{\tau _1}\right)\mathrm{arctan}\sqrt{\frac{\tau _1}{\tau _2\tau _1}}+3\sqrt{\frac{\tau _2\tau _1}{\tau _1}}\right].$$
(67)
To find the velocity correlation function $`\dot{X}(\tau _1)\dot{X}(\tau _2)_R`$ we differentiate with respect to $`\tau _1`$ and $`\tau _2`$. We find

$$\frac{}{\tau _2}X(\tau _1)X(\tau _2)_R=\frac{4}{\pi }\{\begin{array}{cc}\mathrm{arctan}\sqrt{\frac{\tau _1}{\tau _2\tau _1}}+\frac{\sqrt{\tau _1(\tau _2\tau _1)}}{\tau _2}\hfill & \tau _1<\tau _2\hfill \\ 2\mathrm{arctan}\sqrt{\frac{\tau _2}{\tau _1\tau _2}}+2\sqrt{\frac{\tau _1\tau _2}{\tau _2}}\hfill & \tau _1>\tau _2\hfill \end{array}.$$
(68)
For $`\tau _1<\tau _2`$ one has

$$\dot{X}(\tau _1)\dot{X}(\tau _2)_R=\frac{4}{\pi \tau _2}\sqrt{\frac{\tau _2\tau _1}{\tau _1}}$$
(69)
and the same expression with $`\tau _1\tau _2`$ when $`\tau _1>\tau _2`$. At $`\tau _1=\tau _2`$, however, there is a jump in (68) of height $`2`$ which implies that

$$\dot{X}(\tau _1)\dot{X}(\tau _2)_R=2\left[\delta (\tau _2\tau _1)+\frac{2}{\pi \tau _+}\sqrt{\frac{\tau _+\tau _{}}{\tau _{}}}\right].$$
(70)
Thus we conclude our derivation of (47).


format: text

text:
# From Pulsar Spin to Pulsar Wind: The Spectral Connection?

## 1 Introduction

Recent observation of pulsars associated with supernova remnants obtained with the Chandra X-ray observatory (Weisskopf, O’Dell & van Speybroeck 1996) are providing an unprecedentedly detailed view of pulsar wind nebulae. For the first time, emission features involving wisps, co-aligned toroidal structures, and axial jets are fully resolved in X-rays on arc-second scales. These features, similar to those seen from the optical Crab nebula, are now found to be common to young, energetic pulsar in supernova remnants (Gotthelf 2001). Herein we present preliminary spectral analysis of several PWNe observed with Chandra, which, collectively, suggest a fundamental observational relationship between the spectral characteristics of pulsars and their pulsar wind nebulae.

Table 1 presents spectral results from a sample of Chandra pulsars which have characteristic wind nebulae. A summary of these objects along with references and images can be found in Gotthelf (2001). All observations were obtained with the ACIS-CCD camera which is sensitive to X-rays in the 0.2–10 keV band with an energy resolution of $`\mathrm{\Delta }E/E0.1`$ at $`1`$ keV. The on-axis point spread function is slightly undersampled by the CCD pixels ($`0.5^{\prime \prime }`$) allowing us to isolate the pulsar emission from that of the nebula. Except for N157B, which serendipitously fell on ACIS-I0, all data were obtained with ACIS-S3. The data were collected in nominal spectral (“FAINT”) and timing (3.24 s) mode and reduced and analyzed using the latest version of CIAO (CIAO 2.2/CALDB v.2.9). Starting with the Level 1 processed event files we corrected the event data for CTI effects (Townsley et al. 2001) and applied the standard Level 2 filtering criteria, then further rejected time intervals of anomalous background rates.

For each object listed in Table 1, we extracted spectra from the PWN, pulsar core, and background, when available, or obtained spectral parameters from the literature as indicated. For the core spectra, the brightest central pixels were extracted based upon data above 4.0 keV, where the core emission is substantially greater than that of the surrounding nebula. The nebula itself was then extracted from the region above which the background was constant, excluding the core. Large variations in background rates and size of both the nebula and the central source amongst observations precluded the use of a standard extraction aperture. Custom spectral response matrix functions (RMFs) were provided with the CTI correction software, and ancillary response functions (ARFs) were created according to standard CIAO 2.2 procedures, using the QEU calibration files similarly provided. All spectra were grouped to a minimum of 50 counts per spectral bin.

We fit the resultant background-subtracted pulsar and PWN spectra with a power-law model above $`2`$ keV using the latest version of XSPEC (v11.1). The spectral fits to the core included a convolution model to account for pileup effects on the spectra. Table 1 lists the spectral parameters obtained for each object; the values for the Crab and N157B have been obtained from the literature (Willingale et al. (2001), Wang & Gotthelf 1998, respectively). No measurement of the pulsed spectrum for 3C58 is currently available. Absorption values were obtained by fitting the nebula spectra with an absorbed power-law above 0.6 keV, and were subsequently frozen for all future fits.

The best-fit linear relationship between core and PWN photon indices is shown in Figure 1. We expect that the linear relationship seem in Fig. 1 may steepen slightly as we more fully account for pileup effects in the spectra. We note that the pileup corrected photon indices are likely subject to change due to the preliminary nature of the pileup model implementation, and suggest caution until these issues are resolved, but expect the basic result to remain unchanged.

We also note that due to extreme pileup effects and the brightness of its nebulae compared to the brightness of its core, SNR 0540-69 has not been included pending a more detailed analysis of the available data.

As a check to our final spectral results and to look for any systematic effect due to pile-up, we compare our measured power-law indices with those obtained from the literature for the pulsed emission for each source (Figure 2). Although they need not be the same, it is reassuring that they agree to within measurement errors.

## 2 Results

The photon indices obtained from this analysis show a trend that the average PWNe photon indices are consistently steeper than their core’s photon indices by a factor of one-half to one. We suggest that this relationship is linear, with the form $`\mathrm{\Gamma }_{PWN}=0.8\times \mathrm{\Gamma }_{Core}+0.8`$. This trend may be due to a number of factors. It may imply that the synchrotron spectrum frequency break occurs close to the neutron star, that there is a fundamental difference in emission mechanisms between the two regions, or that the nebula is being powered by the core and that the steepening is due to synchrotron burn-off or a comparable aging mechanism.

Though pulsations have not yet been discovered in the nebulae in supernova remnants G21.5-0.9 and G54.1+0.3, we note that their spectral indices are in the range predicted by this formula (Slane et al. 2000, Lu et al. 2002, respectively).

Chandra is the first telescope able to spatially resolve the neutron star from the nebula. Previous spectral observation of pulsars are likely to be strongly contaminated by PWN emission. For example, we see no clear relationship between the power-law spectral index obtained for these pulsars by ASCA with any of those measured by Chandra. As more crab-like pulsars are observed and analyzed, a more complete picture of their morphological and spectral characteristics comes into view, paving the way for future models and simulations.

###### Acknowledgements.

We gratefully acknowledge the support by the Heraeus foundation. This work made possible by NASA LTSA grant NAG 5-7935. C.M.O. also acknowledges the support of the I.I. Rabi Scholars Program at Columbia University.