Family of Probability Distributions Derived from Maximal Entropy Principle with Scale Invariant Restrictions

0.14 MB·English

by  Giorgio Sonnino

Family of Probability Distributions Derived from Maximal Entropy Principle with
Scale Invariant Restrictions
Giorgio Sonnino1∗, Gy¨orgy Steinbrecher2, Alessandro Cardinali3, Alberto Sonnino4, Mustapha Tlidi1
1 Department of Theoretical Physics and Mathematics, Universit´e Libre de Bruxelles (U.L.B.),
Campus de la Plaine C. P. 231 - Boulevard du Triomphe, 1050 Brussels, Belgium
2 Association EURATOM-MEC, Physics Faculty,
University of Craiova, Str.A.I.Cuza 13, 200585 Craiova, Romania
3 EURATOM-ENEA Fusion Association, Via E.Fermi 45, C.P. 65 - 00044 Frascati (Rome), Italy and
3
4 Universit´e Catholique de Louvain (UCL), Ecole Polytechnique de Louvain (EPL),
1
Rue Archim`ede, 1 bte L6.11.01, 1348 Louvain-la-Neuve, Belgium
0
2
Using statistical thermodynamics, we derive a general expression of the stationary probability

n distributionforthermodynamicsystemsdrivenoutofequilibriumbyseveralthermodynamicforces.
a The local equilibrium is defined by imposing the minimum entropy production and the maximum
J
entropy principle under the scale invariance restrictions. The obtained probability distribution

3 presents a singularity that has immediate physical interpretation in terms of the intermittency
2 models. Thederivedreferenceprobabilitydistributionfunctionisinterpretedastimeandensemble
average of the real physical one. A generic family of stochastic processes describing noise-driven
]
intermittency,wherethestationarydensitydistributioncoincidesexactlywiththeoneresultedfrom
h
c entropymaximization, is presented.
e
PACS Numbers: 05.70.Ln 52.25.Dg 05.20.Dd
m
- *Email: [email protected]
t
a
t
s In the Onsagerregion,the minimum entropy production Ent principle on the random variables, under the scale
.
t (MEP) theorem provides a variational characterization invariance restrictions. By definition, the reference DDF
a
m of a stationary states, both for macroscopicsystems and (indicated with F0) is an initial distribution function.
forstochasticmodels. However,tocharacterizesomeun- F0 shoulddepend onlyonthe invariantsofmotion, with
-
d known events with a statistical model, we should choose the property to evolve slowly from the local equilibrium
n the one that has maximum entropy (MaxEnt’s princi- state i.e., it remains confined for sufficiently long time.
o ple). Thismeansthat,outofallprobabilitydistributions Hence, the reference DDF results in a perturbation of
c
consistent with a given set of constraints, we choose the the local equilibrium state. We consider open thermo-
[
one that has maximum uncertainty. Entropy is there- dynamic systems obeying to Prigogine’s statistical ther-

2 fore regarded as a measure of information. The MaxEnt modynamics. We then define the local equilibrium state
v principledevelopedincommunicationandininformation by adopting a minimal number of hypotheses. Finally,
6 ′
technology has recently been found to have a wide rang- we link the density distribution function with particles
2
ingapplicationsinmanyareasofscience. Inevolutionary DDF. The density probability distribution of finding a
3
biology, the MaxEnt principle is used to explain the ob- state in which the values of the fluctuating thermody-
2
. served species abundance distribution. In this case the namic variable, β˜κ, lies between β˜κ and β˜κ+dβ˜κ is
1
constraintmaybesetbythehabitat,whichfixestheaver-
1 agepopulationsizeofthespecies[1]-[7]. Thisprincipleis F =N0exp[−∆IS] (1)
2
1 alsoappliedfordescribingtheactivityinavarietyofneu- whereN0ensuresnormalizationtounity,andwehavein-
: ralnetworks[8]-[12],thestatisticsofaminoacidsubstitu- troduced the dimensionless (density of) entropy produc-
v
tions in protein families [13]-[16] and in material science tion∆ S. The negative signin Eq.(1) is due to the fact
i I
X [17]andinBose-Einsteincondensateinadyemicrocavity that, during the processes, −∆ S ≤0. Indeed, if −∆ S
I I
r [18]. ItisalsoproventhattheMaxEntmodelcanprovide were positive, the transformation β˜κ → β˜κ′ would be a
a agooddescriptionofdatafromseeminglycompletelyun- spontaneous irreversible change and thus be incompati-
related phenomena observed in social contexts [19], [20]. ble with the assumption that the initial state is a stable
For example, it describes very well the equilibrium dis- (local) equilibrium state [22]. We suppose that the sys-
tribution function for the city-population subject to two tem is subject to N˜ thermodynamic forces. The entropy
additional scale-invariant constraints: the normalization production ∆ S is linked to the thermodynamic forces
I
ofthe probabilitydistributionfunctionandthe expected Xκ, the thermodynamic fluctuation β˜ ,andthe thermo-
κ
value of the city-population [21]. dynamic flows J by the following set of equations [23]
κ
The purpose of this paper is to use statisticalthermody-
namicstoderiveareferencedensitydistributionfunction ∂∆ S dβ˜ d S N˜
Xκ = I , J = κ , I = XκJ ≥0 (2)
(DDF) for open thermodynamic systems close to a local ∂β˜ κ dt dt k
equilibrium state. To this end, we assume the validity of κ κX=1
theminimumentropyproductionprincipleandtheMax- Our aim is firstly to derive the reference DDF, F0, by
2
imposing the MEP and the MaxEnt principle, submit- (equilibrium) values of the remaining phase space vari-
ted to the scale invariance constraints. Successively we ables, is expressed by the gamma distribution. Notice
shall determine the class of stochastic processes whose that if we normalize P(w) to unity then it can be in-
stationary PDF includes F0 as a special case. terpreted as the conditional probability density function
Let us consider an open thermodynamic system subject (PDF), ρ(w), conditioned by β = 0 with κ = 1,···N
κ
to N˜ =N+1thermodynamicforces. N thermodynamic i.e., ρ(w)≡P(w |β1 =···=βN =0).
′
forces are linked to N fluctuations of Prigogines type We consider the entropy S[ρ(.)] of probability density
(i.e., the entropy production is expressed in quadratic function (PDF), ρ(w)≥0, given by
form with respect to these fluctuations. For an exact
′ ∞
definition of Prigogines fluctuations refer to [22], [24].
S[ρ(.)]=− ρ(w)log(ρ(w))dw (6)
One thermodynamic force is linked to a fluctuation of a
Z0
different nature. This random variable has the density
distribution function, which satisfies the MaxEnt prin- We investigate the consequences of the MaxEnt princi-
ciple under the scale invariant restrictions. This fluc- ple submittedto the mostgeneralscale-invariantrestric-
tuation will be indicated with w whereas the Prigogine tions. Let us then start to consider the following restric-
fluctuations will be denoted by β , with κ = 1,··· ,,N. tions
κ
Hence, β˜κ = βk for k = 1,··· ,N and β˜N+1 = w. To be ∞
more precise, we identify the local equilibrium state by wαkρ(w)dw =E(wαk)=µ ; k =0,1,...,n (7)
k
imposing the following two conditions. Z0
i) Thelocalequilibriumstatecorrespondstothevalues where, at this stage, αk, and µk are real numbers. How-
of the fluctuations β with κ =1,··· ,N for which ever,forcompleteness,weshouldalsoinclude thefollow-
κ
the entropy production tends to reach an extreme. ing restriction obtained as the limit case of (7)
Underthisassumption,closetothelocalequilibrium,the E(log(w))=ν (8)
entropy production can be brought into the form
where ν is a real number. Indeed, suppose that we have
1 N for some fixed k : αk =ε≪1, then
−∆IS =g0(w)− gκj(w)βκβj +h.o.t. (3)
2
κ,j=1 E(wε)=µ (9)
X k
with g0(w)≡−∆IS |β1···βN=0
From Eqs (7) and by taking into account the normal-
where h.o.t. stands for higher order terms. Coefficients izationconditionon the probability [see the forthcoming
g aredirectlylinkedtothe transportcoefficientsofthe Eq. (11)], we get
κj
system [25]. Therefore, the general expression for F,
given by Eq. (1), becomes the reference DDF, F0, when wε−w0 µk−1
E = (10)
theentropyproductionisprovidedbyEq.(3). TheDDF ε ε
(cid:18) (cid:19)
P(w)≡N0exp[g0(w)]=N0exp[−∆IS |β1···βN=0] (4) If the support of the PDF ρ(w) is concentrated mainly
on the domain where |log(w)| is not too large, then we
related to the variable w, at βκ = 0 (with κ = 1···N), can approximate: wε−w0 /ε ∼= log(w). In this case
is determined by the following condition.
Eq. (10) reduces to Eq.(8). Finally, Eq.(8) can be seen
(cid:0) (cid:1)
as the limit case of (7) and, for the sake of generality,
ii) Attheextremizingvalues β =0withκ=1,··· ,N
κ
it should be taken into accountamong the equations ex-
under the scale invariance restrictions, the system
pressingthescale-invariancerestrictions. Itisworthwhile
tends to evolve towards the maximal entropy con-
mentioning that the equilibrium distribution function,
figurations.
obtainedbyimposingrestriction(8),retrievestheexpres-
Notice that in particular, if the imposed the two scale sioncurrently used for fitting the numericalsteady-state
invariant restrictions, E[w] = const. > 0 and E[ln(w)] = solution of the simulations for the Ion Cyclotron Radia-
const.(whereE[]istheexpectationoperation),together tion Heating (ICRH) FAST-plasmas, and for describing
with normalization, the density distribution for the w various scenarios of tokamak-plasmas [25]. Another ex-
variable is given by a gamma distribution function [26] ample showing the usefulness of constraint (8) can be
found at the end of this report. Let us now consider the
P(w)=N0wγ−1exp(−w/Θ) (5) following simplest physical case
Here, we have introduced the scale parameter Θ and the ∞
shape parameter γ, which could be determined by the ρ(w)dw =1 (11)
fit of experimental data. The goal is now to justify a Z0
∞
posteriori the result found in Ref. [25] : The kinetic en- wρ(w)dw =E(w)=µ1 (12)
ergy dependence of the DDF, P(w), at the most probable Z0
3
Hence, n≥1 and, in particular It is clear that in this simplest and extremal case, given
by the ansatz in the Remark (1), the PDF has the form
µ0 ≡1; α0 =0;α1 =1 (13)
where µ1 = T/m, with T and m are temperature and ρ(w)= ΘΓ1(γ)(w/Θ)γ−1exp(−w/Θ) (16)
particle’s mass, respectively.
TounderstandthemeaningofEqs(7)and(8),forn=1,
where Γ(γ) is the Euler Gamma function. The relations
we find Eq. (5) from the MaxEnt principle submitted to
betweenthe parametersγ,Θ, µ1 =T/mandν aregiven
the constraints (8), (11) and (12). In the following, we
by
shallbeabletogiveonlyapartialanswertothisquestion.
TherestrictionsonthePDFgivenbyEq.(7)andEq.(8) E(w)=T/m=γΘ (17)
are invariant under scale transformations. This means
E(log(w))=ν =Ψ(γ)+log(Θ) (18)
′
thattheeffectofthescaletransformationw→kw onthe
restrictions(7)and(8)isE(w′αk)=µ =µ k−αk and
αk αk where Ψ(γ) is the digamma function. We draw the at-
′ ′
E(log(w ))=ν =ν−log(k),respectively. Sorestrictions
tention to the fact that, in general, the parameters ap-
(7) and (8) remain in the same class specified by the
pearing in the Eqs (7) and (8) are not completely inde-
functions wαk and log(w), respectively. pendent. Indeedsupposethatwehave E(wα1)=µ1 and
Remark (1). It remains to explain the ansatz: n = 1, E(wα2) = µ2 with α1 < α2. Then, by using the Ho¨lder
and Eq. (8) [naturally combined with restrictions (11)
inequality [27], we obtain:
and (12)].
By denoting with λk,0 ≤ k ≤ n+1 the Lagrange mul- µ1/α1 ≤µ1/α2 (19)
tipliers in the problem of maximizing the entropy given 1 2
by Eq. (6), with the restrictions (7) and (8), we get
Inthelimitcase,byusingEqs(8)and(12),togetherwith
n Jensen’s inequality for the exponential function [27], we
log[ρ(w)] =1− λkwαk −λn+1log(w) (14) find
kX=0 µ1/αk ≥exp(E(log(w)))=exp(v);k ≥1 (20)
k
It is convenient to introduce the following parametriza-
tion of the general PDF: The inequalities (19) and (20) are the necessary condi-
tions for the existence and not the necessary conditions
n
ρ(w) = 1wγ−1exp − λ wαk (15) for maximalentropy ofthe PDF. The inequality (19) re-
k
Z sults from the constraint (7) whereas (20) from the con-
" #
k=1
X straints (7) and (8). In particular, for restrictions (11)
withγ−1=−λn+1 and1/Z =exp(1−λ0). The param- and (12) we get µ1 ≥exp(v).
eters Z, γ, λ are fixed by Eqs (7) and (8). We should Due to the very special choice of the restrictions used in
k
have λ > 0 when α is the largest exponent, as well as theMaxEntprinciple(the minimalityofthelog(w)), the
k k
when α is the least negative exponent. resulting PDF is expected to have some special proper-
k
In Eq.(15), despite the condition α <0, with λ >0, is ties. Indeed, the gamma distribution is infinitely divisi-
k k
mathematicallyacceptable,weshouldexcludethispossi- ble and stable in the following sense. Let X1,...,Xn be
bility fromthephysicalpointofview. Becauseinseveral a sequence of n independent random variable having n
situations,like with collisionalplasmas,sucha condition identical gamma distributions (with the same scale pa-
would provide a very reduced population at low energy. rameter and shape parameters γ1,...γn), then the sum
Similarly, also the alternative α > 1, with λ > 0, for n X has the same gamma distribution function,
k k k=1 k
some value of k, should be discarded. This because such with the same scale parameter and with shape param-
a conditionwouldgiveriseto asub-Maxwellian distribu- Peter equal to n γ . An important physical property
k=1 k
tion function when, on the contrary, it is expected that of the distributions (15), in the case 0 < γ < 1, is the
P
thehighenergytailofthePDFshouldbelargerthanthe following. By varying γ in this range of parameters,and
one predicted by the Maxwellian distribution (consider, by keeping constant the other parameters, the shape of
for instance, the case of burning plasmas). the graph of the DDF changes, due to the singularity
Conjecture (2). Physical considerations lead us to se- appearing at w = 0. This singularity can be associated
lect max(α1,··· ,αn =1, and 0),αk <1 if k >2. withtheappearanceofaregimeofintermittentbehavior
Even though we accept the condition α < 0, the case [29]. WeemphasizethatthereferenceDDFisinterpreted
k
α → 0, should anyhow be considered as a limit situa- as time and ensemble averageof the physical, DDF.
k
tionbecausethegraphofρ(w)changesabruptlywiththe We describe a class of stochastic processes admitting
changeofthesignofα . Fromthisreason,therestriction Eq.(15)asstationaryPDFsolutions. NotethatEq.(15)
k
(8) can be considered as an extreme case. Consequently reducestothegammadistributionwhensomecoefficients
the partialjustificationofthe ansatz,inthe Remark(1), are set to zero. To this end, we consider the stochas-
is given for simplicity and extremality reasons. The ex- tic differential equation (SDE) for the random variable
ploration of the case n = 2, with 0< α2 < 1 is ongoing. w(t). In this case w(t) is the energy of an individual
4
charged particle. The SDE includes the simplest soluble By comparing Eq. (28) with Eq. (15), we can link the
cases of the class of intermittency models [28]-[30]. In exponentsα withtheexponentsξ appearinginEq.(27).
k k
the Stratonovich version, the SDE reads Notice that, by setting the parameters entering in the
Eqs (21) and (27) as
dw(t)=(adt+σdB(t))◦w(t)−S[w(t)]dt (21)
2a
where◦standsforthe Stratonovich product. Inaddition, γ = σ2 =γ and αk =ξk (29)
a > 0 is the instability threshold, B(t) is the standard
Brownian motion (Wiener’s process) and σ is the inten- weobtain,byintermittencemechanism,exactlythesame
sity of the multiplicative noise. The function S(w) is as- stationary distribution function derived by the MaxEnt
sociated with the saturation of the instability controlled principle under the scale invariant restrictions. In ad-
by the linear term. Hence, we should have dition, we have λk = Ak/ξk. Let us now consider the
standard Landau type generic form for modeling the in-
S(w)
lim =+∞ (22) stability growth. In the particular case of the simplest
w→∞ w choice S(w) = Aw2, we obtain the PDF from Eq.(5)
Wealsorequirethatthesolutionnearw =0isdominated resulting from the Ansatz in the Remark (1). The high
by the linear term. So the phenomenology described by energytailisinagreementwiththeMaxwelldistribution.
Eq. (21) is still related to the noise-drivenintermittency However, for the case of the Landau term S(w) = Aw3,
if we require that we obtain a distorted Gaussian distribution
S(w) ρ(w)=Cwγ−1exp −Aw2 (30)
lim =0 (23)
w→0+ w
where the high energy tail decrea(cid:0)ses too(cid:1)quickly. This
We will see that, in the particular case S(w)=Aw2, the behaviour is incompatible with many physical situations
stationarysolutionofEq.(21)isthegammadistribution. (such as, for example, the case of collisional plasmas)
Hence,theclassofEqs(21)includesthegenericfamilyof and, for this reason,this possibility should be discarded.
equations describing the instability growth, on the pos- Instead,inaccordancewithConjecture(2),ifinthelead-
itive semi-axes (which corresponds to our case), limited ing terms Aκw1+ξκ, we have ξ1 = 1 and for the rest
by the saturation term. We have slightly modified this 0<ξ <1,the highenergytailwoulddecreasein agood
k
equationbyaddingtherandommultiplicativenoiseterm way. Hence, Conjecture (2) is supported by the inter-
σdB(t). In the Itoˆ formalism we get mittency models. Notice that in the special case n = 0,
′ inEqs.(15, 28) the resulting DDF is non-integrable,but
dw(t)=(adt+σdB(t))w(t)−S[w(t)]dt
self similar, like in the case of self organized criticality
′ σ2 (SOC) models [31], [32]. So restriction given by Eq. (8)
with a =a+ (24)
2 is a remnant SOC, avalanche-like behavior. By replac-
ing the self similarity of the DDF by scale invariance
The stationary Fokker-Planck equation for the density
of the constraints in the MaxEnt principle, we obtain a
distribution ρ(w) reads
class of DDF that better parametrize the intermittent,
∂ σ2 ∂2 avalanche- like behaviour [31], [33] of the physical sys-
[(a′w−S(w))ρ(w)]− w2ρ(w) =0
∂w 2 ∂w2 tems that can be modeled by self organized criticality
models [34].
admitting, up to a normalization const(cid:2)ant, the(cid:3)following
We conclude this brief report by summarizing the main
steady state solution
results. Using statistical thermodynamics, we derive the
ρ(w)=Cwγ−1exp − S(w)dw (25) general expression of the reference PDF, F0, by impos-
w2 ingthe minimumentropyproduction,andthemaximum
(cid:20) Z (cid:21) entropy principle under the scale invariance restrictions.
2a
γ = σ2 >0 (26) Some restrictions on the parameters entering in F0 are
obtained by physical arguments. The obtained PDF ex-
The generalformofS(w), compatible withEqs (22) and hibits a singularity, which has immediate physical inter-
(23), is pretation in terms of the intermittency models. Finally,
m we deriveda family of stochastic processes admitting F0
S(w)= A w1+ξk with ξ >0 (27) as stationary PDF solutions.
k k
This work gives several perspectives. Through the ther-
k=1
X modynamical field theory (TFT) [35] it is possible to es-
with the constraint that at infinity the coefficient of timatethePDFwhenthenonlinearcontributionscannot
the leading term in Eq.(27) should be positive. From be neglected [36]. The next task should be to establish
Eqs (25) and (27) we obtain the relation between the reference PDF herein derived
with the one found by the TFT.The solutionof this dif-
m
A
ρ(w)=Cwγ−1exp − kwξk (28) ficult problem will contribute to provide a link between
" k=1 ξk # a microscopic description and a macroscopic approach
X
5
(TFT). Another problem to be solved is the possibility Universit´e Libre de Bruxelles, for his scientific sugges-
to improve the numerical fit by adding new free param- tions and for his help in the development of this work.
eters according to the principles exposed in this report. M.Tlidi is a Research Associate with the Fonds de la
Recherche Scientifique F.R.S.-FNRS, Belgium.
Acknowledgements
G. Sonnino is very grateful to M.Malek Mansour, of the
[1] B. Shipley, D. Vile, E. Garnier, Science, 314, 812-814 Plastino, Phys. Rev. E,86,066105 (2012).A.Hernando,
(2006). A. Plastino, Eur. Phys. J. B, 85, 293 (2012).
[2] S. Pueyo, F. He, T. Zillio, Ecol. Lett., 10, 1017-1028 [22] I.Prigogine,(1954),ThermodynamicsofIrreversiblepro-
(2007). cesses, (John Wiley & Sons).
[3] J. Harte, T. Zillo, E. Conlisk, A.B. Smith, Ecology, 89, [23] S.R. De Groot and P. Mazur, (1984), Non-Equilibrium
2700-2711 (2008). Thermodynamics, (DoverPublications, Inc.,New York).
[4] J.R. Banavar, A. Maritan, I. Volkov, J. Phys.Condens. [24] I. Prigogine, (1947), Etude Thermodynamique des
Matter, 22, 063101 (2010). Ph`enom`enes Irr´eversibles, Th`ese d’Aggr´egation de
[5] B. Haegeman, R.S. Etienne, Am. Nat., 175, E74-E90 l’EinseignementSup´erieurdel’Universit´eLibredeBrux-
(2010). elles (U.L.B.).
[6] F. He,Oikis, 119, 578-582 (2010). [25] G. Sonnino,A.Cardinali, M.Tlidi, P.Peeters, G.Stein-
[7] S.A.Frank,J. Evol. Biol., 24, 485-496 (2011). brecher, A. Milovanov, General Approach for Deriv-
[8] J. Shlenset Al.,J. Neurosci, 26, 8254-8266 (2006). ing Reference Distribution Functions for Systems out of
[9] S.Yu,D.Huang,W.Singer,D.Nikolic,CerebCortex,18, EquilibriumbyStatisticalThermodynamics,submittedto
2891-2901(2008). Plasma Physics and Controlled Fusion (PPCF), (2012).
[10] A.Tang et Al.,J. Neurosci, 28, 505-518 (2008). [26] A. Papoulis and S.U. Pillai, (2002), Probability, Ran-
[11] I.E. Ohiorhenuan et Al.,Nature, 466, 617-621 (2010). domVariablesandStochastic Processes,(MacGrawHill,
[12] E. Ganmor, R. Segev and E. Schniedman, Proc. Natl. Fourth edition).
Acad. Sci. (USA), 108, 9679-9684 (2011). [27] W.Rudin,(1987),RealandComplexAnalysis,(McGraw-
[13] F. Seno, A. Trovato, J.R. Banavar, A. Maritan, Phys. Hill).
Rev. Lett., 100, 078102 (2008). [28] A.Schenzle,H.Brand,Phys.Rev.A,20, 1628, (1979).
[14] M.WeigtetAl.,Proc.Natl.Acad.Sci.(USA),106,67-72 [29] S. Aumaˆıtre, F. P´etr´elis, K. Mallick, Phys. Rev. Lett.,
(2009). 95, 064101, (2005).
[15] N.Halabi et. Al.,Cell.,138, 774-786 (2009). [30] S.Aumaˆıtre,K.Mallick,F.P´etr´elis,J. Stat. Phys.,123,
[16] J.I. Sulkowskaet Al.,Proc. Natl. Acad. Sci (USA), 109, 909, (2006).
10340-10345 (2012). [31] D. Turcotte, Rep. Progr. Phys., 64 , 1377, (1999).
[17] Y. Nakazawa at Al., J. Phys. Soc. Jpn., 81, 094601 [32] P. Bak, How Nature Works: The Science of Self Orga-
(2012). nized Criticality , New York: Copernicus, (1996).
[18] D.N. Sob’yanin,Phys. Rev. E, 85, 061120 (2012). [33] V. P. Budaev,Phys. Lett., A 373, 856, (2009).
[19] S.K. Baek, S. Bernhardsson and P. Minnhagen, New J. [34] A. Gruzinov, Phys. Rev. Lett., 89, 25501, (2002).
Phys., 13, 043004 (2011). [35] G. Sonnino, Phys. Rev. E.,79, 051126, (5Pt 1),(2009).
[20] A. Hernando, A. Plastino, A.R. Plastino, Eur. Phys. J. [36] G. Sonnino, The European Physical Journal D (EPJD),
B, 85, 147 (2012). 62, Issue 1, 81, (2011).
[21] A.Hernando,R.Hernando,A.Plastino,A.R.Plastino,J.
R. Soc. Interface, 78,20120758 (2013). A.Hernando,A.