Functional-integral based perturbation theory for the Malthus-Verhulst process

Abstract

We apply a functional-integral formalism for Markovian birth and death processes to determine asymptotic corrections to mean-field theory in the Malthus-Verhulst process (MVP). Expanding about the stationary mean-field solution, we identify an expansion parameter that is small in the limit of large mean population, and derive a diagrammatic expansion in powers of this parameter. The series is evaluated to fifth order using computational enumeration of diagrams. Although the MVP has no stationary state, we obtain good agreement with the associated quasi-stationary values for the moments of the population size, provided the mean population size is not small. We compare our results with those of van Kampen's omega-expansion, and apply our method to the MVP with input, for which a stationary state does exist. We also devise a modified Fokker-Planck approach for this case.

Birth-and-death process; Master equation; Path integral; Omega-expansion; Doi formalism

REGULAR ARTICLES

Functional-integral based perturbation theory for the Malthus-Verhulst process

Nicholas R. Moloney^I,II; Ronald Dickman^II

^IInstitute of Theoretical Physics, Eötvös University, Pázmány Péter sétány 1/A, 1117 Budapest, Hungary

^IIDepartamento de Física, ICEx, Universidade Federal de Minas Gerais, 30123-970 Belo Horizonte - Minas Gerais, Brasil

ABSTRACT

We apply a functional-integral formalism for Markovian birth and death processes to determine asymptotic corrections to mean-field theory in the Malthus-Verhulst process (MVP). Expanding about the stationary mean-field solution, we identify an expansion parameter that is small in the limit of large mean population, and derive a diagrammatic expansion in powers of this parameter. The series is evaluated to fifth order using computational enumeration of diagrams. Although the MVP has no stationary state, we obtain good agreement with the associated quasi-stationary values for the moments of the population size, provided the mean population size is not small. We compare our results with those of van Kampen's W-expansion, and apply our method to the MVP with input, for which a stationary state does exist. We also devise a modified Fokker-Planck approach for this case.

Keywords: Birth-and-death process; Master equation; Path integral; Omega-expansion; Doi formalism

I. INTRODUCTION

The need to analyze Markov processes described by a master equation arises frequently in physics and related fields [1,2]. Since such equations do not in general admit an exact solution, approximation methods are of interest. A widely applied approximation scheme is van Kampen's 'W-expansion', which furnishes corrections to the (deterministic) mean-field or macroscopic description in the limit of large effective system size [1]. Another approximation method, based on a path-integral representation for birth-and-death type processes, was proposed by Doi [3] and then by Grassberger and Scheunert [4]. Later Peliti [6] and Goldenfeld [5] revived it. Renewed interest in this type of representation has been stimulated by Cardy and coworkers [7,8]. This approach was recently reviewed and extended [9], and applied to derive a series expansion for the activity in a stochastic sandpile [10], and to study metastability in the contact process [11].

It should be noted that while effectively exact results can be obtained via numerical analysis of the master equation, the calculations become extremely cumbersome for large populations or multivariate processes. A further limitation of numerical analyses is that they do not furnish algebraic expressions that may be required in theoretical developments. For these reasons it is highly desirable to study approximation methods for stochastic processes.

In the present work we apply the path-integral based perturbation approach to a simpler problem, namely, the Malthus-Verhulst process (MVP), a birth-and-death process in which unlimited population growth is prevented by a saturation effect. (The death rate per individual grows linearly with population size.) This is an important, though highly simplified model in population dynamics. Although the master equation for this process is readily solved numerically, the model serves as a useful testing ground for approximation methods. A lattice of coupled MVPs exhibits (in the infinite-size limit) a phase transition belonging to the directed percolation universality class.

In the perturbation approach developed here [6,9], moments of the population size n are expressed as functional integrals over a pair of functions, y(t) and (t), involving an effective action. The latter, obtained from an exact mapping of the original Markov process, generally includes a part that is bilinear in the functions y(t) and (t), whose moments can be determined exactly, and 'interaction' terms of higher order, that are treated in an approximate manner. In the present approach, the interaction terms are analyzed in a perturbative fashion, leading to a diagrammatic series. With increasing order, the number of diagrams grows explosively, so that it becomes convenient to devise a computational algorithm for their enumeration and evaluation. Elaboration of such an algorithm does not, however, require any very sophisticated techniques, and could in fact be applied to a variety of problems. This is the approach that was applied to the stochastic sandpile in Ref. [10]. In the latter case, the evaluation of diagrams involves calculating multidimensional wave-vector integrals. The present example is free of this complication, allowing us to derive a slightly longer series than for the sandpile.

In this work we focus on stationary moments of the MVP. The series expressions are apparently divergent, but nevertheless provide nearly perfect predictions away from the small-population regime, as compared with direct numerical evaluation of quasi-stationary properties. One might suppose that the divergent nature of the perturbation series is due to the MVP not possessing a true stationary state. (The process must eventually become trapped in the absorbing state, although the lifetime grows exponentially with the mean population [12].) Applying our method to the MVP with a steady input, which does possess a stationary state, we find however that the perturbation series continues to be divergent, although again providing excellent predictions over most of parameter space.

In the following section we define the Malthus-Verhlst process, explain the perturbation method and report the series coefficients for the first four moments, up to fifth order in the expansion parameter. In Sec. III we briefly compare these results to those of the W-expansion. Then in Sec. IV we present numerical comparisons of our method (and of the W-expansion) against quasi-stationary properties. We apply our method to the MVP with input in Sec. V, and also discuss an approximation based on the Fokker-Planck equation. We summarize our findings in Sec. VI.

II. PERTURBATION THEORY FOR THE MALTHUS-VERHULST PROCESS

Consider the Malthus-Verhulst process (MVP) n(t), in which each individual has a rate l to reproduce, and a rate of µ + n(n – 1) to die, if the total population is n. By an appropriate choice of time scale we can eliminate one of these parameters; we choose to set µ = 1. Then in what follows we use a dimensionless time variable t' = µt and dimensionless rates l' = l/µ and n' = n/µ. From here on we drop the primes.

The mean-field or rate equation description of the process is

where x º án(t)ñ, with the nontrivial stationary solution = (l – 1)/n for l > 1. We are interested in deriving systematic corrections to this result, and in calculating higher moments of the process.

Our starting point is the expression for the r-th factorial moment of a general process taking non-negative integer values,

where, for simplicity, we have assumed an initial Poisson distribution with parameter p (see Eq. (108) of Ref. [9]), so that the probability generating function at time zero is F₀(z) = e^p(z–1). Here, the kernel is given by the functional integral,

(equivalent to Eq. (106) of [9]), with

containing the bilinear part of the action. In the case of the MVP, the "interaction" part is,

(see Eq. (54) of [9]).

Our goal is to expand each factorial moment about its mean-field value. To this end, consider the shift of variable,

where is a real constant. On performing this shift and letting = (l – 1)/n, which eliminates the term µf⁰, the argument of the exponential (i.e., of the factors [y,]_{z = 1} and exp[-S_I]) in Eq. (3) becomes,

where we have introduced w º l – 1 (equal to -w as defined in [9]). We recognize the exponential of z plus the first term in the integrand as [,f]_{z = 1}; the remaining terms then represent -, the new effective interaction, which will be treated perturbatively. The four terms of - are represented graphically, following the conventions of Ref. [9], in Fig. 1.

We refer to these vertices as source, bifurcation, conjunction and 4-vertex, respectively.

A. Perturbation expansion

The analysis of the MVP now follows the lines of [9]. Let

denote the free expectation of any function of f and . From the discussion of Sec. 4 Ref. [9] we have,

(here we used f(0) = p – , and have already canceled the factor e^z in with the corresponding factor in the inital generating function, F₀(z)),

and the propagator,

Consider now the series for the mean population size án(t)ñ. The zeroth-order term is simply

which converges to the mean-field result as t ® ¥. Similarly, the leading term in án^r (t)ñ_f is ^r, so that the stationary probability distribution, in mean-field approximation, is Poisson with parameter .

Corrections involving are conveniently represented as diagrams with all lines leaving vertices contracted with ingoing lines, either at vertices to the left, or at the "sink" lying to the left of all vertices. (In the calculation of án^r (t)ñ_f the sink is a point with r incoming lines.) The lowest-order diagram in án(t)ñ is the left-most diagram shown in Fig. 2.

This diagram makes the following contribution to án(t)ñ:

The four diagrams (in the series for án(t)ñ) having two vertices are also shown in Fig. 2. Only the first of these four makes a nonzero contribution to án_¥ñ º lim_{t ® ¥}án(t)ñ, namely,

(The combinatorial factor 2 represents the number of ways the lines may be contracted between source and bifurcation.)

If we are only interested in stationary properties, it is convenient to use the Laplace transform. Let f_D be the n-fold integral over time variables in a given n-vertex diagram D. Noting that time-dependent factors are associated with each propagator and each uncontracted incoming line, we see that f_D is of the form:

where b_i is w times the number of uncontracted lines incident on vertex i. (In the first graph of Fig. 2, b₁ = 2w.) The factors a_i are given by w times the number of lines (propagators) running between vertices i and i – 1 (with i = 0 representing the sink), regardless of where these lines originate or terminate. Now consider

Using Eq. (15) we may write

The contribution of diagram D to án_¥ñ is proportional to

by the Final Value Theorem of Laplace transforms. Since all of the a_i are nonzero, _D is zero unless b_i = 0, "i = 1,...,n. In the latter case,

Thus the only diagrams contributing to án_¥ñ (and by extension, to áñ) are those free of uncontracted lines. This is in fact evident on physical grounds: each such line carries a factor p – , whereas stationary properties cannot depend on the initial mean population p. In diagrams contributing to án_¥ñ, the first vertex (i.e., immediately to the right of the sink) must be a conjunction, while the n-the vertex must be a source.

B. Diagrammatic analysis for the stationary mean population

The contribution of a diagram D to án_¥ñ is the product of three factors: _D discussed above; a combinatorial factor counting the number of contractions consistent with the diagram topology; the product of vertex-associated factors shown in Fig. 1.

Certain infinite classes of diagrams may be summed up exactly. Consider the sequence shown in Fig. 3: between the source and conjunction we insert any number of 4-vertices or bifurcation-conjunction pairs.

A short calculation shows that the contribution of this series is given by:

In this way, by multiplying the n = 2 contribution, -, by k º (1 + n/w²)^-1, we have included all diagrams of the form of Fig. 3, i.e., diagrams in which all lines exiting vertex i are contracted on vertex i – 1, "i = 1,...,n. With this simple replacement such reducible diagrams are no longer to be included explicitly.

A similar observation allows us to add arbitrary sequences of 4-vertices or bifurcation-conjunction pairs to any diagram of four or more vertices. (Diagrams outside the series depicted in Fig. 3 have n > 4 vertices.) The diagram consists of a 'body' containing vertices 2,...,n – 1 linked to a conjunction (vertex 1) and a source (vertex n). Call the factor associated with the body A, so that the diagram makes a contribution of -nA to án_¥ñ. A family of diagrams may be constructed by adding sequences, as above, between vertex n and the body; the sum of these contributions is

By the same reasoning we may insert arbitrary sequences between the first vertex and the body, yielding the same set of contributions. As a result, we may multiply the contribution of any diagram (not included in the sequence of Fig. 3) by k², and thereby automatically include all diagrams with the same body but arbitrary linear sequences between the first vertex and the body, and between vertex n and the body. Such diagrams are no longer included explicitly. This procedure will be called "dressing" the conjunction (vertex 1) and the source (vertex n).

At this stage we have essentially exhausted the simplifications (via "dressing" or summing sets of trivial variations to a diagram) that can be realized in the Laplace transform representation. Summarizing, the contributions to án_¥ñ are as follows.

1. The mean-field contribution .

2. The set of diagrams shown in Fig. 3, giving -k/w.

3. Diagrams of four or more vertices, with no dangling lines, and subject to the restrictions mentioned above. Each such contribution is to be multiplied by k².

Including diagrams of up to four vertices we find:

where the higher order terms come from diagrams of five or more vertices. This result suggests that we adopt

as the expansion parameter. Up to diagrams of four vertices we have

The first correction -e µ 1/l, as l ® ¥.

Note that n µ e to lowest order. Thus, qualitatively, the expansion parameter e may also be thought of as describing overcrowding/competition.The lowest order in e at which a given diagram contributes to án_¥ñ, when expressed in the form of Eq. (24), may be found as follows. First observe that the diagram carries a factor n^c+f–s+1 where c, f and s represent the number of conjunctions, four-vertices and sources, respectively, and the contribution of 1 in the exponent is due to the prefactor . Certain constraints exist between the numbers of vertices. Equating the total number of lines emanating from vertices with the total number entering, we find

where b is the number of bifurcations and the 1 on the r.h.s. represents the sink. Thus c = 2s + b – 1. On the other hand the total number of vertices is n = s + b + f + c, and using this we find the power of -n to be c + f – s + 1 = n – c º m.

A simple analysis yields the algebraic factor associated with a given diagram having n > 4 vertices:

This is readily expressed in terms of the parameter e as:

Since we intend to organize the expansion in powers of e, it is of interest to know which values of n correspond to a given m. We begin by noting that in diagrams of four or more vertices, the second vertex (from the left) must be a conjunction, just as the first is. (If it were a bifurcation or a four-vertex the result would be a diagram already included by dressing the first vertex.) Thus c > 2, and there are in fact diagrams with just two conjunctions, for any n > 4. We therefore have m < n – 2 or n > m + 2. To find an upper bound on n, we find an upper bound on c. Recall that c = 2s + b – 1. To maximize c we let b = f = 0 (a diagram consisting of only sources and conjunctions), in which case c + s = n, yielding m = n – c = s = (n + 1)/3, or n < 3m – 1. Thus for m = 2 we have contributions from diagrams of 4 and 5 vertices, while for m = 5 the number of vertices ranges from 7 to 14.

In order to evaluate the contributions from diagrams of five or more vertices we have developed an enumeration code. Reducible diagrams are eliminated by imposing the following rules:

1) The second vertex, like the first, is a conjunction.

2) Vertex n – 1 cannot be a 4-vertex.

3) If vertex n – 1 is a conjunction, vertex n – 2 must be a source.

Using the enumeration code we are able to evaluate contributions up to O(e⁵). At this order there are approximately 8 × 10⁸ diagrams. (There are 9, 1317 and 594 339 diagrams at orders 2, 3 and 4, respectively.) This explosive growth in the number of terms prevents our going to higher order. The fifth-order calculation requires about two days on a fast PC. The series in e begins

where u º 2 – l. Higher order terms may be found in the ^Appendix.

A particularly simple case is l = 2, for which we have

The ratios of successive coefficients in the series are: 1, 5, 8.2, 11.829 and 15.346, suggesting unlimited growth and therefore a divergent series. Numerical results (see below) nevertheless reveal excellent agreement with quasi-stationary properties for l sufficiently large. (The quasi-stationary probability distribution q_n is the limiting t ® ¥ probability of state n, given that the system has not visited the absorbing state up to time t [12].)

C. Higher moments

We turn now to the expansion of higher stationary factorial moments. From Eqs. (3) and (6) we have

Expanding the product, we have first the mean-field contribution

^r, and then a series of terms involving diagrams. The term µ involves (for q < r) diagrams that already appeared in the series for the q-th factorial moment. Thus for r = 2 we have

Note that is simply ₁, the sum of all diagrams with a single line entering the sink, that is, the diagrammatic series for ánñ. Consider now the final term in Eq. (31), the series ₂ of diagrams with two lines entering the sink. Recalling that, in all diagrams in ₁, the first (leftmost) vertex is a conjunction, we see that there is a one-to-one correspondence between ₁ and ₂. For a given diagram, making a contribution of A to ₁, there is a corresponding diagram (with the conjunction and one-line sink replaced by a two-line sink) that contributes -wA/n = -A to ₂, that is, ₂ = -

₁. Thus we have

where in the second line we used ánñ = + ₁. Using Eqs. (65) and (32), we find that

For a Poisson distribution the difference is zero; here it approaches 1/n as l ® ¥.

For r > 3 we do not have a simple relation between D_r and D₁, so the diagrammatic series must be evaluated to the desired order. In general we have

where

_j = , with ₀ = 1. (Note that for r = 3 the j = 1 and j = 2 terms cancel.) As in the case of ánñ, we include an overall factor of ^r, and write

Note that for r > 3 we can dress the leftmost source, as before, but that it is no longer possible to dress the sink. Thus the algebraic factor associated with an arbitrary diagram in (n/w)^r

_r is

An analysis along the lines presented in the r = 1 case leads to m = n – c (as before) for the order in n, and to s – n – r = f, so that

The limits on n, for fixed order m, are m < n < 3m – r. (Note that n = m is possible for r > 3 because there are diagrams with no conjunctions.)

As in the case r = 1, we have constructed an enumeration code to evaluate the corrections due to diagrams. The expansion may be found in the ^Appendix. Numerical comparisons of these series against quasi-stationary properties of the MVP will be discussed in Sec. 4.

III.COMPARISON WITH THE W-EXPANSION

A well known method for obtaining approximate solutions to the master equation is van Kampen's W-expansion [1]. The method depends on the system size (or the expected value of the stochastic variable of interest) being large, so that fluctuations are small compared to the value (t) predicted by the macroscopic equation. The stochastic variable n is then written in the form

where the first, deterministic, term represents the solution to the macroscopic equation, with W denoting the size of the system. The stochastic contribution, represented by the second term, is expected on general grounds to scale as W^1/2. As shown in Ref. [1], z(t) satisfies the macroscopic equation, while the probability density P(x,t) satisfies, to lowest order in W^-1/2, a linear Fokker-Planck equation. In the present case the macroscopic equation is

i.e., the Malthus-Verhulst equation, with nontrivial stationary solution z = (l – 1)/n = . The equation governing the evolution of P is, in the present instance,

This implies that in the stationary state, x is Gaussian with mean zero and variance l/n, so that, setting the formal expansion parameter W to unity in Eq. (38), n is Gaussian with mean and variance l/n. In the perturbation expansion developed in the preceding section, n is, to zeroth order, a Poissonian random variable with mean . Conceptually, a Poisson distribution seems preferable to a Gaussian as the reference distribution (since n is discrete and cannot assume negative values), but in practical terms we expect the difference to be very small.

To first order in e, we find ánñ = (1 – e) and var(n) = (1 + e). Noting that

we see that the difference from the W-expansion result, to first order in W^-1/2, is small for l 1. In the W-expansion, corrections to the moments áx^rñ can be obtained via the procedure detailed in Ref. [1]. In the present case one finds

so that

The series prediction is

so that the two methods agree to first order in e. Extending the W-expansion to include terms of order W^-1 in x, one finds

These results are compared against the e series in the following section.

IV. NUMERICAL COMPARISONS

In this section we compare the e-series predictions with exact (numerical) results for quasi-stationary (QS) properties of the MVP. The latter are obtained via recursion relations leading to the QS distribution as detailed in [12]. QS properties are those obtaining at arbitrarily long times, conditioned on survival (i.e., the process has never visited the absorbing state). Although a condition on survival is not involved in the perturbative analysis of Sec. 2, it seems reasonable to compare its predictions with QS properties, since the true stationary properties are the trivial ones of the absorbing state (population zero, no fluctuations). We note that, as suggested above, the series in e appears to be divergent for any set of parameter values, as reflected in the ratios between coefficients of successive terms, whose ratios appear to grow without limit. For suitably large values of l the ratios are very small, although increasing with order m, so that the series, truncated at fifth order, appears to have "converged". The numerical evidence (limited to m < 5, of course) points nevertheless to a divergent series.

In Fig. 4 we compare the QS mean population size with án_¥ñ as predicted by the series, as a function of l, with n fixed at 0.01. For l less than about 1.2, the series is useless (it yields a negative population size!), and is inferior to mean-field theory. For l > 1.4 on the other hand, we find good agreement with the QS result. For a more detailed comparison we plot, in the inset of Fig. 4, the difference Dn between ánñ (as given by the QS distribution and the e series) and the mean-field prediction ; there is perfect agreement for l > 1.5. (The correction to the mean population size decays µ 1/l for large l. This is readily seen from Eq. (22) if we note that e ~ 1/l for l 1.) Similar behavior is found for n = 0.1 and 0.001. The smaller n, the smaller the value of l required for agreement between series and QS values. (For n = 0.1 agreement sets in for e < 0.024, approximately; in the other cases e < 0.06 is sufficient.) At the point where the series and QS results begin to agree, the mean population is not particularly large; for n = 0.01, l = 1.4 corresponds to ánñ ~ 35.

In order to gauge the importance of successive terms, we plot (Fig. 5) Dn (for n = 0.01) as predicted by the series to order e^m, for m = 1,..., 5. For l > 1.5 very good agreement with the QS mean population is obtained using the result to O(e³). For this value of l the absolute difference between the three- and five-term series is about 0.04 The relative difference is about one part in a thousand; the difference rapidly decreases for larger l.

In light of the poor performance of the e series in the vicinity of l = 1, it is natural to apply a resummation technique such as Padé approximants. We therefore constructed the [2,3], [3,2] and [4,1] approximants to the series 1 + g₁e + ... + g₅e⁵ and to the series for the logarithm of this expression. None of the Padé approximants yielded any improvement over the original series; the approximants are ill-behaved (large and negative) near l = 1 (they typically exhibit a pole in this region) and reproduce the excellent agreement with the QS results whenever the original series does.

Before turning to the behavior of the higher moments, we compare the mean population as predicted by the perturbation series with that predicted by the pseudo-stationary distribution. In [12] the pseudo-stationary distribution (PSD) for the contact process on a complete graph was constructed and found to agree well the exact QS distribution, away from the region of the transition. Formally, the PSD, to be denoted as p_ps,n, is obtained by setting p_ps,0 = 0, adjusting the ratios p_ps,n+1/p_ps,n so that the r.h.s. of the master equation is identically zero (for n > 1), and finally normalizing on the pertinent set of states (n > 1 in the present case). Physically, the resulting distribution represents the stationary solution for the case of a reflecting boundary at the origin, rather than an absorbing one. Thus the PSD should agree well with the quasi-stationary distribution when the latter is very small in the vicinity of n = 0. Otherwise the PSD tends to overestimate the probability near n = 0, due to the reflecting boundary [12].

It is easily seen that the pseudo-stationary distribution for the MVP is:

where C is a normalization constant. Expanding p_ps,n about its maximum, which falls at n* ~ 1 + (l – 1)/n for (l – 1)/n 1, we see that the PSD corresponds, in this limit, to a Gaussian with mean and variance both equal to n*. In Fig. 6 we compare the prediction for Dn derived from the PSD with that of the series and of the quasi-stationary distribution. Evidently the pseudo-stationary prediction converges to the QS value almost as rapidly as the series prediction.

The trends observed for the mean population continue when we compare series and QS predictions for higher moments. In Fig. 7 we plot the difference, Dvar(n), between the variance furnished by these predictions and that given by mean-field theory. (Since the latter yields a Poisson distribution, the variance is simply equal to in mean-field approximation.) For l ~ 1 the series prediction is useless; close agreement with the QS result sets in (for n = 0.01), around l = 1.5. Note that Dvar(n) approaches 1/n as l ® ¥, as predicted by Eq. (33); the absolute value of the difference does not approach zero for large l, as it does for the mean. (The relative difference Dvar(n)/var(n) does of course approach zero in the limit l ® ¥.) Similar comparisons for the third and fourth moments reveal that the series prediction agrees well with the QS result (again, for the case n = 0.01), for l > 2.5.

On the basis of these comparisons we conclude that the perturbation theory developed in the previous section is incapable of describing the vicinity of the transition (l » 1), but agrees extremely well with quasi-stationary properties (including higher moments) above a certain, not very large value of l. In this regime, the perturbation approach yields accurate analytic expressions for QS properties.

In Fig. 8 we compare the lowest order W-expansion prediction for Dn and the corresponding result of the e series (truncated at order e) with the QS result, for n = 0.01. The W-expansion result is seen to be slightly better, although it diverges as l ® 1. The inset of Fig. 8 shows the correction Dvar(n) to the mean-field prediction for the variance for the same parameter values. The W-expansion is again slightly superior. The fifth-order e series is in fact superior to the second order W-expansion result, yielding a nearly exact prediction (for the parameter values of Fig. 8) for l > 1.4.

V. MVP WITH INPUT

In this section we examine the effect of a steady input of individuals at rate g. For g > 0 the n = 0 state is no longer absorbing, and the system approaches a true stationary state as t ® ¥. The forward transition rate is now w_n,n+1 = g + ln; the reverse rate w_n-1,n = n +nn(n – 1) as before. The evolution operator is that of the MV process with the addition of a term g(p – 1), in the notation of [9]. This corresponds to a new term in the action (i.e., the argument of the exponential in Eq. (4)),

We now proceed as in Sec. II, introducing the shift of variable y = + f(t). Equating the coefficient of f⁰ in the action to zero yields the stationary solution of the macroscopic equation, that is,

where w = . After the shift, the action has the same form as for the MV process without a source, but the factors associated with the source and the bifurcation are changed. The factor for the source is now:

while that for the bifurcation is 1 – w. The factors associated with the conjunction and the 4-vertex are -n, as before.

Consider the contribution to án_¥ñ/ due to the lowest order diagram (i.e., the second diagram of Fig. 2); this is readily seen to be -wn/w². As before, we may dress the diagram (see Fig. 3), which again has the effect of multiplying the above result by the factor k = [1 + n/w²]^–1. Letting e = 1 – k, the lowest order contribution, when dressed, is -we. For diagrams with n > 4 both rightmost source and the leftmost conjunction can be dressed, leading to the algebraic factor

where m = n – c and u = 1 – w as before. In terms of the parameter e we have

where

and

With these results the enumerations derived for the original process may be used to develop an e series for the MVP with input. The input rate g enters via the expression for w. Up to third order we find,

Series is compared with the exact value of án_¥ñ (from the stationary solution of the master equation, obtained numerically), and with the mean-field prediction in Fig. 9, for parameters g = n = 0.01. Excellent agreement between series and the exact result is found for l > 1.55 and again (see inset) for l < 0.75. In the regions of agreement e is small (e < 0.14 for l < 0.75; e < 0.032 for l > 1.55), but it becomes large in the intermediate region, attaining a maximum of about 0.96 for l = 1. In this region the series prediction deviates strongly from the true value, and can become negative. Similar results are found for other values of g and n.

A. Fokker-Planck approximation

The failure of the perturbative approach when e is not small, in the case with input, cannot be blamed on the absence of a true stationary state. It appears that in this case mean-field theory does not provide a good first approximation on which to base a perturbative expansion. The true probability distribution, moreover, is far from Poisson-like. The only reliable approximation scheme we have found for this case is one based on the Fokker-Planck equation.

To develop this approach we let the population size assume continuous (non-negative) values x, and denote the probability density by p(x,t). Following the usual procedure [1], the master equation for the process may be transformed to a Fokker-Planck equation,

with

and

The stationary solution to Eq. (55) is [1],

where C is a normalization constant. In the present case the integral can be evaluated in closed form, leading to

where

and

Returning to the discrete case, we normalize the distribution using å_{n> 0}p(n) = 1. The resulting prediction for ánñ follows the qualitative behavior of the exact solution (to the master equation), but is quantitatively accurate only for larger values of l (e.g., for l > 1.5 in case g = n = 0.01).

Truly quantitative agreement is obtained if we impose the condition p(0) = p(1)/g, as implied by the master equation, and use the Fokker-Planck solution for for n > 1. In this case the normalization factor is given by

Figure 10 shows that the resulting prediction agrees very well with the exact result. The probability distribution (Fig. 11) is reasonably smooth for n > 1, but it is clear that the jump from n = 0 to n = 1 cannot be reproduced using a continuum approximation. The Fokker-Planck equation with the added condition p(0) = p(1)/g yields an excellent prediction for the stationary probability distribution.

VI. CONCLUSIONS

We apply the path-integral based perturbation theory for Markovian birth-and-death processes [6,9] to the Malthus-Verhulst process and a variant of this process that includes particle input. At zeroth order the formalism yields a Poisson distribution whose mean is given by the mean-field or macroscopic equation. Computational enumeration of diagrams allows us to derive the series coefficients for the first four stationary moments of the process. The expansion parameter, e, is small when the mean population size ánñ is large, but is of order unity when ánñ is small. In the latter regime the expansion fails, but in the former its predictions are in near-perfect agreement with numerical results, despite evidence that the series is divergent. Truncating the series at third order, we find excellent agreement with numerical results when the mean population ánñ > 35. Our analysis yields asymptotic expressions (valid for large population size) for the moments of MVP.

Our study of the MVP with input shows that the poor behavior of the series, when e is not small, is not due the presence of an absorbing state, since input eliminates such a state from the process. Since the present approach is based on systematic approximations to mean-field theory, one should not expect it to yield useful results where the latter is seriously in error, as is the case for l < 1 in the MVP with or without input. For parameter values such that the mean-field solution is reasonable, the series provides useful corrections to it. Development of a globally accurate perturbation method remains as an open challenge, one we hope to explore in future work. An approach based on the Fokker-Planck equation does yield reliable predictions if we use the master equation to fix the form of the probability distribution at n = 0.

Comparison with van Kampen's W-expansion shows that (for the problems considered here), the latter method, to lowest order, yields results that are very similar, though slightly superior to our first-order expansion. Including higher order terms, the quality of the e series improves. In the present case (and in other problems likely to arise in stochastic modeling), deriving higher-order corrections is simpler using the diagrammatic approach, making our method an attractive alternative to the W-expansion.

Acknowledgments

NRM is grateful to FAPEMIG for financial support, and to the hospitality of the Universidade Federal de Minas Gerais. This work was supported by FAPEMIG and CNPq, Brazil.

[1] N. G. van Kampen, Stochastic Processes in Physics and Chemistry, (North-Holland, Amsterdam, 1992).

[2] C. W. Gardiner, Handbook of Stochastic Methods, (Springer-Verlag, Berlin, 1990).

[3] M. Doi, J. Phys. A9, 1465, 1479 (1976).

[4] P. Grassberger and M. Scheunert, Fort. Phys. 28, 547 (1980).

[5] N. Goldenfeld, J. Phys. A17, 2807 (1984).

[6] L. Peliti, J. Physique 46, 1469 (1985); A. Mecozzi, F. De Pasquale, and L. Peliti, Nuovo Cimento 100B, 733 (1987).

[7] B. P. Lee and J Cardy, J. Stat. Phys. 80, 971 (1995); J. Stat. Phys. 87, 951 (1997).

[8] J. Cardy and U. C. Tauber, Phys. Rev. Lett. 77, 4780 (1996).

[9] R. Dickman and R. Vidigal, Braz. J. Phys. 33, 73 (2003).

[10] R. Vidigal and R. Dickman, J. Stat. Phys. 118, 1 (2005).

[11] C. Deroulers and R. Monasson, Phys. Rev. E 69, 016126 (2004).

[12] R. Dickman and R. Vidigal, J. Phys. A 35, 1145 (2002).

Received on 21 August, 2006

In this appendix we collect the terms of our expansion up to O(e⁵), as described in Section II B.

For higher moments, writing the contribution to án^rñ_f coming from _r in Eq. (35) as ^r å_{m> 1}g_r,me^m, we have, for r = 3,4 and m = 5:

and

Appendix

[1 + å

Writing Eq. (24) in the form án_¥ñ = [1 + å_{m> 1}g_me^m], and using u º 2 – l, we have:

Publication Dates

• Publication in this collection
  21 June 2007
• Date of issue
  Dec 2006

History

• Received
  21 Aug 2006