IOP PUBLISHING JOURNAL OF PHYSICS A: MATHEMATICAL AND THEORETICAL
J. Phys. A: Math. Theor. 41 (2008) 355002 (20pp) doi:10.1088/1751-8113/41/35/355002
The evolution of a single-paired immigration death
process
Colin S Gillespie1 and Eric Renshaw2
1 School of Mathematics and Statistics, University of Newcastle, Newcastle upon Tyne NE1
7RU, UK
2 Department of Statistics and Modelling Science, Livingstone Tower, University of Strathclyde,
26 Richmond Street, Glasgow G1 1XH, UK
Received 30 April 2008, in final form 2 July 2008
Published 28 July 2008
Online at stacks.iop.org/JPhysA/41/355002
Abstract
The general question of whether it is possible to determine the fundamental
structure of a hidden stochastic process purely from counts of escaping
individuals is of immense importance in fields such as quantum optics, where
externally based radiation elucidates the nature of the electromagnetic radiation
process. Although the general probability structure has been derived in an
earlier paper in terms of the joint probability generating function of the (hidden)
population size and (known) counts, its complex nature hides some particularly
intriguing features of the underlying process. Our current objective is therefore
to examine specific immigration regimes in order to highlight the underlying
saw-tooth behaviour of the underlying probability and moment structures. The
paper first explores paired- and triple-immigration schemes, and then introduces
birth in order to show that the technique is equally successful in exposing hidden
multiplicative effects. These analyses uncover novel and highly illuminating
features, and emphasize the potential of this population-counting construct for
expanding into more complex multi-type situations.
PACS numbers: l05.40.−a, 05.49.Df, 89.75.−k
1. Introduction
Over the past 40 years stochastic population processes have figured prominently in the
description of phenomena across a wide variety of fields in Science and Engineering. An area
of application that has been particularly successful has been quantum optics, where population
models have been routinely used to describe the quantum nature of electromagnetic radiation
(Srinivasan 1988). An especially interesting problem concerns the stochastic evolution
of populations of photons within optical cavities. This field has been studied since the
inception of the laser (Shimoda et al 1957), and understanding the counting statistics of photo-
electron pulses registered by detectors such as photomultiplier tubes has been essential for
1751-8113/08/355002+20$30.00 © 2008 IOP Publishing Ltd Printed in the UK 1

J. Phys. A: Math. Theor. 41 (2008) 355002 C S Gillespie and E Renshaw
the interpretation of experimental measurements (Saleh 1978). Jakeman and Shepherd (1984)
and Shepherd (1984) monitor the cavity population of interest via the counting statistics of
emigration, modelled as a simple death process. Properties of the number of emigrants leaving
the population in a fixed time interval correspond precisely to the experimentally measured
photon-counting statistics, and so provide an indirect measure of the evolution of the cavity
population. This approach not only provides additional insight into the quantum formulation
of the problem, but it also enables the interchange of models and techniques with those in the
field of classical population statistics. In this context, Jakeman et al (1995) draw attention to
the unusual properties of an exactly solvable population model which is generic to the recently
developed area of quantum optics known as ‘non-classical light’.
For many years following the invention of the laser, it was generally accepted that an
adequate representation of photodetection was provided by the doubly stochastic Poisson
process. In this classical situation the probability of registering c such pulses during a time
interval of fixed length T (see Mandel (1959), Cox and Lewis (1966)) is given by
pc(t) =
∫
∞
0
I c e−I
c!
f (I) dI, (1.1)
where I is the instantaneous light intensity integrated over the interval (0, t) and f (I) is its
probability density function. It follows from (1.1) that the variance to mean ratio, ρ, of c (i.e.
the Fano factor) must always be greater than or equal to that of the Poisson distribution, for
which ρ = 1. Light with counting statistics which can be represented in this way such as
coherent (laser) light, which is Poisson distributed, or thermal light, which satisfies a geometric
distribution, is termed classical light. However, it is easy to construct discrete models for
the incident photon flux which cannot be derived through the representation (1.1). Whilst
sub-Poissonian models (ρ < 1) clearly cannot be so represented, neither can a wide range of
super-Poissonian models (ρ > 1). Such light is now termed non-classical (for references see
Jakeman et al (1995) and Gillespie and Renshaw (2005)).
Intense activity in the development of experimental methods for generating non-classical
light has produced a range of techniques which provide overwhelming evidence for its
existence. Since most applications have required consideration of both wave and particle
properties, theoretical treatments have generally avoided a classical population statistics
approach. One of the few exceptions is Jakeman et al’s (1995) investigation of light, which
involves the simultaneous emission of pairs of photons at rate α2 and the subsequent death
of each individual single photon at rate µ. This was itself stimulated by the burgeoning
area of non-classical light in the late 1980s (Louden and Knight 1987). Although the
equilibrium distribution is super-Poissonian, odd–even effects ensure the breakdown of (1.1),
so the resulting light is non-classical. Not only does the simplicity of this model mean that
it is analytically amenable, but it reflects one of the earliest mechanisms used to produce
non-classical light through parametric down-conversion in a nonlinear crystal (Burnham and
Weinberg 1970), and so is clearly realizable experimentally.
Gillespie and Renshaw (2005) expand the work of Jakeman et al (1995) by investigating
the general problem where batches of photons enter the population with constant rate αqi ,
where
∑
∞
i=1 qi = 1. This formulation subsumes many previous results:
basic single-immigration-death process q1 = 1 and qi = 0 for i = 1
paired immigrants (Jakeman et al 1995) q2 = 1 and qi = 0 for i = 2
geometric immigration (Matthews et al 2003) qi = (1 − ξ)ξ i−1 for 0 < ξ < 1
power-law process (Hopcraft et al 2002) qi = ν(i − ν)/[(1 − ν)i!]
for 0 < ν < 1.
2

J. Phys. A: Math. Theor. 41 (2008) 355002 C S Gillespie and E Renshaw
Gillespie and Renshaw also examine two further immigration strategies. First, they consider
the case in which each batch comprises exactly k immigrants, so that qk = 1 and qj = 0 for
j = k, and develop attractive expressions for the population size probabilities and cumulants,
together with a particularly neat approximation for the cumulants of the associated counting
process. The population and counting behaviours exhibit similar manifestations, in both their
probability and moment structures, thereby demonstrating that inference can be made on
hidden population behaviour purely by examining the structure of successive counts. They
then tackle a related question of uniqueness. For since an infinite number of probability
measures can give rise to the same set of moments, there is no reason to presuppose that
such similarity will hold universally. So for their second case they exploit a classic example
of this phenomenon due to Schoenberg (1983); see Stoyanov (1988) for a comprehensive
review of such probability classes. Essentially, they retain the individual death rate µ, but now
allow batches of immigrants of size 2k to enter the population. Specifically, for ||
1 the
immigration rate takes the Schoenberg–Poisson form
q2k =
e−22k
k!
{1 + k!(−1)k[(2 − 1)(22 − 1) · · · (2k − 1)]−1} (1.2)
with qi = 0 otherwise. Expressions are then derived to show that although the moment
structure remains independent of  for both the population and counting processes, the
corresponding probability structures do indeed change as  varies.
However, the rather complex nature of these examples hides some particularly intriguing
features of the underlying process, and so our current objective is to use simpler immigration
regimes in order to illuminate the underlying stochastic behaviour.
2. The batch immigration-death-counting process
Suppose that a stochastic process is monitored purely externally, by counting (for example)
the number of individuals emigrating during a fixed time interval (0, t), with emigration being
modelled by a pure death process with rate η. Note the fine distinction between death and
emigration. For as well as allowing µ to denote the individual death rate and η the individual
‘escape’ rate from the system (with all escapees being subsequently counted); we could also
have an individual death rate of µ+η and ‘inefficient counting’, with each individual who dies
having probability η/(µ+η) of having their death recorded. The raison d’eˆtre which underlies
this approach is of fundamental significance, since it provides a means of answering a general
question which is of considerable importance across a wide range of disciplines, namely: ‘If
a stochastic process is developing within a hidden system, with the only information provided
being the event times of ‘escaping’ individuals, can the properties of the hidden process be
inferred purely from knowledge of these counting statistics?’. Though here we present an
algebraically tractable example based on immigration and death (and later on birth as well), the
counting concept is easily extended to all other types of processes. For example, in biological
and epidemiological scenarios we might record death and notifiable illness, respectively.
Let pnc(t) denote the probability that at time t > 0 a population contains n individuals and
that c emigrants have been counted. Moreover, suppose that batches of i immigrants enter the
population at constant rate αqi , where
∑
∞
i=1 qi = 1, and let each individual die at rate µ. The
population is monitored externally by counting the number of individuals who emigrate at rate
η during a fixed time (0, t). Then under the assumption that individuals develop independently
from each other, the joint population size/count probabilities {pnc(t)} are defined through the
3

J. Phys. A: Math. Theor. 41 (2008) 355002 C S Gillespie and E Renshaw
Kolmogorov forward equation
dpnc(t)
dt
= µ(n+ 1)pn+1,c(t)+ η(n+ 1)pn+1,c−1(t) + α
∞
∑
i=1
qipn−i,c(t) − [(µ + η)n + α]pnc(t),
(2.1)
where n, c = 0, 1, 2, . . . and pnc(t) = 0 for n, c = −1,−2, . . . . On defining the joint
probability generating function
Q(z, s; t) ≡
∞
∑
n,c=0
znscpnc(t), (2.2)
and denoting δ = µ + η, equations (2.1) reduce to the single partial differential equation
∂Q
∂t
+ (δz − ηs − µ)
∂Q
∂z
= αQ
∞
∑
i=1
qi(z
i
− 1). (2.3)
Whence on assuming that the population is of size zero at time t = 0, Gillespie and Renshaw
(2005) obtain the general solution
Q(z, s; t) = exp
(
∞
∑
i=1
αt(ηs + µ)iqi
δi
− αt
)
× exp
(
∞
∑
i=1
αqi
δi+1
i
∑
r=1
1
r
(
i
i − r
)
(ηs + µ)i−r (δz − ηs − µ)r(1 − e−rδt )
)
. (2.4)
It should be noted that the above derivations can also be framed in the context of a marked
Poisson process (see Kingman 1993). In this scenario we fix a time τ and mark a particular
batch immigration event formed before τ by its number of descendants at time τ . Each point
in this marked space is a pair (T ,D) denoting the random time, T, of the occurrence of the
cluster and its eventual number of descendants, D. The total number of particles can now be
obtained through Campbell’s formula.
2.1. Paired immigration
Suppose that single and paired immigrants arrive at rate αq1 = α1 and αq2 = α2, respectively,
so qi = 0 (i > 2). Moreover, each individual either dies at rate µ or escapes (i.e. emigrates)
and is subsequently counted at rate η. Then the general p.g.f. solution (2.4) reduces to
Q(z, s; t) = exp
(α1
δ2
(ξ(1 − e−δt ) + δηt (s − 1))
)
exp
( α2
2δ3
(2δt (ηs + µ)2 − 2δ3t)
)
× exp
( α2
2δ3
(ξ 2(1 − e−2δt ) + 4ξ(1 − e−δt )(ηs + µ))
)
, (2.5)
where δ = µ + η and ξ = δz − ηs − µ. Although Gillespie and Renshaw (2005) provide the
general marginal probabilities for this process, these hide the underlying odd–even structure.
To expose this we first note that the generating function representation for Hermite polynomials
is
∞
∑
n=0
Hn(x)z
n/n! ≡ exp(2xz − z2) (2.6)
(Abramowitz and Stegun [A&S] (1970), result 22.9.17), and that we also have the relationships
H2n(x) ≡ (−1)n22nn!L−1/2n (x2) and H2n+1(x) ≡ (−1)n22n+1n!xL1/2n (x2) (2.7)
4

J. Phys. A: Math. Theor. 41 (2008) 355002 C S Gillespie and E Renshaw
(A&S, 22.5.40/41) between Hermite and Laguerre polynomials. Then to recover the even
and odd time-dependent probabilities, pn(t) = pn.(t), from (2.5), extracting the coefficients
of z2n and z2n+1 at s = 1 yields for k = 0, 1,
p2n+k(t) =
n!22n
(2n + k)!
(
α2(1 − e−2µt )
2µ
)n
f1(k) exp
(
(2α1 + α2(3 − e−µt ))(e−µt − 1)
2µ
)
×Lk−1/2n
(
(α1 + α2(1 − e−µt ))2(1 − e−µt )2
2α2µ(e−2µt − 1)
)
, (2.8)
where f1(0) = 1 and f1(1) = (α1 + α2(1 − e−µt ))(1 − e−µt )/µ. Likewise, placing z = 1
recovers the marginal counting probabilities pc2c(t) = p.,2c(t) and pc2c+1(t) = p.,2c+1(t),
namely for k = 0, 1,
pc2c+k(t) =
c!
(2c + k)!
{
2α2η2
δ3
[2δt − (1 − e−δt )(3 − e−δt )]
}c
f2(k) exp
{
α1η
δ2
[1 − δt − e−δt ]
}
× exp
{
α2
2δ3
[η(4µ(1 − e−δt ) + η(1 − e−2δt )) + 2δt (µ2 − δ2)]
}
×Lk−1/2c
{
−
θ
2δ3α2[(1 − e−2δt ) − 4(1 − e−δt ) + 2δt]
}
, (2.9)
where
θ = {α1δ[δt + e−δt − 1] + α2[2(η − µ)(1 − e−δt ) − η(1 − e−2δt ) + 2δµt]}2,
f2(0) = 1,
and
f2(1) =
α2η
δ3
[2(η − µ)(1 − e−δt ) − η(1 − e−2δt ) + 2δµt] + α1η
δ2
(δt + e−δt − 1).
Although at first glance solutions (2.8) and (2.9) look rather opaque, closer inspection soon
exposes the odd–even effect in the probability structure. For p2n(t) and p2n+1(t), and also
pc2c(t) and pc2c+1(t), essentially differ only by the powers of the Laguerre polynomial, namely
L−1/2c (·) and L1/2c (·), respectively. It is this ‘saw-tooth’ effect that causes the generation of
non-classical light.
The factorial moments
{
Mcr (t)
}
of the counting distribution may also be obtained directly,
by successively differentiating Q(z, s; t) with respect to s and forming Q(2r)(1, s; t) and
Q(2r+1)(1, s; t). This procedure yields for k = 0, 1,
Mc2r+k(t) = r!22r
{
α2η2
2δ3
[2δt − (1 − e−δt )(3 − e−δt )]
}r
×Lk−1/2r
{
(α1 + 2α2)2(δt + e−δt − 1)2
2α2δ[(1 − e−δt )(3 − e−δt ) − 2δt]
}
f3(k), (2.10)
where f3(0) = 1 and f3(1) = (η/δ2)(α1 + 2α2)(δt + e−δt − 1). These expressions clearly
highlight the odd–even effect in the same manner as with the counting probabilities; the only
substantive difference between the counting moments Mc2r (t) and Mc2r+1(t) once again lies in
the powers of the Laguerre polynomials, namely L−1/2r (·) and L1/2r (·).
To illustrate how this odd–even structure degrades when single immigrants are allowed,
suppose we take time t = 2 and a total immigration rate of 10, i.e. α1 + 2α2 = 10. Figure 1
shows how changing the balance between α1 and α2 affects the counting probabilities pc(2),
and it is clear that the extreme saw-tooth shape of the distribution is only (visually) present
when α1  0. This can be explained by considering
p˜0(t) = Pr(no deaths from single immigrants by time t)
= exp{(α1/η)(1 − e−ηt − ηt)}. (2.11)
5

J. Phys. A: Math. Theor. 41 (2008) 355002 C S Gillespie and E Renshaw
Cou
nts
0
10
20
30
40
50alpha[1]
0
2
4
6
8
10
Probability
0.00
0.02
0.04
0.06
0.08
0.10
0.12
Figure 1. Count probabilities pcc(t) for η = 300 and µ = 0, with α1 + 2α2 = 10 for varying α1
and fixed time t = 2.
For when α1/η  0, (2.11) simplifies to p˜0(t)  e−α1t , and since this swiftly decays to zero
with increasing α1t the intensity of the odd–even effect must decay with it. Indeed, the time
required to reduce p˜0(t) to (say) 0.5, namely ln(2)/α1, reduces inversely with α1.
The overall effect that α1 induces on the even and odd probabilities is best seen by recalling
that Laguerre and Hermite polynomials are neatly related through (2.7). For on applying
∞
∑
n=0
(−1)n
(2n)!
H2n(x)z
n
≡ ez cos(2x
√
z) (2.12)
∞
∑
n=0
(−1)n
(2n + 1)!
H2n+1(x)z
n
≡ z−1/2 ez sin(2x
√
z) (2.13)
(A&S section 22.9.18/19) to the full probability expressions (2.9) with µ = 0 and η large
compared with (α1 + α2)t , we obtain
∞
∑
c=0
pc2c(t) = e
−α1t sinh(α1t) and
∞
∑
c=0
pc2c+1(t) = e
−α1t cosh(α1t). (2.14)
Since tanh(α1t) ↑ 1 as t → ∞, it follows that when
∑
∞
c=0 p
c
2c(t)  0.5 the odd–even effect is
no longer ‘visible’ as the sums of the even and odd probabilities are equal. Rather surprisingly,
expressions (2.14) do not depend on α2, indicating that the degradation of the odd–even effect
depends solely on α1 and t. Figure 2 illustrates how the probabilities (2.9) change through
time when α1 = 0.075 is substantially smaller than α2 = 1.2. Whilst a strong odd–even
effect is clearly present at times t = 10 and t = 15, as t increases further the saw-tooth shape
becomes much less pronounced. Quantitatively, we see from (2.11) that the probability of
zero single immigrant deaths reduces from p˜0(10) = 0.472 to p˜0(30) = 0.105. So at time
t = 30, the probability of no single-immigrant deaths is about five times smaller than that at
t = 10, resulting in a much reduced odd–even effect.
6

J. Phys. A: Math. Theor. 41 (2008) 355002 C S Gillespie and E Renshaw
0 20 40 60 80 100 120
0.
00
0.
02
0.
04
0.
06
0.
08
Counts
Pr
ob
ab
ilit
y
t=10
t=15
t=20
t=25
t=30
Figure 2. Count probabilities pc(t) for η = 300 and µ = 0, with α1 = 0.075 and α2 = 1.2 for
varying times t = 10, 15, . . . , 30.
To determine the corresponding effects on the second-order moments, first note that the
population variance
Var[n(t)] = {2α2(1 − e−2µt ) + [α1 + α2(1 − e−µt )](1 − e−µt )}/µ (2.15)
is linear in α1 and α2; whilst the autocovariance function
E[n(t)n(t + s)] − E[n(t)]E[n(t + s)] = e−µs(1 − e−µt )[α1 + α2(3 + e−µt )]/µ. (2.16)
Hence as t → ∞, the limiting autocorrelation function
ρ(∞; s) = e−µs, (2.17)
which, perhaps surprisingly, does not depend on either α1 or α2.
To obtain the associated counting statistics we parallel the approach taken by Shepherd
and Jakeman (1987), who develop the correlation properties of the counting process from
the joint probability Pr{c1; (t1, t1 + T ) and c2; (t2, t2 + T )} of counting c1 individuals during
the period (t1, t1 + T ) and c2 individuals during the subsequent, and non-overlapping period,
(t2, t2 + T ). Though our analysis is slightly different, since we use a specific initial condition
(e.g. p00(0) = 1) rather than assuming that the process starts in equilibrium, in order to
examine both transient and persistent behaviour. First note that
Pr{c1; (t1, t1 + T ) and c2; (t2, t2 + T )}
=
∞
∑
N1,N2,N3,N4=0
Pr(size N1 at time t1 | size n0 at time 0)
×Pr(counting c1 in (t1, t1 + T ) and size N2 at time t1 + T | size N1 at time t1)
×Pr(size N3 at time t2 | size N2 at time t1 + T )
×Pr(counting c2 in (t2, t2 + T ) and size N4 at time t2 + T | size N3 at time t2)
=
∞
∑
N1,N2,N3=0
Pr(N1; t1 | n0; 0)Pr(c1; (t1, t1 + T ) and N2; t1 + T | N1; t1)
× Pr(N3; t2 | N2; t1 + T )Pr(c2; (t2, t2 + T ) | N3; t2). (2.18)
Then after some algebra, it may be shown that the counting autocorrelation function for the
single-paired immigration process is given by
7

J. Phys. A: Math. Theor. 41 (2008) 355002 C S Gillespie and E Renshaw
ρc(c1(t1, t1 + T ), c2(t2, t2 + T )) =
α2η2 e−δ(t1+t2)(1 − e−δT )2(eδ(2t1+T ) − 1)
δ3
√
Var(t1, t1 + T )Var(t2, t2 + T )
−
n0η2 e−δ(t1+t2)(1 − e−δT )2
δ2
√
Var(t1, t1 + T )Var(t2, t2 + T )
, (2.19)
where, for i = 1, 2,
Var(ti , ti + T ) =
α1ηT
δ
−
α1η e−δti (1 − e−δT )
δ2
+
2α2ηT (δ + η)
δ2
−
2α2η2(1 − e−δT )
δ3
−
η2 e−2δti (α2 + δn0)(1 − e−δT )2
δ3
+
η e−δti (δn0 − 2α2)(1 − e−δT )
δ2
. (2.20)
Placing t2 = t1 + T + τ in the counting a.c.f. (2.19) and letting t1 → ∞ then yields
ρc(c1(t1, t1 + T ), c2(t1 + T + τ, t1 + 2T + τ))
→
α2η(1 − e−δT )2 e−δτ
2α2(T (2δ2 − µ2 − ηµ) − η(1 − e−δT )) + α1δ2T
. (2.21)
This result reveals four interesting insights into the counting process. First, it shows how
the correlation decreases as τ increases; as the two time periods being compared become further
apart the corresponding (driving) population processes have fewer individuals in common.
Second, it describes how as α1 increases, the correlation decreases; for the process becomes
dominated by single immigrants. Indeed, when α2 = 0 the process reverts back to a Poisson
process for which ρ ≡ 0. Third, when α1 = 0 expression (2.21) is independent of α2. Finally,
as the counting rate η → ∞, we see that ρc → 0, since immigrants entering the system are
immediately removed, and so the non-overlapping time intervals cannot contain individuals
common to both.
2.2. Triple immigration
Having seen that the general solution (2.4) produces a highly definitive saw-tooth effect under
paired immigration, provided α1  α2, questions naturally arise as to the specific form of
the solution under say triple immigration with α3 = αq3 α1, α2, and what happens in a
multiplicative environment in which the complexity of the underlying population process is
increased by including birth. Tackling the former scenario first, we assume, as before, that
individuals leave the system at rate η, die unnoticed at rate µ and are counted at rate η. So the
general forward Kolmogorov equation (2.1) now takes the specific form
dpnc(t)/dt = µ(n + 1)pn+1,c(t) + η(n + 1)pn+1,c−1(t) + α3pn−3,c(t)
− (µn + ηn + α3)pnc(t), (2.22)
for n, c = 0, 1, 2, . . . , where we define p
−1,c(t) = p−2,c(t) = p−3,c(t) = 0. This equation
may be solved in terms of the associated generating function
R(z, s; t) = exp
{
α3
6δ4
18(ηs + µ)2(δz − ηs − µ)(1 − e−δt )
}
× exp
{
α3
6δ4
9(ηs + µ)(δz − ηs − µ)2(1 − e−2δt )
}
× exp
{
α3
6δ4
2(δz − ηs − µ)3(1 − e−3δt )
}
× exp
{
α3
6δ3
6t ((ηs + µ)3 − δ3)
}
, (2.23)
8

J. Phys. A: Math. Theor. 41 (2008) 355002 C S Gillespie and E Renshaw
where δ = µ + η. Whence it follows that the (marginal) population p.g.f. can be expressed as
R(z, 1; t) = exp(ξ0 + ξ1z + ξ2z2 + ξ3z3), (2.24)
where
ξ1 = (α3/µ)(1 − e−µt )3, ξ2 = (α3/2µ)(1 − e−µt )2(1 + 2 e−µt ),
ξ3 =
α3
3µ
(1 − e−3µt ) and ξ0 = −ξ3 − ξ2 − ξ1.
(2.25)
Similarly, the (marginal) counting p.g.f. takes the form
R(1, s; t) = exp(ζ0 + ζ1s + ζ2s2 + ζ3s3), (2.26)
where
ζ1 =
α3η
6δ4
[18µ(−µ + 2η)(1 − e−δt ) + 9η(−2µ + η)(1 − e−2δt )]
−
α3η3(1 − e−3δt )
δ4
+
3α3µ2ηt
δ3
,
ζ2 =
α3η2
6δ4
[18(η − 2µ)(1 − e−δt ) + 9(µ − 2η)(1 − e−2δt ) + 6η(1 − e−3δt )] + 3α3µη
2t
δ3
,
ζ3 =
α3η3
6δ4
[9 e−δt (2 − e−δt ) + 2 e−3δt − 11] + α3η
3t
δ3
, (2.27)
and ζ0 = −ζ3 − ζ2 − ζ1. On transforming the marginal p.g.f.s R(z, 1; t) and R(1, s; t) into
their cumulant generating function (c.g.f.) forms, namely ln[R(eθ , 1; t)] and ln[R(1, eφ; t)],
it immediately follows that the ith cumulants of the population and counting processes are
respectively given by
κi(t) = ξ1 + 2iξ2 + 3iξ3 and κci (t) = ζ1 + 2iζ2 + 3iζ3, (2.28)
highly suggestive of a saw-tooth structure with three teeth. Note that both of these cumulant
forms exhibit the characteristic of κi+1/κi → 3 as i → ∞, in contrast to the value 1 for the
Poisson process generated by the single-immigration-death process. So in the context of light
intensity the triple immigration-death process corresponds to non-classical, and not classical,
light (see the comment following (1.1)).
To obtain the associated marginal population and counting probabilities we need to
extract the coefficients of zn and sc in expressions (2.24) and (2.26), respectively. Inserting
the Hermite expansion (2.6) for the exponential of a quadratic function into the exponential
cubic representation (2.24), and paralleling the derivation for the paired-immigration process,
yields the marginal population size probabilities
pn(t) = e
ξ0
[n/3]
∑
r=0
ξ r3 (−ξ2)
(n−r)/2
r!(n − r)!
Hn−r (ω1). (2.29)
Note that it is the presence of [x] (i.e. the integer part of x) in the upper limit of the summation
that induces the triple effect. Specifically, on splitting pn(t) into p3n(t), p3n+1(t) and p3n+2(t),
we have for i = 0, 1, 2,
p3n+i (t) = e
ξ0
n
∑
r=0
ξ r3 (−ξ2)
3(n−r+i/3)/2
r![3(n − r + i/3)]!H3(n−r+i/3)(ω1). (2.30)
The strong degree of similarity between these three expressions means that on denoting
G(i)m (x1, x2, x3, z) ≡
m
∑
r=0
xr3(−x2)
3(m−r+i/3)/2
r![3(m − r + i/3)]!H3(m−r+i/3)/2(x1)z
3r+i , (2.31)
9

J. Phys. A: Math. Theor. 41 (2008) 355002 C S Gillespie and E Renshaw
they can be written in the much simpler, unified, form
p3n+i (t) = e
ξ0G(i)n (ω1, ξ2, ξ3, 1) where ω1 =
√
α3
√
−2µ
(1 − e−µt )2
√
1 + 2 e−µt
. (2.32)
As with the paired-immigration probabilities (2.8), these triple-immigration probabilities
exhibit a structure that corresponds directly to the number of immigrants that simultaneously
enter the process. Letting t → ∞ yields (relatively) simple expressions for the equilibrium
probabilities, namely
p3n+i (∞) = exp
(
−11α3
6µ
)
G(i)n
(
√
α3
−2µ
,
α3
2µ
,
α3
3µ
)
(2.33)
for i = 0, 1, 2 and n = 0, 1, 2, . . . . A parallel argument for the counting probabilities, this
time based on the p.g.f. (2.26), produces similar triple-effect characteristics.
To calculate the factorial population moments Mr(t) we first substitute z = 1 + z′ into
R(z, 1; t) to form
R(1 + z′, 1; t) = exp{ξ0 + ξ1(1 + z′) + ξ2(1 + z′)2 + ξ3(1 + z′)3}
= exp{(ξ1 + 2ξ2 + 3ξ3)z′ + (ξ2 + 3ξ3)(z′)2 + ξ3(z′)3}, (2.34)
where ξ0, . . . , ξ3 are defined by (2.25). Then on extracting the coefficient of (z′)n we obtain
for i = 0, 1, 2
M3r+i (t) =
r
∑
i=0
ξ i3(−ξ2 − 3ξ3)3(r−i+i/3)/2
i![3(r − i + i/3)]! H3(r−i+i/3)(ω2), (2.35)
where ω2 = (ξ1 + 2ξ2 + 3ξ3)/(2
√
−ξ2 − 3ξ3). As with the population probabilities (2.32), use
of the representation (2.31) yields the much simpler, unified, form,
M3r+i (t) = G
(i)
r (ω2, ξ2 + 3ξ3, ξ1 + 2ξ2 + 3ξ3, 1). (2.36)
Note that this result is a direct parallel of expression (2.32) for the p3n+i (t), and provides a
clear exposition of the triple effect in exactly the same way. The marginal counting factorial
moments
{
Mc3r+i (t)
}
can be calculated in similar manner, based on the counting p.g.f. (2.26)
for R(1, s; t) with s replaced by 1 + s ′.
Figure 3(a) shows the equilibrium population probabilities (2.33) for M1(∞) = 1, . . . , 6.
When the mean is close to unity, the gradient changes clearly illustrate a tripled effect,
which diminishes as M1(∞) increases; this effect parallels that for the paired process, where
an odd–even effect is witnessed only when the mean is close to 1 (Jakeman et al 1995).
Figure 3(b) shows the associated time-dependent p.d.f. at t = 0.1, 1 and 10 for n0 = 1; taking
α3 = 2 and µ = 1 gives an equilibrium mean of 6. Though the gradient changes exhibit a
strong triple effect at t = 0.1, this has already substantially diminished by t = 1, and has
(visually) disappeared by t  10. This swift decline is supported by figure 4(a), which shows
the ratio pi(t)/pi+1(t) for i = 0, 1, . . . , 20. This ratio not only demonstrates that the triple
effect disappears as time increases, but also that the effect itself holds throughout the range
of i-values. Note that a comparison of figures 3(b) and 4(a) shows that this probability ratio
reveals the presence of a triple effect much more strongly than a straight probability plot.
Further comparison with figure 4(b) shows that the power of the factorial population moment
ratio, Mi(t)/Mi+1(t), to highlight the triple effect lies inbetween these two. Here, as with the
probability ratio shown in figure 4(a), the tripled effect is (visually) present only when t is
small. A parallel analysis for the counting probability ratio pci (t)
/
pc+1i (t) (figure 5(a)), and
factorial counting moment ratio Mci (t)
/
Mci+1(t) (figure 5(b)), produces similar results, i.e. a
strong tripled effect is present only for small time values. This close degree of similarity once
again highlights the power of the counting process to capture a considerable amount of the
information generated by the underlying population process.
10

J. Phys. A: Math. Theor. 41 (2008) 355002 C S Gillespie and E Renshaw
and ρc(c1, c2) ≡ ρ(c1; (t1, t1 + T ), c2; (t2, t2 + T )) is evaluated as per expression (2.18). To
examine the equilibrium correlation structure we set t2 = t1 + τ + T , for τ  0, and allow
t1 → ∞, whence (2.39) yields
ρc(c1, c2) =
η e−δτ (1 − e−δT )2
T (2δ2 + η2 − µ2) − 2η(1 − e−δT )
. (2.41)
Thus as the distance, τ , separating the two non-overlapping time intervals (t1, t1 + T ) and
(t2, t2 + T ) increases, the correlation ρc(c1, c2) reduces as e−δτ , as also happens with the
pairs autocorrelation function (2.21). However, whilst both correlation functions have
an e−δτ term, the correlation in the triples process is always greater for T > 0, since
2δ2 + η2 − µ2 < δ2 + η2 + ηµ.
3. Birth–death process
Although the p.g.f. (2.4) of Gillespie and Renshaw (2005) provides the general ‘solution’ to
the (additive) batch immigration-death-counting process, we have shown that detailed analysis
of specific cases (here double and triple immigration) yields far more transparent solutions
that illuminate the saw-tooth behaviour of the underlying probability and moment structures.
Development of specific results for quadruple, etc., immigration could well be rewarding,
but let us now explore the effect of introducing birth, in order to determine whether the
counting process can still provide a high level of information content on the development of
an unobservable (multiplicative) population process.
The simplest and most transparent way of modelling this construct is to employ the
simple birth–death (BD) process (Cox and Miller 1965, Bartlett 1966). For this is one of
the few multiplicative continuous-time processes for which the probability structure can be
expressed in closed form, and suchmathematical tractability suggests that the parallel counting
probabilities will not be too opaque to prevent the extraction of process structure. Moreover, on
including immigration, the resulting birth-immigration-death (BDI) process not only provides
an immediate extension to our earlier immigration-death counting model, but it is also directly
associated with a kind of behaviour commonly encountered in optics (Jakeman et al 2003).
For when the birth and immigration rates are equal, thermal or Bose–Einstein statistics are
predicted for the single-fold population distribution, which is a simple geometric progression.
This feature is found to characterize the photon statistics of a laser below threshold (Shimoda
et al 1957) and also Gaussian speckle noise generated when laser light is scattered by particles
or rough surfaces (Berlolotti 1974). In the case of the laser model, spontaneous and stimulated
emission within a population of photons in a cavity is analogous to immigration and birth,
respectively, and absorption is the analogue of death. Now a key point of Jakeman et al
(2003) is that, in equilibrium, the MID process results in a negative binomial probability
structure that closely resembles that of the BDI process. So their main objective is to devise
and optimize a strategy for distinguishing between these two scenarios when measurement is
restricted to the external monitoring scheme. Here, we shall concentrate on developing the
full time-dependent probability structure for the simpler BD process; that for the BDI model
involves straightforward, albeit more tedious, algebraic development, and in this sense adds
little to our overall understanding.
For a population of size n at time t, in the small time interval (t, t + h) there is a
probability λnh + o(h) that a particular individual gives birth, a probability µnh + o(h) that
it dies unobserved, and a probability ηnh + o(h) that it dies and is counted. So the forward
Kolmogorov equation for this population counting process is
dpnc/dt = λ(n − 1)pn−1,c + µ(n + 1)pn+1,c + η(n + 1)pn+1,c−1 − n(λ + µ + η)pnc, (3.1)
13

J. Phys. A: Math. Theor. 41 (2008) 355002 C S Gillespie and E Renshaw
where t  0 and pn,−c(t) = 0 for c = 1, 2, . . . . This can be solved in terms of the joint
generating function (2.2) to give
R(z, s; t) =
(
θ2z eλ(θ1−θ2)t − θ1θ2 eλ(θ1−θ2)t − θ1z + θ1θ2
z eλ(θ1−θ2)t − θ1 eλ(θ1−θ2)t − z + θ2
)n0
, (3.2)
where
θ1, θ2 =
[
λ + η + µ ∓
√
(λ + η + µ)2 − 4λ(µ + ηs)
]/
(2λ). (3.3)
Note that on setting s = 1 in (3.2) we recover the well-known generating function for the
simple birth–death population process.
3.1. Counting mean and variance
On forming [∂R(1, s; t)/∂s]s=1 it immediately follows that the counting mean
κc1(t) =
{
ηn0(e(λ−µ−η)t − 1)/(λ − µ − η) for λ = µ + η
ηn0t for λ = µ + η.
(3.4)
Whence allowing t → ∞ yields
lim
t→∞
κc1(t) =
{
n0η/(µ + η − λ) for λ < µ + η
∞ for λ  µ + η. (3.5)
So when the population explodes, i.e. λ  µ + η, the counts explode. Conversely, when the
population is certain to become extinct, i.e. λ < µ+η, the counting mean approaches the limit
n0η/(µ + η − λ). Similarly, forming [∂2R(1, s; t)/∂s2]s=1 leads to the counting variance
κc2(t) = n0ηt[3(1 − tη) + 2t2λη]/3
for λ = η + µ, and
κc2(t) = n0η{e
(λ−η−µ)t [(λ − µ)2 − η2 − 4ληt (λ − η − µ)]
+ e2(λ−η−µ)tη(λ + η + µ) − (λ − µ)2 − η(λ + µ)}/(λ − η − µ)3 (3.6)
for λ = η + µ. Whence for large t = T , we see that
κc2(T )  n0η
2 e2(λ−η−µ)T (λ + η + µ)/(λ − η − µ)3 for λ > η + µ
κc2(T )  2n0T
3λη2/3 for λ = η + µ
κc2(T ) = n0η[(λ − µ)2 + η(λ + µ)]/(−λ + η + µ)3 for λ < η + µ.
(3.7)
Thus unlike the population variance which tends to zero as t → ∞ for λ < µ+η, the counting
variance reaches a finite limit. For ultimate extinction is certain to occur, at which point there
are no longer any individuals available for counting.
3.2. Probabilities
To derive the (marginal) counting probability pcc(t) we need to extract the coefficient of sc
from the generating function Q(1, s; t), and to effect this we first need to construct a general
framework. On noting that the modified Bessel function of the first kind is defined (A&S
section 9.6.10) as
Ix(ξ) =
(
ξ
2
)x ∞
∑
k=0
(ξ/2)2k
k!(x + k + 1)
,
14

J. Phys. A: Math. Theor. 41 (2008) 355002 C S Gillespie and E Renshaw
we can write down two particularly useful forms, namely
Ix−1/2(ξ) =
2x+1/2
ξx+1/2
x!
√
π
∞
∑
k=x
(
k
x
)
ξ 2k
(2k)!
I1/2−x(ξ) =
2x+1/2
ξx+1/2
x!
√
π
∞
∑
k=0
(
k + 1/2
x
)
ξ 2k+1
(2k + 1)!
.
(3.8)
It then follows that
eξ2
√
1+s/ξ1
=
∞
∑
n=0
ξn2 (1 + s/ξ1)n/2
n!
=
∞
∑
n=0
ξn2
n!
∞
∑
i=0
(
n/2
i
)
si
ξ i1
=
∞
∑
i=0
si
ξ i1
(
∞
∑
n=i
(
n
i
)
ξ 2n2
(2n)!
+
∞
∑
n=0
(
n + 1/2
i
)
ξ 2n+12
(2n + 1)!
)
=
∞
∑
n=0
ζn(ξ2)s
n (3.9)
where
ζn(ξ2) =
2−1/2−n
√
πξ 1/2+n2 [In−1/2(ξ2) + I1/2−n(ξ2)]
ξn1 n!
. (3.10)
Moreover, a similar expansion leads to
[(1 + eξ2
√
1+s/ξ1)(1 + ξ3
√
1 + s/ξ1) − 2]−1 =
∞
∑
c=0
ϕcs
c, (3.11)
where
ϕc = −
1
2
∞
∑
n=0
1
2n
⎧
⎨
⎩
n
∑
i=0
(
n
i
) n
∑
j=0
(
n
j
)
ξ j3
c
∑
k=0
(
j/2
c − k
)
ζk(iξ2)
ξ c−k1
⎫
⎬
⎭
. (3.12)
Whence on observing that the denominator of the p.g.f. Q(1, s; t) takes the form
ˆQ(1, s; t) =
ξ4
(1 + eξ2
√
1+s/ξ1)(1 + ξ3
√
1 + s/ξ1) − 2
, (3.13)
where ξ1 = [4λµ − (λ + η + µ)2]/[4λη], ξ2 = −t
√
(λ + η + µ)2 − 4λµ, ξ3 = −ξ2ξ4/(2λt)
and ξ4 = 2λ/(λ − µ − η), it follows that the general coefficient of sc in the denominator of
Q(1, s; t) is ξ4ϕc where ϕc is given by (3.12). Finally, on defining ξ5 =
√
(λ + η + µ)2 − 4λµ,
the numerator of Q(1, s; t) can be expressed as
θ2 e
λ(θ1−θ2)t
− θ1θ2 e
λ(θ1−θ2)t
− θ1 + θ1θ2 ≡
∞
∑
c=0
ϑcs
c
=
1
2λ
[eξ2(λ + η + ξ5 − µ) + ξ5 + µ − λ − η]s0
+
η
λξ5
(eξ2λt (λ + η + ξ5 − µ) − λ(1 + eξ2) + ξ5(1 − eξ2))s1
+
∞
∑
n=2
{
ζn(ξ2)(λ + µ + η)(1 − 2µ) − ζn−1(ξ2)2η
− ξ5
(
1/2
n
)
1
ξn1
+
n
∑
i=0
ζi(ξ2)
(
1/2
n − i
)
ξ5
ξn−i1
}
sn
2λ
. (3.14)
15

J. Phys. A: Math. Theor. 41 (2008) 355002 C S Gillespie and E Renshaw
where for i = 1, 2,
Var[(ti , ti + T )] =
n0η eζ1ti [η eζ1ti (λ + η + µ)(1 − eζ1T )2 − eζ1T 4ηλζ1T ]
(λ − η − µ)3
−
n0η eζ1ti (ζ 21 + 4ηλ)(1 − eζ1T )
(λ − η − µ)3
. (3.23)
Figure 7(b) shows (3.22) for t1 = 0, t2 = T , n0 = 1 and µ = 0. The most surprising feature
of this a.c.f. is the presence of negative correlation for 0
T
2, which contrasts with the
population a.c.f. where correlation is always non-negative. The explanation for this is that if
a count is observed in the time period (0, T ) for small T, then this is likely to have resulted
in extinction. Hence, no counts will be observed in the subsequent time interval (t2, t2 + T ).
Note that, as before, the correlation is one or zero for λ larger or smaller than µ.
On setting t2 = t1 + T + τ in (3.22), for τ > 0, and letting T → ∞, it follows that
ρ[c1; (t1, t1 + T ), c2; (t2, t2 + T )] → 1 for λ > µ + η, 3/4 for λ = µ + η and 0 for λ < µ + η.
So we obtain the intuitively reasonable result that when λ is greater or less than µ, the counting
correlation is either one or zero. A surprising feature, however, occurs when λ = µ, for this
results in a correlation of 3/4. This anomaly occurs because although extinction is certain,
the event itself can take an infinitely long time to happen (see figure 7(b)).
4. Summary and discussion
The construction of models is the most basic activity in science, lying at the heart of
understanding, prediction and practical applications, and uncovering the dominantmechanisms
governing the stochastic evolution of dynamical systems in an essential component of this
activity (Jakeman et al 2003). Now in many situations of practical interest, it will not be
possible to make a direct measurement of the population itself, but only of some external
manifestation of its evolution. Here we assume that it is only possible to monitor the rate of
evolution of individuals from the population. So in quantum optics we might count photon
emissions, in biological populations emigrants from the system, and in epidemiology people
either dying, being isolated, or exhibiting specific symptoms.
The time-dependent paired-immigration-death process, which has the merit of being both
easy to formulate and is analytically tractable, produces some neat, and highly illuminating
features, especially the strikingly clear odd–even population effects in both the probability and
moment structures. That these features are similarly present in the counting probabilities and
moments is witness to the strength of the information in the counting process in allowing us to
make inferences on the (unseen) population process. We see that these effects are, nevertheless,
highly susceptible to the presence of single immigrants, only a small infusion of which are
necessary to mask the odd–even phenomenon. A key question relates to the correlation
between counts in non-overlapping time intervals, and whilst the algebraic development is
non-trivial it is perfectly possible to develop relations that are sufficiently transparent to
expose the underlying correlation features. Extending single- and paired-immigration to k-
immigration produces intriguing mathematical results; here we just consider k = 3 in order
to demonstrate that both the approach and the associated conclusions carry over neatly to the
multiple-immigration scenario.
Even under the generalized k-immigration scheme, once individuals enter the population
‘black-box’ they cannot produce yet more individuals, so there is a direct balance between
them and those who are subsequently either counted or die. However, once we introduce
birth the situation changes dramatically. For individuals can be created within the black-
box that are totally unseen from outside, and so the extent to which the population statistics
18

J. Phys. A: Math. Theor. 41 (2008) 355002 C S Gillespie and E Renshaw
can be gleaned from the counting statistics is far less intuitive. Indeed, being multiplicative
by nature, this process has the potential for offering a much greater theoretical challenge.
Although the counting probability structure is algebraically more awkward to derive, in
essence it nevertheless just involves a basic series expansion, which is itself sufficiently
tractable to allow basic inferences to be drawn without further calculation. However, unlike
the immigration processes, we now have to distinguish between the cases λ
µ+η, for which
ultimate extinction is certain, and λ > µ + η, in which case both extinction and unlimited
exponential population/count growth are possible. Intriguingly, the correlation coefficient
(3.22) exhibits an anomalous departure from the values 0 (λ < µ + η) and 1 (λ > µ + η) at
λ = µ + η, since in this case it equals 3/4.
The approach clearly offers considerable future potential for expanding intomore complex
situations. Not only is the analysis of multi-type processes, such as predator–prey or general
epidemic models, a potentially huge field covering applied situations of immediate practical
and intrinsic use, but the associated inferential problems of determining an unknown stochastic
process from known counts offers rich theoretical rewards. Overall the problem of determining
what information is contained within the counting statistics of a stochastic process is complex.
Future research should concentrate on determining methods and techniques that would enable
us to state explicitly the information that the observed process could tell us about the underlying
latent system. Current methods of analysis are based almost exclusively on numerical
techniques centred on Markov chain Monte Carlo procedures (see, for example, Boys
et al 2008 and Golightly and Wilkinson 2008). So the development of parallel theoretical
results would be both novel and timely.
Acknowledgments
We wish to thank Professor Eric Jakeman, Dr Keith Hopcraft and Jonathan Matthews for
valuable discussion. The work was supported by an EPSRC Research Studentship. We also
thank the anonymous referees for their insightful comments that helped to strengthen this
paper.
References
Abramowitz M and Stegun I A 1970 Handbook of Mathematical Functions 9th edn (New York: Dover)
Bartlett M S 1966 An Introduction to Stochastic Processes (Cambridge: Cambridge University Press)
Berlolotti M 1974 Photon statistics Photon Correlation and Light Beating Spectroscopy ed H Z Cummins and
E R Pike (New York: Plenum) pp 41–74
Boys R J, Wilkinson D J and Kirkwood T B L 2008 Bayesian inference for a discretely observed stochastic kinetic
model Stat. Comput. 18 125–35
Burnham D C and Weinberg D L 1970 Observation of simultaneity in parametric production of optical photon pairs
Phys. Rev. Lett. 25 84–7
Cox D R and Lewis P A W 1966 The Statistical Analysis of Series of Events (London: Methuen)
Cox D R and Miller H D 1965 The Theory of Stochastic Processes (London: Methuen)
Gillespie C S and Renshaw E 2005 The evolution of a batch-immigration death process subject to counts Proc. R.
Soc. A 461 1563–81
Golightly A and Wilkinson D J 2008 Bayesian inference for nonlinear multivariate diffusion models observed with
error Comput. Stat. Data Anal. 52 1674–93
Jakeman E, Hopcraft K I and Matthews J O 2003 Distinguishing population processes by external monitoring Proc.
R. Soc. A 459 623–39
Jakeman E, Phayre S and Renshaw E 1995 The evolution and measurement of a population of pairs J. Appl.
Prob. 32 1048–62
Jakeman E and Shepherd T J 1984 Population statistics and the counting process J. Phys. A: Math. Gen. 17 L745–50
Kingman J F C 1993 Poisson Processes (Oxford: Oxford University Press)
19

J. Phys. A: Math. Theor. 41 (2008) 355002 C S Gillespie and E Renshaw
Louden R and Knight P L 1987 Squeezed light J. Mod. Opt. 34 709–59
Mandel L 1959 Fluctuations of photon beams: the distribution of the photo-electrons Proc. Phys. Soc. 74 233–43
Matthews J O, Hopcraft K I and Jakeman E 2003 Generation and monitoring of discrete stable random processes
using multiple immigration population models J. Phys. A: Math. Gen. 36 11585–603
Saleh B E A 1978 Photoelectron Statistics (Berlin: Springer)
Schoenberg I J 1983 Problem 650 (Nieuw Archiff Vor Wiskunde, Series 4 vol 1) pp 377–8
Shepherd T J 1984 Photoelectron counting—semiclassical and population monitoring approaches Opt. Acta 31
1399–407
Shepherd T J and Jakeman E 1987 Statistical analysis of an incoherently coupled, steady state optical amplifier
J. Opt. Soc. Am. B 4 1860–9
Shimoda K, Takahasi H and Townes C H 1957 Fluctuations in amplification of quanta with application to maser
amplifiers J. Phys. Soc. Japan 12 686–700
Srinivasan S K 1988 Point Process Models of Cavity Radiation and Detection (Oxford: Oxford University Press)
Stoyanov J M 1988 Counterexamples in Probability (New York: Wiley)
20