## Abstract

The illumination of single population behaviour subject to the processes of birth, death and immigration has provided a basis for the discussion of the non-Gaussian statistical and temporal correlation properties of scattered radiation. As a first step towards the modelling of its spatial correlations, we consider the populations supported by an infinite chain of discrete sites, each subject to birth, death and immigration and coupled by migration between adjacent sites. To provide some motivation, and illustrate the techniques we will use, the migration process for a single particle on an infinite chain of sites is introduced and its diffusion dynamics derived. A certain continuum limit is identified and its properties studied via asymptotic analysis. This forms the basis of the multi-particle model of a coupled population subject to single site birth, death and immigration processes, in addition to inter-site migration. A discrete rate equation is formulated and its generating function dynamics derived. This facilitates derivation of the equations of motion for the first- and second-order cumulants, thus generalizing the earlier results of Bailey through the incorporation of immigration at each site. We present a novel matrix formalism operating in the time domain that enables solution of these equations yielding the mean occupancy and inter-site variances in the closed form. The results for the first two moments at a single time are used to derive expressions for the asymptotic time-delayed correlation functions, which relates to Glauber’s analysis of an Ising model. The paper concludes with an analysis of the continuum limit of the birth–death–immigration–migration process in terms of a path integral formalism. The continuum rate equation and evolution equation for the generating function are developed, from which the evolution equation of the mean occupancy is derived, in this limit. Its solution is provided in closed form.

## 1. Introduction

The *K*-distribution model of the non-Gaussian single point statistics of optical and microwave radiation scattered from extended objects, such as the sea surface, was introduced some 30 or so years ago (Jakeman & Pusey 1976). This sought to capture the effects of significant large-scale correlated structure within the scatterer that vitiate any appeal to the central limit theorem as a justification for a more conventional Gaussian statistical model. Since then the model has grown to maturity, and the *K*-distribution, whose probability density function (PDF) is expressed in terms of a modified Bessel function of the second kind (Abramowitz & Stegun 1965), and is analytically tractable, adaptable and robust, has found widespread application (Jakeman & Ridley 2006; Ward et al. 2006). Capturing the two-point or correlation properties of the scattered radiation has always been more problematic and invariably based on simplified models (Tough & Ward 2002) which, in the main, exploit the representation of the intensity of the scattered radiation as a Rayleigh process (generated by scattering by the small scale structure on the surface) with a ‘local power’ that itself varies as a gamma process and so models the modulation of the local scattering by large scale structures (Ward 1981; Ward et al. 2006). Jakeman’s (1980) analysis of the dynamics of the *K* process in terms of scattering from a population subject to negative binomial fluctuations (established by competing processes of birth, death and immigration at a single site) provided an early and fruitful example of this approach. Subsequently, Jakeman’s model was reformulated in terms of stochastic differential equations (SDEs) and their associated Fokker Planck equation (FPE; Tough 1987); this SDE formulation has been re-examined and extended considerably in more recent work by Field & Tough (2003*a*,*b*, 2005) and Field (2008). Thus far, only the dynamics of the radiation scattered by a single fluctuating population have been considered in this way; in this paper, we take the first steps in an extension of this analysis to include both spatial and temporal structure in the scatterer and the correlation function of the scattered radiation.

The dynamics of a single population subject to the processes of birth (B), death (D) and immigration (I) lie at the heart of Jakeman’s model and its SDE analogues; not surprisingly their analysis has many parallels in work carried out in biometric and related studies (Renshaw 1991). Here we intend to extend this to an infinite chain of sites, each of whose populations is subject to BDI and migration (M) between nearest neighbours. Once again significant progress in this area has been made in the biometric context which provides us with a useful starting point. Thus Bailey (1968) has considered a stepping stone model, which incorporates B, D and M; in this, populations either decay to zero or grow without limit as time progresses, depending on whether the birth rate is greater or lesser than the death rate. As the multi-site problem is less tractable than its single site analogue, Bailey was only able to evaluate the mean and variance of the population at each site and the inter-site population covariance. This model has been examined in more detail by Renshaw (1974), who also discusses the effect of immigration to a single site. (We note that Glauber’s (1963) classic work on the dynamics of the Ising model has much in common with these analyses.) Here, we wish to consider the rather different case where immigration occurs, at the same rate, at each site; this serves to establish a non-zero mean value for each population. Previously, Bailey (1968) and Renshaw (1974) used techniques based on generating functions and Laplace transforms; we attack the coupled differential equations satisfied by the population cumulants directly and are able to extract more compact and general solutions, as well as incorporating the effects of immigration at every site.

Once we have investigated some of the properties of the solutions obtained, we turn our attention to the limit in which the chain of sites becomes a continuum. In the discrete site case, the generating function provides us with a useful tool; in the continuous limit this is replaced by a generating functional defined as an average over all population profiles, which takes the mathematical form of a path integral. We discuss this in some detail, and show how the equations of motion of the number density and inter-site distribution function can be calculated. Particular attention is paid to the description of the diffusive motion, resulting from the inter-site migration process, within this framework.

## 2. A simple migration model for diffusion

Before we consider the problem of the populations of a set of sites, subject to the processes of birth, death and immigration and coupled by inter-site migration, we will discuss a simpler multi-site problem, of a single particle migrating about an infinite chain with steps made with equal probability to the two neighbours of its current site. This allows us to introduce some of the concepts we will use later, albeit in a simpler context, and should motivate parts of our subsequent analysis. As such, the material in this section is intended to be primarily didactic in nature.

We label the sites by an integer *n*, taking values between plus and minus infinity; the probability that the *n*-th site is occupied at time *t* is denoted by *p*_{n}(*t*). We assume that, at time zero, the particle is sited at the origin, *n*=0. Thus we have
2.1
At any time, the *p*_{n}(*t*) satisfy the rate equation, which encodes the conservation of probability under the inter-site hopping process
2.2
where the parameter *η* characterizes the inter-site transition rate. More precisely, the transition rate in the positive or negative directions is equal to and so the rate of migration in either direction is equal to *η*.

This rate equation is intimately related to the diffusion equation. If we identify the label *n* with a spatial coordinate *x*, divided by the inter-site separation Δ*x*, and develop the right-hand side of equation (2.2) in a Taylor series, we find that, adopting a fairly obvious notation,
2.3
Thus we see that the diffusion equation emerges in the continuum Δ*x*→0 limit, provided that the transition rate tends to infinity in such a way that
2.4
which quantity we identify as the diffusion constant. We remark that this divergence in the transition rate is akin to the property of a diffusion process with an infinite level crossing rate (LCR). Consistently, the above analysis constitutes a construction of the Brownian motion, or continuous valued Wiener process, from a discrete random walk.

The coupled differential equations (2.2) can be solved explicitly, perhaps most straightforwardly by introducing a generating function for the site probabilities
2.5
From equation (2.2) we see that this satisfies
2.6
which we can solve to give
2.7
The probabilities of site population can then be recovered by the Fourier inversion
2.8
Here we have recognized the modified Bessel function *I*_{n} in the form of a well-known integral representation (Abramowitz & Stegun 1965).

We have argued that equation (2.2) reduces to a familiar diffusive form in the continuum limit. It is natural therefore to examine to what extent this is manifest in its solution (2.8). Previously, we noted that, in going to the diffusive continuum limit, it was necessary to allow the transition rate to diverge. This suggests that we can introduce the large argument asymptotic representation of the modified Bessel function (Abramowitz & Stegun 1965) 2.9

As before, we can make the identification 2.10 so that 2.11 When we consider that, in going to the continuum limit, we must impose the relationship (2.4) on the transition rate and the inter-site separation, we see that 2.12 so that equation (2.11) becomes 2.13 Reassuringly, a familiar diffusive result has again emerged, albeit after a more sophisticated analysis.

The approach to the continuum limit can also be discussed in the framework provided by the generating function *c*(*q*,*t*), which then tends to the characteristic function of the PDF of the particle’s position.
2.14
If we invoke equation (2.10) and impose the condition (2.4), we find that the equation of motion (2.6) reduces to
2.15
which, once more, reveals that the PDF *P* satisfies the simple diffusion equation (2.6).

In this section we have discussed a particularly simple multi-site transition problem and highlighted the role in its solution played by the generating function technique and its introduction of the modified Bessel functions. The importance of the constraint (2.4) placed on the vanishing inter-site separation and the diverging transition rate as the continuum limit is approached has also been emphasized. All these features will be encountered in the following sections when we analyse the more complicated problem of multiply (i.e. a given site may be occupied by many individuals at one time) occupied sites, with populations subject to birth, death and immigration, in addition to inter-site migration considered here.

## 3. The analysis of discrete coupled populations subject to birth, death, migration and immigration

The analysis of the temporal behaviour of a single negative-binomially distributed population, established by the competing processes of birth, death and immigration, played a central role in the development of the SDE description of a gamma distributed process. Consequently, as we address the problem of a process displaying spatial and temporal variation, we consider an array of sites, whose populations are, in addition to being subject to birth, death and immigration, coupled by inter-site migration. Bailey (1968) has discussed the simpler problem in which immigration is suppressed, deriving expressions for the temporal evolution of the means and variances of the site occupancies. Renshaw (1974) gives an exhaustive account of this work, and of his subsequent contributions. Here, we will consider a novel mathematical approach to the problem, where immigration takes place at every site. (Of course, assuming that individuals behave independently and identically, the probability generating function (PGF) for the multi-site immigration problem factorizes into the component PGFs for the stepping stone problem considered over an infinite set of sites with a single immigration source at just one site. On the other hand, migration induces spatial correlation and therefore the joint PGF, with respect to different *sites*, does *not* factorize.)

The methods used by Bailey and Renshaw exploit generating function and Laplace transform techniques to solve the coupled differential equations satisfied by these mean values, and lead to various implicit solutions in these terms. We develop an alternative approach to the problem, which operates in the temporal domain throughout, that is sufficiently economical to provide explicit compact solutions of the zero-immigration (or ‘stepping-stone’) case and to be extendible to the more general problem in terms of a novel matrix formalism.

The present approach contains much that complements the intimately related work of Bailey and Renshaw, with the benefit that the methods employed operate entirely in the time domain and are more explicit and compact. Renshaw’s distinguished work in this area has paved the way for us to achieve our results, which are of interest in their own right. Thus, our contribution in this section is to be regarded as a powerful alternative new mathematical approach, which appears to be more tractable.

### (a) Rate equation, generating function and cumulants

The probabilities *P*(**N**,*t*) that the populations **N**≡{*N*_{j}} are established on the sites *j*, , at time *t* obey the coupled rate equations
3.1
which generalize those of Bailey (1968) through the incorporation of terms representing immigration at each site. The derivation of these equations and the issue of ‘double counting’ that arises are dealt with separately in the electronic supplementary material, appendix A. Here **N**[*k*,Δ] represents a population differing from **N** at one site only, i.e. ⋯*N*_{k−2},*N*_{k−1},*N*_{k}+Δ,*N*_{k+1},*N*_{k+2},⋯ . In the case where populations at two sites differ from those specified by **N** we write, for example, **N**[ *j*−1,1;*j*,−1]≡⋯*N*_{j−2},*N*_{j−1}+1,*N*_{j}−1,*N*_{j+1},*N*_{j+2}⋯ . The parameters *λ*,*μ* characterize the processes of birth and death, whose rates are taken to be proportional to the current site population, while *ν* characterizes the immigration, which is taken to occur at a constant rate at each site, independent of its population. The rate of inter-site migration (as a net result of migration in either direction) is taken to be proportional to the population of the exited site and characterized by the parameter *η*. The state at each site is represented by an infinite collection of probabilities {*P*(**N**,*t*)} that can be reduced down to a single generating function that encodes the probabilities of all possible states, and in terms of which the appearance of the underlying rate equation can be simplified. To this end, these probabilities can be combined to form a moment-generating function
3.2
in which 〈·〉 denotes expectation. The time derivative of this generating function includes single site and adjacent site contributions. The birth, death and immigration contributions take the form
3.3
The first of the adjacent site contributions derived from equation (3.1), which represents a sum over migrations from site *j* to site *j*−1, takes the form
3.4
which can be re-written as
3.5
The second adjacent site term in equation (3.1) can be handled similarly; when all the contributions are brought together, we arrive at the multi-site Kolmogorov equation
3.6
This can now be used to derive the equation of motion of the mean population of any site as

We note that, when evaluating this quantity, the ‘second moment’ ∂^{2}*G*(**x**,*t*)/∂*x*_{j}∂*x*_{k} is encountered by virtue of equation (3.6). In the limit **x**=**0**, however, the coefficient of this quantity vanishes, so that the equations of motion of the mean site populations take the form of a set of coupled linear differential equations that do not depend on any higher moments of the site populations. This absence of a ‘hierarchy’ of moment equations renders the problem tractable and is a consequence of the simple dynamics encompassed in equation (3.1).

The evaluation of the site population variances is facilitated by the introduction of the cumulant-generating function
3.7
wherein we have identified the means and variances
3.8
The function *K* in turn satisfies the temporal evolution equation
3.9
from which we can deduce the equations of motion of the first- and second-order cumulants (3.8) as
3.10
and
3.11
These differ from the corresponding results obtained by Bailey (1968) only through the appearance of the immigration rate *ν* in the equations of motion of *m*_{j} and *V* _{jj}.

The equations of motion for the mean occupancies and the variances (3.10) and (3.11) are linear, as are those for the immigration at a single site (cf. §3*e*). Consequently the results obtained for the latter could be superposed to generate solutions of the former. In comparison, the expressions obtained by Renshaw (1974) for the single site immigration do not lend themselves to this procedure. The key point is that our methodology, which eschews the generating function approach, puts forth solutions to all problems of interest that are more compact.

### (b) Properties of the infinite dimensional matrix **M**

If we introduce the infinite dimensional matrix **M** with elements
3.12
we can express the set of equations (3.11) in a compact and transparent form, which is amenable to direct solution. Before we do this, however, we will highlight some of the properties of **M** which we will make use of in our subsequent analysis. This is addressed in detail in the electronic supplementary material, appendix C.

### (c) Time variation of the mean population of an individual site

With these few preliminaries taken care of we turn to the solution of the equations of motion (3.10) and (3.11): the first we can write as
3.13
The vector **n** has unit elements, i.e. *n*_{k}=1,∀*k*. We now introduce an integrating factor
3.14
The identity (see the electronic supplementary material, appendix C, equation (C8)) allows us to evaluate the right-hand side of equation (3.14) in the form
3.15
Straightforward integration leads us to
3.16
This can be written out more explicitly as
3.17
where the modified Bessel functions are identified through equation (C7) in the electronic supplementary material, appendix C. This derivation is altogether more compact than the generating function-based methods similar to those discussed in §2 and applied to this problem by Bailey (1968) and Renshaw (1974). These generating functions are quite distinct from equations (3.2) and (3.7), introduced to facilitate the derivation of the equations of motion (3.10) and (3.11).

### (d) Time variation of the inter-site variance of populations

The equations of motion for the variances (3.11) can be written in the form
3.18
where **V** is the matrix of covariances composed of the elements *V*_{ij}(*t*). The matrix **A** takes the form
3.19
where
3.20

The elements of the matrix are constructed, using the result (3.17), as
3.21
We note that does not commute with **M** and consequently the solution to equation (3.18) must be written as
3.22
Results equivalent to those of Renshaw (1974) can be obtained from this relatively simply. The first term in equation (3.22) separates out the effects of immigration at a constant rate to each site; it vanishes when we set *ν* to zero. The second term captures the transient effects derived from the initial population of the sites and provides a compact representation of the solution of the stepping stone problem discussed by Bailey (1968) and Renshaw (1974). In the absence of this second, ‘symmetry breaking’, term the covariance matrix exhibits a translationally invariant *Toeplitz structure*
3.23
We will now ‘unpack’ the formal solution (3.22). By exploiting the translational invariance of the first term we can recast it in the form
3.24
where
3.25
and (cf. equation (3.11))
When we identify elements of the exponentiated matrix with the modified Bessel functions, as in equation (C7) in the electronic supplementary material, appendix C, we can write equation (3.24) explicitly as
3.26
This result can be simplified by making use of the Bessel function identity (Abramowitz & Stegun 1965)
3.27
which leads us, when we integrate by parts, to
3.28
In the infinite time limit, we find that
3.29
which is seen to be identical with the result
3.30
when use is made of equation (C5) in the electronic supplementary material, appendix C. It is particularly interesting to note that, at the level of this second cumulant, the populations at different sites are uncorrelated, even in the presence of inter-site migration, if the birth process is absent (i.e. *λ*=0).

We now consider the second term in equation (3.22), which we will denote by , in more detail. Introducing the expressions (3.21) for the elements of , we find that the elements of are given by a sum of three contributions 3.31 The second and third of these can be combined to give 3.32

Using equation (3.27) we can recast this into a form
3.33
which facilitates integration by parts. When this is carried through, we find that the transient stepping stone contribution to the variance matrix can be written as
3.34
It is interesting to compare the results we have obtained here with those presented in the literature. Bailey (1968) analysed the stepping stone model, in which the immigration process is absent, using a generating function formalism. His solution for the mean site occupancies is identical to ours; his analysis of the variances encounters some problems. In particular he is able to make progress towards an explicit solution only in the case where the birth rate *λ* is set to zero; he then finds that the variance of the population of the *i*-th site at time *t*, when the site at the origin alone is occupied at time zero, is given by
3.35
Our solution (3.34) reduces to this special case when we set *λ*=0,*N*_{p}(0)=*δ*_{p,0}. The second, integral, contribution to equation (3.33) vanishes when *λ*=0 and presumably captures those terms that Bailey is unable to accommodate in his analysis when *λ*≠0. It is interesting to note that Bailey claims that in the *λ*≠0 case ‘to find *V* _{i,j}… appears to be a matter of some complexity, and involves evaluating the appropriate complex integral for the relevant term in a two-dimensional Laurent expansion’. Appendix B (electronic supplementary material) discusses how the integral contribution in equation (3.34) can be cast in a form that makes contact with these remarks.

Renshaw (1974) has also considered the stepping stone model, again adopting a generating-function approach. His analysis is carried through for *λ*≠0. The results obtained are rather cumbersome and do not make direct contact with Bailey’s simple result (3.35).

### (e) Immigration at a single site

In addition, Renshaw (1974) has discussed a model incorporating birth, death, inter-site migration and immigration to just a single site, again making use of generating functions. The methods we have developed here for the analysis of the case where immigration takes place at each site can be applied to this problem; we will now outline briefly how this can be done (cf. remarks concerning component PGFs in the opening paragraph to this section). This is addressed in detail in the electronic supplementary material, appendix D.

### (f) Inter-site temporal correlation functions

In the previous section, we analysed the variances of the site populations, evaluated at some time *t*, and were able to separate out their transient parts from those that remain in the long time limit. Here we consider the long time limit of the correlation function between site populations with a time delay, adopting a method similar to that used by Glauber (1963) in his analysis of the dynamics of an Ising model. (Glauber’s model corresponds more closely to the stepping stone model, and does not include the terms that we would ascribe here to immigration.)

Thus far we have considered the mean and variance of the site occupancies of the form
3.36
We now study the time-delayed correlation functions of the form
3.37
in the limit where *τ* becomes very large. Following Glauber, we can put these together ‘by hand’ from the preceding results. Thus we consider
3.38
Here the angular bracket signifies an average over the stochastic evolution of *N*_{j}(*t*), conditional on a given initial configuration {*N*_{i}(0)} as we have discussed previously. The square bracket signifies an average over the initial configuration {*N*_{i}(0)}, using the asymptotic equilibrium results we have just derived. Thus, Glauber’s procedure to calculate the asymptotic correlation function (3.37) is to calculate the same for the finite time *τ*=0 and then average over initial equilibrium configurations that conform to the asymptotic distribution. According to this idea, the representation (3.38) is a consequence of the ‘tower law’ or iterated expectation rule.

Following this procedure, 3.39

From equation (3.17) we see that 3.40

Carrying out the square bracket averaging using equation (3.40) then leads us to 3.41

This gives us the asymptotic time-delayed two-site correlation function.

In figure 1, we show the behaviour of *r*_{i,j}(*t*), with the values *μ*=2.0, *λ*=1.0, *ν*=1.5, *η*=1.0 and |*i*−*j*|=0,1,2.

We note that the correlation function of the populations at different sites can increase in value at short times. This behaviour becomes more pronounced as the inter-site migration parameter increases and is illustrated in figure 2; we have taken |*i*−*j*|=5, *μ*=2.0, *λ*=1.0, *ν*=1.5, and *η*=1,3,5 and 10.

We remark, in a more applied setting, that such results could form the basis of comparison in experimental and simulation studies of spatially correlated radar cross-section (RCS) data, for example.

## 4. A path integral formulation of the continuum limit

Thus far we have considered the dynamics of an infinite set of populations, supported on a chain of discrete sites. We will now consider how our formalism might be modified to accommodate the merging of these sites into a continuous line. To do this, we generalize the generating function encountered already to a generating functional, whose construction we then discuss in terms of path or functional integration. At the outset, we shall focus on the pure migration (diffusive) component of the process, and incorporate terms due to birth, death and immigration subsequently. For pure migration, the rate equation simplifies to 4.1

As we have seen previously, the corresponding moment-generating function obeys the evolution equation
4.2
To address the continuum limit of this process, i.e. for which the inter-site separation distance tends to zero, we introduce the number *density* (i.e. per unit length) *n*(*z*), analogous to *N*_{j}, and the test function *l*(*z*), corresponding to the site labelling vector **x**. Corresponding to equation (3.2), the generating functional takes the form
4.3
where the angular brackets represent an average over an ensemble of density profiles. Thus the continuous variable *z* labels space, and in going from the discrete to continuous space models we have made the following identifications:
4.4
As we shall discuss in detail, this functional can be considered as a *path integral* somewhat akin to that encountered in quantum field theory (QFT).

The mean number density can be derived from the generating function by functional differentiation as follows: 4.5

Here, we have identified a delta function through 4.6 This enables us to extract the mean number density through 4.7 Higher order moments of the number density can be generated in much the same way, by repeated functional differentiation.

In order to determine the dynamics of the mean number density, we now require the equation of motion of the generating function. To begin with we retain a positive inter-site separation Δ*z*, which we shall then let tend to zero, while the transition (migration) rate tends to infinity in such a way as to yield a finite diffusion constant. In this respect, the requirements of the limiting behaviour for the continuum are consistent with the discussion in §2 surrounding equation (2.4). Thus, as we shall see, the limiting behaviour that we require is
4.8
We observe, in this limit, the vanishing properties
4.9
which will be required in the calculations that follow.

In accordance with these remarks, we can generalize equation (7.2) to the continuous case as
4.10
where, in the last equation, we have used the Taylor series identity .^{1}

Now the hyperbolic differential operator in equation (4.10) can be expanded as 4.11 and in the continuum limit, according to equation (4.9), only the first term of the right hand side survives. Using the second derivative identity 4.12 the result of the spatially continuous limit of equation (4.10) is 4.13 To determine the time evolution of the mean number density, first observe, via equation (4.7) and interchange of functional and partial derivatives, that 4.14 Using this we find, most agreeably, that the equation of motion of the mean number density is 4.15

Here we have used the Leibnitz property for the functional derivative of products. As we might have expected, we obtain the diffusion equation here.

It is illuminating, and perhaps indeed more fundamental, to achieve this same result via the continuum limit of the rate equation itself, as opposed to that of the evolution equation for the generating function. Before we put this together, we make a few remarks. The analogue of *N*_{j} in the continuous case is the number density *n*(*z*); the number of point particles found in the range between *z*−d*z*/2 and *z*+d*z*/2 is *n*(*z*) d*z*. In the discrete case, we encountered unit changes in the necessarily integer site populations; the continuous analogues of these, as we shall see, are delta function valued changes in the number density *n*(*z*).

To achieve this aim, we need first to introduce some appropriate probability measure *P* on the space of continuous population profiles, which we shall denote as *P*[*n*(*z*)]. Such a measure should satisfy the continuum form of the discrete-rate equation (4.1):
4.16

In the above rate equation, the derivatives of the delta functions represent the changes in the population density due to migration, as the inter-site separation Δ*z*→0. The precise way in which the delta functions arise rests on the observation that, for a migration towards the origin from the positive direction, the resulting change in density can be considered as −*δ*(*z*−Δ*z*)+*δ*(*z*)∼*δ*′(*z*−Δ*z*/2)Δ*z* which, as must be the case for a density occurring in the argument of the measure *P*, has dimensions of the reciprocal distance.^{2}

Now the expression for the generating function can be expressed in terms of the measure *P* as a path integral as in the following:
4.17
which is the continuum analogue of equation (3.2). This path integral expression enables one to calculate the evolution equation for the generating function directly from the above continuum rate equation, as follows:
4.18

It is readily observed that the final term in the above equation contributes . For the remaining terms, first observe that occurs over all profiles, so that we can transfer the derivative of the delta function occurring in the probability measure to the exponential, just like in the discrete case. Thus, for instance, in the first such term, we introduce the change of the density profile variable
4.19
so that *n*(*z*_{0})+Δ*n*=*n*⌣(*z*_{0}). This change of variable obviates any difficulties that could have arisen in assigning a numerical value to the quantity Δ*n*, which is distributional-valued in nature. As a result of this procedure, we obtain, for the first such term, taking careful note of the change of variable,
4.20
In this expression, the integral involving the delta function equates to −*l*′(*z*_{0}−Δ*z*/2)Δ*z*. On taking the path integral, we are thus left with
4.21
Similarly, for the second such term, we make the change of variable
4.22
so that, like before, . This yields the contribution
4.23
Thus, altogether we find, direct from the continuum form of the rate equation itself, the evolution equation for the generating function is
4.24

Writing the Taylor series expansions of the above test functions in the familiar exponential form
4.25
we find, for their sum occurring in the integrand,
4.26
Expressed in terms of hyperbolic functions this is equal to
4.27
Taking all the derivatives explicitly, we determine the following evolution equation for the generating function, for Δ*z*>0:
4.28

Observe that all terms in the integrand are regular in the continuum limit represented by equation (4.8), and moreover that only the first two terms in brackets […] survive, which equate to ((*l*′)^{2}+*l*′′)Δ*z*^{2}, whose expression was encountered earlier in equation (4.13). Thus, the resulting dynamics of the generating function, and therefore also the dynamics for the mean occupancy, are identical to those obtained via the previous method wherein the continuum limit of the discrete form of the evolution equation for the generating function was tackled directly.

It is satisfying and indeed intriguing that, despite the obvious differences between the evolution equations (4.10) and (4.24) (which are considered for positive inter-site separation Δ*z*), in the continuum limit Δ*z*→0 both approaches yield the same result, as they must for consistency and the limit to be well defined.

This analysis generalizes naturally to the case where the single site processes of birth, death and immigration also occur. Translating equation (3.6) to the continuum, we obtain the following equation of motion for the generating function: 4.29

Now applying *δ*/*δl*(*z*_{0})|_{l=0} to the above equation of motion and using equation (4.14), the extra contributions to ∂〈*n*(*z*_{0},*t*)〉/∂*t* involving the birth, death and immigration rates are as follows:
4.30
Since *G*(*l*,*t*)|_{l=0}=1 and *δG*(*l*,*t*)/*δl*(*z*)|_{l=0}=〈*n*(*z*)〉, upon integration against the delta function, we deduce that the mean number density for the continuum limit of the BDIM model satisfies the equation of motion
4.31
In this equation, for an equilibrium solution to exist (cf. the asymptotic solution for the single site BDI population model), we require *μ*−*λ*>0 and thus (below) set this equal to a quantity *κ*^{2}.

This result makes a direct contact with the situation for the discrete case as represented by equation (3.13). Moreover, an exact solution of the above equation can be obtained, as follows. The PDE 4.32 (which can be interpreted in terms of a diffusion supplemented by a uniform source, and depleted by some local decay process) can be solved in free space (i.e. without boundary conditions) by the Fourier (space) and the Laplace (time) transformation.

Thus we find that 4.33 in which the tilde and hat represent the Fourier and the Laplace transforms with respect to space and time according to the scheme 4.34 This can be re-arranged to give 4.35 Fourier–Laplace inversion then yields 4.36 This continuum solution for the mean occupancy is the analogue to the solution of the discrete case, and emerges as the continuum limit of equation (3.17), in a similar fashion to the simple diffusion case studied in §2.

Observe, as is well known, that such a diffusive model exhibits acausal propagation. We refer the interested reader to Field (2008, ch. 1), where related ideas are discussed in a historical context.

It is interesting to compare our application of path integral methods with those exploited by Dashen (1979) in his analysis of multiple scattering of waves in random media. The methodologies of the two approaches have relatively little in common and are, in some sense, complementary. Thus Dashen’s path integral analysis explicitly captures propagation through a random medium, describing it with the parabolic wave approximation to a Helmholtz equation that incorporates a fluctuating refractive index. As the parabolic wave equation is very similar in structure to the two-dimensional non-relativistic Schrödinger equation (distance along the axis of propagation in the former plays a role that is formally identical to time in the latter) Dashen is able to co-opt much of the formal machinery developed by Feynman & Hibbs (1965) and apply this to the calculation of statistical properties of the propagated field. Furthermore Dashen ascribes Gaussian statistics to the refractive index fluctuations, or assumes they are very small. Our analysis focuses on the properties of the scattering cross section itself, and makes no explicit reference to the propagation process. In particular, our path integral is not taken over propagation paths, as in Dashen’s work, but over an ensemble of density profiles; the underlying equation of motion does not contain an imaginary term and so is more diffusive than Schrödinger-like. It is expected that results of further applications of our formalism will characterize the statistics of the scattering cross section of a continuous medium; a full description of the scattering process might then be achieved by fusing elements of our analysis and that of Dashen in a path integral formulation that captures all aspects of the problem.

## 5. Discussion and conclusions

Initially the work we have presented here was motivated by its potential usefulness in the modelling of the spatial dependence of the non-Gaussian statistics of scattered light. Populations subject to birth, death and inter-site migration had already been discussed in the biometric literature; to establish a uniform lattice of stable non-zero populations, however, immigration at each site must be incorporated into the model. As these non-zero populations are to be identified with local scattering powers in the generalization of Jakeman’s (1980) earlier single site model, the analysis of this more general BDMI model is a necessary precursor to any scattering calculation. The methods used previously to characterize the dynamics of the stepping stone and BDM become rather intractable when they are modified to take account of immigration at each site. We have proposed an alternative approach, which eschews the generating function formalism adopted by Bailey (1968) and Renshaw (1974). This is able to capture the effect of immigration, and has produced a detailed solution of the problem, in terms of a novel matrix formalism, that operates entirely in the time domain and yields explicit closed form solutions.

The approximation of an underlying discrete lattice of populations of scatterers can be criticized on the grounds that it is not realistic; as Bailey (1968) has found, its generalization to more than one spatial dimension is also quite intractable. It is hoped that these difficulties might be alleviated by a spatially continuous description of the problem. We have developed such a formulation and applied it in a preliminary investigation of the dynamics of the population density. The simple results we have obtained make direct contact with those obtained in the discrete site case. The generalization of our formalism to more than one spatial dimension also appears to carry through straightforwardly. However it should be stressed that, in this model, it is merely the spatial support that is continuous; we have not considered the limit that allows us to proceed to the stochastic differential equation description of scattering by a single population with a very large mean occupancy (Tough 1987; Field & Tough 2003*a*,*b*). So, while some progress has been made towards our goal of capturing the spatial correlation characteristics of non-Gaussian scattered light, we expect that a significant quantity of work remains to be done in this area. It seems reasonable to assert that, in a more applied setting, this approach has significant potential for generating appropriate spatially correlated RCS patterns in simulation models for radar backscatter. The focus of this work is largely theoretical, and accordingly simulation studies will be discussed in a separate paper.

The path integral formulation has proved to be very illuminating and challenging. It may be worthwhile, in light of these results, to consider whether the same types of techniques could be applied successfully to QFT.

In the present context, the existence and finiteness of an underlying discrete interacting population serves to place the emergent path integral formalism on a much more rigorous footing. It should be emphasized here that the continuum description is indeed an *emergent* and apparent one, which is generated from the discrete structure.

Identification of such a rigorous and finite underlying discrete structure could well provide some insights in QFT in connection with some of the issues there to do with infinities, renormalization and the precise meaning of the path integral formalism. Such issues will be contemplated elsewhere in the literature.

## Footnotes

- Received February 1, 2010.
- Accepted March 25, 2010.

- © 2010 The Royal Society