The fractional kinetic equation and thermonuclear functions

arXiv:astro-ph/0001279v1 16 Jan 2000 THE FRACTIONAL KINETIC EQUATION AND THERMONUCLEAR FUNCTIONS H.J. Haubold Outer Space Office, United Nations, Vienna International Centre, P.O.Box 500, 1400 Vienna, Austria and A.M. Mathai Department of Mathematics and Statistics, McGill University, 805 Sherbrooke Street West, Montreal, Quebec, Canada H3A 2K6 Abstract. The paper discusses the solution of a simple kinetic equation of the type used for the computation of the change of the chemical composition in stars like the Sun. Starting from the standard form of the kinetic equation it is generalized to a fractional kinetic equation and its solutions in terms of H-functions are obtained. The role of thermonuclear functions, which are also represented in terms of G- and H-functions, in such a fractional kinetic equation is emphasized. Results contained in this paper are related to recent investigations of possible astrophysical solutions of the solar neutrino problem. 1 1. Introduction A spherically symmetric, non-rotating, non-magnetic, self-gravitating model of a star like the Sun is assumed to be in thermal equilibrium and hydrostatic equilibrium, with a non-uniform chemical composition throughout. The star is characterized by its mass, luminosity, effective surface temperature, radius, central density, and central temperature. For a given mass, four of these variables are independent and are governed by four simultaneous, nonlinear, ordinary differential equations of the first order and four boundary conditions. Since there are four equations but more than four unknowns, additional information must be provided through the equation of state, nuclear energy generation rate, and the opacity (constitutive equations). The assumptions of thermal equilibrium and hydrostatic equilibrium imply that there is no time dependence in the equations describing the internal structure of the star (Kourganoff, 1973, Perdang, 1976, Clayton, 1983). The evolution of a star like the Sun is governed by a second system of differential equations, the kinetic equations, describing the rate of change of chemical composition of the star for each species in terms of the reaction rates for destruction and production of that species (Kourganoff, 1973, Perdang, 1976, Clayton, 1983). Methods for modeling processes of destruction and production have been developed for bio-chemical reactions and their unstable equilibrium states (Murray, 1989) and for chemical reaction networks with unstable states, oscillations, and hysteresis (Nicolis and Prigogine, 1977). Stability investigations of thermonuclear reactions of stellar interest have not yet been worked out in detail. However, the potentiality of instabilities in thermonuclear chains may not be overlooked, since, as was pointed out once by Eddington, “what is possible in the (Cavendish) Laboratory may not be too difficult in the Sun” (Perdang, 1976, Mestel, 1999). Consider an arbitrary reaction characterized by a time dependent quantity N = N(t). It is possible to equate the rate of change dN/dt to a balance between the destruction rate d and the production rate p of N, that is dN/dt = −d + p. In general, through feedback or other interaction mechanisms, destruction and production depend on the quantity N itself: d = d(N) or p = p(N). This dependence is complicated since the destruction or production at time t depends not only on N(t) but also on the past history 2 N(τ ), τ < t, of the variable N. This may be formally represented by dN/dt = −d(Nt ) + p(Nt ), (1) where Nt denotes the function defined by Nt (t∗ ) = N(t − t∗ ), t∗ > 0. Here d and p are functionals and eq. (1) represents a functional-differential equation. In the following we study a special case of this equation, namely the equation dN/dt = −αN(t) (2) with a constant α > 0. Eq. (2) implies that spatial fluctuations or inhomogenieties in the quantity N(t) are neglected. The standard solution of the differential equation (2) will be briefly discussed in Section 2 and the generalization to a fractional differential equation and its solution will be derived in Section 3. Conclusions will be drawn in Section 4. 2. Standard Kinetic Equation The production and destruction of species is described by kinetic equations governing the change of the number density Ni of species i over time, that is, X X dNi Ni Nj < συ >ij + Nk Nl < συ >kl , (3) =− dt j k,l6=i where < συ >mn denotes the reaction probability for an interaction involving species m and n, and the summation is taken over all reactions which either produce or destroy the species i (Haubold and Mathai, 1995). For a gas of mass density ρ, the number density Ni of the species i is expressed in terms of its abundance Xi , by the relation Ni = ρNA Xi /Ai , where NA is Avogadro’s number and Ai is the mass of i in mass units. The mean lifetime τj (i) of species i for destruction by species j is given by the relation λj (i) = 1 Xj = Nj < συ >ij = ρNA < συ >ij , τj (i) Aj (4) where λj (i) is the decay rate of i for interactions with j. Eq. (4) reveals the physical importance of < συ >ij for the kinetic equation (3). 3 In the case of a nondegenerate, nonrelativistic gas, if a nonresonant charged nuclear reaction proceeds at low energies dominated by Coulombbarrier penetration, the reaction probability < συ >mn takes the form (Clayton, 1983, Bergstroem et al., 1999) < συ >mn = ( 2 1 8 1/2 X S (ν) (0) ) I2 (ν, a, z, ρ), −ν+1/2 πµ ν! ν=0 (kT ) (5) where I2 represents a thermonuclear function given by I2 (ν − 1, a, z, ρ) = Z ∞ 0 dyy ν−1e−ay−zy , ν > 0, a > 0, z > 0, ρ > 0 −ρ 1 a−ν 2,0 = H0,2 az ρ ρ  (ν,1),(0, ρ1 )  , (6) where H denotes Fox’s H-function (Mathai, 1993, Haubold and Mathai, 1998), which was introduced into mathematical analysis to unify and extend existing results on symmetrical Fourier kernels. Fox’s H-function (Mathai and Saxena, 1973, 1978) is defined in terms of a Mellin-Barnes type integral H(z) = where h(s) = m,n Hp,q (a ,α ),...,(a ,α ) z|(b11,β11),...,(bqp,βqp) h i 1 Z = h(s)z −s ds, 2πi L n {Πm j=1 Γ(bj + βj s)}{Πj=1 Γ(1 − aj − αj s)} {Πqj=m+1 Γ(1 − bj − βj s)}{Πpj=n+1Γ(aj + αj s)} (7) (8) √ i = −1, and L is a suitable path which will be briefly described in the following. An empty product is interpreted as unity and it is assumed that the poles of Γ(bj +βj s), j = 1, . . . , m, are separated from the poles of Γ(1−aj −αj s), j = 1, . . . , n; a1 , . . . , ap , b1 , . . . , bq are complex numbers; α1 , . . . , αp , β1 , . . . , βq are positive real numbers. The poles of Γ(bj + βj s), j = 1, . . . , m, are at the points s = −(bj + ν)/βj , j = 1, . . . , m, ν = 0, 1, . . . and the poles of Γ(1 − aj − αj s), j = 1, . . . , n, are at s = (1 − ak + λ)/αk , k = 1, 2, . . . , n, λ = 0, 1, . . . The condition of separability of these two sets of poles imposes that there be a strip in the complex s-plane where the H-function has no poles. There are 4 three types of paths L possible. These correspond to paths 1,2, and 3 shown in Fig. 1. Fig. 1. For all practical problems where H-functions are to be applied one mainly requires paths 2 and 3. Hence the conditions coming from these two paths are given in the following. It is to be pointed out that when more than one path L makes sense then it can be shown that they lead to the same function and thus there will be no ambiguity. Let p q µ= X j=1 and βj − X αj (9) j=1 −βj β = Πpj=1 αjαj Πqj=1βj . (10) The H-function exists for the following cases: Case Case Case Case (i) (ii) (iii) (iv) q ≥ 1, µ > 0, H(z) q ≥ 1, µ = 0, H(z) p ≥ 1, µ < 0, H(z) p ≥ 1, µ = 0, H(z) exists exists exists exists for all z, z 6= 0. for |z| < β −1 . for all z, z 6= 0. for |z| > β −1 . In the above cases it is assumed that the basic condition is satisfied that the poles of Γ(bj + βj s), j = 1, . . . , m and Γ(1 − aj − αj s), j = 1, . . . , n are separated. Note that in cases (i) and (ii) the H-function is evaluated as the sum of the residues at the poles of Γ(bj + βj s), j = 1, . . . , m and in cases (iii) and (iv) the H-function is evaluated as the sum of the residues at the poles of Γ(1 − aj − αj s), j = 1, . . . , n. When α1 = . . . = αp = β1 = . . . = βq = 1, the H-function reduces to Meijers’s G-function. That is, h (a ,1),...,(a ,1) i h a ,...,a i p 1 m,n Hp,q z|(b11,1),...,(bqp,1) = Gm,n p,q z|b1 ,...,bq . (11) When ρ in eq. (6) is real and rational, then the H-function can be reduced to a G-function by using the multiplication formula for gamma functions.   1−m 1 j Γ(mz) = (2π) 2 mmz− 2 Πm−1 , m = 1, 2, . . . (12) Γ z + j=0 m 5 In the case of Coulomb-barrier penetration (Gamow factor), ρ = 12 , I2 reduces to # " 1 a−(ν+1) 3,0 az 2 1 , (13) I2 (ν, a, z, ) = G0,3 2 π 1/2 4 ν+1,a, 2 where G denotes Meijer’s G-function, which was introduced into mathematical analysis in attempts to give meaning to the generalized hypergeometric function p Fq in the case p > q + 1. Proceeding with eq. (3), the first sum in eq. (3) can also be written as − X j Ni Nj < συ >ij = −Ni ( X Nj < συ >ij ) = Ni ai , (14) j where ai is the statistically expected number of reactions per unit volume per unit time destroying the species i. It is also a measure of the speed in which the reaction proceeds. In the following we are assuming that there are Nj (j = 1, . . . , i, . . .) species j per unit volume and that for a fixed Ni the number of other reacting species that interact with the i-th species is constant in a unit volume. Following the same argument we have for the second sum in eq. (3) accordingly, + X Nk Nl < συ >kl = +Ni bi , (15) k,l6=i where Ni bi is the statistically expected number of the i-th species produced per unit volume per unit time for a fixed Ni . Note that the number density of species i, Ni = Ni (t), is a function of time while the < συ >mn , containing the thermonuclear functions (see eqs. (5) and (6)), are assumed to depend only on the temperature of the gas but not on the time t and number densities Ni . Then eq. (1) implies that dNi (t) = −(ai − bi )Ni (t). dt (16) For eq. (16) we have three distinct cases, ci = ai − bi > 0, ci < 0, and ci = 0, of which the last case says that Ni does not vary over time, which means that the forward and reverse reactions involving species i are in equilibrium; such a value for Ni is called a fixed point and corresponds to a steady-state behavior. The first two cases exhibit that either the destruction (ci > 0) of species i or production (ci < 0) of species i dominates. 6 For the case ci > 0 we have dNi (t) = −ci Ni (t), dt (17) with the initial condition that Ni (t = 0) = N0 is the number density of species i at time t = 0, and it follows that Ni (t)dt = N0 e−ci t dt. (18) The exponential function in eq. (18) represents the solution of the linear one-dimensional differential equation (17) in which the rate of destruction of the variable is proportional to the value of the variable. Eq. (17) does not exhibit instabilities, oscillations, or chaotic dynamics, in striking contrast to its cousin, the logistic finite-difference equation (Perdang, 1976, Haubold and Mathai, 1995). A thorough discussion of eq. (17) and its standard solution in eq. (18) is given in Kourganoff (1973). 3. Fractional Kinetic Equation In the following, for the sake of brevity, the index i in eq. (17) will be dropped. The standard kinetic equation (17) can be integrated N(t) − N0 = −c 0 Dt−1 N(t), (19) where 0 Dt−1 is the standard Riemann integral operator. The generalization of this operator to the fractional integral of order p > 0 is denoted by a Dt−p and is defined, following Riemann-Liouville, based on the Cauchy formula, by 1 Zt −p dτ f (τ )(t − τ )p−1 , p > 0 (20) a Dt f (t) = Γ(p) a with 0 a Dt f (t) = f (t) (Oldham and Spanier, 1974, Miller and Ross, 1993). The most general fractional integral operator of the type (20) contains Fox’s H-function (7) as the 7 kernel function. If f (t) is continuous for t ≥ a, then integration of arbitrary real order has the property −q −p a Dt (a Dt f (t)) =a Dt−p−q f (t). Replacing the Riemann integral operator by the fractional Riemann-Liouville operator 0 Dt−ν in eq. (19), we obtain a fractional integral equation corresponding to eq. (17) N(t) − N0 = −cν 0 Dt−ν N(t). (21) For dimensional reasons, the coefficient c in eq. (19), containing the probabilities of the reaction under consideration, had to be replaced by cν accordingly. The Laplace transform of the Riemann-Liouville fractional integral is L where n o −p −p 0 Dt f (t); p = p F (p), F (p) = Γ(p) Z ∞ τ =0 (22) dτ e−pτ f (τ ). In order to solve eq. (21), the integral equation is exposed to a Laplace transformation leading to N(p) = L{N(t); p} = N0 p−1 . 1 + ( pc )−ν (23) To arrive at a representation of eq. (23) in terms of Fox’s H-function, Mathai and Saxena’s (1978) result can be used, β β zβ 1,1 −α α ( α ,1) = a az H , β 1,1 (α ,1) 1 + az α " !ν # 1 1,1 c ( ν1 ,1) . N(p) = N0 H1,1 ( ν1 ,1) c p   (24) To prepare eq.(24) for an inverse Laplace transform, the following two fundamental properties of an H-function can be used, i h 1 m,n h (a1 ,α1 ),...,(ap ,αp ) i (a ,kα ),...,(a ,kα ) m,n z k |(b11,kβ11),...,(bqp,kβqp) , k > 0, Hp,q z|(b1 ,β1 ),...,(bq ,βq ) = Hp,q k 8 (25) h (a ,α ),...,(a ,α ) m,n Hp,q z|(b11,β11),...,(bqp,βqp) leading to i n,m = Hq,p  1 (1−b1 ,β1),...,(1−bq ,βq ) , | z (1−a1 ,α1 ),...,(1−ap ,αp )  (26) 1 1,1 p (1− ν1 , ν1 ) , (27) H 1 1 cν 1,1 c (1− ν , ν ) where the H-function is defined in eq. (7). The Laplace transform of Fox’s H-function (7) is given in terms of another H-function by   N(p) = N0 1 n+1,m h L {H(z); p} = Hq,p+1 p p (1−bq ,βq ) (1,1),(1−ap ,αp ) i (28) for 0 ≤ µ ≤ 1 in (9), and " 1 m,n+1 1 L {H(z); p} = Hp+1,q p p (0,1),(ap ,αp ) (bq ,βq ) # (29) for µ ≥ 1 in (9), respectively. Further, having H(p), the inverse Laplace transform of this H-function is given by i 1 n,m h (1−bq ,βq ) (30) z (1−ap ,αp ),(1,1) H(z) = L−1 {H(p), z} = Hq,p+1 z for 0 ≤ µ ≤ 1 in (9) and 1 m,n h H(z) = L−1 {H(p), z} = Hp+1,q z z for µ ≥ 1 in (9), respectively. The above four transforms hold for  Laplace  max1≤j≤n ℜ a−1 p αp < min1≤j≤m ℜ  bq βq  (ap ,αp ),(0,1) (bq ,βq ) i (31) . Applying an inverse Laplace transform to the H-function in eq. (27) gives 1 1,1 N(t) = N0 H1,2 ct ν  (0, ν1 ) (0, ν1 ),(0,1)  , (32) which is the solution of the fractional kinetic equation (21). For the Hfunction in eq. (7) with (32), the following computable representation can be derived (Mathai and Saxena, 1978). When the poles of Πm j=1 Γ(bj − βj s) are simple, that is, βh (bj + λ) 6= βj (bh + ν) 9 for j 6= h; j, h = 1, . . . , m; λ, ν = 0, 1, 2, . . . . Then one obtains the following expansion for the H-function, m,n Hp,q (z) = ∞ m X X h=1 ν=0 × n n n n (bh +ν) Πm j=1,j6=h Γ bj − βj βh n Πqj=m+1 Γ 1 − bj + βj (bhβ+ν) h n Πnj=1 Γ 1 − aj + αj (bhβ+ν) h n oo n Πpj=n+1 Γ aj − αj (bhβ+ν) h oo oo oo (−1)ν z (bh +ν)/βh , (ν)!βh (33) which exists for all z 6= 0 if µ > 0 and for 0 < |z| < β −1 if µ = 0, where µ and β are given in eqs. (9) and (10). Comparing (32) with eq. (33), one obtains the series expansion N(t) = N0 ∞ X (−1)k (ct)νk . Γ(νk + 1) k=0 (34) For ν = 1, the exponential solution of the standard kinetic equation (18) is recovered. 4. Conclusions The gravitationally stabilized solar fusion reactor, because of its density and temperature, is not a weakly interacting or high temperature plasma. Recently, Kaniadakis et al. (1997) and Coraddu et al. (1998) have explored the possibility that the electrical microfield distribution influences the particle dynamics, collisions are not two-body processes and retain memory beyond single scattering events, and that the velocity correlation function has longtime memory arising from the coupling of collective and individual degrees of freedom. In this connection particle diffusion is a non-Markovian process (anomalous diffusion) and diffusion and frictional coefficients are energy dependent. Based on these results, they conclude that the equilibrium statistical distribution function differs from the Maxwell-Boltzmannian one and is governed by generalized Boltzmann-Gibbs statistics. Even if the deviations from the Maxwell-Boltzmann distribution, that are compatible with the current knowledge of the solar core behavior, are small, they are sufficient to sensibly modify the sub-barrier particle reaction rates and subsequently solar neutrino fluxes. The above authors have also shown 10 that the respective modifications of the reaction rates do not affect bulk properties of the gravitationally stabilized solar core such as sound of speed or hydrostatic equilibrium which depend on the mean values obtained by averaging over the Maxwell-Boltzmann distribution function (Degli’Innocenti et al. 1998). Haubold and Mathai (1998) have derived general closed-form representations of particle reaction rates that are suitable to incorporate changes of the Maxwell-Boltzmann distribution function, mindfully taking into account the fact that the whole distribution has physical meaning, contrary to the case of, say, the Gaussian law of errors. In this paper we proceeded one step further in generalizing the standard kinetic equation (17) to a fractional kinetic equation (21) and derived solutions of a fractional kinetic equation that contains the particle reaction rate (or thermonuclear function) (5) as a time constant, and provided the analytic technique to further investigate possible modifications of the reaction rate through a kinetic equation. The Riemann-Liouville operator in the fractional kinetic equation introduces a convolution integral with slowly-decaying power-law kernel, which is typical for memory effects referred to in Kaniadakis et al. (1997) and Coraddu et al. (1998). This technique may open an avenue to accommodate changes in standard solar model core physics (Dar and Shaviv, 1999) as proposed by Shaviv and Shaviv (1999) and Schatzman (1999). In the solution of the fractional kinetic equation (21), given in eqs. (32) and (34), the standard exponential decay is recovered for ν = 1. However, eqs. (32) and (34), for 0 < ν < 1, show power law behavior for t → ∞ and are constant (initial value N0 ) for t → 0. In the investigations in this paper, (∞) we used, as an example, the standard thermonuclear function I2 = I1 ; the (d) same computations can be done for I2 , I3 , or I4 with additional parameters showing modifications of the Maxwell-Boltzmann distribution function (Haubold and Mathai, 1998). All analytic results in this paper have been achieved by the application of the theory of generalized hypergeometric functions, particularly Meijer’s G-function and Fox’s H-function, which seem to be the natural means for tackling the problems referred to above. 11 Acknowledgments The investigation of formation of structure in a non-equilibrium solar fusion plasma was discussed in-depth on 12 September 1984 with Hans-Juergen Treder at the Einstein Laboratory for Theoretical Physics in Caputh. Shorttime instabilities in solar nuclear reaction kinetics have been explored briefly in a letter of Jean Perdang dated 8 January 1991. The corresponding author of this paper (HJH) is grateful to both colleagues and wishes to put this on record even at such a late date. The same author acknowledges inspiring discussions on non-Markovian processes in a solar nuclear plasma with Piero Quarati at the occasion of the sixth UN/ESA Workshop on Basic Space Science at the Max-Planck-Institute for Radioastronomy in Bonn, 9-13 September 1996. References Bergstroem, L., Iguri, S., and Rubinstein, H.: 1999, Constraints on the variation of the fine structure constant from big bang nucleosynthesis, Phys. Rev. D60, 045005-1 - 045005-9. Clayton, D.D.: 1983, Principles of Stellar Evolution and Nucleosynthesis, Second Edition, The University of Chicago Press, Chicago and London. Coraddu, M., Kaniadakis, G., Lavagno, A., Lissia, M., Mezzorani, G., and Quarati, P.: 1998, Thermal distributions in stellar plasma, nuclear reactions and solar neutrinos, http://xxx.lanl.gov/abs/nucl-th/981108. Dar, A. and Shaviv, G.: 1999, The solar neutrino problem - an update, Physics Reports 311, 115-141. Degli’Innocenti, S., Fiorentini, G., Lissia, M., Quarati, P., and Ricci, B.: 1998, Helioseismology can test the Maxwell-Boltzmann distribution, http://xxx.lanl.gov/abs/astro-ph/9807078. Haubold, H.J. and Mathai, A.M.: 1995, A heuristic remark on the periodic variation in the number of solar neutrinos detected on Earth, Astrophys. 12 Space Sci. 228, 113-134. Haubold, H.J. and Mathai, A.M.: 1998, On thermonuclear reaction rates, Astrophys. Space Sci. 258, 185-199. Kaniadakis, G., Lavagno, A., Lissia, M., and Quarati, P.: 1998, Anomalous diffusion modifies solar neutrino fluxes, Physica A261, 359-373, http://xxx.lanl.gov/abs/astro-ph/9710173. Kourganoff, V.: 1973, Introduction to the Physics of Stellar Interiors, D. Reidel Publishing Company, Dordrecht. Mathai, A.M.: 1993, A Handbook of Generalized Special Functions for Statistical and Physical Sciences, Clarendon Press, Oxford. Mathai, A.M. and Saxena, R.K.: 1973, Generalized Hypergeometric Functions with Applications in Statistics and Physical Sciences, SpringerVerlag, Lecture Notes in Mathematics Vol. 348, Berlin Heidelberg New York. Mathai, A.M. and Saxena, R.K.: 1978, The H-function with Applications in Statistics and Other Disciplines, John Wiley and Sons, New Delhi. Mestel, L.: 1999, The early days of solar structure theory, Physics Reports 311, 295-305. Miller, K.S. and Ross, B.: 1993, An Introduction to the Fractional Calculus and Fractional Differential Equations, John Wiley and Sons, New York. Murray, J.D.: 1989, Mathematical Biology, Biomathematics Texts Vol. 19, Springer-Verlag, Berlin. Nicolis, G. and Prigogine, I.: 1977, Self-Organization in Nonequilibrium Systems - From Dissipative Structures to Order Through Fluctuations, John Wiley and Sons, New York. Oldham, K.B. and Spanier, J.: 1974, The Fractional Calculus, Academic Press, New York. Perdang, J.: 1976, Lecture Notes in Stellar Stability, Part I and II, Instituto 13 di Astronomia, Padova. Schatzman, E.: 1999, Role of gravity waves in the solar neutrino problem, Physics Reports 311, 143-150. Shaviv, G. and Shaviv, N.J.: 1999, Is there a dynamic effect in the screening of nuclear reactions in stellar plasma?, Physics Reports 311, 99-114. 14