International Journal of Scientific & Engineering Research, Volume 4, Issue 11, November-2013 156

ISSN 2229-5518

A Note on Integro-Differential Explicit Method of the Master Equation


Ahmed Abdul-Razzaq Selman

Department of Astronomy and Space, College of Science, University of Baghdad



A new theoretical approach is presented for solving the master equation of nuclear
preequilibrium states, based on an explicitly integrated method. Earlier ad hoc hypothesis was proposed to follow mathematical derivation of the master equation reasonably, thus providing an explicit method to derive the dependence of probability distributions of exciton levels occupation. Few theoretical comparisons were made, based on approximation formulae of the present solution.

1. Introduction

The preequilibrium models represent a group of models that try to describe the intermediate nuclear reactions based on the contribution weight of many statistical states in each reaction channel. A model with significance is the exciton model [1-4], which depicts intermediate stages of nuclear reactions based on excitations occurring due to development of residual two-body interaction. The exciton model implies that a hypothetically and temporarily created entity, the exciton, a notion for particle-hole (p,h) pair, is responsible of this development. The idea was first described by Griffin [4]. Intermediate states of interaction will ensure that the number of excitons n
(=p+h) distinguishing each stage in the equilibration process will change by an even
integer, that is

n = ±2

or zero. At every stage in this course, equilibration is specified
by an n and excitation energy E, both of which will affect a certain decay probability

W . Transition between adjacent stages is characterized by the transition rate,

a b

between stages a and b as given by Fermi rule, depending on the matrix element of the interaction [2].
During exciton model calculations, the master equation ME, will describe the time evolution of occupation probability of each stage based on a first order differential equation. It expresses various states development at a given time. In [1,5]
a comprehensive review and solution methods of ME description was given. ME in

IJSER © 2013

International Journal of Scientific & Engineering Research, Volume 4, Issue 11, November-2013 157

ISSN 2229-5518

terms of one-component system, when protons and neutrons are approximated to have the same behavioral characteristics, is given as [6]

dP(n, E, t ) = λ

(n, E) P(n + 2, E, t ) + λ+

(n, E) P(n − 2, E, t ) − P(n, E, t )


dt n + 2

n − 2

τ n (n, E)

where P(n,E,t) is the occupation probability of the nth stage with E at t, and

λ+ , λare

adjacent transition rates with resp. The mean lifetime of this stage τ n is

− +


τ n (n, E) =  λn (n, E) + λn (n, E) + W (n, E) 


 

and W(n,E) is the exciton rate of decay to the continuum. The total lifetime for each
state is given by

T (n, E) =


P(n, E, t )dt


For two-component system of reaction where neutrons are distinguished from protons, the ME is given as [5]

dP(N , h

, t )

π =


λ+ + (E, N − 1, h

) + λ+ + (E, N − 1, h

)]P( N − 1, h

− 1, t )

v π π π π π π

+ [λ+ 0 (E, N − 1, h

) + λ+ 0 (E, N − 1, h

)]P( N − 1, h

, t )

π v π v v π π

+ [λ−0 (E, N + 1, h

) + λ−0 (E, N + 1, h

)]P( N + 1, h

, t )

v v π π v π π

+ [λ− − (E, N + 1, h

+ 1) + λ− − (E, N + 1, h

+ 1)]P( N + 1, h

+ 1, t )

π π π π v π π

+ [λ0+ (E, N , h

− 1)]P( N , h

− 1, t )

v π π π

+ [λ0− (E, N , h

+ 1)]P( N , h

+ 1, t )

π v π π

[λ+ 0 (E, N , h

) + λ+ 0 (E, N , h

) + λ+ + (E, N , h )

π v π

+ λ+ + (E, N , h

v v π

) + λ−0 (E, N , h

v π π

) + λ−0 (E, N , h )

π π π v v

π π v π

+ λ− − (E, N , h

) + λ− − (E, N , h

) + λ0+ (E, N , h )

π v π π π

π v π π

+ λ0− (E, N , h

) + W (E, N , h

) ]× P( N , h

, t )


π v π π π


nπ , nν are the exciton numbers for protons and neutrons, respectively,

n = nπ + nν


P( N , hπ , t ) of

nπ , nν number of excitons. N is a function of

nπ , nν

number as:

n = pπ

+ hπ

+ pv + hv

= 2(N + 1) + no , and no is the initial exciton number. h

IJSER © 2013

International Journal of Scientific & Engineering Research, Volume 4, Issue 11, November-2013 158

ISSN 2229-5518

and p represent hole and particle numbers, and the subscripts

π ,ν represent proton and

neutron types, respectively. The decay rate W ( N , hπ ) is defined as [7],

W ( N ,

h ) = W ( N ,

π β


h ) dε




W ( N , hπ )

is the decay rate of a particle of type β to the continuum with energy ε
from the state described by N and

hπ . This problem is of major significance in

preequilibrium models, specially the exciton model [8-10].
In the present work we provide a note on the solution of eq.(4), which is applied for a ne solution of the two-component ME.

2. Theory: The Integro-Differential Master Equation

If one examines eq.(4), it can be shown that it is equivalent to,

dP /

(t ) 10

= λ


P (t ) −  λ/

+ W P / (t )


dt j =1 j j


j =1

where an assumption was made that various population probabilities are,

P( N , h

) = P / = P / = P / = ... = P /

= P /

= P / , and the prime (magnified acute accent) sign is

π 1 2 3

10 11

added in order to distinguish these (successive) probabilities from (initial) P, the
prime here does not mean a derivative operation. Eq.(6) is ME rewritten in a compact
form. Noting that
λ/ 's are not functions of t, neither is W, then we'll further define the


B(t ) = λ j

j =1

P (t )




A = λ/ + W


≠ ( A(t ))


j =1

then eq.(6) can simply be turned to,

dP /

(t ) + A P / (t ) = B(t )



IJSER © 2013

International Journal of Scientific & Engineering Research, Volume 4, Issue 11, November-2013 159

ISSN 2229-5518

Let's define, P / (t ) = u v , where u and v are partial functions of P, and let us assume a
general form of

v = exp(− A t ) , this definition of v implies that

dv + Av = 0, thus,


Pj (t = 0) = Po j

= k = constant

and one may write,
after few steps,


B(t = 0) = B = λ k

o j



P / (t ) = exp(− A t )


t / =0

exp( A t / ) B(t / ) dt / + P / exp(− A t )




/ = P / (t = 0) . This is the integrated form of ME. The exponential function was

added here to describe a general behavior of nuclear decay. It is consistent with the
results of the simpler one-component ME system found in [1].

2.1. The Solution:

The importance of eq.(12) is seen here because it was given in a form that can be solved analytically. In previous methods of ME solution, no direct expression for the probability distribution P/(t) was given, but instead, only the integral form of P/(t) is usually reached that is similar to eq.(12) -see [1] for detailed survey on solution methods of ME.
Let us try to rewrite eq.(6) again as,

dP S S j

k =

dt j =1

λ P

j , k j

j =1

λ P

k , j k


k k


λ j, k is transition rate from the state j to the state k and Sj in the summation
limit is the number of states surrounding the jth state, without including the jth state itself, because there is no transition from the jth state to itself. Note that here the index
k has a similar meaning given by Herman et al.[6] for index representation,

(n n

)(n n

+ 2)

k = o o + h

+ 1 , and it leads to the scheme shown in Fig.(1).

8 π

IJSER © 2013

International Journal of Scientific & Engineering Research, Volume 4, Issue 11, November-2013 160

ISSN 2229-5518

k =

k = 2 3

k = 4 5 6

. . . etc 7

8 9 10

Fig.(1). A scheme representing the transition from various states.

To put eq.(13) in general form let us define

jk as,


∆ = 1 − δ

jk jk

j = k

j k


where δ jk is the usual Kronecker delta function. Obviously,

∆ = ∆

jk kj


Then eq.(13) can be safely written as follows,

dP S

S j

k =

dt j =1

λ P

jk j , k j

− 

 j =1

λ + W

jk k , j k


 k


As before we put the definitions of Bk and Ak as follows,



Bk =

j =1

λ P

jk j ,k j

(17 − a),



A =


j =1

λ + W

kj k , j k

(17 − b),

IJSER © 2013

International Journal of Scientific & Engineering Research, Volume 4, Issue 11, November-2013 161

ISSN 2229-5518

here Sj represents all the states surrounding the jth state with no restrictions, because
we added the analytical function


which forces this conditions as a selection rule
kj = 0 when k = j. Eq.(13) is now equivalent to,


k = B



and eq.(12) becomes,

dt k k k

P = exp(− A t )



t / =0

exp( A t / ) B


dt /

+ P exp(− A t )

o k k


Once again, the functions Ak is not a function of time. t/ also is used for integration purpose only. The solution of the explicitly-integrated form of the master equation is given below. To integrate eq.(19) let's define part of the integral as,

I =

Then after few steps,



j =1


jk j ,k


t / =0

exp( A t / )P (t / ) dt / ,



P = j k

λ P

jk j ,k j


+ P exp(− A t )

ok k


jk λ




j ,k j

× T T


A exp(− A t )

k k

exp( A t ) T


T )dt 



j k k

W exp(− A t )

P exp( A t ) dt


which represents the explicit solution of ME. The



stat is presented to differ from



the jth one, from defining B


 =1


λ, j P


A j =

 =1


λ j, 

+ W j

. Eq.(21) was
derived without approximation. Apparently, the first difficulty faced is the integrations of T’s in the second and third terms.

IJSER © 2013

International Journal of Scientific & Engineering Research, Volume 4, Issue 11, November-2013 162

ISSN 2229-5518

2.2. Approximations

However, one may apply some approximations at this point because of the significant and specific shape of eq.(21). As a first approximation let the equilibration time teq. be long enough so that the exponentials with negative exponents fall fast with time, and W are not significant at each stage, and let the mean lifetimes of adjacent stages


be almost equal, i.e.,



W j exp(− Ak t ) ≈ 0 . Then

P = j + P

k A ok

exp(− A t ),




Pk Pj . This means that the occupation probability of the state k relates to all

occupation probabilities of the surrounding jth states because every kth state is actually a newly born daughter and a mother of the surrounding states at the same time which is the idea of detailed balance behind ME. The proportionality above is therefore an image taken from eq.(5).
At small values of time, then if one accepts the approximation

 

 

jk λ j, k

jk λ j, k + W j  ≈ 1 which means that W is much smaller (in

 

magnitude) than the sum of transition rates then one can have,


+ P (t = 0) exp(− A t ) .

k j k k

3. Discussions and Conclusion

The physical interpretation of eq.(13) ensures that for stage with k=1 (first stage), the
following approximation holds,

exp(− A t )


t /

exp( A t / )


B dt /


≅ 0 because the first stage

mainly gives to other stages, its chance to receive from other stages is almost zero.

Thus, P P

1 o1

exp(− A t ) which is very similar to the well-known nuclear decay


N = No exp(−λ t ) , noting that in here A = λ j . Therefore we choose the expression


Pok exp(− Ak t ) rather than


only. One should be careful how to apply and use the

IJSER © 2013

International Journal of Scientific & Engineering Research, Volume 4, Issue 11, November-2013 163

ISSN 2229-5518

approximation above because we have,
where β B .Thus we'll have,


exp( A1t )

j =1

j1 λ j,1 Pj

= constant =

1 exp(− A1t )


j =1

λ P

j1 j ,1 j

= β exp(− A t )

1 1


Now if we take the time derivative for both sides of eq.(22) then,

− 1 3 dP

β exp(− A t ) = ∆ λ j

(23) ,

and from this we have,

1 1 A


j =1

j1 j ,1 dt



 1


j ,1 A


j + P  = 0



j =1

 1 dt 

and clearly this equation holds if ( 1




+ Pj )=0,

λ j,1 = 0 , or j=1; all of which will
logically lead to, P


= P exp(− A t ),

oj 1

( j = 2 and 3 only) which is the exact form of

eq.(24) for j=1. Eq.(21) will lead to the solution of all stages with the same form, i.e.,
will give a solution of the type, P


= P exp(− A t ) , for all values of j and also it reveals

oj 1


P / (t ) = u v

actually reduced to

P / (t ) = v , and this isn't quite efficient as a full

solution, but it is an approximation. Therefore eq.(21) must be handled with care as being for the case with j=1 only because is not valid for all stages of the preequilibrium reaction, but the more reasonable form, eq.(19) is always valid which is a useful relation that describes the population probability of the kth state is the same of the jth state plus the initial population of the same kth state, which decreases
exponentially with time t. Therefore, the asymptotic behavior of the kth occupation
probability will directly be as Pk
equilibration condition.

= Pj and this is the situation that describes the

A third point seen from the above eq.(21) is that, if one assumed that initially Pj
= 0, which means that all surrounding states are empty. Then

P = P (t = 0) exp(− A t ) ,

k k k

IJSER © 2013

International Journal of Scientific & Engineering Research, Volume 4, Issue 11, November-2013 164

ISSN 2229-5518

and this again is the well-known decay law of any excited equilibrated, or compound nucleus. In the case of Pj = 0 then we will have no jth states, therefore there will be no
transition rates from or to these states and we must set all λ j
to zero, i.e., we will
have Ak = Wk only. This leads immediately to have a solution form as
P = P (t = 0) exp(−W t ) , simply, if we labeled W by λ, and Pk (t=0) by Po, then we write

k k k

P = P0 exp(−λ t ) , the natural radioactive decay law.

To compare with earlier ad hoc hypothesis, Griffin [4] put the following formula -the
same symbols are used as found in Ref.4 for convenience,

ρ (u ) ρ

n −1 n

(E )

P ( E ) dE = n dE ,


and Blann [10, 11] put the following,


n −1


nˆ ρ

(U ) 

λ (E) 

P (E) dE =

f n −1 c D

x n n =n


x 

ρ n (E)   λc

(E) + λ

n+ 2

(E)  n

ρn (E) is the density of states with exciton number n and energy E,

n f x

is the
emission ratio of a particle of type x,

c (E) is the decay rate to the continuum, i.e.,

λ = W , and


Dn =


(1 − P / ).

n −2

In this paper we've tried to find an analytical

n =n +2


reason behind this assumption. Comparing Blann's [11] result in particular

λ    


P with our approximation

λ λ

+ W  ≈ 1

λc + λ

n + 2 n /

(1 −

n / − 2 ) 

jk j, k

jk j, k j

a general analogy is seen, as both have ratio of transition rates multiplied by probability function, and the definition of compensates for the term inside the multiplication procedure.


IJSER © 2013

International Journal of Scientific & Engineering Research, Volume 4, Issue 11, November-2013 165

ISSN 2229-5518

[1] A. A. Selman, "Neutron Induced Preequilibrium Nuclear Reactions Using the
Exciton Model", Ph.D. Thesis, University of Baghdad (2009).
[2] A. A. Selman, "Preequilibrium Nuclear Reactions Using the Exciton Model
Review and Applications", 1st Ed., LAP Lambert Academic Pub. GmbH, (2012).
[3] A. A. Selman, "The Exciton Model Part I. Basic Calculations of Nuclear Level
Density", 1st Ed., LAP Lambert Academic Pub. GmbH, (2012).
[4] J. J. Griffin, “Statistical Model of Intermediate Structure”, Phys. Rev. Lett.,


[5] M. H. Jasim, S. S. Sahik and A. A. Selman, " A Suggested Numerical Solution for the Master Equation of Nuclear Reaction", J. Kerbala Uni., 7,(2009) 271-282.
[6] M. Herman, G. Reffo, and C. Costa, "Early Stage Equilibrium Dynamics in a
Two-Component Nuclear System", Phys. Rev. C39 (1989) 1269.
[7] J. Dobeš and E. Bĕták, "Fast and Precise Exciton Model Calculations of the
Nuclear Reactions", Z. Phys. A288 (1978)175.
[8] F. J. Luider, “Note on the Solution of the Master Equation in the Exciton Model of
Preequilibrium Theory”, Z. Phys. A284 (1987)187.
[9] J. Dobeš and E. Bĕták, “Two-Component Exciton Model”, Z. Phys. A310 (1983)329.
[10] M. Blann, “Extension of Griffin’s model for medium-energy nuclear reactions”, Phys. Rev. Lett. 18(1968)1357,
[11] M. Blann, “Hybrid Mode for Pre-Equilibrium Decay in Nuclear Reactions”, Phys. Rev. Lett. 27(1971)337

IJSER © 2013