text:
# Invariance of density correlations with charge density in polyelectrolyte solutions

Present Address: \]Sandia National Laboratories, Albuquerque, NM 87185

## Abstract

We present a theory for the equilibrium structure of polyelectrolyte solutions. The main element is a simple, new optimization scheme that allows theories such as the random phase approximation (RPA) to handle the harsh repulsive forces present in such systems. Comparison is made with data from recent neutron scattering experiments of randomly charged, hydrophilic polymers in salt-free, semi-dilute solution at various charge densities. Models with varying degrees of realism are examined. The usual explanation of the invariance observed at high charge density has been counterion condensation. However, when polymer-polymer correlations are treated properly, we find that modeling polymer-counterion correlations at the level of Debye-Hückel theory is sufficient.

preprint: APS/123-QED

Polyelectrolytes are polymers with ionizable groups that dissociate in polar solvent leaving ions bound to the chain and counterions free in solution. The study of polyelectrolytes has traditionally been driven by their omnipresence in living things moleccellbook , but there has been rapidly growing use of and interest in synthetic ones for technological applications, typically at the nano-scale level applications . The theoretical study of them is one of the most challenging in the field of complex fluids, however, because of the strong repulsions and attractions present in such systems forster95 .

In this Letter, we present a scheme that allows a potentially large class of theories to deal properly with these strong long-range interactions. We demonstrate the utility of this approach by comparing with some recent intriguing experiments of Nishida, Kaji and Kanaya (NKK) nishida95 , and Essafi, Lafuma and Williams(ELW) essafi99 . Counterion condensation oosawamanning  commonly is invoked as the cause of the invariance at high charge density observed in these experiments. While there is some ambiguity in the literature about what precisely is counterion condensation, it is generally regarded that the mechanism is not captured in any theory in which ion-counterion correlations are included at the level of Debye-Hückel theory oosawamanning ; ccrefs . However, we show in the theory described below that Debye-Hückel modeling of ion-counterion correlations can itself produce an invariance at high charge density.

There are few theories that can handle the harsh repulsions and attractions in polyelectrolyte solutions deshkovski01 . Further, there is currently no theory that provides reliable predictions for liquid structure in the semi-dilute regime—which is the region of most experimental interest. The theory of Laria, Wu and Chandler laria91  has been applied successfully to examine density correlations there as a function of number density shew96 , but performs poorly at large, experimentally relevant, interaction energies and charge densities dhwunpub . The theory of Donley, Rajasekaran and Liu (DRL) seems to perform properly at large interaction energies and charge densities, but presently only for purely repulsive systems donley98 . Despite its simple form, the random phase approximation (RPA) can in some instances give reasonable predictions for the free energy of polyelectrolyte solutions rpafreeenergy . Here, we assume that theories such as the RPA are also adequate for density correlations given a suitable reinterpretion of the interaction potential that appears in them andersen71 . Let

$$S_{MkM^{}k^{}}(r)=\left[\widehat{\rho }_{Mk}(𝐫)\rho _{Mk}\right]\left[\widehat{\rho }_{M^{}k^{}}(0)\rho _{M^{}k^{}}\right]$$
(1)
be the density-density correlation function between a $`k`$ type monomer (or site) on a molecule of type $`M`$ and a $`k^{}`$ type monomer on a molecule of type $`M^{}`$ a distance $`r=|𝐫|`$ apart. Here, $`\widehat{\rho }_{Mk}(𝐫)`$ is the microscopic density of $`k`$-type monomers on molecules of type $`M`$ at position $`𝐫`$ with average value $`\rho _{Mk}`$. The brackets denote a thermodynamic average. The Fourier transform of $`S_{MkM^{}k^{}}(r)`$ is the partial structure factor and links the theory to scattering experiments. The correlation function can be separated into $`\mathrm{𝑖𝑛𝑡𝑟𝑎}`$ and $`\mathrm{𝑖𝑛𝑡𝑒𝑟}`$molecular pieces, $`\mathrm{\Omega }_{Mkk^{}}(r)`$ and $`H_{MkM^{}k^{}}(r)`$:

$$S_{MkM^{}k^{}}(r)=\mathrm{\Omega }_{Mkk^{}}(r)\delta _{MM^{}}+H_{MkM^{}k^{}}(r),$$
(2)
where $`\delta _{MM^{}}`$ is the Kronecker delta. In this work, the molecular structure function $`\mathrm{\Omega }_{Mkk^{}}(r)`$ is assumed known. The intermolecular correlation function $`H_{MkM^{}k^{}}(r)=\rho _{Mk}\rho _{M^{}k^{}}\left[g_{MkM^{}k^{}}(r)1\right]`$, where $`g_{MkM^{}k^{}}(r)`$ is the radial distribution function. Thus, to compare with experiment we need a form for this function in terms of the potentials and molecular structures. Consider for simplicity a system where all the interactions can be decomposed as pair potentials $`u_{MkM^{}k^{}}(r)`$ between sites. Here, we take these potentials to be hard-core for distances $`r`$ less than some range $`\sigma _{MkM^{}k^{}}`$, and Coulombic outside. Consider also a class of theories for which $`g_{MkM^{}k^{}}(r)`$ is a function only of the interactions $`\{u\}`$ and the total molecular structure functions $`\{\mathrm{\Omega }\}`$. Theories in this class include the RPA and possible higher loop improvements, and the polymer version melenkevitz97  of the theory of Chandler, Silbey and Ladanyi chandler82 . One failure of theories such as the RPA is that they predict a negative $`g_{MkM^{}k^{}}(r)`$ for short distances if the interaction $`u_{MkM^{}k^{}}(r)`$ is strongly repulsive even though the radial distribution function is intrinsically non-negative. This shortcoming is thought to be one reason that RPA underestimates the strength of the dependence of $`q_{max}`$ with $`\rho _m`$ in the semi-dilute regime.

For short-range repulsive, e.g., hard-core, potentials one remedy to this deficiency was proposed by Andersen and Chandler(AC) andersen71 . Their idea was to replace the true potential $`u`$ by a pseudo or optimized one, $`\stackrel{~}{u}_{MkM^{}k^{}}(r)`$. For distances $`r<\sigma _{MkM^{}k^{}}`$, $`\stackrel{~}{u}_{MkM^{}k^{}}(r)`$ takes a value such that the hard-core condition is satisfied, i.e., $`g_{MkM^{}k^{}}(r)=0`$. Outside the core, $`\stackrel{~}{u}_{MkM^{}k^{}}(r)`$ can be set equal to the original potential $`u`$ which is long-ranged and presumably weak. In this way, higher order diagrams neglected in the theory are summed in an approximate manner to enforce the core condition. When applied to the RPA, this optimization scheme is a generalization of the MSA closure of the Orstein-Zernike equation lebowitz66  and a closely related scheme has been shown to be diagramatically proper for CSL theory melenkevitz97 . This optimized RPA or ORPA, usually with long-range potentials absent, is more popularly known as RISM for small molecules chandler82b  and PRISM for polymers curro87 .

Unfortunately this AC optimization is not very useful for polyelectrolytes. One reason is that the contact energy between ion and counterion seems to be less important than the interaction energy of the counterion with the whole chain. Hence, enforcing the hard-core behavior between opposite charges matters less than modeling correctly the chain structure, e.g., determining the bond length $`b`$. A more important reason is that the repulsion between like-charged polymers is usually very large. This causes the radial distribution function to be effectively zero out to a distance that scales with the Debye-Hückel screening length for semi-dilute densities which is usually much larger than the hard-core distance shew96 ; donley98 ; canessa91 . Hence, enforcing the core condition produces little improvement in the theory for long polyelectrolytes.

An alternative scheme is to optimize the $`\mathrm{𝑟𝑎𝑛𝑔𝑒}`$ of the pseudo-hard-core portion of $`\stackrel{~}{u}`$, and not just its amplitude. By this we mean that if $`g(r)`$ is nearly zero out to some distance, then it makes no difference whether this exclusion zone is caused by a hard-core potential or a Coulombic one. This effective hard-core diameter $`\sigma _{MkM^{}k^{}}^{eff}`$ is determined by requiring that $`g_{MkM^{}k^{}}(r)`$ be non-negative everywhere (not just zero inside the core). The closure to the theory then is:

$$\begin{array}{ccc}g_{MkM^{}k^{}}(r)=& 0,& r<\sigma _{MkM^{}k^{}}^{eff},\\ \stackrel{~}{u}_{MkM^{}k^{}}(r)=& u_{MkM^{}k^{}}(r),& r>\sigma _{MkM^{}k^{}}^{eff},\end{array}$$
(3)
and the range $`\sigma _{MkM^{}k^{}}^{eff}`$ of the hard-core interaction is chosen to have the smallest value such that $`g_{MkM^{}k^{}}(r)>0`$ for $`r>\sigma _{MkM^{}k^{}}^{eff}`$ subject to the constraint that $`\sigma _{MkM^{}k^{}}^{eff}\sigma _{MkM^{}k^{}}`$. With guesses for the $`\sigma _{MkM^{}k^{}}^{eff}`$’s, one can solve the RPA (or similar theory) numerically in a manner very similar to what is done in PRISM. One then adjusts the $`\sigma _{MkM^{}k^{}}^{eff}`$ by iteration. Eq. (3) is the primary result of this paper.

For simplicity, we show results for the range optimized version of RPA, denoted as RO-RPA. We compare with scattering data of NKK and ELW. In these experiments, the molecules were hydrophilic, sulfonated vinyl polymers that were linear and randomly charged, with monovalent ions and counterions. Small angle X-ray (SAXS) and neutron (SANS) scattering methods were used to examine the equlibrium structure of polyelectrolytes and counterions in their salt-free solution, water. Previous experiments on polyelectrolytes have examined how density correlations in the liquid change with polymer monomer number density $`\rho _m`$. Here, though, $`\rho _m`$ was held constant at a semi-dilute or near semi-dilute density and the per chain average fraction of charged monomers $`f`$ was varied. Batches of polymers with different values of $`f`$ were created and the scattered intensity $`I(q)`$ as a function of wavevector $`q`$ was measured for each. NKK found that at small $`f`$, the peak position of $`I(q)`$, $`q_{max}`$, obeyed a power law, but at larger $`f0.4`$ appeared to reach an asymptote. ELW extended these measurements to higher $`f`$ and found that $`q_{max}`$ was effectively constant for $`f>0.4`$. SANS measurements also allowed them to extract the polymer monomer-monomer structure factor $`\widehat{S}_{mm}(q)`$. They found that at least for the range $`0.55<f<0.81`$, $`\widehat{S}_{mm}(q)`$ was invariant for the wavevectors measured. A picture that ELW offer is that as one increases $`f`$ above $`0.4`$, the system acts as if the charges on the chain are renormalized such that the effective $`f`$ is constant. Counterion condensation is usually given as the explanation.

The neutron scattering data of ELW on poly-AMAMPS is sufficiently well characterized to allow quantitative comparison via $`\widehat{S}_{mm}(q)`$ essafi99 . We show results for $`\widehat{S}_{mm}(q)`$ for three models, the “minimal”, “primitive” and “two-state”. In all models the polymer chains are linear and rod-like. The chains are uniformly charged and all monomers are identical so that the monomer valency $`Z_m=f`$, which can be varied continuously from $`0`$ to $`1`$. The strength of the system interactions is characterized partly by the ratio of the Bjerrum length $`l_B=e^2/(ϵk_BT)`$ to the polymer bond length $`b`$. Here, $`e`$, $`ϵ`$ and $`k_BT`$ are the electron charge, solvent (water) dielectric constant and thermal energy, respectively. In the models here, the effect of the salt-free solvent enters only through $`ϵ`$. Both monomer and counterion are spherical and we set their diameters, $`\sigma `$, equal to each other and equal to the bond length $`b`$. For ELW, $`ł_B7.1\AA `$, $`b2.5\AA `$ and $`\rho _m\sigma ^33\times 10^3`$. The counterion density $`\rho _c`$ was determined by charge neutrality, $`Z_m\rho _m+Z_c\rho _c=0`$ and counterion valency $`Z_c=1`$. The chain length N was set to 500. In the minimal model the system consists only of linear polymers. Monomers on the same chain or different chains interact with a potential which is hard-core for $`r<\sigma `$. Outside the core the potential has a screened, Debye-Hückel form: $`u_{mm}(r)=\frac{Z_m^2l_B}{r}e^{\kappa r}`$, where the inverse screening length $`\kappa =\sqrt{4\pi Z_c^2l_B\rho _c}`$. In the primitive model, the system consists of linear polymers and counterions. Polymer monomers and counterions interact with a potential which is also hard-core for $`r<\sigma `$. Outside the core the potential has a bare Coulomb form: $`u_{ij}(r)=\frac{Z_iZ_jl_B}{r}`$, where $`i`$ or $`j`$ denote monomers or counterions.

Comparison of the predictions of RO-RPA in the primitive model with simulation data stevens95  shows that RO-RPA tends to underestimate the strength of polymer-counterion correlations dhwunpub . The two-state model is a common way to improve these correlations oosawamanning ; rpafreeenergy . In this model counterions are divided into two species, free and condensed. In one version of the model, free ones are counterions as in the primitive model, and condensed ones are bound to the surface of the chain, but able to translate along its length. Each chain has the same number of condensed counterions, $`N_c`$. An average $`N_c`$ is determined by constructing a free energy for the system and then minimizing it with respect to $`N_c`$. For a given $`N_c`$ we computed the intramolecular structure functions $`\mathrm{\Omega }_{Mkk^{}}(r)`$ for the polymer-condensed counterion “molecule” by a single chain Monte Carlo simulation. Here, the intramolecular effective interaction between chain monomers and condensed counterions was taken to have a screened Debye-Hückel form similar to above, but the inverse screening length $`\kappa =\sqrt{4\pi l_B(Z_m^2\rho _m+Z_c^2\rho _c)}`$ was due to both ions and counterions. We calculated the free energy using the “charging” formula chandler82b , assuming that the molecular structure was constant during the charging. This constancy is not correct obviously; however, the goal here is not a value for the free energy, but the position of its minimum with respect to $`N_c`$. We expect for the minimum that this approximation is a reasonable one.

Figures 1 and 2 show results for the peak position $`q_{max}`$ and peak height $`\widehat{S}_{mm}(q_{max})`$ of the monomer-monomer partial structure factor. The solid, dashed and dashed-dotted lines correspond to results of RO-RPA using the minimal, primitive and two-state models, respectively. Results from the RPA are also shown as dotted lines. Diamonds are data of ELW. At small $`f`$, both quantities seem to obey power laws: $`q_{max}f^\alpha `$ and $`\widehat{S}_{mm}(q_{max})f^\eta `$. All RO-RPA models predict $`\alpha 2/3`$ and $`\eta 0.9`$. For flexible hydrophilic polyelectrolytes $`\alpha `$ is expected to be $`1/3`$ nishida95 ; pfeuty78 . The difference is due probably to our modeling the chains as rods while they are predicted to be rod-like only for a restricted range of lengthscales pfeuty78 . We know of no data for $`\eta `$ for polyelectrolytes. At large $`f0.4`$, it is clear that range optimization completely changes the character of the RPA theory. In contrast to the RPA, RO-RPA in the minimal model predicts that both quantities become pretty much constant, in agreement with NKK nishida95  and ELW essafi99 . Note that adding explicit counterions and condensed counterions produces only moderate improvement in the theoretical trends comment1 . As a check on the validity of the RO-RPA theory, we also show results (dot-dot-dashed lines) of the DRL theory donley98  in the minimal model. DRL is an approximation to the “two-chain” equation for $`g(r)`$ laria91 ; donley94 . As can be seen, this theory also produces an invariance at large $`f`$ and agrees almost quantitatively with RO-RPA at all values of $`f`$. Agreement between RO-RPA and experiment for $`\widehat{S}_{mm}(q_{max})`$ is very good with the former predicting values less than $`5\%`$ lower than the latter. Agreement for $`q_{max}`$ is less satisfactory with RO-RPA predicting values about $`35\%`$ higher than experiment. The most probable cause of this discrepancy is our neglect of chain flexibility since the Debye-Hückel screening length $`\xi 5\sigma `$ at large $`f`$ implying a flexible chain on scales larger than that. If the chains were partially collapsed, then the repulsion and average distance between monomers on different chains would be larger. Thus $`q_{max}`$ would decrease and $`\widehat{S}_{mm}(q_{max})`$ would increase. On the other hand, if the solvent—which is “good” for poly-AMAMPS—were included explicitly then $`\widehat{S}_{mm}(q_{max})`$ would decrease. An explanation for this invariance with $`f`$ seen by the theory is due to screening: as $`f`$ increases, the repulsion between polymer chains increases; on the other hand, charge neutrality forces more and more counterions into the solution which increases the screening between polymers. At small $`f`$ polymer and counterion correlations are weak, but at large $`f`$ they increase such that one reaches a balance between polymer-polymer repulsion and counterion screening.

In conclusion, we have presented a new optimization scheme that appears to improve substantially the predictive power of the RPA theory for the structure of polyelectrolyte solutions. We showed that this range-optimized RPA in the minimal model yields all trends exhibited in the experimental data of ELW at large $`f`$ and is in moderate agreement with that data for the peak position and height. The crudeness of our chain model, however, precludes us from drawing any firm conclusions about the origins of the effects seen in ELW. On the other hand, since the minimal model $`\mathrm{𝑑𝑜𝑒𝑠}`$ predict an invariance and the comparison is reasonable enough, it is apparent that a number of questions need to be answered before the cause(s) of this phenomenon is determined. But it is interesting that a model of screened Coulomb chains can produce effects that would normally be associated with more complex mechanisms such as counterion condensation.

###### Acknowledgements.

We thank John G. Curro, Andrea J. Liu, John D. McCoy, Monica Olvera de la Cruz and Craig E. Pryor for helpful conversations and correspondence.


format: text

text:
# Stress-Energy Tensor for the Massless Spin 1/2 Field in Static Black Hole Spacetimes

## Abstract

The stress-energy tensor for the massless spin 1/2 field is numerically computed outside and on the event horizons of both charged and uncharged static non-rotating black holes, corresponding to the Schwarzschild, Reissner-Nordström, and extreme Reissner-Nordström solutions of Einstein’s equations. The field is assumed to be in a thermal state at the black hole temperature. Comparison is made between the numerical results and previous analytic approximations for the stress-energy tensor in these spacetimes. For the Schwarzschild (charge zero) solution, it is shown that the stress-energy differs even in sign from the analytic approximation. For the Reissner-Nordström and extreme Reissner-Nordström solutions, divergences predicted by the analytic approximations are shown not to exist.

Classically, the most general static black hole solution to Einstein’s equations is spherically symmetric and is described by just two parameters, the mass ($`M`$) and charge ($`Q`$). Quantum effects may be taken into account by solving the semi-classical backreaction equations ,

$$G_{\mu \nu }=8\pi T_{\mu \nu },$$
(1)
where $`T_{\mu \nu }`$ is the expectation value of the stress-energy tensor operator, which we shall hereafter refer to as the stress-energy tensor for the quantized fields. It includes the effects both of physical particles and vacuum polarization. A non-zero stress-energy tensor will alter the black hole’s geometry. Even if Eq. (1) is not explicitly solved, finding the stress-energy tensor in a given background spacetime allows one to estimate the size of quantum effects in that spacetime. For example, if the stress-energy diverges at the event horizon, quantum effects may prevent the existence of the black hole. Such a divergence is predicted in one component of the stress-energy tensor by analytic approximations for massless spin 0, 1/2, and 1 fields in Reissner-Nordström ($`0<|Q|<M`$) and extreme Reissner-Nordström ($`|Q|=M`$) spacetimes  and, in one case, by a two dimensional calculation of the full renormalized stress-energy tensor for a massless spin 0 field .

Numerical computations of the stress-energy tensor for quantized fields in black hole spacetimes provide important tests of both analytic approximations and two dimensional calculations. Previous numerical computations of the stress-energy of quantized fields in black hole spacetimes have primarily been for scalar fields in Schwarzschild ($`Q=0`$), Reissner-Nordström, and extreme Reissner-Nordström spacetimes . While computations relating to scalar fields are important and provide some insight into the effects of quantized fields, it is important to also make computations with more realistic fields of nonzero spin. Significant differences can exist between fields of zero versus nonzero spin . Because of difficulties in working with fields of higher spin in curved space, the only previous numerical calculation of the full renormalized stress-energy tensor in a black hole spacetime for a field of nonzero spin was for the massless spin $`1`$ field in Schwarzschild spacetime .

In this paper we present the results of numerical computations of the stress-energy tensor for the massless spin 1/2 field in Schwarzschild, Reissner-Nordström, and extreme Reissner-Nordström spacetimes in the Hartle-Hawking-Israel state , a thermal state at the black hole temperature. Our results address the validity of the analytic approximations that have been derived for these spacetimes  as well as the question of whether the stress-energy is finite on the horizon. We find that in Schwarzschild spacetime the analytic approximation for the field is very poor, not even giving the correct sign for the energy density on the horizon. This is in contrast to the cases of massless scalar fields  and the spin $`1`$ field in Schwarzschild spacetime  where it was found that analytic approximations were good both on and everywhere outside of the event horizon. We also find that for Reissner-Nordström and extreme Reissner-Nordström spacetimes all components of the stress-energy tensor are finite on the event horizon. This is in contrast to the analytic approximation of Frolov and Zelnikov  which predicts a divergence in one component on the horizon. It is also in contrast to the result found by Trivedi  that one component of the stress-energy tensor for the conformally invariant field in two dimensions is divergent on the horizon of an extreme Reissner-Nordström black hole. However, it is in agreement with the numerical results found previously for massless scalar fields with arbitrary curvature couplings .

In Ref.  a method was given to compute the stress-energy tensor for the massless spin 1/2 field in a general static spherically symmetric spacetime. The method involves writing the stress-energy tensor in terms of the Euclidean Green’s function for the field and renormalizing using the method of point splitting . An explicit expression for the Euclidean Green’s function in a static spherically symmetric spacetime was derived in terms of mode solutions to the Euclidean Dirac equation. The time dependent part of the modes and the angular dependent part were both given in terms of known functions. The radial part was given in terms of solutions to a coupled set of ordinary differential equations. For most spacetimes, including black hole spacetimes, these equations must be solved numerically. The final expression for the stress-energy tensor consists of two parts both of which are finite and separately conserved. For one part the mode sums have been computed analytically for a general static spherically symmetric spacetime. It has a trace equal to the trace anomaly and thus can serve as an analytic approximation for the entire stress-energy tensor. It is equivalent to previous analytic approximations for the massless spin 1/2 field . The other part must usually be computed numerically.

To simplify the numerical computation of the stress-energy tensor it is helpful to use the WKB approximation for the high frequency and/or large angular momentum modes . The higher the order of the WKB approximation used the fewer modes that must be numerically computed to give an answer to a given accuracy. However the number of terms in the WKB approximation grows rapidly with order so it is not practical to use too high an order.

The metric for a Reissner-Nordström spacetime can be written

$$ds^2=f(r)dt^2+f(r)^1dr^2+r^2d\mathrm{\Omega }^2,$$
(2)
where

$$f(r)=1\frac{2M}{r}+\frac{Q^2}{r^2}.$$
(3)
The temperature of the black hole is

$$2\pi T=\frac{\sqrt{M^2Q^2}}{2M^2Q^2+2M\sqrt{M^2Q^2}}.$$
(4)
Recall that in the Hartle-Hawking-Israel state, the field is assumed to be in a thermal state at this temperature. Note that the extreme Reissner-Nordström black hole is at zero temperature.

In Figs. 1 and 2 our results for Schwarzschild spacetime are displayed. As in the case of conformally invariant spin $`0`$  and $`1`$  fields, the weak energy condition is violated near the horizon***With our conventions the energy density is $`\rho =T_{t}^{}{}_{}{}^{t}`$.. The analytic approximation is very poor near the horizon (although it is accurate far from it), and even gets the sign of the energy density wrong. This is in contrast to the spin $`0`$ and $`1`$ cases where the analytic approximation is good near the horizon .

As a check on this surprising result, we also computed the energy density at the horizon using radial point splitting with one point fixed on the horizon. This method was developed by Candelas  who used it to compute the stress-energy tensor at the horizon for the conformally invariant scalar field in Schwarzschild spacetime. It has the advantage that modes of only the two smallest allowed frequencies contribute to the mode sum. The method in  involves splitting the points in the time but not the space direction. The resulting expression for the stress-energy tensor contains a sum over modes with arbitrarily large frequencies. Thus the two methods, while not completely independent, are significantly different. The renormalization counterterms are also different for these two methods. The energy density of the spin 1/2 field on the horizon in the radial point splitting calculation was found to agree to better than four digits with that obtained in the time point splitting calculation.

Some of our results for Reissner-Nordström spacetimes are displayed in Fig. 3. For non-zero charge $`Q`$, there is a well-known ambiguity in the analytic portion of the stress-energy tensor caused by renormalization . This ambiguity does not arise for the Schwarzschild metric, nor does it exist on the horizon for $`T_{}^{t}{}_{t}{}^{}`$, $`T_{}^{r}{}_{r}{}^{}`$, or $`T_{}^{\theta }{}_{\theta }{}^{}`$ in Reissner-Nordström or extreme Reissner-Nordström spacetimes. For definiteness, we have chosen the renormalization parameter $`\mu `$ so that $`\mu M=1`$. From Fig. 3, it is clear that the energy density on the horizon is negative and its magnitude increases with increasing $`|Q|/M`$. Thus the weak energy condition is violated in these cases as well.

One of the predictions of the analytic approximation for the massless spin 1/2 field is that one component of the stress-energy diverges on the event horizon for all Reissner-Nordström black holes. The component is related to the energy density seen by a freely falling observer passing through the event horizon and is explicitly

$$\frac{T_{}^{r}{}_{r}{}^{}T_{}^{t}{}_{t}{}^{}}{f(r)}.$$
(5)
For $`|Q|<M`$ the divergence in this component of the analytic approximation is logarithmic and of the form

$$\frac{Q^2}{20\pi ^2r^6}\mathrm{log}\left[\mu ^2M^2f(r)\right],$$
(6)
where $`\mu `$ is the renormalization parameter discussed above. This divergence also occurs in the case $`|Q|=M`$. However in that case there is a term in the analytic approximation for the component (5) with a stronger divergence. It is

$$\frac{M^2(32M23r)}{180\pi ^2r^6(rM)}.$$
(7)
This is also the type of divergence found by Trivedi  for the same component of the full renormalized stress-energy tensor of a conformally invariant scalar field in a two dimensional extreme Reissner-Nordström spacetime.

However, when we numerically computed the component (5) for the massless spin 1/2 field in Reissner-Nordström and extreme Reissner-Nordström spacetimes, we found no evidence that this component of the stress-energy tensor diverges on the event horizon. This is in agreement with results found for massless scalar fields with arbitrary curvature couplings in these spacetimes . Some of our results are shown in Fig. 4.

In summary the full renormalized stress-energy tensor for the massless spin 1/2 field has been numerically computed in Schwarzschild, Reissner-Nordström, and extreme Reissner-Nordström spacetimes in the case that the field is in the Hartle-Hawking-Israel state. The analytic approximation for the stress-energy tensor in Schwarzschild spacetime is found to be very poor near the event horizon. Divergences in the stress-energy, predicted by the analytic approximations for Reissner-Nordström and extreme Reissner-Nordström, are found not to occur. Thus it is clear that the analytic approximations for the massless spin 1/2 field in static black hole spacetimes are of extremely limited use and must usually be supplemented by numerical calculations of the full renormalized stress-energy tensor.

###### Acknowledgements.

This work has been supported in part by grant number PHY-0070981 from the National Science Foundation.


text:
# References

THE UNITED THEORY OF THE TWO FIELDS OF THE ELECTRIC AND MAGNETIC NATURE

Rasulkhozha S. Sharafiddinov

Institute of Nuclear Physics, Uzbekistan Academy of Sciences, Tashkent, 702132 Ulugbek, Uzbekistan

Each of electrically charged particles testifies in favor of the existence of a kind of the magnetically charged monoparticle. As a consequence, only the corresponding mononeutrinos answer for quantization of the electric charges of all the neutrinos. Therefore, to understand the nature of matter at the fundamental level one must use the electromagnetic field as the field of the unified system of the photon and monophoton where the electric and magnetic forces of the nature are united. Some logical and laboratory confirmations of the availability of compound structure of gauge invariance have been listed which say also about the connection between the states of elementary particles and monoparticles.

A study of the behavior of electrons and their neutrinos in a nucleus Coulomb field shows clearly that between the mass of a particle and its charge there exists a sharp dependence . For the light Dirac neutrino $`(\nu =\nu _e)`$ it has the form

$$e_\nu ^E=\frac{3eG_F\left(m_\nu ^E\right)^2}{4\pi ^2\sqrt{2}},e=\left|e\right|.$$
(1)
Here $`m_\nu ^E`$ and $`e_\nu ^E`$ are the neutrino electric rest mass and charge.

Such an intimate connection saying about the gravitational structure of a Coulomb interaction and reflects the fact that the mass and charge of a particle correspond to the most diverse form of the same regularity of the nature of this field .

On the other hand, it is known that in the framework of the hypothesis of field mass based on the classical theory of an extensive electron , a particle all the mass is purely electric. However, according to such a point of view, the charge distribution of the neutrino will not be steady.

Our conclusion is that each of existing types of charges come forwards in the system as the source of a kind of the inertial mass . Herewith all the mass of the neutrino coincides with its united rest mass $`m_\nu ^U`$ which includes in self the electric $`m_\nu ^E,`$ weak $`m_\nu ^W,`$ strong $`m_\nu ^S`$ and some other components:

$$m_\nu ^U=m_\nu ^E+m_\nu ^W+m_\nu ^S+\mathrm{}.$$
(2)
Such a compound structure giving the possibility to directly look at the nature of matter is observed as a consequence of the availability of the united charge $`e_\nu ^U`$ in neutrino equal to its all the charge :

$$e_\nu ^U=e_\nu ^E+e_\nu ^W+e_\nu ^S+\mathrm{},$$
(3)
where the indices $`E,`$ $`W`$ and $`S`$ correspond to the electric, weak and strong parts of charge.

Thus, it follows that if the neutrino possesses a Dirac mass then all possible components of this mass establish the intraneutrino interratio between the forces of the electric and unelectric nature. Therefore, the charge distribution of the neutrino must be steady .

The purpose of the present article is to elucidate what the above - mentioned duality of matter say about structures of electromagnetic field, gauge invariance and charge quantization law.

In the framework of the loop phenomena, the neutrino must be electrically neutral at the conservation of gauge invariance . It appears that here on the basis of (1) one can will decide a question about the equality of the neutrino electric mass to zero.

But if $`m_\nu ^E=0,`$ the neutrino is strictly Majorana particle . At the same time the existence of massive Dirac neutrinos is by no means excluded experimentally. Insofar as their neutrality in the loop approximation is concerned, it reflects just a latent regularity of general picture of gauge invariance which requires the study of the structure of electromagnetic field.

At first sight it appears that a connection between $`e_\nu ^E`$ and $`m_\nu ^E`$ is incompatible with charge quantization. Such a regularity, however, takes place owing to the magnetic monopoles. It is also known that the standard $`SU(2)_L`$$``$$`U(1)`$ theory, by itself, does not require  that the magnetic monopoles be exist in the form as was suggested by Dirac . Therefore, from point of view of a massive four - component neutrino, it should be expected that each of electrically charged particles testifies in favor of the existence of a kind of a ¡monoparticle¿ namely a particle having the magnetic mass and charge. Then it is possible, for example, only the corresponding mononeutrinos are responsible for quantization of the electric charges of all the neutrinos.

There exists a range of other phenomena, in which a fundamental part is said of monoparticles. One of such systems may serve as the source of electromagnetic field. But unlike earlier known, it must not be usual gauge boson. The point is that the same photon does not lead to the appearance of electromagnetic field. It possesses the electric mass  and charge .

From point of view of a Dirac neutrino, this gauge state will indicate to the existence of the monophoton with the magnetic mass and charge. In other words, the photon $`\gamma _E`$ and monophoton $`\gamma _H`$ must be source of the electric $`\stackrel{}{E}`$ and magnetic $`\stackrel{}{H}`$ fields respectively. Under shuch circumstances the electromagnetic field $`(\stackrel{}{E},\stackrel{}{H})`$ appears as the field of the unified system of the photon and monophoton $`(\gamma _E,\gamma _H)`$ where the electric $`\stackrel{}{F}_E`$ and magnetic $`\stackrel{}{F}_H`$ forces of the nature are united.

To elucidate the compound structure of its gauge invariance one must apply once more to the question whether the same particle has simultaneously both electric and magnetic charges. There exist, however, many uncertainties both in the nature and in the behavior of these types of charges. Nevertheless, if we assume  that the possibility of the existence of such Dirac fermions is not excluded, the neutrino interaction with an electromagnetic field of emission arises at the expense of exchange simultaneously both by the photon and by the monophoton. The latter would lead us to the implication that the invariance of these types of gauge bosons concerning C, P and T, and also their combinations CP and CPT are not different.

Finally, insofar as the appearance of magnetic field in usual laboratory conditions is concerned, to this one must apply as to one of the available experimental data confirming the existence of monoparticles and saying that any electrically charged particle can be converted into the corresponding monoparticle and vice versa. Of course, these transitions and many other aspects of compound structure of each of existing types of charges and masses open of principle possibilities for creation at the fundamental level of truly physical theories of naturally united gauge fields.


text:
# Power Law Distributions of Offspring and Generation Numbers in Branching Models of Earthquake Triggering

## I Introduction

All large earthquakes are followed by an increase of seismic activity known as “aftershocks.” Aftershock sequences of small earthquakes are less obvious because the aftershock productivity is weaker, but can be clearly observed when stacking many sequences alpha . It is thus natural to assume that each earthquake can trigger its own aftershock sequence, and that observed aftershock sequences result from the cascade of direct aftershocks (triggered directly by the mainshock) and indirect aftershocks (triggered by a previous aftershock of the mainshock). This assumption is the basis of the Epidemic Type Aftershock Sequence model (ETAS) of seismicity Ogata88 ; Ogata99 ; KK81 ; KK87 ; helmsor1 , which describes earthquake triggering as a branching process. In addition, the ETAS model includes the Omori law decay $`1/(t+c)^p`$ of the number of direct aftershocks as a function of the time $`t`$ since the mainshock, where $`c`$ is a small constant and $`p`$ is an exponent close to but larger than 1. Previous works on this model have shown that the ETAS model provides a better fit to aftershock sequences than a single Omori law (no secondary aftershocks) GuoO97  and that a significant fraction of aftershocks, of the order of 80%, are secondary aftershocks Felzer1 ; Felzer2 . The ETAS model has been used in many studies KK81 ; Felzer1  to describe or predict the spatio-temporal distribution of seismicity and reproduces many properties of real seismicity, including a renormalization of the Omori exponent from the local Omori law (direct aftershocks) to the global Omori law (observed rate of aftershocks including secondary aftershocks) Ogata99 ; SorSor ; helmsor1 , Båth’s law Bathlawpap , a diffusion of aftershocks and realistic foreshock properties foreshockp ; Forexp . In this work, we present an analytical derivation of the distribution of the total number of aftershocks, summed over all generations of the cascade of aftershock triggering, and of the distribution of the total number of generations of aftershocks before extinction.

There are two well known statistical laws that describe the scale-invariance of earthquake physics with respect to magnitudes. First, the (complementary cumulative) Gutenberg-Richter (GR) distribution of earthquake magnitudes gives the probability

$$p_m(m)10^{bm}$$
(1)
that an earthquake has magnitude equal to or larger than $`m`$. This magnitude distribution $`p_m(m)`$ is not dependent on the magnitude of the triggering earthquake, i.e., a large earthquake can be triggered by a smaller one alpha ; Forexp .

Second, the average number of aftershocks triggered directly by an earthquake of magnitude $`m`$ is found to increase with $`m`$ as

$$N_m=K10^{\alpha m}$$
(2)
where $`K`$ is a numerical factor independent of the magnitude. The number of direct aftershocks cannot be measured, because what is observed is the total number of direct and secondary aftershocks. If earthquake triggering can be described by a branching process such as the ETAS model, then it can be shown that the scaling of the total number of aftershocks with the mainshock magnitude also obeys the law (2), but with a larger factor $`K^{}`$ which accounts for the cascades of secondary aftershocks helmsor1 . The exponent $`\alpha `$ in (2) can thus be measured from a fit of the total number of aftershocks with the mainshock magnitude. Fits of the total number of aftershocks as a function of the mainshock magnitude in individual sequences support (2) with an exponent $`\alpha `$ in the range 0.75-1 6 ; 9 ; 11 ; 13 . However, the precision of these studies is limited by the narrow range of mainshock magnitudes considered and the large scatter of the number of aftershocks per mainshock. The value of $`\alpha `$ estimated in these studies may also be biased by the arbitrary constraint that aftershocks must be smaller than the mainshock, by the incompleteness of the catalog just after the mainshock, and by the background seismicity at large times after the mainshock. Other studies have measured the exponent $`\alpha `$ in (2) by calibrating the ETAS model to real data (individual aftershock sequences or complete catalogs) using maximum likelihood methods Ogata88 ; K91 ; GuoO97 ; Zhuangetal . These studies found a large scatter of the $`\alpha `$ exponent in the range $`0.22`$. It is not yet clear if this range of values reflects a real variability of $`\alpha `$ or the inaccuracy of the estimation of $`\alpha `$. One of us used a stacking method to estimate the average rate of earthquakes triggered (directly or indirectly) by a previous earthquake as a function of the magnitude of the triggering earthquake alpha , without constrain on the aftershock magnitude. For the catalog of the Southern California Data Center for Southern California, using the time period 1975-2003 and $`m3`$ earthquakes, $`\alpha `$ is found equal to $`0.8\pm 0.1`$, smaller than $`b=1\pm 0.1`$. Small earthquakes are thus collectively more important than larger earthquakes for earthquake triggering if $`\alpha b`$, because they are much more numerous than larger ones.

Let us combine the two laws (1) and (2) to get the unconditional probability density for the number $`N_m`$ of events triggered directly by any event whose magnitude $`m`$ is drawn at random from the GR law. For this, we note that

$$\mathrm{Pr}(\mathrm{fertility}N_m)=\mathrm{Pr}(K10^{\alpha m}N_m)=\mathrm{Pr}(m\frac{1}{\alpha }\mathrm{log}_{10}\frac{N_m}{K})10^{\frac{b}{\alpha }\mathrm{log}_{10}\frac{N_m}{K}}\frac{1}{N_m^{b/\alpha }},$$
(3)
where $`Pr`$ means “Probability”. The first equality makes use of (2) and the third equivalence makes use of (1). Hence, the probability density of the fertilities $`N_m`$, for a magnitude $`m`$ drawn at random in the GR law, is

$$p_{N_m}(N_m)\frac{1}{N_{m}^{}{}_{}{}^{1+\gamma }},$$
(4)
with $`\gamma =b/\alpha `$. Typically, $`b=1`$ and $`\alpha 0.8`$ alpha  leading to $`\gamma 1.25`$. Because $`\gamma <2`$, the variance of the fertility, for an earthquake of arbitrary magnitude drawn from the GR law, is mathematically infinite (if we assume the absence of a roll-off in the GR distribution, see below).

Beyond earthquakes, a multitude of phenomena can also be described by branching processes with power-law distributions of fertility. Stochastic branching processes indeed describe well a multitude of phenomena Athreya ; Sankaranarayanan  from chain reactions in nuclear and particle physics, material rupture, fragmentation and earthquake processes, critical percolation cluster sizes and population growth models, to population and biological dynamics, epidemics, economic and social cascades and so on. Branching processes are also of particular interest because deep connections have been established with critical phenomena Vere76 ; Vere77 . Epidemic transmission of diseases, and more generally transmission processes involving avalanches spreading on networks such as the World Wide Web, cellular metabolic network, ecological food webs, social networks, and so on exhibit such heavy-tail probability density functions (PDF) given by (5) below, as a consequence of the well-documented power law distribution of connectivities among nodes. Our results are thus relevant to systems in which the number of offsprings may be large due to long-range interactions, long-memory effects or large deviation processes. Goh et al. Goh  actually derive results that overlap with ours in the context of avalanches in social networks.

In branching processes with a finite variance of the number of daughters per mothers, various quantities exhibit power law distributions with universal exponents at criticality (in statistical physics, the term “universal” refers to the independence of the critical exponents on the microscopic details of the physics). This includes the distributions of cluster sizes, of the number of generations before extinction and of durations which are mean field Athreya ; Sankaranarayanan  (“mean field” refers to the branching approximation which leads to a lack of dependence on the space dimension). In the case of earthquakes and for other systems mentioned above, the distribution of fertilities has an infinite variance, leading to an anomalous scaling of offsprings and generation numbers. While the number of direct aftershocks per mainshock for a fixed mainshock magnitude has a finite variance, usually modeled by a Poisson distribution, the effect of multiple cascades of triggering and the variability of the fertility of each earthquake lead to a much larger power-law distribution for the total number of aftershocks summed over all generations. As a consequence, there are huge fluctuations of the total number of aftershocks from one sequence to another one, for the same mainshock magnitude.

The goal of this paper is to provide a general exact derivation of the distribution of the total aftershock productivity of any earthquake, summed over all generations of the cascade of aftershock triggering. We also derive the exact distribution of the total number of generations of aftershocks before extinction in the case of a power-law distribution of fertility relevant for earthquakes and for many other systems. Beyond earthquakes our results apply to any branching process with a power-law distribution of fertilities.

## II Model and main results

We consider a general branching process in which each progenitor or mother is characterized by its average number $`N_m\kappa \mu (m)`$ of children (first generation offsprings), where $`\mu (m)=10^{\alpha (mm_0)}`$ is a mark associated with an earthquake of magnitude $`mm_0`$, $`\kappa `$ is a constant factor and $`m_0`$ is the minimum magnitude of earthquakes capable of triggering other earthquakes. We note $`\mu `$ the mark of an earthquake which has an arbitrary magnitude $`m`$ drawn according to the GR law. According to (3), the mark $`\mu `$ is distributed according to

$$p_\mu (\mu )=\frac{\gamma }{\mu ^{1+\gamma }},1\mu <+\mathrm{},\gamma =b/\alpha .$$
(5)
Note that $`p_\mu (\mu )`$ is normalized: $`_1^+\mathrm{}𝑑\mu p_\mu (\mu )=1`$. The relation $`N_m\kappa \mu (m)`$ together with (5) thus ensures that the fertility obeys the law (4). For a fixed $`\gamma `$, the coefficient $`\kappa `$ then controls the value of the average number $`n`$ of children of first generation per mother:

$$n=N_m=\kappa \mu =\kappa \frac{\gamma }{\gamma 1},$$
(6)
where the average $`N_m`$ is taken over all mothers’ magnitudes drawn in the GR law. In the terminology of branching processes, $`n`$ is called the branching ratio. For $`n<1`$, there are less than one child per mother: this corresponds to transient (sub-critical) branching processes with finite lifetimes with probability one. For $`n>1`$, there are more than one child per mother: this corresponds to explosive (super-critical) branching processes with a number of events growing exponentially with time. The value $`n=1`$ of exactly one child per mother on average is the critical point separating the two regimes.

The realized number of children of an earthquake of fixed magnitude $`m`$ can be deterministic or may result from a Poisson or other more general distribution with mean $`\kappa \mu (m)`$. Any given child $`i`$ may then generate an average number $`\kappa \mu _i`$ of children, where the mark $`\mu _i`$ is specific to the child $`i`$ and is drawn from the PDF (5). These grand-children of the initial progenitor in turn generate new children with marks drawn from (5), and so on. The process cascades down along generations.

Here, we focus on global quantities $`p_r(r)`$ (total number $`r`$ of offsprings until extinction following a mainshock) and $`p_g(g)`$ (total number $`g`$ of generations until extinction following a mainshock). Therefore, arbitrary time-dependent branching processes can be considered and our results thus apply to arbitrary stochastic marked point-processes in discrete or continuous time. In particular, our results apply to the ETAS model. We have previously shown that, for $`n<1`$, the infinite variance of the number of first-generation events leads to anomalous global direct and inverse Omori law with apparent exponents varying continuously with $`\gamma `$ foreshockp . Our results also apply to variations of the ETAS models with different arbitrary time dependences.

In this work, we present an analytical derivation of the distribution of the total number of aftershocks, summed over all generations of the cascade of aftershock triggering, and of the distribution of the total number of generations of aftershocks before extinction. Specifically, we uncover a novel regime with continuously varying exponents for the probability density function (PDF) $`p_r(r)`$ of the total number $`r`$ of progenies and $`p_g(g)`$ of their total number $`g`$ of generations of aftershocks of a mainshock before extinction. This regime appears in the regime of large deviations $`\gamma <2`$ relevant for earthquakes, when the distribution of the number of first-generation offsprings from any mother has a power law tail with infinite variance.

We study the sub-critical and critical regimes for which the number $`n`$ of children per mother, averaged over all possible numbers of children per mother, is less than or equal to $`1`$. This condition $`n1`$ ensures that, with probability $`1`$, the cascade ends after a finite number of generations with a finite total number of offsprings Athreya ; Sankaranarayanan . Figure 1 presents comparisons between numerical simulations of branching processes for different values of $`\gamma `$ and our main result (at criticality $`n=1`$) derived below:

$$p_r(r)C_r/r^{1+\frac{1}{\gamma }};p_g(g)C_g/g^{1+\frac{1}{\gamma 1}},$$
(7)
for $`1<\gamma 2`$. $`C_r`$ and $`C_g`$ are two constants independent of $`r`$ and $`g`$, respectively, such that $`p_r(r)`$ and $`p_g(g)`$ are normalized to 1. For $`n<1`$, the intermediate asymptotics (7) holds for $`r<R_21/(1n)^{\gamma /(\gamma 1)}`$ beyond which there is another power-law asymptotic

$$p_r(r)1/r^{1+\gamma }.$$
(8)
For $`\gamma 2`$ and $`n=1`$, the standard mean field scaling

$$p_r(r)1/r^{3/2}\mathrm{and}p_g(g)1/g^2$$
(9)
are retrieved. The regime $`0<\gamma 1`$ gives rise to explosions and, in continuous time, to stochastic finite-time singularities sorhelmsing . If the GR law is truncated or exhibits a deviation from its standard form in the large magnitude range, our results will hold for intermediate values of $`r`$ and $`g`$ and will be truncated at large $`r`$’s and large $`g`$’s.

The following section gives the technical derivation of the results (7,8,9).

## III General formulation

### III.1 Formal solution in terms of Probability Generating Functions (PGF)

The most general relations of branching theory are expressed via probability generating functions (PGF) defined by

$$z^R=\underset{r=0}{\overset{\mathrm{}}{}}p_r(r)z^r$$
(10)
where $`\mathrm{}`$ means statistical averaging. By definition, $`p_r(r)=\text{Pr}(R=r)`$ is the probability that the random variable $`R`$ takes the value $`r`$. It can be obtained from its PGF by the relation

$$p_r(r)=\text{Pr}(R=r)=\frac{1}{r!}\frac{d^rz^R}{dz^r}|_{z=0}.$$
(11)
The key property of PGFs is that the PGF of the sum of statistically independent summands is equal to the product of the summands PGF’s. This is useful for branching processes for which the number of daughters triggered by different mothers are statistically independent. We introduce four PGF’s: (i) $`G_{\mu _0}(z)`$ is the PGF of the number $`R_1`$ of daughters generated from a given mother with fertility $`\kappa \mu `$ at the first generation; (ii) $`G(z)`$ is the same as $`G_\mu (z)`$ but averaged over all mothers’ fertility $`\kappa \mu `$; (iii) $`\mathrm{\Theta }_\mu (z)`$ is the PGF of the total number of daughters generated from a given mother with fertility $`\kappa \mu `$ summed over all generations; (iv) $`\mathrm{\Theta }(z)`$ is the same as $`\mathrm{\Theta }_\mu (z)`$ but for a mother of arbitrary fertility. Then, general branching theory Athreya ; Sankaranarayanan  gives the functional equations

$$\mathrm{\Theta }(z)=G(z\mathrm{\Theta }(z));\mathrm{\Theta }_\mu (z)=G_\mu (z\mathrm{\Theta }(z)).$$
(12)
Using $`R=d\mathrm{\Theta }(z)/dz|_{z=1}`$ and similarly for the other PGF, the formulas (12) lead to $`R=n/(1n)`$ for the number of daughters summed over all generations and averaged over all possible mother’s fertilities, where $`n`$ defined by (6) is the average number of first-generation daughters from mothers of arbitrary fertility.

Using Lagrange series, we transform the implicit equation (12) in $`\mathrm{\Theta }_\mu `$ into an explicit equation. Recall that a Lagrange series is the Taylor expansion of the function $`G_\mu (y(z,a))`$ with respect to $`z`$, where $`y(z,a)`$ is solution of the implicit equation $`y=a+zG(y)`$. For infinitely differentiable functions $`G_\mu (y)`$ and $`G(y)`$ in the neighborhood of $`z=a`$, we have the following Taylor series

$$G_\mu (y)=G_\mu (a)+\underset{k=1}{\overset{\mathrm{}}{}}\frac{z^k}{k!}\frac{d^{k1}}{da^{k1}}\left[G^k(a)\frac{dG_\mu (a)}{da}\right].$$
(13)
In (12), we have $`y(z,a=0)=z\mathrm{\Theta }(z)`$, and using the identity $`p_{r,\mu }(r)=(1/r!)d^r\mathrm{\Theta }_\mu (z)/dz^r|_{z=0}`$ allowing to recover the PDF from its generating function, we obtain the explicit formula for the PDF $`p_{r,\mu }(r)`$ of the total number of daughters born from a mother with fertility $`\mu `$:

$`p_{r,\mu }(0)`$ $`=`$ $`G_\mu (0)`$ (14)
$`p_{r,\mu }(r)`$ $`=`$ $`{\displaystyle \frac{1}{r!}}{\displaystyle \frac{d^{r1}}{dz^{r1}}}\left[G^r(z){\displaystyle \frac{dG_\mu (z)}{dz}}\right]|_{z=0},r>0.`$ (15)
Let us now specialize to the Poisson distribution

$$p_{r_1,\mu }(r_1)=\frac{(\mu \kappa )^{r_1}}{r_1!}e^{\mu \kappa },$$
(16)
giving the PDF of the number $`R_1`$ of daughters of first generation born from a mother with fertility $`\mu `$. By construction, the average of $`R_1`$ at fixed mother fertility $`\mu `$ over an ensemble of Poisson realizations is $`R_1_\mu =\mu \kappa `$. The PGF associated with (16) is

$$G_\mu (z)=e^{\mu \kappa (z1)}.$$
(17)
From (17), we obtain the expression of $`G(z)`$ by averaging over all possible fertilities $`\mu `$ according to the PDF (5). Using explicitly the normalized PDF $`p_\mu (\mu )=\gamma /\mu ^{1+\gamma }`$ for $`\mu 1`$ and $`0`$ otherwise, we obtain

$$G(z)=\gamma \kappa ^\gamma (1z)^\gamma \mathrm{\Gamma }(\gamma ,\kappa (1z)).$$
(18)
Expression (18) can be expanded as

$$G(z)1a_1ya_\gamma y^\gamma +a_2y^2+\mathrm{}$$
(19)
where $`y=\kappa (1z)`$, $`a_1=\frac{\gamma }{\gamma 1}`$, $`a_\gamma =\mathrm{\Gamma }(1\gamma )`$ and $`a_2=\frac{\gamma }{2(\gamma 2)}`$. This expansion (19), valid for any $`\gamma `$, applies beyond the specific Poisson process (16). It is solely based on the asymptotic power law PDF (5) of the average number of daughters at the first generation. Thus, our asymptotic results (7) derived below hold under very general conditions.

### III.2 Probabilistic interpretation

Let us provide a probabilistic interpretation of formula (15). In addition to be of intuitive appeal, it will be very useful for the derivation of our results.

Let us consider the random integer $`R`$ with PDF $`P(k)=\text{Pr}(R=k)`$ and its PGF as

$$S(z)=\underset{k=0}{\overset{\mathrm{}}{}}P(k)z^k,$$
(20)
and reciprocally

$$P(k)=\text{Pr}(R=k)=\frac{1}{k!}\frac{d^kS(z)}{dz^k}|_{z=0}.$$
(21)
Let us decompose $`R=X+Y`$, where $`X`$ and $`Y`$ are statistically independent random integers with PGFs $`Q(z)`$ and $`G(z)`$ respectively. Due to the statistical independence of $`X`$ and $`Y`$, we have $`S(z)=Q(z)G(z)`$ and relation (21) takes the form

$$P(k)=\text{Pr}(X+Y=k)=\frac{1}{k!}\frac{d^k}{dz^k}\left[Q(z)G(z)\right]|_{z=0}.$$
(22)
Let in turn

$$Y=Y(r)=Y_1+Y_2+\mathrm{}+Y_r$$
(23)
where $`\{Y_1,\mathrm{},Y_r\}`$ are $`r`$ statistically independent random integers with the same PGF $`G(z)`$. Then, formula (22) takes the form

$$P(k)=\text{Pr}(X+Y(r)=k)=\frac{1}{k!}\frac{d^k}{dz^k}\left[Q(z)G^r(z)\right]|_{z=0}.$$
(24)
Let us now rewrite relation (15) in a form similar to (24). For this, let us introduce the auxiliary PGF

$$Q_\mu (z)=\frac{1}{r_1(\mu )}\frac{dG_\mu (z)}{dz}.$$
(25)
It is easy to prove rigorously that, for an arbitrary PGF $`G_\mu (z)`$ possessing a finite mean value

$$r_1(\mu )=\frac{dG_\mu (z)}{dz}|_{z=1}<\mathrm{},$$
(26)
then the auxiliary function $`Q_\mu (z)`$ is indeed the PGF of some random variable. Moreover, in the case under consideration of a Poissonian PGF, we have

$$G_\mu (z)=e^{\mu \kappa (z1)}Q_\mu (z)G_\mu (z)=e^{\mu \kappa (z1)}$$
(27)
and

$$r_1(\mu )=\mu \kappa .$$
(28)
One can thus rewrite relation (15) in the form

$$p_{r,\mu }(r)=\frac{\mu \kappa }{r!}\frac{d^{r1}}{dz^{r1}}\left[G^r(z)G_\mu (z)\right]|_{z=0}.$$
(29)
Taking into account relation (24) and introducing the new random integer $`X(\kappa \mu )`$ possessing the Poissonian PGF (27), we can rewrite (29) in the form

$$p_{r,\mu }(r)=\frac{\mu \kappa }{r}\text{Pr}\left[X(\kappa \mu )+Y(r)=r1\right].$$
(30)
This is the key formula that we will use for deriving the asymptotic relations for the PDF $`p_{r,\mu }(r)`$ using limit theorems.

## IV Distribution of the total number of aftershocks

### IV.1 Case $`\gamma \mathrm{}`$

Let us first consider the limiting case $`\gamma \mathrm{}`$. The formal limit of a power law with an exponent going to infinity can be assimilated to an exponential in the following sense. Writing the tail of the complementary distribution as $`A^\gamma /(A+\mu )^\gamma `$, and assuming that $`A`$ also grows with $`\gamma `$ as $`A=\gamma /a`$, then $`A^\gamma /(A+\mu )^\gamma =1/(1+a\mu /\gamma )^\gamma _{\gamma +\mathrm{}}\mathrm{exp}[a\mu ]`$. This corresponds to an exponential tail for the distribution of the number of first generation daughters for arbitrary mother’s fertilities, as in (16) but with the same fertility for all mothers. In this case, $`n=\kappa `$ as can be seen from (6). This amounts to replacing expression (18) by $`G(z)=e^{\kappa (z1)}`$ and we also have $`G_\mu (z)=e^{\mu \kappa (z1)}`$. Thus,

$$G^r(z)\frac{dG_\mu (z)}{dz}=\mu \kappa e^{\kappa (z1)(r+\mu )}.$$

Correspondingly,

$$\frac{d^{r1}}{dz^{r1}}\left[G^r(z)\frac{dG_\mu (z)}{dz}\right]=\mu \kappa ^r(r+\mu )^{r1}e^{\kappa (z1)(r+\mu )}.$$
(31)
By substituting (31) into the right-hand-side of (15) gives exactly

$$p_\mu (r)=\mu \kappa ^re^{\kappa \mu }\frac{\left[(\mu +r)e^\kappa \right]^r}{r!(\mu +r)}.$$
(32)
There are two ways of deriving asymptotic formulas for $`p_{r,\mu }(r)`$ in the case $`\gamma =\mathrm{}`$. The first one relies on Stirling’s formula $`r!\sqrt{2\pi r}r^re^r`$ applied to expression (32). The second method relies on the probabilistic interpretation of the general formula (15) given in section III.2 together with the central limit theorem. The formulas obtained by these two methods differ in details, but are actually very close quantitatively over the most part of the tail, and are easily checked to converge to the same result asymptotically. The approach using the Stirling formula works only for $`\gamma =\mathrm{}`$ and does not allow to derive more general results for $`1<\gamma <\mathrm{}`$. The second approach is much more powerful and elegant. In particular, all asymptotics and crossovers in the case $`1<\gamma <\mathrm{}`$ are obtained by the second approach.

Let us first examine the first method. Using the above mentioned Stirling’s formula, one can rewrite (32) in the form

$$p_{r,\mu }(r)=\frac{1}{\mu }f(x)x=\frac{r}{\mu },$$
(33)
where

$$f(x)\frac{C}{\sqrt{2\pi x}(1+x)}\left(1+\frac{1}{x}\right)^{\mu x}e^{\mu \delta x}$$
(34)
and

$$C=\sqrt{\mu }e^{\kappa \mu },\delta =\kappa \mathrm{ln}\kappa 1.$$
(35)
Because Stirling’s formula is very precise even for $`r1`$ (for example $`1!=1`$ while Stirling’s formula gives $`1!0.922`$, $`2!=2`$ while Stirling’s gives $`2!1.919`$ and so on), formula (34) actually works well even for not too large $`\mu `$ and intermediate $`\kappa <1`$. The advantage of (34) is that it gives a clear understanding of the structure of the PDF $`f(x)`$. It consists of (i) a characteristic power tail

$$f(x)\frac{1}{\sqrt{2\pi x}(1+x)}x^{3/2},$$
(36)
(ii) an exponential decaying factor

$$f(x)e^{\mu \delta x}$$
(37)
which disappears when $`\kappa =1`$, i.e., $`\delta =0`$, and (iii) an algebraic factor

$$f(x)\left(1+\frac{1}{x}\right)^{\mu x}e^\mu \mathrm{exp}\left(\frac{\mu }{x}\right)$$
(38)
which possesses a lower cut-off at $`x\mu `$.

The shortcoming of (34) is that we can derive it only when an explicit expression such as (32) is available. We have not such luxury in the more general case $`1<\gamma <\mathrm{}`$ for which we need the second approach. Let us first quote the asymptotic formula (39) given below for $`\gamma +\mathrm{}`$ obtained from the application of the second probabilistic method. We shall show how to derive (39) as a special case of the next section for $`2<\gamma `$.

When the variance $`\sigma _1^2`$ of the number $`r_1`$ of daughters of the first generation from mothers of arbitrary fertilities is much greater than $`1\kappa `$, we find (see next section) that, due to central limit theorem, expression (32) reduces to

$$p_\mu (r)\frac{\mu \kappa }{r\sqrt{2\pi \kappa (r+\mu )}}\mathrm{exp}\left(\frac{[(1\kappa )r\mu \kappa 1]^2}{2\kappa (r+\mu )}\right).$$
(39)
At criticality, $`n=\kappa =1`$, this expression becomes

$$p_\mu (r)\frac{\mu }{r\sqrt{2\pi r}}\mathrm{exp}\left(\frac{\mu ^2}{2r}\right),$$
(40)
which retrieves the announced well-known mean field asymptotics $`p_\mu (r)1/r^{3/2}`$. The exponential term $`\mathrm{exp}(\mu ^2/(2r))`$ describes the roll-off of the number of aftershocks for small numbers close to the characteristic value $`r\mu ^2`$.

In the subcritical regime $`n<1`$, expression (32 ) gives an exponential decay of the number of aftershocks for large $`r`$.

These results can be checked explicitly as follows. The case of $`\gamma +\mathrm{}`$ corresponds either to $`b\mathrm{}`$ (all events have the same magnitude) or to $`\alpha =0`$ (all the events trigger the same expected number of offsprings independently of their magnitude). In both cases, $`\mu `$ tends to a constant independent of the mainshock magnitude. In order to check the above results, let us take this constant equal to $`1`$. This choice will modify the constants but not the functional dependences. Assuming $`\mu =1`$ transforms (32) into Taking $`\mu =1`$ transforms (32) into

$$p_\mu (r)=\kappa ^r\frac{(r+1)^r}{(r+1)!}e^{\kappa (r+1)}.$$
(41)
By using the Stirling formula, it has the approximation

$$p_\mu (r)\frac{e^{r\mathrm{ln}\kappa +(r+1)(1\kappa )}}{\sqrt{2\pi }(r+1)^{3/2}}.$$
(42)
When $`\kappa =n=1`$ (critical case), we retrieve $`p_\mu (r)r^{3/2}`$. When $`\kappa =n<1`$ is close to $`1`$, $`p_\mu (r)r^{3/2}e^{\delta r}`$ where $`\delta `$ is defined in (35).

### IV.2 Case of finite variance $`2<\gamma <\mathrm{}`$

In this case, the summands in the sum (23) have a finite mean $`Y_1`$ and finite variance $`\sigma _1^2`$. Thus, the Central Limit Theorem (CLT) holds (see for instance Sorbook  for a pedagogical exposition of the CLT). Let us recall the explicit expressions of the mean and variance:

$$Y_1=n=\frac{\kappa \gamma }{\gamma 1},\sigma _1^2=\frac{n^2}{\gamma (\gamma 2)}+n,$$
(43)
If $`r1`$ (in practice, it is sufficient that $`r6`$), the sum $`Y(r)`$ in (23) converges to a Gaussian variable with the following mean and variance

$$Y(r)=rn,\sigma ^2(r)=r\sigma _1^2.$$
(44)
Let $`\mu \kappa `$ be integer (just for illustrative purpose). Then, one can imagine $`X(\mu \kappa )`$ as a sum

$$X(\mu \kappa )=\underset{m=1}{\overset{\mu \kappa }{}}X_m$$
(45)
of independent summands $`\{X_1,X_2,\mathrm{},X_{\mu \kappa }\}`$ with the following PGF, mean and variance:

$$Q(z)=e^{z1},X_1=\sigma ^2=1.$$
(46)
If the number of summands $`\mu \kappa 1`$, then the sum (45) is asymptotically Gaussian as well, with mean and variance equal to $`\mu \kappa `$.

Thus if $`r1`$ and $`\mu \kappa 1`$, then the sum $`X(\mu \kappa )+Y(r)`$ is asymptotically Gaussian and we can use the following asymptotic formula

$$\text{Pr}\left[X(\kappa \mu )+Y(r)=m\right]=\frac{1}{\sqrt{2\pi \sigma _\mu ^2(r)}}\mathrm{exp}\left[\frac{(ma_\mu (r))^2}{2\sigma _\mu ^2(r)}\right].$$
(47)
Here $`a_\mu (r)`$ and $`\sigma _\mu ^2(r)`$ are respectively the mean and variance of the sum $`X(\mu \kappa )+Y(r)`$:

$$a_\mu (r)=rn+\mu \kappa ,\sigma _\mu ^2(r)=r\sigma _1^2+\mu \kappa .$$
(48)
Substituting (48) into (47) and (47) into (30), we obtain

$$p_{r,\mu }(r)\frac{\mu \kappa }{r\sqrt{2\pi (\sigma _1^2r+\mu \kappa )}}\mathrm{exp}\left(\frac{[(1n)r\mu \kappa 1]^2}{2(\sigma _1^2r+\mu \kappa )}\right).$$
(49)
For $`\gamma =\mathrm{}`$, for which we have the mean and variance defined in (43) reduce to $`\sigma _1^2=n=\kappa `$, relation (49) transforms into the previously announced result (39). In this case, we can test the quality of the general asymptotic relation (49) by comparing its particular application (39) with the exact expression (32) and its Stirling asymptotics (34). The exact and Stirling’s approximation essentially coincide while the “Gaussian” approximation is also excellent and goes closer and closer to the exact formula the larger $`\mu `$ is and the closer $`n=\kappa `$ is to $`1`$.

For $`\sigma _1^2r>\mu \kappa `$, expression (49) can be simplified into

$$p_\mu (r)\frac{\mu \kappa }{r\sqrt{2\pi \sigma _1^2r}}\mathrm{exp}\left(\frac{[(1n)r\mu \kappa 1]^2}{2\sigma _1^2r}\right),$$
(50)
which reduces again to the standard mean field asymptotics (40) at the critical point $`n=1`$ and for $`r\mu ^2\kappa ^2`$. Expression (50) is obtained using the following approximation

$$\sigma _1^2r+\mu \kappa \sigma _1^2r,$$
(51)
that is, by neglecting the term $`\mu \kappa `$. From a probabilistic point of view, this corresponds to using the Law of Large Numbers for the sum (45). Namely, if $`\mu \kappa 1`$, then one can replace (45) by its mean field limit

$$X(\mu \kappa )X(\mu \kappa )=\mu \kappa ,$$
(52)
which provides the truncated equality (51). Another reason for using the mean field limit (52) is more technical: in order to be able to effectively use the same technique in the case $`1<\gamma <2`$, we need to implement such a “mean field approximation.” Otherwise, the formulas would include very complicated convolutions of Gaussian and Levy stable laws.

We have checked the validity of the mean field approximation by comparing numerically the general asymptotics (49) and of its truncated mean field limit version (50) for various values of $`\gamma `$, $`n`$ and $`\mu \kappa `$. We find always an excellent convergence in the tail. For instance, consider the simplest case $`\gamma =\mathrm{}`$, for which the general asymptotic formula (49) transforms into (39), while its truncated mean field version is

$$p_{r,\mu }(r)\frac{\mu \kappa }{r\sqrt{2\pi \kappa r}}\mathrm{exp}\left(\frac{[(1\kappa )r\mu \kappa 1]^2}{2\kappa r}\right).$$
(53)
The agreement between all these expressions is excellent in the tails for all values of the parameters. The truncated mean field approximation describes satisfactory the body of distribution and becomes very precise in the tail of the distribution, even for moderate values of $`\mu 25`$ and $`n=\kappa <0.9`$. Even better agreement is obtained for larger $`\mu `$ and closer to the critical case $`n=1`$.

Let us complement this section by a few additional useful formulas. The expansion (19) with $`r_1=dG(z)/dz|_{z=1}`$ shows that the average number $`r_1`$ of first-generation daughters from mothers of arbitrary fertilities is $`n=\mu \kappa =\kappa \gamma /(\gamma 1)`$ and that, for $`1<\gamma <2`$, the variance of $`r_1`$ is infinite. For $`\gamma >2`$, the variance of $`r_1`$ is $`\sigma _1^2=n+n^2/[\gamma (\gamma 2)]`$. The form of (19) is essentially controlled by the expression of the characteristic function of the PDF (5) and it is thus not a surprise that the PDF of the number of daughters of first generation born from a mother with arbitrary fertility has the same asymptotic form as (5). For a mother with fixed fertility $`\mu `$, the average and variance of the total number of daughters are respectively

$$r_\mu =\frac{\mu \kappa }{1n};\sigma _\mu ^2=\frac{\mu \kappa }{(1n)^3}\left[1+\frac{n^2}{\gamma (\gamma 2)}\right],\mathrm{for}\gamma >2.$$
(54)
### IV.3 Case of infinite variance $`1<\gamma <2`$

We now turn to the novel regime $`1<\gamma <2`$. In this case, it is convenient to rewrite relation (15) in a more transparent probabilistic form

$$p_\mu (r)=\frac{\mu \kappa }{r}\text{Pr}\left[X+Y(r)=r1\right],Y(r)=\underset{k=1}{\overset{r}{}}Y_k,$$
(55)
where $`\{X,Y_1,\mathrm{},Y_r\}`$ are mutually independent random integers such that the PGF of $`X`$ is $`G_\mu (z)`$ given by (17) while the PGF’s of the remaining integers are equal to $`G(z)`$ given by (18). For $`1<\gamma <2`$, the variance of each variable $`Y_k`$ is infinite, for the same reason that $`\sigma _1^2`$ is infinite. From the generalized limit theorem Major , $`\stackrel{~}{Y}(r)=Y(r)nr`$ converges in distribution to a PDF $`f_y(s)`$ which is proportional to a stable infinitely-divisible PDF $`\phi _\gamma (x)`$ with exponent $`\gamma `$

$$f_y(s)=\frac{1}{\nu r^{1/\gamma }}\phi _\gamma \left(\frac{s}{\nu r^{1/\gamma }}\right),$$
(56)
where

$$\nu =\kappa [\gamma \mathrm{\Gamma }(\gamma )]^{1/\gamma },$$
(57)
and

$$\phi _\gamma (x)=\underset{0}{\overset{\mathrm{}}{}}\mathrm{exp}\left(u^\gamma \mathrm{cos}\left(\frac{\pi \gamma }{2}\right)\right)\mathrm{cos}\left(u^\gamma \mathrm{sin}(\frac{\pi \gamma }{2})+ux\right)𝑑u$$
(58)
with $`\phi _\gamma (0)=1/\mathrm{\Gamma }\left(1/\gamma \right)`$. Thus we obtain the distribution of aftershock number from (55) and (56)

$$p_\mu (r)\frac{\mu \kappa }{\nu r^{1+\frac{1}{\gamma }}}\phi _\gamma \left(\frac{(1n)r\mu \kappa 1}{\nu r^{1/\gamma }}\right).$$
(59)
For large number of aftershocks $`r>>R_1`$ where $`R_1`$ is defined by

$$R_1=\left(\frac{\mu \kappa +1}{\nu }\right)^\gamma \frac{\mu ^\gamma }{\gamma \mathrm{\Gamma }(\gamma )},$$
(60)
we can use the following asymptotic of $`\phi _\gamma (x)`$ for $`x\mathrm{}`$

$$\phi _\gamma (x)\frac{x^{1\gamma }}{\mathrm{\Gamma }(\gamma )}$$
(61)
Putting (61) in (59) retrieves at criticality ($`n=1`$) our main result announced in (7) with $`C_r=\mu \phi _\gamma (0)/[\gamma \mathrm{\Gamma }(\gamma )]^{1/\gamma }`$.

In the subcritical case ($`n<1`$), there is another power asymptotics. For $`rR_1`$, one can rewrite (59) in the form

$$p_\mu (r)\frac{\mu }{[\gamma \mathrm{\Gamma }(\gamma )]^{1/\gamma }r^\beta }\phi _\gamma \left(\frac{1n}{\nu }r^{\frac{\gamma 1}{\gamma }}\right).$$
(62)
If additionally $`rR_2`$, where $`R_2`$ is defined by

$$R_2=\left(\frac{\nu }{1n}\right)^{\gamma /(\gamma 1)},$$
(63)
then, using the asymptotic (61) of $`\phi _\gamma (x)`$, we obtain

$$p_\mu (r)\mu \gamma \left(\frac{n(\gamma 1)}{(1n)\gamma }\right)^{\gamma +1}\frac{1}{r^{1+\gamma }},$$
(64)
for $`r\text{max}\{R_1,R_2\}`$. If additionally $`R_2R_1`$, which holds if $`\frac{1n}{n}(\gamma 1)\mathrm{\Gamma }(\gamma )\mu ^{1\gamma }`$, then for $`R_1rR_2`$, there is an intermediate asymptotics following expression (7). These results are checked in Figure 2.

## V Distribution of the total number of generations

We now turn to the determination of the PDF of the total number of branching generations, in other words, to the probability that the branching process terminates at a given generation number. Let $`U_\mu (z;g)`$ (respectively $`U(z;g)`$) denote the PGF corresponding to the number of daughters born from a mother with fertility $`\mu `$ (respectively with arbitrary fertility) at the $`g`$-th generation. A standard result of branching theory is Athreya

$$U_\mu (z;g+1)=G_\mu (U(z;g));U(z;g+1)=G(U(z;g)),$$
(65)
where $`U(z;0)=z`$. The probability that the branching process survives at the $`g`$-th generation is

$$p_\mu (g)=1U_\mu (0,g).$$
(66)
The probability of termination of the branching process at the $`g`$-th generation is then given by

$$p_{g,\mu }(g)=p_\mu (g)p_\mu (g+1).$$
(67)
From (65,66), the probability for the branching process to survive at the $`g`$-th generation is

$$p_\mu (g+1)=1G_\mu (1p(g)).$$
(68)
$`p(g)`$ is the surviving probability at the $`g`$-th generation triggered by a mother of arbitrary fertility $`\mu `$, which obeys the recurrence

$$p(g+1)=1G(1p(g)),p(0)=1.$$
(69)
Equations (68) and (69) are easily solved numerically. We now extract the asymptotic power law behavior of $`p_{g,\mu }(g)`$. For small $`p(g)`$, the right-hand-side of (69) can be expanded as

$$1G(1p)np\frac{C}{\gamma 1}p^\gamma Dp^2,$$
(70)
where $`C=n^\gamma \mathrm{\Gamma }(2\gamma )[(\gamma 1)/\gamma ]^\gamma `$ and $`D=n^2(\gamma 1)^2/[2\gamma (\gamma 2)]`$. For $`\gamma >2`$, it is enough to take the leading behavior of (70) up to second order in powers of $`p`$: $`1G(1p)npDp^2`$. With (69), this gives $`p(g+1)p(g)=(1n)p(g)Dp^2`$. Close to criticality $`n1`$ for which a large number $`g`$ of generations occur, the leading behavior of the survival probability can be obtained by taking a continuous approximation to (69), giving $`\frac{dp(g)}{dg}=(1n)pDp^2`$. Its solution with initial condition $`p(g_0)=p_0`$ has the form

$$p(g)=\frac{p_0(1n)}{(1n+p_0D)e^{(1n)(gg_0)}p_0D}.$$
(71)
Expression (71) gives a power law $`p(g)1/(Dg)`$ for $`1g<1/(1n)`$ and crosses over to an exponential law for $`g1/(1n)`$. Knowing $`p(g)`$, $`p_\mu (g)`$ given by (68) can be obtained. For instance, in the case of a Poisson distribution (16) giving (17), we obtain $`p_\mu (g+1)=1\mathrm{exp}\left[\mu \kappa p(g)\right]`$. For $`\mu \kappa p(g)1`$, $`p_\mu (g+1)`$ becomes $`p_\mu (g+1)\mu \kappa p(g)`$. The sought PDF of $`p_{g,\mu }(g)`$ defined by (67) is given by

$$p_{g,\mu }(g)\mu \kappa \frac{dp(g)}{dg}\frac{\mu \kappa }{Dg^2}.$$
(72)
The last equality holds for $`g1/(1n)`$. This retrieves the standard mean field asymptotics $`p_{g,\mu }(g)1/g^2`$.

For $`1<\gamma <2`$, it is sufficient to keep only the following terms $`1G(1p)np\frac{C}{\gamma 1}p^\gamma `$ in the expansion (70). Taking again a continuous approximation of (69) gives $`\frac{dp(g)}{dg}=(1n)p\frac{C}{\gamma 1}p^\gamma `$, whose solution, for $`g1/[(\gamma 1)(1n)]`$, reads

$$p(g)\left[p_0^{1\gamma }+C(gg_0)\right]^{1/(\gamma 1)}\frac{C}{(gg_0)}^{1/(\gamma 1)},$$
(73)
Correspondingly, the PDF $`p_{g,\mu }(g)`$ that the number of generations is exactly $`g`$ is given by (7) with $`C_g=\frac{\mu n}{\gamma }C^{1/(\gamma 1)}`$, valid for $`\gamma `$ not too close to $`2`$.

## VI Concluding remarks

We have shown that the existence of cascades of triggered seismicity produces huge fluctuations of the total number of aftershocks and of the total number of generations from one sequence to another one for the same mainshock magnitude, characterized by power-law asymptotics. If the distribution of offsprings per mothers (of direct aftershocks per mainshock for seismicity) has a finite variance, we recover the well known mean-field results (9). In the case of an infinite variance relevant for earthquakes, we have discovered a new regime with exponents that varies continuously with the exponent $`\gamma =b/\alpha `$. The anomalous scaling reported here for $`1<\gamma <2`$ gives rise to less wild fluctuations in the total number of daughters from one mother to the next, compared with the mean field regime $`\gamma 2`$. For instance, for earthquakes, we have probably $`\gamma 1.25`$ for the preferred values $`b1`$ and $`\alpha 0.8`$ alpha , which leads to $`p_r(r)1/r^{1+0.8}`$ and $`p_g(g)1/g^{1+4}`$ at $`n=1`$, compared with $`p_r(r)1/r^{1+0.5}`$ and $`p_g(g)1/g^{1+1}`$ for $`\gamma 2`$. The reason for this behavior lies in the important role exerted by the mothers with the largest fertilities $`\mu `$ on the rate of daughter births. In particular, this role explains why, away from criticality but close to it ($`n1`$), $`p_r(r)`$ crosses over from $`1/r^{1+\frac{1}{\gamma }}`$ to $`1/r^{1+\gamma }`$ as $`r`$ increases: the latter asymptotic regime is nothing but the tail of the PDF (5) which controls the number of daughters born in the first generation from a mother with arbitrary fertility.

How relevant is the infinite variance regime $`\gamma <2`$ in view of the uncertainties in the exponents $`\alpha `$ and $`b`$? Indeed, the value of $`\alpha `$ is not only approximately known but its constancy (as a function of time) and uniformity (as a function of space) remains to be ascertained. In addition, the $`b`$-value also fluctuates and is typically in the range $`0.751.25`$ Frohlich . If we take the largest $`b`$-value of this interval ($`b=1.25`$), we remain in the infinite variance regime $`\gamma <2`$ as long as $`\alpha >0.63`$. Most present estimations discussed in section 2 give values above this lower bound. It thus seems that the results presented here should remain pertinent as more precise estimations of $`\alpha `$ and $`b`$ become available.

All these discussions assume untruncated power law distributions. However, it is well-known that the Gutenberg-Richter law exhibits an upper magnitude cut-off Kagan11 ; Pisor22 . By the logic leading to (4), this automatically implies also an upper fertility cut-off. This in turn removes the divergence of the variance for $`\gamma <2`$. How does this affect the distributions $`p_r(r)`$ and $`p_g(g)`$? This question is standard in the theory of power laws (see for instance Sorbook  and references therein). Since the upper cutoff is very large so that the actual range of events is very large, our results obtained without truncation hold for a large range of values of $`r`$ and of $`g`$ respectively. However, the predicted power laws (7) have to cross-over to the mean-field ones (9) at threshold values beyond which the finiteness of the variance of fertilities become flagrant. In contrast, the presence of a lower cutoff magnitude benzion  does not modify our results on the power law tail and has an influence only in shaping the bulk (small values) of the distributions.

The direct validation of our predictions (7) and (8) on earthquake catalogs is not feasible, due to the impossibility to distinguish aftershocks from uncorrelated events at large times after the mainshock, and due to the limited number of large aftershock sequences. Our results have however some consequences for the statistical properties of aftershock sequences. We have indeed shown in Bathlawpap  that the existence of large fluctuations of the number of aftershocks per mainshock (summed over all generations) induces a non-trivial scaling of the difference in magnitude between a mainshock and its largest aftershock, so that Båth law can be recovered in the regime $`b\alpha `$. These results also suggest that the common use of a Poisson distribution in seismicity forecasts to model the distribution of the number of events within a finite space-time window is questionable, since the simple physics of cascades of earthquake triggering gives a very different (power law) distribution. Recent observations show non-Poisson power law distributions of seismic rates (see PG96  and work in progress). We suggest that these observations and our results could be used to improve earthquake forecasting by providing a more realistic distribution of the number of events.

Acknowledgments: We are grateful to two unknown referees and to Y. Ben-Zion as the editor, for useful remarks that helped improve the manuscript. This work is partially supported by NSF-EAR02-30429, by the Southern California Earthquake Center (SCEC) and by the James S. Mc Donnell Foundation 21st century scientist award/studying complex system. SCEC is funded by NSF Cooperative Agreement EAR-0106924 and USGS Cooperative Agreement 02HQAG0008. The SCEC contribution number for this paper is 741.