## Abstract

We analyse a number of stochastic processes that give rise to first-order number statistics governed by Laguerre distributions with properties that lie between Poisson and geometric random variables. These distributions have hitherto been used to characterize the photon statistics of a coherent mixture of thermal and laser light. Here, we explore a number of discrete population processes that can lead to the same first-order statistics and correlation functions and highlight the distinguishing features that may be used to identify the relevant model from data.

## 1. Introduction

In the semi-classical theory of photoelectric detection, the time-varying electromagnetic field falling on a photo-detector (for example, photomultiplier or photo-diode) is assumed to be a classical wave field. The detector responds to this field by emitting photons in proportion to the square of its envelope or *intensity* *I*(*t*). The random train of events generated in this way is a doubly stochastic Poisson process governed by the so-called Mandel formula (Mandel 1959) for the probability of finding *n* events in time *T*:
1.1
(see Jakeman & Pike (1968) and references therein). In equation (1.1), *I* is the instantaneous intensity *i*(*t*) averaged over a response or integration time *T*
1.2
*η* is the quantum efficiency of the detector and *P*(*I*) is the probability density of *I*. Note that if *T* is much shorter than the fluctuation time of the intensity, then the intensity is almost constant over the integration time and (1.2) reduces approximately to *I* (*t*;*T*)=*i*(*t*)*T*. The relationship (1.1) implies that the *factorial* moments of the photo-count number distribution are proportional to the *ordinary* moments of the integrated intensity fluctuation distribution
1.3
Indeed, the Laplace transform of the latter quantity is identical with a generating function for the photon-counting fluctuations, apart from a scaling factor arising from the detector efficiency:
1.4
Relation (1.1) implies that the *constant* intensity or brightness *I*_{0} characterizing light from an ideal amplitude-stabilized laser produces a random Poisson train of events, since if we assume *P*(*I*)=*δ*(*I*−*I*_{0}) it predicts that
1.5
Here, the mean photo-count number is given by *n*_{0}=*ηI*_{0}=*ηi*_{0}*T*.

Another well-known result is obtained in the case of the so-called Gaussian light with intensity fluctuations governed by a negative exponential distribution. When the integration time is sufficiently short, *P*(*I*) is also approximately a negative exponential leading from (1.1) to a *geometric* or *thermal* distribution of photo-counts:
1.6
In this formula, the mean number is given by where and are the mean intensity and the average-integrated intensity, respectively. The generating functions corresponding to (1.5) and (1.6) take particularly simple forms:
1.7a
and
1.7b
In many, early photon-counting experiments *coherent mixtures* of laser light and Gaussian light were encountered, either by design or by accident. This configuration corresponds to the classical problem of a constant amplitude signal embedded in Gaussian noise and is generally termed *heterodyning*. A number of results for the related photo-electron statistics have been obtained for the special case in which the carrier frequency of the signal and noise is the same, usually referred to as *homodyning*. In §2, this problem will be briefly reviewed, but suffice it to say here that the resulting photon statistics for this case can be derived from the generating function (see Jakeman & Pike 1969).
1.8
Here *n*_{0} and are the average number of laser and ‘thermal’ counts that would be measured if photons from the different kinds of light could be distinguished, again assuming that the integration time *T* is sufficiently short to resolve the fluctuations. The result (1.8) is closely related to the generating function for the Laguerre polynomials and it can easily be shown that the corresponding distributions are given by
1.9
where the L_{n}(*x*) are Laguerre polynomials of order *n* (Gradshteyn & Ryzhik 1994).

A gauge for the size of the fluctuations is the Fano factor, and for Poisson fluctuations, this is unity, whereas for thermal fluctuations . For the Laguerre fluctuations, with and these distributions interpolate between the thermal and Poisson results (1.7). The Laguerre distributions provide a potentially useful model for the analysis of series of events in scenarios other than photon-counting. The results of Gaussian noise theory can also be invoked to derive fluctuation properties so that the *correlation* of events can be calculated as well as their single interval number statistics.

We have demonstrated in previous papers (Matthews *et al*. 2003; French *et al*. 2008) that the statistical properties of the time series of events generated when individuals leave an evolving population reflect the properties of the population itself (apart from a scaling of the mean) provided that any fluctuations in the population take place on a time scale much longer than the integration time *T*. This suggests that the statistical properties attributed to the train of photon detection events described above may also be generated by an appropriate discrete population model. For example, the equilibrium distribution of the number of individuals in a population where there are only deaths and immigrations is Poisson (1.5), while the thermal distribution (1.6) is found to characterize the number of individuals in a population governed by a birth–death–immigration process in which the birth and immigration rates are equal.

In this paper, we examine a number of discrete population models that generate homodyne first-order statistics of the type (1.9). Just as in the case of a doubly stochastic process driven by a signal in Gaussian noise, such population processes will also give rise to fluctuations. However, in what follows, we shall only consider first-order Markov processes so that the spectral characteristics of the population fluctuations will be Lorentzian in character. This may be contrasted with the Gaussian model, where the spectral characteristics need not have a particular prescribed form.

In §2, we present a review of the doubly stochastic model based on results from Gaussian noise theory, deriving some new results for the joint-probability of event numbers at different times. In §3, we analyse a simple non-stationary population model proposed by Schell & Barakat (1973) that exhibits fluctuations governed by the first-order probability distribution (1.9) and show that it is characterized by a single time constant. Section 4 examines the potential of multiple-immigration models (Jakeman *et al*. 2003) for generating Laguerre statistics. It is shown that these models can be tailored to give correlation properties more closely resembling those found in the case of the homodyne problem. In §5, the predictions of the various models are compared and distinguishing features that might be used to develop the most appropriate model for a given dataset are identified. Conclusions and suggestions for further work are summarized in §6.

## 2. Photo-electron statistics in homodyne configurations

In this section, we shall review the way in which Laguerre number fluctuations arise in the context of optical detection of Gaussian fields. The problem of a constant amplitude signal embedded in Gaussian noise is of interest in many areas of science and engineering and received a great deal of attention early in the development of radio frequency communications, radar and sonar. It is indeed a ubiquitous problem of signal processing and many technical results have appeared in the literature following the fundamental papers of Rice in the 1940s (see Wax 1954). We shall not restate the fundamentals of Gaussian noise theory here, but remind the reader of well-known relevant results as they are required in the present context.

For a narrowband process with a spectral bandwidth of amplitude fluctuations much less than the frequency of the carrier wave, a single polarized component of the electromagnetic field can be represented by the scalar quantity (Middleton 1996)
2.1
Here, *w* is the carrier frequency and the square of the envelope of the field is
2.2
In the case of an un-modulated carrier wave, *a* and *b* do not vary with time so that *I*=*I*_{0} is also constant and formula (1.1) for a doubly stochastic Poisson process predicts result (1.5) as indicated in §1. In the case of narrowband Gaussian noise, *a*(*t*) and *b*(*t*) are statistically independent with Gaussian single interval probability densities. We shall assume henceforth that these densities have zero mean and equal variances, i.e. are the components of a zero-mean, circular, complex Gaussian process. This ensures that the intensity (2.2) has a negative exponential probability density and leads, for a doubly stochastic Poisson process, to result (1.6).

The case of homodyning is a little more complicated. In this case, we note that (2.1) can be written without loss of generality in the form:
2.3
where is the amplitude of the constant component. Assuming that *a* and *b* are independent zero-mean, equal-variance Gaussian variables, a calculation of the Laplace transform of the probability density of intensity fluctuations in the short sample time limit requires evaluation of the following
2.4
The integral with respect to *θ* gives a modified Bessel function of zero-order, and the integral with respect to *I* can then be evaluated exactly to obtain
2.5
Apart from scaling factors, this is seen to be of the form (1.8) and through identity (1.4) generates the Laguerre distribution (1.9) in the case of a doubly stochastic Poisson process. Note that the probability density of the *continuous* variable *I* is obtained by inverse Laplace transformation of (2.5) and takes the well-known form (Middleton 1996)
2.6
Here *I*_{0}(⋅) denotes a zero-order-modified Bessel function of the first kind (Gradshteyn & Ryzhik 1994). Substituting a suitably scaled version of this *Rice* distribution into (1.1) confirms (1.9):
As indicated in §1, the normalized factorial moments of this distribution are identical with the normalized ordinary moments of (2.6) because of the form of relation (1.1):
2.7
So far our results invoke only the properties of Gaussian *variables* although we have introduced at several points the notion of time evolution. When Gaussian variables change continuously in time their behaviour constitutes a process, Gaussian noise, and in addition to the multi-variate Gaussian probability density, this is characterized by correlation properties. In fact, only knowledge of the first-order correlation function of the field *e*(*t*) (or equivalently its Fourier transform, the noise spectrum) is required to completely characterize a *zero-mean circular complex Gaussian process*. Assuming that the process is stationary with a symmetric spectrum, the joint-probability density of intensity fluctuations may be written (Middleton 1996)
2.8
In this formula, the prime indicates the value at a second time and the quantity *g*(*t*) is the correlation function
2.9
This quantity approaches zero as *t* becomes large since *a* and *b* are then independent from their values at later times and have zero mean. At zero-time separation, evidently *g*(0)=1.

An important result that follows from (2.8) is the Siegert relation (Siegert 1943) or Reed factorization theorem (Reed 1962)
2.10
The corresponding *discrete* result is obtained from a simple generalization of formula (1.1) to the joint-probability of finding *n* events at one time and *m* events at a different time, assuming that the efficiency is the same for each measurement and that the integration times are non-overlapping and equal in duration:
2.11
Equation (2.11) preserves the relationship between moments of the joint probability of the number of events and those of the joint probability density of intensity fluctuations and leads to the simple normalized event-number correlation function
2.12
This formula figures prominently in the measurement and analysis of weak optical fields that require the use of photo-electron counting techniques (for early studies on this topic, see Cummins & Pike 1974).

A result for the joint-probability density of intensity fluctuations for the case of a signal in Gaussian noise has been given by Middleton (1996). For the homodyne case considered here, this can be expressed as a sum over modified Bessel functions of the first kind
2.13
Here *ε*_{0}=1, *ε*_{n}=2 for *n*≥1. Fortunately, the double Laplace transform of this result takes a simpler form that is more convenient for calculating correlation properties and for comparison with results based on other models:
2.14
where *τ* refers to the time dependence of *g*. Note that when *s*=0 or *s*′=0, this reduces to the correct result of the form (2.5). When *g*=1 and *g*=0 (corresponding to small and large time separations) it also correctly reduces, respectively, to
2.15
and
2.16
The first of these results corresponds to the two intensities being fully correlated. The second reflects the property of statistical independence possessed by uncorrelated variables.

Note the linear dependence on *g* in the numerator of the exponent in (2.14). This embodies the presence of *interference* between the signal and noise components in this model. The structure of the denominator of the formula and of the exponent reflects the fluctuation behaviour of the Gaussian noise:
2.17
where
2.18
is the joint-generating function of Gaussian noise derived from (2.8) with correlation function *g*.

Formula (2.11) implies that the correlation function of *events* corresponding to these results is
2.19
where the total mean intensity is with *η* the quantum efficiency.

Result (2.14) provides a rigorous benchmark for comparison with alternative models generating Laguerre statistics. In particular for the case of a first-order Gaussian Markov process, it embodies *all* the essential properties associated with the doubly stochastic model. On the other hand, the correlation function (2.19) is more commonly measured and is easier to mimic with alternative models so in what follows our first objective will be to generate Laguerre fluctuations with this kind of correlation function.

## 3. The Schell–Barakat population model

A fundamentally discrete model that can lead to Laguerre fluctuations without recourse to Gaussian noise theory is the birth–death–immigration population model. In the basic birth–death process, individuals are born and die randomly in time at constant rates that are proportional to the extant population size and this leads to a population that either grows without limit or becomes extinct according to whether the birth rate is, respectively, greater or less than the death rate. Immigration stabilizes this critical and non-stationary behaviour so that provided the death-rate is greater than the birth rate, an equilibrium state is reached where the number of individuals in the population simply exhibits fluctuations about a non-vanishing finite mean number. The statistical properties of the population are then stationary (i.e. time-independent) and a wide range of quantities of practical interest can be predicted.

In addition to its role in characterizing the evolution of biological populations (Renshaw 2011), the birth–death–immigration process has been used to model the interaction of species in chemical reactions (Gillespie 1977). It also plays a fundamental role in the theory of photon-counting statistics where it was first used to describe fluctuations of the number of photons in a laser cavity operating below threshold (Shimoda *et al*. 1957). In this case, births and deaths provide analogues of stimulated photon emission and absorption by atoms in the cavity, while immigration is the analogue of ‘spontaneous’ photon creation. The thermal model of a laser below threshold is obtained when the spontaneous (immigration) and stimulated emission (birth) rates are equal as required by quantum mechanics. This special case of the birth–death–immigration process is identical to the Bose–Einstein model encountered in thermodynamics, and the distribution of photons in the cavity is then thermal or geometric as given by equation (1.6). Schell & Barakat (1973) explored the way in which an initially Poisson-distributed population approached this thermal equilibrium state and found that the *transient* state of the population was governed by the Laguerre distribution (1.9) that is of interest here. In this section, we give a brief overview of the thermal model and reproduce the results of Schell & Barakat. This provides a useful introduction to the population models that we shall explore in the next section as well as highlighting differences with the Gaussian noise model that we have already solved.

The birth–death–immigration model is a first-order Markov process that can be expressed entirely in terms of transitions in which the population changes by unity. When individuals appear and die *randomly* in time, a stochastic representation of the problem is required and this is provided by consideration of the way in which the probability of finding a given number in the population evolves. Thus, births and immigration in a population of *n*−1 individuals and deaths in a population of *n*+1 individuals will all increase the probability of finding a population of size *n*. On the other hand, if any one of these takes place in a population of *n* individuals, then the probability of finding *n* individuals will be reduced. Thus, we can write down a rate equation for the evolution of the probability of finding *n* individuals in the population at time *t* in the form (Bartlett 1966):
3.1
Here *μ*, *λ*, *ν* are the positive proportionality constants corresponding, respectively, to the rates at which deaths, births and immigration are taking place. Solving this equation can be accomplished most easily by first transforming (3.1) into an equation for the generating function *Q*(*s*;*t*)=〈(1−*s*)^{n}〉
3.2
The solution to this linear, first-order partial differential equation must satisfy the initial condition at *t*=0 and the also condition *Q* (0;*t*)=1, which ensures that the corresponding probability distribution is normalized to unity at all times. The result is well known (Bartlett 1966) and may be expressed in the form:
3.3
Here and , assuming that *t*>0 and *μ*>*λ*. In the case of equal birth and immigration rates, this reduces to
3.4
In this last formula, the mean number of individuals in the population is . This clearly satisfies the initial boundary condition and also the normalization requirement. Schell and Barakat considered the case when the initial state of the population was Poisson distributed (or in their terminology ‘a single mode of the radiation field’). According to equation (1.7), the generating function is then a negative exponential function of *s* so that, on setting *λ*/*ν*=1, (3.4) can be written
3.5
In this equation, we have written
where *n*_{0} is the mean of the initial Poisson distribution. Result (3.5) is evidently of the Laguerre form (1.9) except that the parameters of the distribution are evolving in time. Indeed, the distribution changes from the initial Poisson distribution to a thermal distribution as a result of the population processes taking place. However, note that during the transient phase, the two population types do not simply add as this would lead to a product of the two generating functions (1.7).

Although the Schell–Barakat model is non-stationary, the coherence function,
3.6
with *τ* the delay time and *t* the running time, can be separated. This quantity can be calculated by expressing the correlation coefficient in terms of a sum over the conditional distribution that is related to the solution (3.4):
3.7
3.8
In this formula, the prime denotes the second time and *p*(*n*′|*n*) is the probability of finding *n*′ individuals in the population at time *t*+*τ* given that there were *n* present at the earlier time *t*. The generating function in (3.7) is result (3.4) evaluated at the initial condition corresponding to exactly *n* individuals in the population, i.e. *Q* (*s*;0)=(1−*s*)^{n}. The resulting coherence function for the Schell–Barakat model is thus
3.9
This may be compared with the corresponding Gaussian-noise result that can be obtained from (2.19):
3.10
For a meaningful comparison, we must assume that with *τ*>0 since we are comparing the doubly stochastic Gaussian-noise result with that predicted by a first-order Markov process. There are two significant differences between (3.9) and (3.10). The less important one is the presence in (3.9) of additional terms linear in the mean values, *n*_{0} and . We have shown in recent papers (Matthews *et al*. 2003; French *et al*. 2008) that these terms are a common feature of the correlation properties of populations, reflecting the discrete nature of the models. However, the linear terms are eliminated when these events are monitored via individuals that leave the population during fixed time intervals that are sufficiently short and close to each other. This is because detecting an individual diminishes the population size from *N* to *N*−1, and the measurement interval is so short that the population does not have time to change during the next measurement interval. Consequently, the correlation function of the monitored individuals is proportional to the population size at these two measurement times, i.e. 〈*n*(*n*−1)〉 rather than 〈*n*^{2}〉 which produces the linear term, as explained in French *et al*. (2008). The doubly stochastic Poisson process, on the other hand, directly models the properties of this measured train of events and the linear terms do not appear. In making comparisons between the predictions of the doubly stochastic Poisson model and population models, we can, therefore, neglect linear terms of this kind. The second and more significant difference between (3.9) and (3.10) is the presence of *two* time scales in (3.10) but only *one* in (3.9). It is not clear how the two models can be reconciled in this respect and in view of the non-stationary nature of the Schell–Barakat model, we shall not investigate it any further in the present paper.

## 4. Multiple-immigration models

This section describes a family of Markovian population models that are based on the processes described in the last section, but include the possibility of *multiple* immigrations, i.e. simultaneous input of one, two or more individuals into the population, irrespective of the number already present. This type of model has been described earlier (Jakeman *et al*. 2003) and allows certain statistical properties of the population to be tailored to requirements. In particular, the first-order distribution and the correlation function of the number fluctuations can be determined by the choice of model for the immigrant community. In the earlier work, the rates of immigration for the different numbers of individuals were drawn from a statistical distribution that could describe another type of population evolving by different processes and characterized by different fluctuation times. This class of multiple-immigration model, therefore, predicts the properties of a population that interacts with an *ensemble* of identical populations (which may be of a different type) in different states of their evolution. It is analogous to a grand-canonical ensemble, where a population is influenced by those forming its environment, but without that population affecting its environment.

### (a) Model 1

This consists of a thermal model (1) driving another thermal model (2) by means of additional multiple immigrations. The parameters of the thermal populations are a scaled version of each other. This is taken into account here by scaling the time constants. The generating functions for the basic thermal model and the driven population satisfy the equations (where *Q*=*Q*(*s*;*t*))
4.1
and
4.2
The first of these equations is the same as (3.2) apart from the factor scaling the time derivative on the left-hand side. The final term on the right-hand side of the second equation is the multiple-immigration term that couples the population of interest to the ensemble of populations characterized by (4.1). As discussed in the previous paper (Jakeman *et al*. 2003), this arises from the contribution of multiple immigration to the evolution of the population distribution:
4.3
The rates, *ε*_{r}, are interpreted as probabilities corresponding to the population process (1) with *p*_{r} being the probability of finding *r* individuals in population (2). Multiplying (4.3) by (1−*s*)^{n} and summing over *n* leads to the last term on the right-hand side of (4.2) with the constant rate *ε* being the strength of the coupling.

The stationary solution of the first equation is (1.7b): 4.4 and it is not difficult to show that the stationary solution of the second equation is then 4.5 which reduces to (1.8) upon making the identification .

The correlation function of number fluctuations in population (2) can be obtained by first solving the full time-dependent equations for the mean as a function of fixed population number at *t*=0. Multiplying the number at time *t* by the number at *t*=0 and averaging over the equilibrium distribution then gives the correlation function. Using the fact that at any given instant of time, the population number in *n*_{1} is independent of those in *n*_{2} by virtue of (4.1), enables the following result to be obtained:
4.6
Note that this result is obtained assuming that the *Q*_{1} appearing in the multiple-immigration term appearing in equation (4.2) is the time-dependent solution of equation (4.1). The coherence function corresponding to the result (4.6) is identical to that obtained for the Schell–Barakat model and differs from that predicted for a doubly stochastic Poisson process driven by the square of a zero-mean complex Gaussian Markov process. In particular although (4.6) reduces to the expected results in all the relevant limiting cases, it is not of the form (2.19) by virtue of the time-dependent terms exhibiting different decay rates. In the present model, this is because the equilibrium correlation function of population (2) is independent of the fluctuations in population (1). Further examination of this model has failed to reveal any reasonable assumptions that would achieve two time scales in the correlation function despite the fact that two appear in the combined model if *α*≠*β*.

### (b) Model 2

It would seem likely that obtaining a correlation function with two decay times requires a model for the *sum* of two populations. This can be achieved by dropping the single immigration term in equation (4.2) and assuming that we require the properties of the *combined* populations. Thus, equation (4.2) is replaced by
4.7
where again *Q*=*Q* (*s*;*t*). The stationary solution of this equation may be written
4.8
The corresponding mean and second factorial moment are
4.9
and
4.10
Combining equations (4.4) and (4.8) leads to the generating function for the Laguerre distributions (1.8).

Calculation of the correlation function of the combined populations can be carried out as before but now
4.11
It is not difficult to show that
4.12
4.13
4.14
4.15
The difference between results (4.14) and (4.15) arises from the fact that the number in population (1) at time *t*=*τ* is not related to the number in population (2) at *t*=0, while the number in population (2) at time *τ* is correlated with the number at earlier times in population (1). It is clear that the correlation function of the total population contains two time scales unlike model 1. However, these cannot be expressed in the form (2.19).

### (c) Model 3

A very simple modification of model 2 does generate a correlation function of the form (2.19). We shall simply assume that the multiple-immigration rates in (4.7) are the *stationary* solutions of equation (4.1) rather than being time-dependent. This means that the populations (1) and (2) are statistically independent. Result (4.8) for the equilibrium-generating function remains valid but now
4.16
From results (4.12) and (4.13), which also remain valid for this model, we obtain the correlation function
4.17
This is of the form (2.19), when *α*=*β*/2 apart from the terms linear in *n*_{0}, *n*_{G} that do not occur when making measurements on a train of events (cf. earlier remarks on the suppression of these terms, as described in French *et al*. (2008)).

It is interesting to examine some other properties of population (2) for model 3. For example by the use of Laplace transforms, the full time-dependent generating function can be shown to be of the form:
4.18
Here, the first factor on the right-hand side is the initial generating function of the distribution and
4.19
The equilibrium probability distribution and its factorial moments can be obtained by expanding (4.8) in powers of (1−*s*) and *s*, respectively:
4.20
4.21
and
4.22
Using *L*_{1}(*x*)=1−*x*, *L*_{0}(*x*)=1 (Gradshteyn & Rhyzik 1994) immediately obtains the correct result (4.9) for the mean.

Note that there exists a continuous probability density that has a Laplace transform of the form (4.8). This may be expressed in terms of *I*_{1}(⋅), a modified Bessel function of the first kind:
4.23
A continuous variable with probability density of this form would lead to the discrete distribution (4.21) through the doubly stochastic Poisson process (1.1). The distribution (4.23) arises in the context of the diffusion approximation of a branching process (Cox & Miller 1995), and in scattering of electromagnetic radiation from the sea surface (Clements & Yertsever 1997; Ward *et al*. 2006).

It is important to note that in the basic equation (4.7), the parameter *ε* is defined to be so that the multiple-immigration coupling term can be written explicitly as
4.24
This means that the limit that corresponds to *λ*→0 results in a simple death–immigration process for population 2.

## 5. Discussion

For Markovian systems, the lowest order correlation function and joint statistics at two successive times will be sufficient to specify the entire statistical model. We have seen that, although non-stationary, the Schell–Barakat population model can generate first-order Laguerre statistics that are the same as those for a doubly stochastic Poisson process driven by the intensity fluctuations of a signal in circular complex Gaussian noise, but that the correlation coefficient is characterized by a single time constant, lacking the term normally associated with ‘interference’ between signal and noise in the correlation of optical intensities. Similarly, the first of the stationary multiple-immigration models can reproduce the first-order Laguerre statistics, but its correlation function is again characterized by a single exponential decay constant. Although multiple-immigration model 2 is generally characterized by first-order Laguerre statistics and *two* correlation times, the structure of its correlation function is different from (2.19).

By contrast, the correlation function of a population governed by model 3 is characterized by two time-constants and these *can* be chosen to have the relationship found in result (2.19) so that measurements of the first-order statistics and the correlation function cannot be used to distinguish the models. However, it is still possible that the higher order correlation properties of the two models may be different and this requires knowledge of the full joint-distribution of population numbers at two different times. In the case of the doubly stochastic Poisson process driven by the intensity of a signal in complex Gaussian noise, we have shown that this is given by equations (2.13) or (2.14) and it is interesting to calculate a comparable quantity for the multiple-immigration model 3, as follows.

Assuming that the two discrete populations involved are statistically independent, their two joint-generating functions can be simply multiplied together. For the thermal population: we use the well-known result 5.1 with and .

For the case of the second population governed by rate equation (4.7) with the equilibrium form for *Q*_{1}, we find
5.2
The full result for the joint-generating function for the discrete model 3 is the product of these two results. It is immediately apparent from considerations of symmetry with respect to interchange of *s* and *s*′ that this product differs from (2.14). Moreover for model 3, the decay rates appearing in *θ*_{1} and *θ*_{2} must differ by a factor of 2 to obtain the correlation function (2.19) and this renders the denominator in (5.1) different from the denominator of the exponent in (5.2), unlike the structure of (2.14). Thus, although the first-order statistics and second-order correlation function predicted by model 3 are the same as those predicted for the doubly stochastic model driven by a complex Gaussian noise process, the higher order correlation properties of the models are manifestly different. The simplest quantity that might be measured to reveal this difference is the higher order correlation property
5.3
On neglecting terms *O*(*n*^{2}) in accord with the remarks made earlier concerning detection, we obtain for the discrete model
5.4a
and
5.4b
whereas
5.5
Noting that the second-order correlations coincide if *θ*_{2}=*g*, *θ*_{1}=*g*^{2}, inspection of the above shows that with this choice (5.4b) measure is identical with the Gaussian result (5.5), and that differences between this and the other correlation function occur only in the mixing term that is proportional to *n*_{0}*n*^{2}_{G}. The difference between (5.5) and (5.4a) is
which at its greatest discrepancy is ∼*n*_{0}/2*n*_{G} if *n*_{0}/*n*_{G}≪1, or ∼(*n*_{0}/*n*_{G})^{−2}/2 if *n*_{0}/*n*_{G}≫1. Indeed, the maximum value of the difference is 2/27∼7%, which occurs when *g*=*n*_{0}/*n*_{G}=1/2.

## 6. Concluding remarks

In this paper, we have investigated a number of stochastic models that can lead to discrete Laguerre statistics that interpolate between Poisson and thermal fluctuations. We have shown that although, in general, the models are characterized by different kinds of temporal evolution, examples can be found where the second-order correlation functions of two different models are identical and where higher order coherence properties may be almost indistinguishable. Generalizations of these results can be envisaged. For example, the thermal model is but one member of the class of birth–death–immigration processes and the analogous continuous process that may be used instead of *Circular* Gaussian statistics used here to drive a doubly stochastic Poisson model would be a gamma process based on *multi-dimensional* Gaussian statistics (Jakeman 1991).

Finally, we note that the work reported here forms part of an ongoing investigation of coherent interactions between discrete populations. A study of measurements on populations using detection processes that cannot distinguish the origin of individuals being measured is currently in progress.

## Acknowledgements

We would like to thank the referees for their suggestions for improving the clarity of the paper, and for drawing our attention to the presence of the delta-function appearing in equation (4.23), and the interpretation of the distribution.

- Received September 13, 2011.
- Accepted January 30, 2012.

- This journal is © 2012 The Royal Society