Approximate time-dependent solution of a master equation with full linear birth-death rates

There are growing interests on dynamics of phase-singularities (PSs) in complex systems such as ventricular fibrillation, defect in fluids and liquid crystals, living creatures, quantum vortex and so on. A master equation approach on the number of PS for studying birth-death dynamics of PSs is invented first by Gil, Lega and Meunier. Although their approach is applied to various complex systems including non-linear birth-death rates, time-dependent solution of related master equation is obtained only rarely. Even a master equation with full linear birth-death rates, time-dependent solution is not also given due to the analytical complexity and the existence of singularity in the probability generating function. In this paper, an approximate time-dependent solution of the master equation and the associated waiting time distribution are obtained explicitly with the aid of the method of the Poisson transform. Numerical evaluation of the obtained approximate solution teaches us that there exists the universal scaling law in the waiting time distribution.


Introduction
There are growing interests on dynamics of phase-singularities (PSs) in complex systems such as ventricular fibrillation [1,2,[4][5][6], defect dynamics in fluids [7,8] and liquid crystals [9], living creatures [10], quantum vortex dynamics [11] and so on. A master equation approach on the number n of PS for studying birth-death dynamics of PSs in a 2D Complex Ginzburg-Landau equation is invented first by Gil, Lega and Meunier [12]. Although the method is applied in various real systems [7,8,13], it can be utilized to get only the equilibrium distribution of PSs, P s (n). To catch the true dynamical features of PSs, it is necessary to get the information on (i) the PS number distribution, (ii) the waiting time (lifetime) distribution, (iii) the velocity distribution of PS. For the case of 1D CGLE in a region of the amplitude turbulence, we could describe the related stochastic processes having a long-memory through the three distributions (i), (ii) and (iii) by taking a birth-death process (λ n =ν and μ n =μn) of the PS number n [14], and a generalized Cauchy process of the velocity of PS [15].
For the cases of 3D CDLE and Barkley model, Davidsen et al [1] and Reid et al [2] reported a complicated birth-death dynamics on filaments of PSs. Interestingly, they reported that a Markovian master equation with the birth (μ n =ò+μn ) and death (λ n =ν+λn) rate of filaments in a state of the negative line tension (NLT). Their results on the birth and death rates of filaments indicate that (a) the event of pair-extinction of filaments is not dominant in a state of the NLT. Is it possible to get satisfactory description of (i), (ii) and (iii) in a 3D dynamical model on the NLT dynamics of filaments. For the full linear birth-death rates, λ n =λn+ν and the death μ n =μn+ò, there are various problems due to the existence of ò (the details will be shown in later sections): (A) the probability mass (PM) of n particles at time t, P(n, t) under the initial condition P n, 0 n n , 0 d = ( ) is not given explicitly; (B) the extinction probability and the waiting-time (lifetime) distribution associated with P(0, t) are not given explicitly; (C) Method of estimating the parameters (λ, μ, ν, ò) in the model is not known due to the complex expression of the steady state distribution.
The paper is organized as follows. Section 2 describes our model and the exact solution of the generating function. It is explained that there are problems for obtaining the expressions of the PM, the extinction probability and the waiting time distribution (WTD) from the obtained generating function. Section 3 presents an approximate solution obtained by the method of the Poisson transform. Approximate analytic expressions of the PM and the WTD are given explicitly. Section 4 presents numerical example of them for the full linear birthdeath process. Section 5 is devoted to summary and remarks.

Model
Consider a birth-death master equation for the probability mass P(n, t) of n particles at time t with the rate constants with λ n and μ n as d dt P n t P n t P n t P n t P n t , 1 , , The paper focuses on the full linear rates: where λ, μ, ν, ò are constants. The initial condition is assumed as Applications with the rates in equation (2) without ò can be found in [16][17][18][19][20]. The general theory of birthdeath stochastic processes in equation (1) has been formulated by Kerlin and McGregor [21] based on the orthogonal polynomials. For λ n =λ and μ n =μ, the probability mass of n particles at time t P n t n t , , , 0 0 ( | ) under the initial condition P n t , is expressed by the Chevishev polynomials. For λ n =ν and μ n =μn, it is given by the Charier polynomials. For λ n =λn+ν and μ n =μn [22], it is given by the Meixner polynomials. For other nonlinear birth-death rates, one can see related solutions in references by Ismail [23], Sasaki [24] and Schoutens [25]).
Throughout the paper, we will focus on the master equation (1) in the birth-death rates in equation (2) with μ>λ that there exists the steady state P s (n). Since the detailed balance is satisfied, where P 0 is the normalization constant. The distribution can be rewritten as where the GF is defined by , . 9 From the initial condition in equation (3), g z z , 0 . The method of characteristics provides the exact solution of the generating function. The characteristic equation is given by The solution g(z, t) is readily obtained in the form: A derivation of the solution is described in appendix A. As far as we know, this is the first time to show explicitly the generating function (GF) for the full linear birth-death rates. In the famous literatures such as Kerlin and McGregor [21,22] and Goel and Dyn [17], the result in equation (11) is not exhibited.
To demonstrate the validity of equation (11), let us expand the equation in terms of the Meixner polynomials of the first kind, M n (x; β, c), which is defined by is the hypergeometric function. Applying the mathematical formulae on the Meixner polynomials to equation (11) by choosing ; , 1 exp . 15 The functional forms of the Meixner polynomials are Equations (15) and (16) provides the expressions of the time-dependent mean n t g z t , z 1 á ñ = ¢ = ( ) ( ( )| ) and variance t g z t g z t g z t , , , ) (g¢ and g denote as the first and the second derivative with respect to z, respectively) as n t n e e 1 1 8 These time-dependent mean and variance reduce to the steady state values in equation (7).
On the other hand, the GF in equation (11), it is difficult to evaluate the probability mass (PM), 1 , , 20 the extinction probability, P(0, t), and the waiting time distribution, f P 0, ), in the model due to the singularity at z=0 in equation (11). Therefore, it is worth to derive approximate analytic functions of the extinction probability and the waiting time distribution as will be shown in the next section.

Poisson transform
There are various methods [16] for approximating the master equation in equation (1). One of the most popular method for approximating the master equation is the system size Ω expansion method (SSEM) proposed by van Kampen [26]. The SSEM utilizes an expansion of a random state vector X as X x y = W + W . So, one must consider coupled equations of a scaled mean vector x and a scaled fluctuation vector y at the same time to treat a problem.
As far as a single-variable master equation is concerned, the Poisson transform provides more tractable mathematical tool than the SSEM. Let us now utilize the method of the Poisson transform (PT) [16], P n t n n n y y y t dy , e x p , , 2 1 to get an approximate solution of the master equation in equation (1). The PT involves the parameter n á ñ, which is equivalent to the mean (the variance) of the Poisson distribution. Applying the inverse Poisson transform to equation (1), one obtains an approximate Fokker-Planck equation (FPE) [18], where the regression and diffusion term are defined by where α, β, γ and δ are given by n n n , , a n d .
The FPE in equation (23) is an approximate one for the master equation in equation (1) with the birth-death rates in equation (2), since the higher order derivatives n y y t n y y t , , .
are neglected by accounting up to the second order derivative y t ,  (25) and (26)where ò and n á ñ are involved, one can see that (i) the accuracy of the FP approximation with the PT depends ò and n ; á ñ (ii) the parameter n á ñ plays a role of the SSEM (see [26] in details).
Interestingly, the Langevin equation associated with the FPE in equation (23) is written in the form with the use of a single noise dW as

Steady state solution
Rewriting the FPE into a simple form, let us take an affine transform z y c y .
Consequently, the regression and diffusion term in equation (23) and D (z)=δz, respectively. The FPE in the z coordinate takes the same form as the gamma process of the first kind [16]. In mathematical economics, they call equation (23) (or equation (27)) in the z-coordinate as the CIR model [19].
When the parameters α, β, γ and δ take positive values, it is easy to see that the steady state exist as t  ¥. The steady state solution of the FPE in the z coordinate is given by the gamma distribution is the incomplete Gamma function of the second kind. Note that y 0, ]since z is defined by equation (28). So, the normalization factor is Γ(a, bc), not Γ(a). As far as a small value of c>0 is concerned, Γ (a, bc)≈Γ(a). The parameters (α, β, γ, δ) involved in the gamma distribution in equation (29) are expressed in terms of the birth-death rates (λ, μ, ν, ò) and the parameter of the Poisson transform n á ñ as defined in equation (21).
By rewiring equation (29) in the y coordinate, one obtains and a, b and c are expressed by the rate of events (λ, μ, ν, ò) as A steady-state exists provided that the parameters (λ, μ, ν, ò) are positive constants, μ>λ and ν>ò. The PM corresponding to the distribution ψ s (y) with (a, b, c) is given by a e e y y c d y e n a n a n bc , , 35 This is a modified Pascal (or a modified negative polynomial) distribution. The parameter ò is contained in the expression in equation (36) (cf the exact steady state solution in equation (6)). It is clear that when ò=0 (in this case, a b n c , , and n a n bc n a 1, 1 ; G + -á ñ +  G + ( ( )) ( )), the above PM, P FP, s (n), reduces to the Pascal distribution (or the negative binomial distribution) as where L x where I q (x) is the modified Bessel function of the first kind as The PM under the FP approximation in y-coordinate is given by Putting n=0 in equation (44), the extinction probability is obtained as  (44) and (45)give the PM and the extinction probability under the FP approximation. Using the formula f FP (τ)=−dP FP (0, τ)/dτ, one can evaluate numerically the waiting time distribution (WTD) f FP (τ). When ò=0 (viz., c = 0), the Poisson transform can be performed explicitly, which are given below. Namely, the PM for ò=0 is obtained as

Direct numerical evaluation of equations
By applying the integral formulae described in equation (C6) in appendix C to equation (42), we obtain In the special case of ò=0, the extinction probability with the initial condition P n, 0

Steady-state solution
The steady-state PM, P s (n), in equation (6) Figure 2 shows the waiting time distribution (WTD) f (τ) for the set of parameters (λ, μ, ν, ò): (i) (1, 2, 3, 0) (solid line), (ii) (1, 2, 12, 0) (dashed line) and (iii) (1, 2, 20, 0) (dotted line). As the value of a increases, the simple exponential decay tend to change into the double-exponential decay as seen in the figure. The WTD curve with direct numerical evaluation with for the case (i) agrees completely with the exact WTD in equation (49).  . In any case, the slope in each figure takes almost the same value β. This means that the slope for large τ is universal constant independent of the value of ò. If this is true, one can easily estimate the four parameters (λ, μ, ν, ò) with the use of the set of data, the PM and the WTD, as will be shown later.

Parameter estimation
From the viewpoint of real experimental (or numerical) data analysis, estimation of parameters is important. There are various methods for estimating parameters: (i) the maximum likelihood estimator (MLE); (ii) the moments; (iii) the least square fitting. In estimating parameters of the hypergeometric distribution in equation (6), the method of MLE becomes difficult since the derivatives of the hypergeometric function are contained in the expressions of MLEs. Here we will propose a simple method based on the method of moments below.  In the analytical expression of the PM in equation (6), the three parameters , a n d 5 1 q q q q q q q q q q q q q = - Therefore, the three parameters θ 1 , θ 2 and θ 3 are estimated by the observed values of M, V and K as , a n d 2 .

Relation to numerical experiments on 3D scroll waves
Davidsen et al [1], Reid et al [2] and Clayton [4,5] studied 3D scroll waves in the state of negative linear tension and of turbulence in the 3D Barkley model [3], the 3D complex Ginzburg-Landau equation [28] and 3D 3V-SIM model [29]. In the case of the negative linear tension, they reported that the birth-death rates are identified well by  values of parameters reported by Reid et al [2] are λ=0.018, μ=0.08, ν=1 and ò=−2.8 in the state of negative line tension. As far as the analytical approach based on the generation function (GF) is concerned, (i) the mean and the variance can be estimated by the GF in equation (11); (ii) the steady-state PM in equation (6) can reproduce the probability distributions obtained in the numerical experiments. However, it is difficult to obtain (iii) the time dependent solution of P(n, t), the extinction probability P(0, t) and the waiting time distribution f (τ) for the process, due to the singularity located at z=0 in the GF in equation (11). When the value ò<0 is small enough to give the parameter a>0 in equation (33) where the state variable z in equation (29) takes positive value, the method of our approximation works. However, the present method fails to get timedependent solutions for large negative value of ò. Concerning the steady-state properties, it is worth to mention that (i) the measured PM by Reid et al [2] can be evaluated by the theoretical PM in equation (6), and that (ii) the parameter estimation in equations (51)-(55)works to evaluate them even for a large negative values of ò (λ=0.018, μ=0.08, ν=1 and ò=−2.8).

Summary and remarks
The master equation in equation (1) with the full linear birth-death rates in equation (2) is studied. Although the GF of the master equation is obtained in equation (11), due to the singularity located at z=0, the PM P(n, t), the extinction probability P(0, t) and the WTD f P FP (n, t) is given under the constraint that the initial value n 0 is distributed by the Poisson law. As far as we know, this is the first time to obtain the explicit time-dependent solution for the full linear birth-death rates in equation (2). These results are useful for various birth-death processes wherein the four contributions of events with the four parameters (λ, μ, ν, ò), λ n =λn+ν and μ n =μn+ò, are relevant in the systems considered (cf [2,[16][17][18]20]). From the PM and the WTD, one can infer the four parameters as was shown in section 4.4.
In the special case of c=ò=0, the results of the FP equation in equation (23) are expressed by the Laguerre polynomials as shown in equation (41). Kerlin and McGregor [22] obtained the polynomial solution of P(n, t) in terms of the Meixner polynomials for the special case ò=0 of the birth-death rates. We have obtained an approximate solution of the PM under the initial condition of Poisson distribution of n 0 , at least, for positive values of parameters (λ, μ, ν, ò). Therefore, analysis of initial transients can be possible under the Poisson distribution of initial values when the initial condition is uncontrollable in real and/or numerical experiments.
Eliminating z by substituting (A3) into the second equation (A2) for g(z, t), one obtains the time dependent coefficient V(t) and W(t), and that of time-independent constant y. Inserting (A3)-(A4) into (A5), we have , On account of the relation (A7), the solution of (A6) is given by ,0e x p , A 8 where the integrals I 1 and I 2 are defined as I V s y W s ds A9

Appendix B. Relation between the master equation and the Fokker-Planck equation
The one-to-one correspondence of each term of the master equation to the Fokker-Planck equation is summarized in the table below. There are exact one to one correspondent exist for the items (1), (3) and (4). On the other hand, the correspondence is not exact for the item (2). Namely, the higher order derivatives are neglected (see text). where I α (x) is the modified Bessel function of the second kind.
The following integral formula for Laguerre Polynomials is also well known: where I q (z) is the fractional order modified Bessel function of the first kind and F a b z , ; 1 1 ( )is the Wittacker function (cf [27]).