format: text

text:
# Electromagnetic Catastrophe in Ultrarelativistic Shocks and the Prompt Emission of Gamma-Ray Bursts

## 1 Introduction

Astrophysical objects where we observe (active galactic nuclei) or infer (gamma-ray bursts, soft gamma-repeaters) a relativistic bulk motion demonstrate surprisingly high gamma-ray luminosity. In some cases the energy release in gamma-rays can be comparable to the total kinetic energy of the bulk motion in the source. This implies a very efficient mechanism of the energy conversion of the bulk motion into the hard radiation. It is clear that in the case of explosive phenomena like GRBs the radiation originates from ultrarelativistic shocks (Cavallo & Rees 1978; Rees & Mészáros 1992). This is less certain in the case of AGN jets, however there exist a large number of works considering shocks as the main source of AGN jet radiation (for a review see Kirk, 1997).

One can describe the energy stored in the shock in two ways: as a total kinetic energy of ejecta and shocked matter or as their internal energy. There are known two different ways of shock energy release: the first order Fermi acceleration and different mechanisms of internal energy dissipation. Fermi acceleration does not require the dissipation as it is fed directly by the kinetic energy of the bulk motion. On the other hand, Fermi acceleration has a number of limiting factors in available power and radiative efficiency (see, e.g., Bednarz & Ostrovsky 1999) . With this reason Fermi acceleration is usually considered as a source of high energy cosmic rays while the intensive radiation is traditionally explained in terms of the internal energy dissipation.

This work suggests a highly radiative analog of Fermi acceleration where the main role is played by photons rather than by charged particles. A similar idea about important role of neutral particles in shocks has been suggested independently in a very recent paper by Derishev et al. (2003). For a simple illustration the shock and the external medium can be represented as two mirrors moving with an ultrarelativistic velocity respectively to each other. Let $`\xi _s`$ and $`\xi _e`$ be their reflection coefficients (depending on the kind and the energy of a particle) and $`\mathrm{\Gamma }`$ – the Lorentz factor of the shock. A photon being elastically reflected head-on from a moving mirror gains a factor $`\mathrm{\Gamma }^2`$ in its energy. Thus the total energy of photons participating in back and forth motion between the mirrors can be roughly described by the equation

$$dE/dt=(\xi _s\xi _e\mathrm{\Gamma }^2S)E/t_c$$
(1)
where $`t_c`$ is a characteristic time of reflection cycle and $`S`$ is the probability of the particle escape. Actually equation (1) represents not more than an instructive toy model since a real situation is essentially non-local ($`t_c`$ is uncertain) and require an integration through particle spectra. If $`\xi _s\xi _e\mathrm{\Gamma }^2>S`$, the total energy of particles will grow exponentially until mirrors decelerate.

In the case of the first order Fermi acceleration the reflection results from the diffusion of supra-thermal charged particles in magnetic field. How this two-mirror approach will look when applied to photons? If a soft photon (e.g. of a synchrotron nature) interacts only with electrons, the reflection coefficients $`\xi _s`$ and $`\xi _e`$ are of order of the Thomson optical depth of each medium ($`\xi _s\xi _e\tau _T`$, see section 2). This case does not differ essentially from the well known thermal Comptonization in relativistic regime at electron temperature $`\mathrm{\Gamma }1`$. The photon boosting terminates as soon as the photon gains a sufficient energy to interact in Klein-Nishina regime when reflection coefficients decrease. The Compton boosting is probably not very important by itself, however it could provide photons above pair production cutoff and give rise to a new much more powerful effect.

The importance of pair production in GRB shocks and AGN jets has been emphasized in many works, e.g., by Thompson (1997), Ramires-Ruis et al. (2000) and Mészáros, Ramires-Ruis & Rees (2001). Thompson (1997) and Ramires-Ruis et al. (2000) stated that the pair production can enhance the radiation efficiency of the shock providing new supra-thermal particles as seeds for Fermi acceleration. In this work we show that the pair production plays a crucial role initiating a runaway energy release.

A necessary condition to trigger this process is the existence of a seed radiation field in the vicinity of the shock. This radiation should be dense in terms of the photon number: $`n_pR\sigma _T1`$ where $`R`$ is a characteristic scale, both $`n_p`$ and $`R`$ correspond to the comoving frame of the shock where photons are approximately isotropic. Such photons can appear e.g. as the result of an internal energy dissipation and undergo the Compton boosting.

Now let a photon of energy $`\epsilon 1`$ (hereafter $`\epsilon `$ is the photon energy in electron mass units) moves upstream across the shock front. Then, interacting with a soft photon, it produces an e+e- pair in unshocked matter ($`\xi _e1`$ in this case). The pair is bound to the rest frame of the external matter if the e+ and e- Larmor radii do not exceed the scale characterizing the problem. Particles are advected downstream across the shock front, losing a fraction of their energy to Comptonization. Some of downstream Comptonized photons can produce new pairs behind the shock front which in turn upscatter soft photons upstream. So we deal with multiple electromagnetic interactions and the whole process can be characterized as an electromagnetic cascade fed by kinetic energy of two converging media or cascade boosting. The energy involved in the cascade can grow exponentially and eventually the process can reach a nonlinear stage when the boosting rate increases.

All above conclusions have been made in assumption that the shock front separating unshocked and shocked media is very sharp and can be characterized as a viscous jump. It is known from studies of nonrelativistic shocks that the viscous jump can disappear transforming into a smother transition due to radiation pressure ahead the shock. Thompson & Madau (2000) have studied preacceleration of external matter for the case of a GRB external shock and demonstrated that the pair production (pair loading) is extremely important in this respect and feeds back the radiation efficiency. Beloborodov (2002) has considered a more extreme case when the radiation of GRB ejecta sweeps up the external medium far ahead the shock. In that case the energy should be supplied from internal shocks. Our case is closer to that of Thompson & Madau (2000). The preacceleration of external medium transforms the viscous jump into a smoother velocity gradient. This transformation removes the radiative divergence: now an upstream photon most probably interacts in preaccelerated medium and the energy gain in the interaction becomes much less than $`\mathrm{\Gamma }^2`$. We should expect the formation of some steady state regime where the fluid velocity pattern is adjusted in a way to satisfy the energy conservation law. The luminosity at this stage should be close to its ultimate value $`L\mathrm{\Gamma }^2dm/dt`$ where $`m`$ is the mass of the external matter swept up by the shock.

Summarizing, we can expect the following evolution: (i) a seed radiation due to dissipation of the internal energy; (ii) Compton boosting and the cascade boosting of the electromagnetic component in linear exponential regime; (iii) a sharper nonlinear evolution when the boosted high energy component contributes to the soft target photon field; (vi) transformation of the shock front and setting up a highly luminous quasi-steady state regime. The entire event is so dramatic that it deserves the name “electromagnetic catastrophe”.

This qualitative illustration is confirmed below at a quantitative level. The main tool in this study is the large particle nonlinear Monte-Carlo (LPMC) code described in Stern et al. (1995).

While the mechanism can work in different astrophysical phenomena this study is focused on its application to GRBs. Correspondingly, the explored region of parameters represents the prevailing view on GRBs (large Lorentz factors, moderate densities).

In Section 2 we describe the general formulation of the problem, the accepted assumptions and the range of parameters. Section 3 briefly describes the approach to the numerical simulations of the effect. Section 4 presents the results of simulation runs demonstrating the evolution of the system from a low luminosity initial state through the catastrophe to a highly radiative post-catastrophe regime. Section 5 extends the range of the phenomenon to lower densities and outlines the threshold where the catastrophe is possible at certain assumptions. Section 6 describes a phenomenon of GeV bursts which should arise when the compactness is insufficient to provide optically thick pair loading. Possible associations with existing data are discussed. In section 7 we discuss various GRB scenario associated with electromagnetic catastrophe and the issue of the GRBs time variability.

## 2 Formulation of the problem

### 2.1 The simplified description of the shock

For convenience we consider the shock in the comoving system which is preferable for its higher symmetry: the density contrast is moderate (factor 4 in an idealized case), upstream and downstream energy fluxes can be comparable. More explicitly, we bound our reference system to the contact discontinuity, its Lorentz factor hereafter is $`\mathrm{\Gamma }`$. The shock front in this frame moves upstream with semirelativistic velocity $`\beta _s`$ (hereafter, $`\beta `$ is the velocity in units of the velocity of light). We assume that the shocked external matter between the shock front and the contact discontinuity is of a constant density and in rest relatively to the contact discontinuity. Actually this region can be subject to a semirelativistic bulk motion (perhaps with a turbulent component) depending on the distribution of external matter and shock dynamics. Neglecting such motion we remove an additional source of the free energy and since we are interested in an explosive energy release this is a conservative assumption.

We do not know and do not consider the state of ejecta behind the contact discontinuity, prescribing to it a role of a heavy piston, carrying the main fraction of the kinetic energy. It is assumed that ejecta does not emit or reflect radiation. Actually ejecta can contribute the seed radiation and the downstream-upstream reflectivity.

We assume a spherical symmetry at least within a cone of opening angle $`1/\mathrm{\Gamma }`$ respectively to the line of sight. Due to relativistic effects we almost do not observe the emission of the shock beyond this cone. Events in the shock separated by angle more than $`1/\mathrm{\Gamma }`$ are not causally connected since the acceleration stage. Most probably the ejecta is strongly beamed within a solid angle $`\mathrm{\Omega }1`$, we exclude $`\mathrm{\Omega }`$ from the consideration using the isotropic equivalents for global energy parameters: the kinetic energy of ejecta, $`E`$, and the luminosity, $`L`$. The actual energetics is then described by $`\mathrm{\Omega }E`$ and $`\mathrm{\Omega }L`$.

The spherical shock viewed from the comoving system is parabolically curved and the matter moves outward from the line of sight, relativistically at the angular displacement $`1/\mathrm{\Gamma }`$ from the line of sight. For simplicity we adopt one-dimensional slab geometry with no bulk motion in the axial direction. This is not a conservative simplification because in this way we violate kinematic and causality constraints on interactions of photons emitted from distant regions of the shock. In order to avoid this problem, the consideration was restricted to a comparatively small time (less than a half deceleration time) when the above constraints are still not important.

The unshocked interstellar matter (ISM) in the comoving system is a cold fluid with bulk Lorentz factor $`\mathrm{\Gamma }`$. The ordered particle motion dissipates into chaotic motion with the same average thermal Lorentz factor $`\mathrm{\Gamma }`$ when the fluid crosses the shock front. The energy density contrast across the front and the front velocity in the comoving system, $`\beta _s`$, result from the energy conservation:

$$U_s\beta _s=U_\mathrm{o}(\beta +\beta _s),$$
(2)
where $`U_s`$ is the energy density of the shocked matter,$`\beta `$ is the external fluid velocity, $`U_\mathrm{o}=\mathrm{\Gamma }^2\rho m_pc^2`$, $`\rho `$ is the rest frame ISM number density, and from the pressure – momentum balance:

$$P_s=\frac{1}{3}U_s=U_\mathrm{o}\beta (\beta +\beta _s).$$
(3)
In the ultrarelativistic limit ($`\beta =1`$) equations yield $`\beta _s=1/3`$ $`U_s=4U_\mathrm{o}`$. For a more detailed description of a relativistic shock see Blandford & McKee (1976).

The thickness of the shock (between the shock front and the contact discontinuity) at the constant external density and spherical geometry is $`\mathrm{\Delta }X=\frac{1}{12}R/\mathrm{\Gamma }`$ in the comoving system.

### 2.2 Main parameters and dimensionless units

External magnetic field $`H_e`$ is different for cases of the stellar wind dominated and the ISM dominated environments. In the first case, the field is predominantly radial and processes leading to the catastrophe are inhibited ($`\xi _e`$ is very low) unless the field is affected by radiation front ahead the shock. Most probably it must be affected by the two stream instability generating a chaotic component. Moreover, a spherical asymmetry of the environment and ejecta together with the charge asymmetry (Compton scattering on electrons) can induce a large scale transversal field similar to that induced in an atmospheric nuclear explosion.

In the ISM dominated case, the field originally should have a large scale geometry and a substantial transversal component. In this work only the ISM case is studied quantitatively, note that the interaction of the shock with a slow shell of matter ejected by a GRB progenitor essentially does not differs from the ISM case. ISM magnetic field depends on the external density $`\rho `$ and the reasonable assumption is this case is $`H_e10^5\sqrt{\rho }`$ G, where $`\rho `$ is in cm-3 units. Hereafter we treat $`H_e`$ as the transversal component.

The lower limit for the magnetic field in the shock $`H`$ follows from Lorentz transformation and the compression by factor 4: $`H>H_g=4H_e\mathrm{\Gamma }`$. However, this field can be much stronger being generated at the shock. At this step we can rely on numerical study of the two stream instability by Medvedev & Loeb (1999). According to their simulations the ensured value of energy density of generated field, $`U_B`$, is $`0.010.1`$ of electron energy density in the shock, $`U_e`$, and, tentatively, the field can be amplified up to $`U_B0.010.1`$ of proton energy density.

The parameters of the problem are: initial Lorentz factor $`\mathrm{\Gamma }`$, total kinetic energy $`E`$, external number density $`\rho `$, characteristic rest frame size $`R`$, seed luminosity of soft photons $`L_s`$, their spectrum, values of magnetic field in both regions, $`H`$ (shock) and $`H_e`$. As a reasonable hypothesis for the seed spectrum we accept a power law $`dN/d\epsilon \epsilon ^\alpha `$, in the range $`\epsilon _1<\epsilon <\epsilon _2`$ (see §2.3).

In the case of GRBs the isotropic kinetic energy of ejecta can vary within two orders of magnitude or more. Maximal observed hard X-ray – gamma-ray GRB energy release is $`510^{54}`$erg. We accept $`E=10^{54}`$ erg as a baseline. In this work we study the case of a constant external density which is a reasonable starting point. The scale of the problem can be characterized by the deceleration radius (Rees & Mészáros 1992):

$$R_d=\left[E/(m_pc^2\frac{4}{3}\pi \rho \mathrm{\Gamma }^2)\right]^{1/3}=0.2410^{16}E_{54}^{1/3}\rho _6^{1/3}\mathrm{\Gamma }_2^{2/3}\mathrm{cm}$$
(4)
In the spherical geometry, the scale changes as the shock propagates, therefore we adopt $`R=0.5R_d`$ which approximately describes the shock propagation between $`0.5R_d`$ and $`R_d`$. Using the rest frame scale $`R`$ we can define dimensionless distance, $`x`$, and time, $`t`$, to work within the comoving system: $`x=X/(R/\mathrm{\Gamma })`$, $`t=T/(R/\mathrm{\Gamma }c)`$, where $`X`$ and $`T`$ are comoving distance and time. Note that the allowed working range is $`x<1`$ and $`t<1`$. The time in observer frame is $`T_{\mathrm{obs}}=tR/(2\mathrm{\Gamma }^2c)`$.

It is convenient to describe the energy budget in terms of dimensionless compactness. The latter defines the importance of pair production and, more generally, the level of nonlinearity of the system. The simplest way to introduce this parameter is via the energy column density through the system normalized to the electron mass and multiplied by Thomson cross section:

$$\omega =\frac{U}{m_ec^2}\sigma _T\frac{R}{\mathrm{\Gamma }},$$
(5)
where $`U`$ is an energy density in the comoving system. Particularly, the magnetic compactness is:

$$\omega _B=\frac{H^2}{8\pi m_ec^2}\sigma _T\frac{R}{\mathrm{\Gamma }}.$$
(6)
This is the compactness describing the energy content of the system. Traditionally the compactness is expressed through the luminosity of the system $`L`$:

$$\mathrm{}=\frac{L\sigma _T}{Rm_ec^3}.$$
(7)
The rest (observer) frame compactness for the ultrarelativistic shock is meaningless since the radiation is strongly collimated. The parameter has the stated above physical meaning only in the shock comoving system where particles have a wide angular distribution. Assuming the comoving size of the system $`R/\mathrm{\Gamma }`$ (with the transversal size of the same order, see §2.1) and substituting $`L`$ by the luminosity flux: $`LF(R/\mathrm{\Gamma })^2`$, with the accuracy of the order of $`\pi `$ one gets:

$$\mathrm{}=\frac{F}{m_ec^3}\sigma _T\frac{R}{\mathrm{\Gamma }}.$$
(8)
The total energy budget in the comoving frame constitutes the kinetic energy flux of the fluid crossing the shock front. The corresponding dimensionless expression, kinetic compactness, is:

$$\mathrm{}_\mathrm{o}=\rho \mathrm{\Gamma }^2m_pc^3\frac{\sigma _T}{m_ec^3}\frac{R}{\mathrm{\Gamma }}=\tau _T\frac{m_p}{m_e}\mathrm{\Gamma }.$$
(9)
Hereafter we denote energy compactness as $`\omega `$ and the luminosity (power) compactness as $`\mathrm{}`$, and $`\mathrm{}=d\omega /dt`$.

### 2.3 Seed radiation

The seed radiation in the shock can be contributed by (i) the synchrotron – self-Compton radiation of shocked electrons which have the thermal Lorentz factor $`\mathrm{\Gamma }`$; (ii) other channels of internal energy dissipation (e.g. dissipation of magnetic field, plasma waves, turbulence); (iii) external (partially side-scattered) photons from the fireball and a progenitor star; (iv) high energy seed photons associated, e.g, with Fermi acceleration and photo-meson production.

The most confident source is (i) unless synchrotron radiation of shocked electrons is strongly self-absorbed. Its fraction in the total energy budget is limited by the factor $`m_e/m_p`$. At some parameters, this seed radiation is sufficient to start the evolution leading to the catastrophe (see Section 3). In such case we deal with the most conservative “minimal hypothesis”. At other parameters the minimal hypothesis is insufficient for the catastrophic evolution. Then we have have to consider additional channels of internal energy dissipation (ii) or external photons (iii) and (iv). We would like to avoid a serious consideration of corresponding mechanisms in this work and use a very moderate ad hoc hypothesis of the seed soft photon field.

As a baseline, we accept the compactness of such emission $`\mathrm{}_s\tau _T\mathrm{\Gamma }`$ which is only $`m_e/m_p`$ fraction of the total energy budget. The model spectrum of this seed radiation has been taken as a power law fragment with $`2.5<\alpha <1.5`$ (a free parameter) in the range between $`\epsilon _1`$ (typically $`\epsilon _110^8`$ which is close to the self-absorption cutoff at $`H100`$ G) and $`\epsilon _2`$ (a free parameter). Such weak seed radiation is sufficient to give rise to the catastrophe in the most of trials described in this work. The situations when the effect implies a stronger seed radiation are discussed in section 7.

### 2.4 Reflection of the energy flux

Reflection of a particle flux from the external environment and from the shock is of a primary importance for the cascade boosting (see Eq. 1). Here we consider the reflection of a high energy component of the electromagnetic cascade where the direct Compton reflection of photons is inefficient. In this case the reflectivity of the medium is due to photon-photon generated pairs gyrating in the magnetic field. The requirement that an electron can change its direction to the opposite without a substantial energy loss is a serious limiting factor for the boosting rate.

In the case of the first order Fermi acceleration the main process responsible for the particle scattering upstream across the shock is the diffusion in a chaotic magnetic field (for a recent review see Gollant 2002). In our case the charged particle does not need to cross the shock, it is sufficient if it just changes its direction upstream. Then photons Comptonized by the particle moving upstream can cross the shock. The deflection of a particle from downstream to upstream hemispheres is a much more probable event than its diffusion upstream the shock, required for Fermi acceleration. Indeed, in the comoving frame of the shocked matter the shock front moves ahead semirelativistically, with $`\beta _s=1/3`$, and a charge particle has to diffuse upstream faster.

The downstream – upstream reflection of a particle can proceed in the following ways:

(i). If the Larmor radius is less than the correlation length of the magnetic field, the particle turns at a half Larmor orbit. This can be the case for particles of a moderate energy.

(ii). If the Larmor radius exceeds the correlation length, the particle diffuses in the chaotic field.

(iii). Except the chaotic field $`H`$ there could exist a large scale component, resulting from a large scale external magnetic field due to the flux conservation: $`H_g4H_e\mathrm{\Gamma }`$. Then a high energy particle can turn around at a half Larmor orbit in the field $`H_g`$. Processes (ii) and (iii) are competing and the leading one depends on concrete parameters. In this work only the effect of the global component has been taken into account since the latter should exist in the adopted model of the external environment. The effect of the diffusion is the matter for future studies.

The reflection efficiency depends on the fraction of energy which the particle loses before it turns around upstream. The energy loss rate is described by

$$d\gamma /dX=\left(\frac{H}{8\pi }+U_{\mathrm{ph}}\right)\frac{\sigma _T}{m_ec^2}\gamma ^2=C\gamma ^2,$$
(10)
where $`X`$ is the comoving distance. The equation yields

$$\gamma (X)=\gamma _\mathrm{o}/(1+\gamma _\mathrm{o}XC),$$
(11)
where

$$C=H^2(1+U_{\mathrm{ph}}/U_B)3.310^{20}\mathrm{cm}^1.$$
(12)
where $`U_B`$ and $`U_{\mathrm{ph}}`$ are energy densities of magnetic field and soft radiation field respectively. The term $`U_{\mathrm{ph}}/U_B`$ describes the relative contribution of Compton losses. From equation (11) one estimates the maximum energy of downstream – upstream reflected particle $`\gamma _{\mathrm{max}}1/XC`$. Substituting X by the half Larmor orbit

$$\pi R_L=5.310^3\frac{\gamma _{\mathrm{max}}}{H_g}cm,$$
(13)
one obtains the final expression for the maximal energy of a reflected particle:

$$\gamma _{\mathrm{max}}=0.7510^8\frac{H_g^{1/2}}{H}(1+U_{\mathrm{ph}}/U_B)^{1/2}.$$
(14)
The downstream – upstream reflection is inefficient if $`U_{\mathrm{ph}}/U_B1`$. Then the main energy of an electron is emitted as a relatively soft synchrotron radiation rather than a hard Comptonized photons. The limiting energy of synchrotron photons emitted by an electron after half Larmor orbit is $`\epsilon 200`$ at $`H_gH`$. This energy is still sufficient for pair production, however the cascade boosting in this case will be substantially slower and the compactness threshold for the electromagnetic catastrophe will be higher. In this work we study the case of low $`H`$ and weak dissipation. If the generated field is much stronger but its dissipation into electromagnetic component is fast, then we have $`U_{\mathrm{ph}}/U_B1`$ again. The boosting is inhibited if the field is strong and its dissipation is slow ($`U_{\mathrm{ph}}U_B`$) .

In the case of the external environment (upstream-downstream reflection) the electron synchrotron losses are small. Indeed, in the rest frame equation (14) with $`H_g=H=10^5\sqrt{\rho }`$ yields

$$\gamma _{r,\mathrm{max}}=1.310^{10}\rho _6^{1/2}(1+U_{\mathrm{ph}}/U_B)^{1/2}$$
(15)
where the subscript $`r`$ refers to the rest frame and $`\gamma _{\mathrm{max}}\mathrm{\Gamma }\gamma _{r,\mathrm{max}}`$. The term $`U_{\mathrm{ph}}/U_B`$ is large in the external environment, (i.e. Compton losses are much larger). However, Comptonized photons have predominantly downstream direction and contribute to the cascade boosting.

The main limiting factor in upstream – downstream reflection is the requirement that the Larmor radius should be smaller than the characteristic size of the problem. Actually, the upstream particle can be advected back across the shock front when it deflects bu just the angle $`1/\mathrm{\Gamma }`$. In this case it has a much less energy gain than $`\mathrm{\Gamma }^2`$. Both effects were accounted for in the numerical simulation. From the limit on the gyroradius $`R_L<R`$ equations (4) and (13) yield:

$$\gamma _{r,\mathrm{max}}<1.410^{12}H_eE_{54}^{1/3}\rho _6^{1/3}\mathrm{\Gamma }_2^{2/3}.$$
(16)
Using the hypothesis $`H_e=10^5\sqrt{\rho }`$ and transforming to the comoving frame one gets:

$$\gamma _{\mathrm{max}}<1.410^{11}E_{54}^{1/3}\rho _6^{1/6}\mathrm{\Gamma }_2^{1/3}.$$
(17)
Note that the limits (14) and (16) have an opposite dependence on the magnetic field. A faster boosting takes place in the case of a stronger external field and a weaker field in the shock. The above estimates are given for the orientation while the numerical simulation reproduces directly the particle trajectory and the energy loss along it. The quantitative effect of these parameters on the boosting rate and examples of upstream and downstream particle spectra are demonstrated in Section 5.

## 3 Monte-Carlo simulations

The entire problem outlined above is too difficult for an analytic treatment. In this work it is studied numerically using a Large Particle Monte-Carlo code (LPMC) developed by Stern (1985) and Stern et al. (1995). The code is essentially nonlinear: the simulated particles constitute at the same time a target medium for other particles. Large particle (LP) method in this application means that each simulated particle represents $`w`$ real particles (for GRBs $`w>10^{50}`$). The most reasonable weighting scheme except some special cases is to set $`w\epsilon =const`$, where $`\epsilon `$ is the energy of the particle and the constant can vary since the total energy involved in the simulation can change. The number of LPs was $`2^{17}=131072`$.

The version of LPMC used here treats Compton scattering, synchrotron radiation, photon-photon pair production, and pair annihilation. Synchrotron self-absorption was neglected as it consumes too much computing power and is not very important in this application. All these processes are reproduced without any simplifications at the microphysics level. On the other hand, a number of serious macroscopic simplifications has been done.

Since the number of large particles which is limited by the computing power available at the moment did not allow full three-dimensional simulation, the problem was reduced to one dimension, $`x`$, along the shock propagation. While locations and momenta of LPs were three dimensional, the target density was averaged over slab layers and the fluid representation was one-dimensional.

The trajectories of electron and positrons in the magnetic field were simulated directly assuming transversal geometry of the field $`H_g`$ in the shock and $`H_e`$ in the external matter. Thus the problem of the reflection qualitatively discussed in §2.4 was implemented in the numerical simulation. Protons were assumed to be cold in the fluid frame unless they have crossed the dissipative shock front.

Hydrodynamic part of the problem becomes very difficult since the energy and momentum transferred by LPs to the fluid fluctuate. Attempts to simulate the fluid with internal pressure have led to rising instabilities at semirelativistic velocities. Therefore a dust approximation (where the internal fluid pressure is neglected) has been adopted and the fluid velocity behind the shock front was artificially fixed to zero.

The dust approximation works satisfactory until the fluid has a moderate temperature and energy density in its comoving frame (i.e. when it is ultrarelativistic in the shock comoving frame). This circumstance was used by Beloborodov (2002) and Thompson & Madau (2000) who also described the preacceleration of the external medium by the radiation front in terms of dust approximation. The approximation does not work at all in the shock where the pressure at a relativistic temperature is comparable to the radiation pressure. The adopted hydrodynamic anzats allows to describe the effect in general and fails to account for possible interesting effects associated with semirelativistic motion behind the radiation front. The approach can be considered as a “zero approximation”.

## 4 Evolution of the system through the catastrophe

To illustrate qualitative arguments discussed in Introduction with numerical simulations, we present results of two runs with different parameters. First, for Lorentz factor $`\mathrm{\Gamma }=100`$ and the density of external environment $`\rho =310^5`$ cm-3. Second, for $`\mathrm{\Gamma }=30`$, $`\rho =10^6`$cm-3. Other parameters for these cases are $`R=1.910^{15}`$ cm, $`\mathrm{}_\mathrm{o}=63`$, $`T_\mathrm{o}=3.1`$s and $`R=2.810^{15}`$cm, $`\mathrm{}_\mathrm{o}=94`$, $`T_\mathrm{o}=52`$s. The first case is closer to the traditionally assumed value of $`\mathrm{\Gamma }`$ for GRBs, however it is more difficult for numerical treatment since statistical fluctuations associated with finite number of large particles being amplified by factor $`\mathrm{\Gamma }^2`$. With this reason the first run tracks the evolution of the system only till the beginning of the post-catastrophe stage.

Initial state is empty of photons and consists only of fluid of cold protons and electrons with bulk Lorentz factor $`\mathrm{\Gamma }`$, representing external medium viewed from the comoving system. Seed photons appear as the (partially self-absorbed) synchrotron radiation of shocked electrons. In this case, the self-absorption is not very strong and can be ignored. The Comptonization of the synchrotron peak gives second and third peaks at higher energies (further Comptonization is inefficient because it proceeds in Klein-Nishina regime). These peaks are visible in Figure 1 at $`t=0.3`$. The Comptonized photons are produced both by shocked and external electrons. The latter cause an anisotropy: the second peak is stronger in downstream photons (Figure 1, lower panel), while the third peak is slightly more intensive upstream (Figure 1, upper panel).

At $`t=0.4`$ a high energy tail appears due to pair production by photons of the third peak in the external medium (which provides comoving energy gain by factor $`\mathrm{\Gamma }^2`$). Then, around $`t=0.5`$, the evolution accelerates dramatically, the total energy of photons rises by more than 3 orders of magnitude (the power supplied to photons jumps by factor $`m_p/m_e`$), and the spectrum transforms from harder to softer one. The evolution from $`t=0.49`$ to $`t=0.52`$ will take about 0.1 second from the point of view of an external observer. The simulation has been terminated at $`t=0.56`$ because numerical problems associated with the high Lorentz factor become further serious. These problems are not principal and can be overcome just with the increase in the computing power. The evolution of the system with time is shown in Figure 3.

This example does not require any additional source of energy except the bulk kinetic energy of the fluid and any external photons. In this sense the case can be interpreted as a “minimal model”, however one still has to assume a non-radial external magnetic field.

The case of $`\mathrm{\Gamma }=30`$, $`\rho =10^6`$cm-3 does require additional assumptions. First of all, the synchrotron radiation of electrons with $`\epsilon =30`$ is deeply self-absorbed. Therefore an internal energy source accelerating a fraction of shocked electrons up to $`\epsilon 200300`$ is required. The required power is, nevertheless, much less than the total energy budget: $`\mathrm{}_sm_e/m_p\mathrm{}_\mathrm{o}`$. To avoid extra details we just ignored the self-absorption since the nature and energy of seed photons is not very important. Another problem is that the energy of the second Comptonization peak is smaller than in the previous example and the catastrophe develops too slowly. In section 5 it is shown that a much faster evolution can be induced by seed high energy photons which population undergoes an exponential breeding. In this run a portion of upstream high energy photons with total energy compactness $`\omega _s10^5`$ (which is $`10^7`$ of the total energy budget) has been injected at $`t=0.25`$. Then the catastrophe has followed at $`t0.3`$.

Figure 2 shows post-catastrophe photon spectra which can be treated as a continuation of the spectra set presented in Figure 1 to later ages. Figure 4 shows the spatial distribution of main quantities: the energy density of upstream photons, the pair number distribution and the bulk Lorentz factor. At the moment of the catastrophe the main fraction of photon energy is tightly concentrated to the fluid velocity jump which has a step-like character (see Figure 4, upper and middle panels, curve at $`t=029`$). The reason is that the evolution at that moment is so fast, that abruptly released photons have no time to disperse wider. Then photons disperse and a leading fraction of them form a radiation front moving with the speed of light. The distribution of the bulk Lorentz factor smooths and adjusts to the radiation front (Fig. 4, middle panel) moving ahead with a near-luminal velocity. The absence of a viscous jump at the post-catastrophe stage means that new protons join the shock being cold: they are smoothly decelerated (or accelerated in the rest frame) by radiative reaction via pairs and magnetic field.

The fluid behind the main radiation front (where it has a semirelativistic velocity) is not pressure balanced in this approximation. The thermal pressure of protons at $`x<0.15`$ swept up before the catastrophe exceeds the radiation pressure of downstream photons. This disbalance should affect the velocity distribution and a real pattern can be more complicated. One can expect a semirelativistic forward bulk motion at the left from the broad peak of the photon distribution, which will compress and push pair-photon gas ahead.

As far as the model is very simplified and works progressively worse in attempts to trace the evolution further, we terminate the simulation at $`t=0.67`$. At the present level without a careful hydrodynamic treatment we are not able to study shock deceleration stage and cannot reproduce the final time profile of the radiation. We also cannot say when and how the highly radiative post-catastrophe state terminates.

In order to represent results in terms of observed gamma-ray energy release we have to consider different photon components in the comoving system:

(a). The downstream component which takes the main fraction of the momentum and the energy of the decelerating external medium, after the viscous shock front disappears. Its energy is related to the mass of the external medium $`m`$ swept up after the catastrophe as $`E_dmc^2\mathrm{\Gamma }`$.

(b). The upstream component (which is actually isotropic in the comoving frame) results from the absorption and quasi-isotropic reemission of the downstream flux: $`E_u\eta E_d`$ where $`\eta `$ is the average opacity of the system for the downstream photon flux. After the transformation of the upstream component to the rest (observer) frame we have:

$$E_{u,o}\eta mc^2\mathrm{\Gamma }^2$$
(18)
The opacity rises during the evolution and approaches unity in the case of the high compactness. Indeed, all soft photons scatter in optically thick layer of pairs, while hard photons are absorbed due to pair production. The deceleration of the shock can be roughly described as

$$\frac{d\mathrm{\Gamma }}{dm}=\eta \frac{\mathrm{\Gamma }^2}{M}$$
(19)
Where $`M`$ is the invariant mass (comoving energy) of the shock. Note that at $`\eta <1`$ the deceleration is slower than in adiabatic case which is described in the same way as (19) at $`\eta =1`$.

It is remarkable that the upstream spectrum evolves towards a GRB-like spectral shape, however still does not reach it. The peak energy of the latest spectrum in Fig. 2 is around 70 keV. When blueshifting the spectrum to the observer frame and assuming a cosmological redshift $`z=1`$ the peak energy comes out around 1 MeV. There exist GRBs with such peak energy, however the typical peak energy is less - between 200 and 300 keV. The low energy slope in this example is $`\alpha 1.5`$ while for a typical GRB $`\alpha 1`$. The values of the peak energy and $`\alpha `$ resulting from this run are quite robust: the peak energy comes out from transition from Thomson to Klein-Nishina regime in Compton scattering on nonrelativistic or semirelativistic electrons, $`\alpha =1.5`$ corresponds to cooling spectrum of relativistic electrons (see e.g. Ghisellini, Celotti, & Lazatti 2000).

Summarizing, one should conclude that at the present level the model does not reproduce GRB spectra. Nevertheless, the phenomenon provides a perspective to address the problem due to efficient pair production and generation of optically thick layer of pairs behind the radiation front (see Fig. 3 and Fig. 4, lower panels). Once we have an optically thick pair plasma, we can expect that it will produce a thermal Comptonization peak which may describe GRB spectra (Thompson 1997, Ghisellini & Celotti 1999, Liang 1999). In this simulation such peak cannot be reproduced since cooled pairs are not Maxwellized: their energy distribution is in Compton equilibrium with the photon spectrum. This is a consequence of the simplified approach omitting various processes of potential importance. Particularly, the spectra could be modified in a proper direction by thermal bremsstrahlung and double Compton scattering (for a review see Thompson 1997).

Actually, cool pairs should be subject to a variety of collective phenomena, i.e. dissipation of plasma waves or turbulence. As the result a prominent Comptonization peak can appear. Note, that the asymptotic case of the Comptonization peak is Wien distribution which is much harder in low energy part ($`\alpha =+2`$) than GRB spectra. A partially developed Wien peak could be a promising solution of the problem.

The contradiction in the location of the spectral break $`\epsilon _p`$ can be relaxed if the initial Lorentz factor of the shock is 20 - 40 rather than 100 – 300 and the peak luminosity being emitted by a partially decelerated shock.

## 5 Exponential Boosting of a Seed High Energy Component

The first example in Section 4 presents a full nonlinear simulation of the system without any initial seed photons. Here we show that the existence of a very small amount of high energy photons at the start extends the possible range of the catastrophic evolution down to a much lower external density than it was probed in Section 4. The most suitable source of such photons is Fermi acceleration of protons with photo-meson production and possible $`pn`$ conversion suggested by Derishev et al. (2003).

Because of a high computing cost of the full nonlinear simulation in this section we explore the parameter space simulating the cascade boosting in a linear regime using linearized version of LPMC code. This means that the target LP field is not affected by the simulated high energy component. A high boosting rate ensures that the system evolution eventually reaches the catastrophe stage even at a small seed amount of high energy photons.

During the boosting, the spectrum of the high energy component evolves until it reaches some steady state shape that we do not know a priori. Therefore we have to perform a full simulation of the high energy component tracking its growth to a sufficiently long time. Then, in order to handle an exponential branching of the cascade tree we have to use statistical weights: to cut some branches compensating this by an increase of the statistical weight of remaining branches. The main technical problem in this simulation is large fluctuations in cascade histories which have a very slow convergence when we accumulate the average history. Since the seed target photon field has an upper cutoff at $`\epsilon _2`$ the high energy component was simulated down to the lower cutoff $`1/\epsilon _2`$: photons below this cutoff do not produce pairs.

Magnetic field and the seed (target) radiation field were taken according to §§2.2 and 2.3. The parameters of the seed spectrum were frozen at $`\alpha =2`$ and $`\epsilon _2=10^4`$. With this assumptions we reduce the parameter space to two dimensions: $`\mathrm{\Gamma }`$ and $`\rho `$. Because the assumptions are not obvious, it would be reasonable to accept the following strategy: to explore $`\mathrm{\Gamma },\rho `$ space with the standard set of assumptions and, at one point, to defreeze other parameters and to estimate the effect of their variation.

The average breeding time profiles for $`\mathrm{\Gamma }=50`$$`\rho =10^4`$ cm-3 corresponding to the first three lines in Table 1 are shown in Fig. 5. One can see a confident exponential rise in the number of hard photons crossing the shock upstream. The results on the breading rate for $`\mathrm{\Gamma }=50`$, $`\rho =10^4`$ with different sets of other parameters are summarized in Table 1.

|  | “Nonstandard” parameters | $`S_{14}`$ | $`R/S`$ |
| --- | --- | --- | --- |
| 1 | Standard set | 5.5 | 21.0 |
| 2 | $`H_e=0.2H_e^\mathrm{o}`$ | 31. | 4.8 |
| 3 | $`H_e=5H_e^\mathrm{o}`$ | 2.4 | 48.2 |
| 4 | $`\omega _B=10\omega _B^\mathrm{o}`$ | 37. | 3.1 |
| 5 | $`\omega _B=0.1\omega _B^\mathrm{o}`$ | 2.0 | 58.8 |
| 6 | $`\alpha =1.75`$ | 2.9 | 40.3 |
| 7 | $`\alpha =1.5`$ | 5.1 | 23.0 |
| 8 | $`\alpha =2.25`$ | 21.0 | 2.5 |
| 9 | $`10^8<\epsilon <10^3`$ | 4.7 | 25.0 |
| 10 | $`10^8<\epsilon <10^5`$ | 15.3 | 7.6 |

Table1 The boosting rate of the high energy component for $`\mathrm{\Gamma }=50`$, $`\rho =10^4`$ ($`\mathrm{}_\mathrm{o}=7.3`$) depending on other parameters. Second column gives those parameters which differ from the standard set (see text). Third column shows the boosting rate in terms of the distance of the shock propagation $`S`$ where the high energy component grows by an order of magnitude$`S`$. Fourth column gives the total growth of the high energy component (orders of magnitude) when the shock propagates the distance $`R=0.5R_d=1.1710^{16}`$cm, where the deceleration distance $`R_d`$ implies the initial isotropic kinetic energy $`10^{54}`$ erg. The standard parameters are: $`\mathrm{}_e=3.510^3`$, $`H^\mathrm{o}=14.7`$ G, $`\omega _B^\mathrm{o}=1.410^3`$, $`H_e^\mathrm{o}=10^3`$G, $`\epsilon _2=10^{04}`$, $`\alpha =2`$

First of all it should be noticed that the breeding rate at moderate parameters $`\mathrm{\Gamma }=50,\rho =10^4`$cm-3 is in a large excess over required to trigger the catastrophe in the most of trials. One can see a strong dependence on magnetic field. The effects of $`H`$ and $`H_e`$ variation are of opposite signs as it was qualitatively shown in §2.4.

Figure 6 demonstrates spectra of particles crossing the shock (energy of all particles is in the comoving system). One can see that the main fraction of energy is transferred downstream across the shock by high energy Comptonized photons rather then by advected pairs. The turnover in spectra of upstream photons qualitatively agrees with equation (14) yielding $`\gamma _{\mathrm{max}}310^6`$ for $`\rho =10^4`$ and $`\gamma _{\mathrm{max}}\mathrm{0.7\; 10}^7`$ for $`\rho =10^3`$. A larger difference in these two cases visible in Fig. 6 appears due to a higher pair production opacity for high energy photons in the first case.

Fig. 7 demonstrates the map of $`\mathrm{\Gamma }\rho `$ space when other parameters are frozen as described above. Levels of the constant compactness and the constant time scale are given for total isotropic kinetic energy $`E=10^{54}`$ erg. Both compactness and the time scale vary as $`E^{1/3}`$ in the spherical geometry and at the constant density. The broken solid line shows the catastrophe threshold which is conventionally defined as the condition that the cascade is boosted at least by 10 orders of magnitude when the shock propagates the distance $`R=0.5R_d`$. Perhaps one could use a weaker requirement: we prefer to use this conservative criterion since the real intensity of the seed high energy radiation is unknown.

The threshold line in Fig. 7 depends on the frozen parameters, i.e. it is model dependent. As it was mentioned in §2.4 the strong slowly dissipating magnetic field in the shock would substantially rise the threshold. On the other hand if the field dissipates very fast leaving a dominating large scale component, the threshold would be lower. The range of the catastrophe would be extended if an external field stronger than assumed is induced due to spherically asymmetric Compton drag as discussed in §2.2.

## 6 GeV bursts

Two runs described in Section 4 correspond to a comparatively high compactness parameter when the main energy release occurs in the soft gamma range. In Section 5 it is shown that the catastrophe is possible at a much less compactness (see Fig. 7). In such case we can expect to observe a much harder radiation which would approximately correspond to early spectra in Fig.1. In order to demonstrate this directly we performed a run for $`\mathrm{\Gamma }=50`$, $`\rho =10^3`$cm-3. The compactness is $`\mathrm{}_\mathrm{o}=1.1`$, “standard” magnetic field in the shock $`H=4.4`$ G, characteristic radius $`R=2.010^{16}`$cm, observer time scale $`T_\mathrm{o}=133`$ s. The “minimal model” does not work at such parameters as the primary synchrotron radiation is strongly self-absorbed, therefore the same standard ad hoc hypothesis of the seed radiation as in Section 5 has been used: $`\mathrm{}_s=\mathrm{}_\mathrm{o}m_e/m_p`$, $`10^8<\epsilon <10^4`$, $`\alpha =2`$. The catastrophe at such parameters develops due to the breeding of external high energy photons. It was assumed that the high energy component due to preceding breeding has gained the luminosity $`\omega _h=0.005\omega _s2.510^6\mathrm{}_\mathrm{o}`$. The corresponding amount of upstream high energy photons with $`\epsilon =10^5`$, has been injected at the start of the simulation.

The resulting sequence of instantaneous upstream photon spectra is shown in Fig. 8. Generally, at a low compactness we can expect a maximum of $`\epsilon ^2dN/d\epsilon `$ distribution in high energies and an approximate power law spectrum in a wide range with the index between $`\alpha =1.5`$ (fast cooling spectrum) and $`\alpha =2`$ (saturated pair cascade, Svensson, 1987 ). In this case the maximum energy release is at the rest frame energy $`5`$ GeV (the spectrum in Fig 8 is not final, however, since some hardest photons will be absorbed on the way to the observer).

According to Fig. 7, the catastrophe threshold at $`\mathrm{\Gamma }50`$ is $`\mathrm{}_\mathrm{o}1`$ when the main energy is still released in hard gamma-rays. This value is model dependent, the threshold can be higher for less favorable combination of parameters, e.g. for higher $`\omega _B/\omega _s`$ ratio. Nevertheless, if the environment of GRB progenitors is diverse, we can expect that GeV bursts should be as usual phenomenon as classic GRBs. On the other hand, they are more difficult for observations. Such events should be weaker in soft gamma range (soft gamma fluence is one – two orders of magnitude less than the high energy fluence). And, more importantly, they should be longer and smoother (see Figs. 7 and 9). The example of this section corresponds to characteristic time scale 133 s (the total duration can be longer). Such events are more difficult for a triggering and an off-line search than classic GRBs, unless an event is strong enough. The long duration also makes difficult a time-location cluster analysis for high energy photons in existing data.

Has something like this been observed? A candidate for such a GeV burst is GRB 970417a detected simultaneously by BATSE and Milagrito (Atkins et al. 2000) in the range above 50 GeV. The inferred fluence above 50 GeV is at least an order of magnitude greater than the fluence in the BATSE range which corresponds to the type of the spectrum shown in Fig. 8 if the maximum is at an order of magnitude higher energy. However, it was a comparatively short burst with $`T_{90}=7.9`$ s. Then following Fig. 7 one has to assume a high Lorenz factor ($`\mathrm{\Gamma }200`$), or a less kinetic energy at a lower $`\mathrm{\Gamma }`$. Note that there should be a strong selection effect in favor of short GeV bursts.

Another hint at GeV bursts provide EGRET data. The famous GRB 940217 (Hurley et al. 1994) where EGRET has detected a delayed GeV emission during an hour after BATSE trigger hardly can match the concept of a GeV burst. It rather looks like a superposition of a classic GRB and a GeV burst. In principle this is possible when the shock meets localized dense clouds while propagating in a less dense environment. Better candidates for GeV bursts can be several EGRET detected GRBs reported by Catelly, Dingus & Schneid (1997) and Dingus, Catelly & Schneid (1997). Some of them have high energy spectral indices $`\alpha 2`$ up to GeV range. Eventually the issue of GeV bursts will be clarified by GLAST mission.

## 7 Discussion

While the effect of the electromagnetic catastrophe can, in principle, be applied to different astrophysical phenomena, its application to GRBs is of the highest interest. Below we discuss the main issues where the effect can be relevant.

### 7.1 The radiative efficiency and the GRB scenario

The considered effect straightforwardly provides the ultimate ($``$ 100% radiation efficiency and, in this respect, it probably has no competing mechanisms. The phenomenon of GRBs does require a high radiation efficiency. The total isotropic soft gamma ray energetics of GRBs with measured redshift varies between $`10^{52}`$ and $`510^{54}`$ erg. A low efficiency would imply much higher values for the isotropic kinetic energy of the ejecta which probably contradicts the afterglow data.

The catastrophe is a model independent effect which, however, has a set of threshold conditions. If GRBs were associated with coalescence of neuron star binaries, the effect would never work since the expected density of the environment is too low: 0.1 – 1 cm-3. However a wealth of recent data support the collapsar scenario where we can expect a much higher density and the phenomenon may occur (see Mészáros 2002 for a review). Is it possible that the threshold conditions are still never satisfied despite a sufficient density of the environment? Formally, it could be if the external magnetic field is always purely radial. But as is was discussed in §2.2 the radial field structure can hardly survive the early stage of a GRB in any scenario. Other inhibiting factor could be a strong domination of magnetic energy density over radiation energy density. However, we do observe GRBs implying a high radiation density in the source. Whatever produces this radiation at the start should end up with the electromagnetic catastrophe.

The effect is associated mainly with external shock scenario as internal shocks have too low relative Lorentz factor. On the other hand, internal shocks can radiate only a relatively small fraction of the total kinetic energy because of kinematic constraints. A competition with a much more efficient mechanism presents a new difficulty for the internal shock scenario. For example, let us consider the typical internal shock scenario: prompt GRB emission from $`R10^{13}10^{14}`$ cm due to collisions of multiple shocks and the following afterglow due to an external shock at larger radii. An important modification to this scenario introduces the effect discovered by Beloborodov (2002): the radiation emitted by internal shocks sweeps out the external matter up to radius $`R10^{15}10^{16}`$cm leaving an empty cavity. Now let us consider the situation when the shock reaches the edge of the cavity. It meets the environment with enhanced density and optical depth (since the swept out matter is compressed and loaded by a large amount of pairs), with a non-radial magnetic field (the radial field structure can hardly survive such event), with intense seed radiation (side scattered photons). These are perfect conditions to give rise to the electromagnetic catastrophe with a proper GRB time scale. If so, the radiative energy release in such collision would outshine a burst from internal shocks by 1.5 – 2 orders of magnitude (the main fraction of total kinetic energy of the ejecta at the full radiative efficiency versus a few per cent of the kinetic energy at unknown radiative efficiency). This would look as the main event, rather than the afterglow. Even if the main energy release is in the GeV range (low compactness, see Section 6), such events can hardly escape systematic observations.

### 7.2 Issue of the time variability

The internal shock model is motivated mainly by difficulties with description of the GRB time variability arising in the external shock scenario (Fenimore 1996). Dermer (2000), on the other hand, suggests a number of possible solutions for external shocks. The effect of the electromagnetic catastrophe may extend the list of possible solutions.

1. The catastrophe is a nonlinear effect which can amplify any fluctuations in external density and seed radiation. In addition it is sensitive to the geometry of magnetic field. One dimensional treatment performed in this work loses many important features which could appear in a 3D simulation with realistic hydrodynamic treatment. The catastrophe most probably will develop locally in many spots within $`1/\mathrm{\Gamma }`$ cone at different time. The spherical symmetry of the shock should be distorted due to a strong feedback between the radiation and the fluid pattern. As a result one can observe a complex time structure even at an approximate symmetry in the initial state.

2. The post-catastrophe evolution was interpreted here as a quasi-steady state. Actually, especially in the case of decreasing density, it could be unstable and recurrent. This could be, e.g., due to a version of the Beloborodov effect: the radiation of external shock passing through the dense matter sweeps out a less dense matter ahead, radiation decays until the shock catches up the swept matter. Then the shock regenerates producing a new emission episode. In this way one can explain episodes separated by long quiescent intervals.

3. The spherical symmetry of ejecta can be completely broken prior the emission stage up to fragmentation into a bunch of separate droplets of different transversal size $`X_T<R/\mathrm{\Gamma }`$ (this implies a hydrodynamic transversal confinement of droplets which needs a separate study). This require a strong instability (most probably of Rayleigh-Taylor type) at some stage. A model representing a GRB ejecta as a shower of blobs was proposed by Heinz & Begelman (1999). This is the most attractive possibility because in addition to the complex time structure it can produce a chain reaction: the catastrophe in one droplet supply seed photons to a neighboring droplet and induce the next catastrophe. It was shown by Stern & Svensson (1996) that the chain reaction gives a statistical and qualitative description of GRBs temporal diversity.

All these suggestions require 3D numerical studies.

### 7.3 Problem of the spectral break

The observed spectral break in the sub-MeV range with a very hard low energy slope is still a serious problem. It seems that this problem is common for all GRB models assuming a high Lorentz factor and a small optical depth. Two alternative explanations of the spectral break for optically thin emitting medium were suggested: optically thin synchrotron emission in slow cooling regime and synchrotron self-absorption break. For the criticism of the former see Ghisellini, Celotti & Lazatti (2000). The latter implies too high magnetic field (in our case the self-absorption break should appear at $`\epsilon 10^610^5`$ in observer frame).

As it was discussed in Section 4 the electromagnetic catastrophe at a high compactness provides an optically thick pair loading and thus gives a perspective to explain the break with thermal Comptonization peak. However there remains the issue whether the thermal Comptonization is sufficiently fast and efficient. If it works, then we should find an anticorrelation between the sharpness of the break (implies high compactness) and the high energy emission (implies low compactness). If the break is still too persistent then we have to conclude that our consideration misses some important detail. For example, a possible role of ejecta behind the contact discontinuity was ignored while actually it could be an efficient photon reflector.

The location of the break hints at a comparatively low Lorentz factor and therefore a high density and small deceleration radius. It could be that the range of parameters studied in this work is actually far from the typical GRB “working regime”.

### 7.4 Summary

The effect of the electromagnetic catastrophe outlined here in the “zero approximation” can explain the high radiative efficiency of ultrarelativistic shocks while the explanation of the GRB time variability and spectra still require a lot of work. The study of the time variability should certainly rely on a detailed three dimensional hydrodynamics. The understanding of the spectra requires a more complete description of the shock structure and physics of particle interactions for nonrelativistic and semirelativistic pairs.

## ACKNOWLEDGMENTS

The author thanks K. Hurley, A. Illarionov, P. Ivanov, J. Poutanen for useful discussions and Y. Tikhomirova for assistance. This work was supported by the Russian Foundation for Basic Research (grant 00-02-16135). A part of this work was done during a visit to the University of Oulu with the support of the grant from the Jenny and Antti Wihuri foundation.


text:
# Wave nucleation rate in excitable systems in the low noise limit

## Abstract

Motivated by recent experiments on intracellular calcium dynamics, we study the general issue of fluctuation-induced nucleation of waves in excitable media. We utilize a stochastic Fitzhugh-Nagumo model for this study, a spatially-extended non-potential pair of equations driven by thermal (i.e. white) noise. The nucleation rate is determined by finding the most probable escape path via minimization of an action related to the deviation of the fields from their deterministic trajectories. Our results pave the way both for studies of more realistic models of calcium dynamics as well as of nucleation phenomena in other non-equilibrium pattern-forming processes.

One very important class of non-equilibrium spatially-extended systems is that of excitable media. In these, a quiescent state is linearly stable but nonlinear waves can nonetheless propagate without decaying. These waves can be generated by above-threshold local perturbations and they also can become self-sustaining in the form of rotating spiralsWinfree (1972). Examples of excitable media include many biological systems such as the cAMP waves seen in Dictyostelium amoebae aggregationLoomis (1975), electrical waves in cardiac and neural tissueWinfree (1987) and, primary for our focus here, intracellular calcium wavesBerridge et al. (2000).

Most excitable systems are sufficiently macroscopic as to render unimportant the role of thermodynamic fluctuations and allow for a description in terms of deterministic pde models. For these cases, noise effects can still be studied via the imposition of external variation in time and/or space (e.g. by varying illumination in a light-sensitive BZ reaction Sendi a-Nadal et al. (2001); Jung et al. (1998)); however, there is no need to include noise in a description of the “natural” version of these systems.

This is manifestly not the case for intracellular calcium dynamics; since the excitability here arises through the opening and closing of a small number of ion channels (allowing/preventing calcium efflux from calcium storesShuai and Jung (2002)), fluctuations are inherently non-negligible. Indeed, experiments show direct evidence of noise effects in the form of abortive waves and spontaneous wave nucleationLlano et al. (2000); Marchant et al. (1999).

In this paper, we study the process of spontaneous wave nucleation for a 1d generic excitable system modeled by the Fitzhugh-Nagumo equations

$`_tu`$ $`=`$ $`(3uu^3v)/ϵ+^2u+\eta _u`$ (1)
$`_tv`$ $`=`$ $`u+(1+\gamma )+\eta _v`$ (2)
Here $`\eta _v`$ and $`\eta _u`$ are small independent white noise terms modeling fluctuation effects with covariance equal to:

$$<\eta _i(x,t)\eta _j(x^{},t^{})>=\beta \delta _{ij}\delta (xx^{})\delta (tt^{})$$
(3)
At $`\beta =0`$ and for positive values of $`\gamma `$, this system is excitable with a single stable equilibrium point $`u_0=(1+\gamma ),v_0=3u_0u_0^3`$. As already mentioned, a wave of excitation can propagate through the system (see fig.1); a counter-propagating pair of such waves will be generated if a local perturbation above a threshold value is applied. For negative values of $`\gamma `$, the system becomes oscillatory. This model is not meant to be a realistic approximation for any specific physical or biological process; instead we use the model to understand the generic features of wave nucleation due to noise.

Specifically, at finite $`\beta `$ the noise will allow the birth of pairs of counter-propagating wave at a rate that will depend on the amplitude of the noise. That rate can be determined in the case of relatively high noise using direct numerical simulations. However, in the case of low noise, such a method becomes computationally prohibitive. Here we present a computation of the transition rate using a most probable escape path (MPEP) approach that is based on solving the Fokker-Planck equation Kramers (1940); Freidlin and Wentzell (1984). This method has been successfully applied to dynamical systemsGraham and Tél (1985); Maier and Stein (1993) and to a few cases of spatially extended systems, namely the transition from creep to fractureMarder (1995) and magnetic domain reversalWeinan et al. . To derive the equation for the MPEP, we use the fact that the solution of the Fokker-Planck equation with given initial and final field configurations for time interval $`[0,T]`$ can be written in terms of the path integral

$$P(T)=𝒟(u,v)\mathrm{exp}\left(\frac{1}{\beta }𝑑t𝑑xS(x,t)\right)$$
(4)
with the “action density” $`S`$ given by the sum of the squared deviations of the time derivatives of the fields from their deterministic values: $`S\delta _u^2+\delta _v^2`$

$`\delta _u`$ $`=`$ $`_tu(3uu^3v)/ϵ\mathrm{\Delta }u`$ (5)
$`\delta _v`$ $`=`$ $`_tvu(1+\gamma )`$ (6)
The functional integral is taken over all paths that begin at $`t=0`$ in equilibrium and end with a given final counter-propagating wave state $`(u_f,v_f)`$ at $`t=T`$. Since we take $`\beta `$ close to 0, the r.h.s. of eq. (4) is dominated by the path (called most probable escape path or MPEP) that maximizes the integrand over all paths; thus, the transition rate is found to be proportional to:

$$\mathrm{exp}(E/\beta )$$
(7)
where E is the minimum over all paths of $`𝑑x𝑑tS`$

Therefore, in order to compute the transition rate between the rest state and a pair of counter propagating waves, one has only to compute the minimum of the quantity in eq. (4). This minimum can be expressed using a variational principle as the solution of a PDEMarder (1996); however, using that approach to finding the actual MPEP between two different states involves the use of a shooting method with numerous parameters, which turns out to be numerically quite difficult. Instead, we used an alternative method based on discretizing the above path integral on a space-time grid and directly using a quasi-NewtonNocedal (1980) method to find the minimum. One difficulty with this approach is that the Fitzhugh Nagumo model has a slow recovery time compared to the timescale associated with a pulse (width of pulse/speed). This then necessitates having a very large spatial domain, if one attempts to completely encompass the region over which the nucleated wave configuration differs from the quiescent fixed point. To get around this difficulty, we use a grid moving with the pulse. That is, the grid moves in the same direction as the pulse in order to keep the wave front at a fixed distance from the boundary. Thus, the boundaries of the domain used to compute the MPEP path are no longer $`[0,L],[0,T]`$, but $`[x_{min}(t)=0+\text{max}(x_f(t)x_{off},0),x_{max}(t)=L+\text{max}(x_f(t)x_{off},0)],[0,T]`$, where $`x_f(t)`$ is the position of the wave front defined as the point where $`u`$ goes above 0 and $`x_{off}`$ is an arbitrary value lower than $`L`$ and significantly bigger than the wave front width. The following boundary conditions are applied:

$$\{\begin{array}{c}\text{at }x=x_{max}(t)\text{}u=u_0,v=v_0\hfill \\ \text{at }x=x_{min}(t)\text{}\{\begin{array}{ccc}_xu(0)=0\hfill & \text{ if }& x_{min}(t)=0\hfill \\ \delta _u(0)=0\hfill & \text{ if }& x_{min}(t)0\hfill \end{array}\hfill \end{array}$$
(8)
The condition $`\delta _u=0`$ implies that the recovery past $`x_{min}`$ is purely deterministic and hence does not contribute to $`E`$. A check on our procedure is afforded by the fact that as long as the distance the wave had traveled (in the final state) is significantly bigger than the region over which the wave initiates, the results obtained are independent both of the distance the wave had traveled and of the width of the space window used ($`L`$). We fix $`T`$ to be large enough that the deviation from deterministic dynamics for very small time is negligible. Once $`T`$ is fixed and the specific final wave state chosen, there is no time-translation invariance in the MPEP.

We now describe the results obtained using this methodology. In the small $`ϵ`$ limit ($`ϵ<0.1`$), the shape of the MPEP is quantitatively independent of $`ϵ`$ and $`\gamma `$ and even for high values of $`ϵ`$, there is no significant qualitative difference. In fig. 2, we present such a typical escape path. As shown on the contour plots, wave nucleation is found to be a very localized event. Essentially, the noise acts to create a local dip in value of the $`v`$ field which is then followed by a large positive excursion for the $`u`$ field as it goes into the excited phase. To describe this mechanism more fully, we present in fig. 3 the phase-plane trajectory at the center of nucleation as well as several snapshots of the spatial form of the fields during the nucleation process. One can see then that the escape path consists of the center of nucleation being driven by noise below the minimum of the $`f`$ nullcline; this then quickly drives the $`u`$ field positive and leads after, $`v`$-field driven relaxation, to the pair of counter-propagating pulses. There is a significant fluctuation contribution to the nucleation event in a small region around the nucleation point(see fig. 4). Results using different values of $`\gamma `$ and $`ϵ`$ show that the width of that small region is proportional to the width of a front, that is $`\sqrt{ϵ}`$ (see fig. 5). Note that the other simple possibility, that of nucleating an excited region for the $`u`$ field at a fixed value of $`v`$Pumir (1995), is not observed.

In accord with the shape of the trajectory, our calculations show that the main contribution to $`E`$ during the MPEP come from the $`\delta _v^2`$, at least for small values of $`ϵ`$. Thus for $`\gamma =0.2`$ the value of the ratio $`\delta _v^2/\delta _u^2`$ never went below 30 for $`ϵ=0.1`$, and went up to 1000 for the minimal value of $`ϵ`$ used ($`ϵ=0.001`$). For higher values of $`ϵ`$, the ratio was significantly lower (3 for $`ϵ=0.8`$, close to the limit of propagation for this value of $`\gamma `$ ). Furthermore for small $`ϵ`$, $`E`$ scales like $`\sqrt{ϵ}`$ (see fig. 5). For higher values of $`ϵ`$, this simple scaling is no longer valid. The fact that for small values of $`ϵ`$, $`E`$ scales like $`\sqrt{ϵ}`$ and that a wave is therefore much easier to nucleate can be explained by a simple argument.

As already mentioned, the phase-plane trajectory of the MPEP is mainly independent of the value of $`ϵ`$. The only dependence is then due to the spatial scale appearing in the integral which is $`\sqrt{ϵ}`$, giving rise to the observed result

We now describe results obtained when varying $`\gamma `$ with $`ϵ`$ held constant. This means that we consider the nucleation rates for differing excitabilities but with the same time-scale ratio between the $`u`$ and $`v`$ field dynamics. For a wide range of values of $`\gamma `$, we find that $`E`$ scaled like $`\gamma ^2`$ (see fig. 5). The more excitable the system, the more likely it is to nucleate a wave. This scaling is that same as one would obtain analytically in the much simpler zero-dimensional version of this problem. Here the analog of wave nucleation is the noise-induced creation of a large transient excursion away from the fixed point. If $`ϵ`$ is small, the escape path will follow the $`f`$ nullcline down to its minimum and then follow the noiseless dynamics to reach the other stable branch. In such a situation, one can compute the transition rate analytically in a piece-wise linear version of the FH model,

$`\dot{u}`$ $`=`$ $`\{\begin{array}{cc}\frac{1}{ϵ}(v+u)\hfill & \text{if }u<0,\hfill \\ \frac{1}{ϵ}(1vu)\hfill & \text{if }u0,\hfill \end{array}`$ (11)
$`\dot{v}`$ $`=`$ $`u\gamma .`$ (12)
After elimination of $`u`$, it can easily be shown that the variational equation for the MPEP has the form

$$\frac{d^2v}{dt^2}=v+\gamma $$
(13)
with $`v(0)=v_0=\gamma `$ and $`v(T)=0`$. A straightforward calculation shows that the action is equal to

$$\gamma ^2\frac{1}{2(1\mathrm{exp}(2T))^2}$$
(14)
and that its minimum, reached for $`T=\mathrm{}`$, is equal to $`\gamma ^2/2`$. This result shows that the $`\gamma ^2`$ scaling of $`E`$ found in our computations can be interpreted as being due to $`\gamma `$ equaling the “distance” of the stable equilibrium point from the border of its basin of attraction. Putting it all together, our data yield a log transition rate proportional to $`\frac{\gamma ^2\sqrt{ϵ}}{\beta }`$. This result could be tested experimentally, perhaps by adding illumination noise to the light-sensitive BZ reaction.

It is worth mentioning that there are other potential applications of the MPEP approach to nucleation in spatially-extended non-equilibrium systems. One example concerns the thermal generations of localized patches of traveling rolls in electro-convection, as studied recentlyBisang and Ahlers (1998). Also, the method used here is not limited to white noise. A simple generalization of the derivation allows for the incorporation of multiplicative noise via dividing the $`\delta _u^2`$ and $`\delta _v^2`$ terms in the action density by the corresponding (possibly field dependent) variances of the noises added to the $`u`$ and $`v`$ equations respectively. This will clearly be necessary for the study of realistic calcium models. Finally, there is a similarity between the MPEP method and what must be done to consider quantum tunneling in spatially extended systemsFreire et al. (1997), where one also must find the entire space-time path tunneling trajectory in order to find a leading estimate of the rate.

It is a pleasure to acknowledge useful discussions with D. Kessler. This research is supported by the National Science Foundation through Grant No. DMR-0101793.


text:
# Acoustic Emission from Paper Fracture

## Abstract

We report tensile failure experiments on paper sheets. The acoustic emission energy and the waiting times between acoustic events follow power-law distributions. This remains true while the strain rate is varied by more than two orders of magnitude. The energy statistics has the exponent $`\beta 1.25\pm 0.10`$ and the waiting times the exponent $`\tau 1.0\pm 0.1`$, in particular for the energy roughly independent of the strain rate. These results do not compare well with fracture models, for (brittle) disordered media, which as such exhibit criticality. One reason may be residual stresses, neglected in most theories.

A simple, common-day statistical physics experiment is to tear a piece of paper into two parts. During this noise is heard sethrev , evidently originating from damage to the paper or the propagation of a crack. It is easy to make the fracture surfaces “rough” as told by the naked eye without any more sophisticated analysis, and the experience is that it is difficult to do the tearing so that the cut is at all close to straight. Meanwhile, one can confirm the existence of disorder in the “sample” by just looking through the sheet against any reasonable light source. The paper looks cloudy, which comes to a large degree from the local fluctuations of areal mass density. This tells that elastic and fracture properties vary locally. Such a test poses the related questions: how does a crack develop in an inhomogeneous material, and what kind of properties does the acoustic emission have?

In fracture there is evidence of scaling properties familiar from statistical physics. The roughness of cracks, as measured by e.g. the root mean square height fluctuations vs. length scale, can often be described by self-affine fractal scaling Mandel ; Bourev ; DNBC . Another such quantity is the distribution of strength, and its average versus sample size. In particular for (quasi)-brittle materials, the average strength can follow scaling laws, if the sample size is varied, which derive from the presence of disorder in the material Dux ; Curtin ; Korte .

Investigating such questions combines materials science and statistical mechanics. Two fundamental problems are, what is the role of disorder in fracture? And, how do cracks develop and interact in real materials? These join forces in the failure of a notched three-dimensional sample, visualized by the advancement of an one-dimensional crack front or line in the presence of disorder HaH95 ; Fisher . It has an equation of motion, in an environment with varying properties as local strength and elastic modulus. The trail left by the line-like front is the two-dimensional crack surface, which can become rough (self-affine) if the propagating line-like object develops critical fluctuations. However, it is difficult to propose such a model that would match with the experimental data. For one, simple models tell that the disorder and microcracks interact in crack dynamics, so that sometimes it is not correct to consider an isolated object like the front or a crack tip. From the materials research viewpoint the dynamics is interesting since it deals with the engineering quantities of strength and toughness, or the resistance to intrinsic flaws. For instance in fiber and two-phase ceramic composites one can control these two by varying the mixing and constituents.

We study fracture in ordinary paper, as a two-dimensional disordered material, via acoustic emission (AE) analysis. The release of acoustic energy is related in paper to irreversible deformation, microcracks, and, perhaps, to plasticity. A large-scale analogy is earthquakes. Their energies are described by a power-law probability distribution, the Gutenberg-Richter-law stars  with an energy exponent $`\beta `$. Cracks in paper may be self-affine Kertesz ; Rosti , with a roughness exponent close to that in a scalar two-dimensional fracture model, random fuse networks (RFN’s) Maloy ; letteri . Both these are near 2/3, the one for surfaces of minimal energy (HaH95 ; letteri ; Rosti ). Paper provides an example of crack growth in the presence of disorder. The probability distributions of the released energy and the temporal statistics, describing dynamical aspects, can be compared to models for two-dimensional failure, and to other work on AE.

Here the strain rate $`\dot{ϵ}`$ is varied in the tensile test, to study possibly time-dependent effects in the fracture process. Our first main conclusion is that in strain-controlled tests power-laws are found in the statistics, with either no or at most a weak dependence on the strain rate. The power-laws do not follow the predictions of theoretical models outlined below, as fuse networks or mean-field ones. It is also apparent that an eventual localization of the crack Ciliberto ; penn  does not change the microscopic properties of the crackling noise, as long as the crack propagation is stable. This implies that paper failure is not a “phase transition” with a diverging correlation length, but leaves open the origins of the observed scalings.

To describe the interactions in an elastic medium with randomness incorporated one starts from mean-field stress-sharing. In fiber bundle models (FBM) the applied external force $`F(t)`$ is shared by all the $`N(t)`$ fibers, and as the (random) failure threshold of the currently weakest one is reached, the stress $`\sigma (t)=F(t)/N(t)`$ is instantaneously distributed among the remaining fibers. As a sign of criticality there is a divergent scale, the size of a typical avalanche which is made up of all the fibers that break down due to a single, slow-time scale stress-increment. The total avalanche size ($`N`$) distribution follows $`P(N)N^{5/2}`$ hansen ; sornette . Thus global load sharing fracture resembles a second-order phase transition.

Local stress enhancements can be added to FBM’s e.g. by considering fiber chains in which the stress from microcracks, of adjacent missing fibers, is transfered to the ones neighboring the crack. A more catastrophic growth results (with an exponent much larger than 5/2), resembling a first-order phase transition since the elastic modulus has a finite drop at $`\sigma _c`$/$`ϵ_c`$ zhang ; hansen ; phoenix . A finite-dimensional, but scalar, approximation is given by random fuse networks review . RFN’s reproduce with strong enough disorder many of the features of mean-field FBM’s, in spite of including the stress-enhancements and shielding from other microcracks Zapperi ; Z2 . The dynamics of crack growth is not understood well even in this simple model, including why the cracks seem to become self-affine, with a roughness exponent close to 2/3 letteri . The models’ event statistics can be compared with experiments by considering the energies. For the RFN, in an event the lost energy is $`E_i\mathrm{\Delta }G_iϵ_i^2`$, where $`\mathrm{\Delta }G_iN_{broken}`$ is the change in the “elastic modulus” due to the AE event, at $`ϵ_i`$. Simulations on RFN’s for strong disorder that reproduces the FBM exponent (5/2) for the avalanche sizes yield the exponent $`\beta _{RFN}1.8`$ and the FBM value is $`\beta _{FBM}2`$ Rostia .

Experiments on, mostly, three-dimensional systems Ciliberto ; Petri2 ; Lockner ; geofyysikoita ; penn  have yielded power-laws for the AE energy release; the typical exponent $`\beta `$ is 1 $`\mathrm{}`$ 1.5. One idea is that approaching final failure resembles a phase transition: in the Lyon group’s stress-controlled experiments indications were seen of a critical energy release rate Ciliberto . This means that the AE energy would diverge like a power-law in the proximity of a critical point, the sample strength $`\sigma _c`$. Another problem is the nature of the disorder in the material: whether quenched disorder is able to give such power-law statistics for the energy Petri2 ; Zapperi .

Normal newsprint paper samples were tested in the machine direction md/cd  on a mode I laboratory testing machine of type MTS 400/M. Due to the lack of constraints the samples could have out-of-plane deformations, and none of the three fracture modes (I, II, III) is excluded on the microscopic level. The deformation rates $`\dot{ϵ}`$ varied between 0.1 %/min and 100 %/min. The AE system consisted of a piezoelectric transducer, a rectifying amplifier and continuous data-acquisition. The time-resolution of the measurements was 10 $`\mu `$s and the data-acquisition was free of deadtime. The stress was measured simultaneously to AE with a time resolution of 0.01 s. We made 20 identical repetitions for statistics. The 100 by 100 mm samples had initial notches (size 15 mm) to achieve stable crack growth. Typical sheet thickness is about 100 $`\mu `$m. The fracture statistics are not affected by such a notch Rostia . For these strain rates the sound velocity timescale is much faster than that implied by $`1/\dot{ϵ}`$. Each individual test contributes, at most, 1000 - 2000 events, so we can only look at integrated probability distributions and not in detail at local averages of e.g the acoustic event size vs. $`ϵ`$.

Fig. (1) shows an example of two tests under strain-control. Stress-strain curves have typically three parts: pre-failure (almost) linearly elastic one, the regime close to the maximum stress where the the final crack starts to propagate or is formed, and a tail that arises due to the cohesive properties of paper which allow for stable crack propagation. The faster the strain rate the less is the role of the tail. For the smallest strain rate most of the AE originates from tail (more than 90 %) while for the highest the situation is the opposite. Quantities of interest are the statistical properties before the maximum, after it, and the integrated totals, in particular the energy distributions. The time series of events allows to make qualitative observations of the correlations between subsequent events and to draw conclusions about the event properties as such. For this rather brittle paper grade, and the strain rates used, the elastic modulus is independent of $`\dot{ϵ}`$, and by AE we are able to detect a constant fraction of the elastic energy Rostia .

Fig. (2) shows the scaling of the energy for a fixed strain rate. The behavior is power-law -like, with several orders of magnitude of scaling. The same exponent fits all three different cases: the pre and post maximum stress cases, and the sum distribution. If the strain rate is varied the same conclusion holds and the exponent only fluctuates at random Rostia . We have $`\beta =1.25\pm 0.10`$, in disagreement with the fiber bundle ones and with that from the fuse networks, though $`\beta _{RFN}`$ is closer than $`\beta _{FBM}`$. The practical implication is, that the material can withstand more damage than expected since $`\beta <\beta _{models}`$. In paper, the post-maximum events correlate with the advancement of the final crack. The remaining ligament length/width contracts, thus the elastic modulus drops should (assuming a constant stress state) here be qualitatively related to AE event energies (see also balankin ).

In an experiment with a varying strain rate $`\dot{ϵ}`$ both the event durations and the waiting time between any two events may (consider Fig. (1)) depend on $`\dot{ϵ}`$. If the crack dynamics becomes “fast” then events take place on timescales set by the sound velocity. This establishes a timescale $`t_s=\mathrm{\Delta }x/v_s`$ where $`\mathrm{\Delta }x`$ is the spatial separation of two events and the sound velocity $`v_s`$ 2 $`\times `$ 103 m/s. In a sample of linear size 0.1 m this results in a maximum $`t_s10^4s`$. In the failure of brittle carbon foams the eventual critical crack growth is dictated by such fast events penn .

In quasi-static fracture models the dynamics of cracks is assumed such that the stress field is equilibrated infinitely fast during microfracture events, between further adiabatic increases of strain or stress. Thus the only time-dependence of any temporal statistics is in the average time interval, between AE events, proportional to $`1/\dot{ϵ}`$. For the FBM model there are no correlations between the waiting times and event size or durations, except for the trend that the average waiting time decreases as $`\sigma _c`$ is approached. It is ’critical’ and thus the energy release rate follows a power-law close to $`\sigma _c`$. Waiting time results do not generally exist for more complicated models herrmann . One can compare to the experimental signatures of the intervals between AE events and the integrated energy release rate integral $`𝑑ϵE_{event}`$. As noted above there is some evidence - from stress-controlled experiments - for a critical behavior for this quantity Ciliberto .

Fig. (3) demonstrates waiting time $`\tau `$ distributions for different strain rates. There is a clear power-law, whose exponent ($`\tau `$) remains roughly the same for all the strain rates cut-off . Importantly, this is true for the post-events regardless of the origin (before/after $`\sigma _c`$) of the majority of the AE energy. For the pre-events there might be some evidence of the exponent increasing with $`\dot{ϵ}`$. In the time-series of events, those with long durations are separated, on the average by shorter intervals from the neighboring events before/after. Fig. (4) depicts the waiting times prior to an event with two different data analysis methods for distinguishing between possibly overlapping events Rostia . The interval separations are similar for both the post/pre-phases, implying that the microscopic failure dynamics does not differ, and that $`\sigma _c`$ or $`ϵ_c`$ can not be inferred from the event characteristics. It may be so that events which are relatively long are precluded by longer waiting times. The durations of the events $`\delta t`$ and the sizes are roughly power-law related as $`E(\delta t)^3`$.

The integral $`𝑑ϵE_{event}`$ demonstrates a rapid exponential growth above a typical strain of $`ϵ0.5`$ %. This originates from increasing event sizes, not from an increasing density vs. strain. For $`ϵ>0.6`$ % the samples start to fail. In the regime where the exponential growth takes place the samples develop plastic, irreversible strain $`ϵ_{pl}`$, with a roughly exponential dependence on strain Rostia . This does not imply that the AE measures plastic deformation work, mostly, since the rate of increase of $`ϵ_{pl}`$ is much less than that of the energy integral. The exponential strain-dependence of the AE implies a typical lengthscale (as should be true for the development of plasticity), which in turn should be related to crack localization.

Concluding, failure of ordinary paper shows several features associated with critical phenomena. Our take on the experiment is that it show that i) there is no clear sign of a “critical point”, or a phase transition, in spite of the fact that the material is close to linearly elastic. In particular, while the event intervals and the event energies follow power-law -like statistics, not all the quantities do so. Also, ii) the temporal behavior hints of complicated time-dependent phenomena not directly related to the fast relaxation of stress. A possible candidate is the viscoelastic nature of the wood fibers in paper, but note that the macroscopic stress-strain behavior remains almost linearly elastic while AE events already occur. We suspect the gradual release of internal stresses plays a role, perhaps due to frictional pull-out of fibers from the network. Finally, iii) the power-laws as obtained are off those predicted by simple fracture models. It remains to be seen whether these models can be tailored closer to such tensile experiments visco . One suggestion is that the dynamics of energy release during the events follows a different course from the model rules, in particular finite-rate dynamics allows for stress overshoots fiprl ; recent . The relation of the acoustic emission to why cracks get rough may be indirect. The latter could relate to the development of pre-failure plastic deformation.

We acknowledge support by the Academy of Finland’s Center of Excellence program, and discussions with K. Niskanen, J. Weiss, and S. Zapperi.