text:
# Interplay between the repulsive and attractive interaction and the spacial dimensionality of an excess electron in a simple fluid

## I Introduction

The behavior of an excess eletron in a wide variety of fluids has been an interesting topic for many years ref1 ; ref2 ; ref3 ; ref4 ; ref5 ; ref6 ; ref6 ; ref8 ; ref9 ; ref10 In the gas phase or the dilute liquid phase, the electron behaves almost like a free particle. As the solvent density increases, the electron exihibits different properties depending upon the nature of the solvent and electron-solvent interactions. At liquid density, the electron may become self-trapped in a cavity of solvent particles or remain quasifree depending upon the nature of electron-solvent interaction. Observed properties such as the electron mobilityref11  and the absorption spectra ref12  probe the nature of the electronic states in the fluid and phenomena of localizationref13 .

There are a broad range of problems in condensed matter physics that are intimately related to the problem of excess electrons in deformable medium. These include charge transfer kinetics in biological reactions, metal-insulator transitions in fluids, polarons, phonon-assisted hopping of charge carriers in semiconductors and insulators, quantum-tunneling, etc. While the excess electron problem belongs to the general problem of electrons in disordered materials, the liquid environment is in many ways different from the solid medium. In liquids, the constituent particles can diffuse, and local environment around the solute electron can be substantially different from that in solid.

When an electron is solvated in a polar liquid such as water or ammonia, the strong anisotropic electron-solvent interaction causes significant local modification of the equilibrium fluid structure ref4 ; ref10 ; ref14. The electron becomes localized in a small cavity because molecules in a solvation shell orient to create a potential minimum. Even simple fluids are found to exhibit electron mobilities that change by many orders of magnitude as the density of the fluid is altered slightly. In super critical helium, for example, the electron mobility drops by over 4 orders of magnitude as fluid density increased by a factor of 2 in low density regimeref2 . The reason for this behavior is strong repulsion between electron and solvent atom. This causes the the depletion of the solvent atoms from the region of the electron and forms a highly localized state of the electron.

In many other nonpolar fluids such as Ar, CH4, etc. the electron always remains in a state of high mobilityref15  comparable to many semiconducting materials. An interesting density dependence of the mobility has been observed in them. It shows a minimum near critical fluid density and a maximum at liquid density.

Electronic states in reduced dimension are of considerable interest. For example, for a system less than two spatial dimension, electrons are localized with an infinitesimal amount of disorder ref16 . The interest in problems of electron or phonon propagation in a one dimensional random potential stems from the discovery and extensive experimental study of a certain class of organic or metallo-organic materials ref17 . These matarials exhibit strongly anisotropic, quasi-one-dimensional behavior attributed to the fact that they consist of long chains, weakly interacting with each other. In many of these, the presense of a random potential has been proposed in order to explain their behavior. Electronic surface states play an important role in a wide variety of physical problems. For example, surface electrons on liquid helium has shown many interesting properties and led to important theoretical advances such as the spectrum of bound electronic states, electron transport on the He surface, effects due to deformation of the He surface, and the possibility that the electrons may form a two-dimensional crystalref18  in the field of low-dimensionality physics.

The theory for the excess electrons in fluids developed by Chandler, Singh, and Richardsonref1  (CSR) is based on the path integral formulation of quantum theory which maps the behavior of the electron on to that of a classical isomorphic polymerref19 . The solvent-induced potential surface for the self-interaction of the isomorphic polymer is evaluated using an integral equation (e.g., reference interaction site model). With known potential surface, the polymer statistics is solved using variational approachref20  that allows the determination of electronic properties and the structure of the liquid near the electron. The input of the theory is the pure solvent structure factor and the electron-solvent particle interaction potential. The CSR theory in its formulation is applicable to an adiabatic solvent (i.e., solvent particles are treated classically), but has been extended to treat the effect of the quantum mechanical charge density fluctuation in the solvent particles ref21 . The calculated electron-absorption line shape and mobility are in good agreement with the simulations22 and experiments ref6 . The predictions of CSR theory were verified by computer simulations ref23 ; ref24 .

For a one-dimensional system we ref25  have shown recently that the repulsive and attractive parts of electron-solvent interaction potential lead separately to localization of electron, respectively, by creating cavity or forming a cluster of the fluids atoms around it. In two dimensional system we ref26  have shown that dispersion of the wavepacket associated with the solvated electron is very sensitive to the interaction between the electron and the fluid atoms, and exhibits complicated behavior in its density dependence. CSR theory has been extended to calculate the effective mass ref27  as well as the density matrix of the excess electron in fluid ref28 .

The CSR theory involves three or more characteristic lengths, depending upon the nature of electron-solvent interaction. These lengths are the thermal wavelength of excess electron $`\lambda _e=(\frac{\beta \mathrm{}^2}{m})^{\frac{1}{2}}`$ (where $`\beta `$ is the inverse of temperature in the unit of the Boltzmann constant $`k_B`$ , $`m`$ is the mass of a bare electron, and $`\mathrm{}`$ is the Planck’s constant divided by 2$`\pi `$), characteristic length associated with the electron-solvent pseudopotential, and a length associated with the mean volume occupied by each solvent atom, which is related to $`\rho ^{^{\frac{1}{D}}}`$, where $`\rho ^{}=\rho _s\sigma ^D,\rho _s`$ being the number density of the solvent and $`\sigma `$ is the diameter of the solvent atom, and $`D`$ is the spatial dimensionality of the system. The behavior of the excess electron is expected to depend sensitively on these lengths. Laria and Chandler ref29  have attempted to explain the contrasting behaviors of the electron in super critical helium and xenon on the basis of different ranges of the electron-solvent repulsive interactions.

In the present work we examine in detail the role played by different lengths and the spatial dimensionality of the system (in which we have considered the same solvent-solvent and electron-solvent model potential) to study the self-trapping behavior of the electron. The organization of the rest of the paper is as follows. In Sec. II we briefly review the CSR theory. In Sec. III we have presented the results and their discussions. Finally, Section IV presents concluding remarks. Appendix A provide some mathematical material for $`D`$\- dimensional integration

## II Theory

The system we consider is a single electron dissolved in a single component classical solvent. In the CSR theory, an excess electron is mapped, using a discretized version of the path integral formulation of quantum mechanics, onto a polymer of P interaction sites or beads. Under this isomorphism ref19 , the electron can be viewed as a classical ring polymer.

The total potential energy can be written as

$`U=U_{es}(𝐫,\{𝐑_i\})+U_{ss}(\{𝐑_i\})`$ (1)
with

$`U_{ss}(\{𝐑_i\})={\displaystyle \underset{i>j=1}{\overset{N}{}}}u_{ss}(|𝐑_i𝐑_j|),`$ (2)
and

$`U_{es}(𝐫,\{𝐑_i\})={\displaystyle \underset{i=1}{\overset{N}{}}}u_{es}(|𝐫𝐑_i|)`$ (3)
Here, r denotes the position of the excess electron, $`𝐑_i`$ is the collection of the coordinates for a solvent atom, and N is the number of solvent atoms and $`u_{es}(r)`$ and $`u_{ss}(r)`$ are, respectively, electron-solvent atom and solvent atom-solvent atom interaction potentian. We consider a $`D`$-dimensional system of spheres of diameter $`\sigma `$ in which the pair interaction between the solvent atoms is taken to be

$`u_{ss}(|𝐫|)=\{\begin{array}{cc}\mathrm{}\hfill & \text{for }|𝐫|\sigma \text{;}\hfill \\ 0\hfill & \text{for }|𝐫|>\sigma \text{.}\hfill \end{array}`$ (4)
where $`|𝐫|`$ is the $`D`$-dimensional distance (for notational convenience the $`D`$ dependence will not always be explicitly indicated). The electron-solvent atom interaction in a real system consists of a strong repulsion at short distance due to orthogonality requirements between wavefunctions of core electrons in the solvent particle and that of the excess electron and attraction at large distances due to dispersion interaction. However, in a system of neutral atoms the electronic states are determined primarily by the short range repulsive interaction or excluded-volume effect. The attractive interaction becomes important only at low densities. The interaction between the electron and solvent atom is taken to be

$`u_{es}(|𝐫|)=\{\begin{array}{cc}\mathrm{}\hfill & \text{for }rd\text{;}\hfill \\ ϵ\frac{\mathrm{exp}(\alpha r)}{\alpha r}\hfill & \text{for }r>d\text{.}\hfill \end{array}`$ (5)
Here, $`d`$ is the distance of closest approach between electron and solvent atom.

In CSR theoryref1  the partition function Z for an electron in a bath of classical particles is written as the functional integral

$`Z={\displaystyle D𝐫(u)d\{𝐑_i\}\mathrm{exp}\left[\frac{1}{\mathrm{}}_0^\beta \mathrm{}𝑑u\left\{\frac{1}{2}m|\dot{𝐫}(u)|^2+U_{es}(𝐫,\{𝐑_i\})\right\}\beta U_{ss}(𝐑_i)\right]}`$ (6)
where r(u) is the electronic path in imaginary time which is periodic in time interval $`0u\beta \mathrm{}`$, i.e., r(0)=r($`\beta \mathrm{}`$). To concentrate our attention on the electron degrees of freedom, the partition function given by Eq.(6) can be written as

$`Z=Z_s{\displaystyle }D𝐫(u)\mathrm{exp}\{\beta S_{}[𝐫(u)\beta \mathrm{\Delta }\mu [𝐫(u)]\}`$ (7)
where $`Z_s`$ denotes the partition function of the solvent, $`\mathrm{\Delta }\mu [𝐫(u)]`$ is the excess chemical potential for the fixed electronic path,

$`\text{with},\beta S_{}[𝐫(u)]={\displaystyle \frac{1}{\mathrm{}}}{\displaystyle _0^\beta \mathrm{}}𝑑u{\displaystyle \frac{1}{2}}m|\dot{𝐫}(u)|^2`$ (8)
and $`\mathrm{exp}\{\beta \mathrm{\Delta }\mu [𝐫(u)]\}`$ is called the influence functional which represents the solvent effects on the electon. In the continuum limit ref1 ,

$`\mathrm{\Delta }\mu [𝐫(u)]=\rho _s\widehat{c}_{es}(0)+{\displaystyle \frac{1}{2}}(\beta \mathrm{})^2{\displaystyle _0^\beta \mathrm{}}𝑑u{\displaystyle _0^\beta \mathrm{}}𝑑u^{}v(|𝐫(u)𝐫(u^{})|),`$ (9)
where $`\widehat{c}_{es}(0)`$ is the k = 0 spatial Fourier transform of $`c_{es}(r)`$.

$`v(|𝐫(u)𝐫(u^{})|={\displaystyle }d𝐫^{}{\displaystyle }d𝐫^{\prime \prime }c_{es}(𝐫^{},u)\chi _{ss}(|𝐫^{}𝐫^{\prime \prime }|)c_{es}(|𝐫^{\prime \prime }𝐫|,u^{})`$ (10)
Here r and $`u`$ appear as independent coordinates, r is the distance between two sites and u measures the length along the contour of the polymer, and

$`\chi _{ss}(|𝐫𝐫^{}|)=\delta \rho _s(𝐫)\delta \rho _s(𝐫^{})`$ (11)
is the density-density correlation function of the unpurturbed bath. In Eqs. (9) and (10), $`c_{es}`$ is the direct correlation function. Its value is determined from the equation ref1

$`\rho _sh(𝐫)={\displaystyle 𝑑𝐫^{}𝑑𝐫^{\prime \prime }\omega (|𝐫𝐫^{}|)c_{es}(|𝐫^{}𝐫^{\prime \prime }|)\chi _{ss}(𝐫^{\prime \prime })}`$ (12)
where

$`\omega (|𝐫𝐫^{}|)=(\beta \mathrm{})^1{\displaystyle d(uu^{})\omega (|𝐫𝐫^{}|;uu^{})}.`$ (13)
is the intrapolymer correlation function.

Eq. (12) is solved for $`c_{es}`$ and $`h`$ using suitable closure relation ref25 . Since all sites of a ring polymer on the average are equivalent, the site dependence disappears from Eq.(12) and only the zero-frequency component $`\omega (|𝐫|)`$ of the equilibrium response function is required in Eq.(12).

To complete the evaluation of excess chemical potential, the electronic path integral still has to be performed. Following Feynmam ref20  and Chandler et. al., ref1  the excess chemical potential for the fixed electronic path is mimicked by a Gaussian functional,

$`\beta \mathrm{\Delta }\mu _{\mathrm{ref}}[𝐫(u)]=\mathrm{\Gamma }_{}+{\displaystyle \frac{1}{2}}(\beta \mathrm{})^2{\displaystyle _0^\beta \mathrm{}}𝑑u{\displaystyle _0^\beta \mathrm{}}𝑑u^{}\mathrm{\Gamma }(uu^{})\times |𝐫(u)𝐫(u^{})|^2`$ (14)
where $`\mathrm{\Gamma }(uu^{})`$ is a solvent-induced force constant between different sites on the electron polymer and $`\mathrm{\Gamma }_{}`$ merely determines the zero of energy. The Bogoliubov inequality provides a upper bound for the excess chemical potential,

$`\mathrm{\Delta }\mu \beta ^1\mathrm{ln}Z_{\mathrm{ref}}+<\mathrm{\Delta }\mu [𝐫(u)]\mathrm{\Delta }\mu _{\mathrm{ref}}[𝐫(u)]>_{\mathrm{ref}}`$ (15)
Here, $`Z_{\mathrm{ref}}`$ is the electronic partion function for the Gaussian reference system and $`<>_{\mathrm{ref}}`$ means the average over the reference system weight determined by $`S_{}+\mathrm{\Delta }\mu [𝐫(u)]`$. Minimizing the right hand side of Eq.(14) provides the optimal Gaussian reference system. This procedure leads to the following equations. The correlation function for the intrapolymer correlation in k-space:

$`\widehat{\omega }(k,\tau )=\mathrm{exp}[{\displaystyle \frac{k^2R^2}{2D}}]`$ (16)
where

$`R^2(\tau )=|𝐫(u)𝐫(u^{})|^2=4D{\displaystyle \underset{n1}{\overset{\mathrm{}}{}}}A_n[1\mathrm{cos}(\mathrm{\Omega }_n\tau )]`$ (17)
is the mean square displacement between two points on the electron path separated by a imaginary time increament $`0uu^{}\beta \mathrm{}`$ with

$`A_n=(\beta m\mathrm{\Omega }_n^2+\gamma _n)^1`$ (18)
where $`\mathrm{\Omega }_n=\frac{2\pi n}{\beta \mathrm{}}`$, and

$`\gamma _n=(D\beta \mathrm{})^1{\displaystyle _0^\beta \mathrm{}}𝑑u[1\mathrm{cos}(\mathrm{\Omega }_nu)]{\displaystyle \frac{d^Dkk^2}{(2\pi )^D}v(k)\mathrm{exp}(k^2R^2(u)/2D)}.`$ (19)
We solve Eq.(12) for $`h`$ and $`c_{es}`$ using closure relation

$`g=0\text{for}rd`$ (20)
$`c_{es}=\beta ϵ{\displaystyle \frac{\mathrm{exp}(\alpha r)}{\alpha r}}\text{for}r>d`$ (21)
We can express Eq.(12) and the closure (20) in the variational form

$`{\displaystyle \frac{\delta I_{RISM}}{\delta c_{es}}}=0`$ (22)
where

$`I_{RISM}=\rho _s\widehat{c}_{es}(0)+{\displaystyle \frac{1}{2}}{\displaystyle \frac{d^Dk}{(2\pi )^D}\widehat{c}_{es}^2(k)\widehat{\chi }_{ss}(k)\widehat{\omega }(k)}`$ (23)
where $`\widehat{c}_{es}(0)`$ is the k = 0 spatial Fourier transform of $`c_{es}(r)`$, $`D`$ is the dimensionality of the space, $`d^Dk`$ is the dimensionality dependence volume element, $`u`$ labels the beads in the polymer ring $`(0u\beta \mathrm{})`$, $`m`$ is the bare electron mass, and $`v(k)=\widehat{c}_{es}^2(k)\widehat{\chi }_{ss}(k)`$ is the Fourier transform of the potential between beads, which is found in Eq.(10). The information about the electron-solvent atom interation is contained in the closure relation.

In Eq.(11), $`\omega `$(r) is the intra polymer distribution function averaged over all beads of the ring polymer. In writing Eq.(11) it has been assumed that for each polymer configuration, the solvent sees only average polymer rather than individual beads. The intra polymer distribution function $`\omega (r,\tau )`$ is determined in the polaron approximation ref1 ,$`\omega (k,u)`$ is the Fourier transform of $`\omega (r,u)`$ and is given by

$`\widehat{\omega }(k,u)=\mathrm{exp}[k^2R^2(u)/2D].`$ (24)
Eq.(22) is solved self-consistently ref30  for a given solvent and model potential representing the electron-fluid particle interaction. This solution gives information about $`v(𝐫)`$, $`\omega (𝐫)`$, $`R(u)`$, $`\gamma _n`$ and $`g`$(r) \[ = 1+$`h`$(r)\]. Note that the quantity R(u) is the root mean square (RMS) value of the displacement between two points on the electron path separated by a time increment $`0u\beta \mathrm{}`$. The characteristic size or breadth of the polymer is measured by $`R(\beta \mathrm{}/2)`$. This is a measure of the spread of the wave packet associated with the particle. Since in the CSR theory a periodic boundary condition , $`𝐫(0)=𝐫(\beta \mathrm{})`$, has been imposed on the path of the electron, $`R(u)`$ is found to be symmetric about $`u=\frac{1}{2}\beta \mathrm{}`$, i.e. it starts from zero at u = 0, attain a maximum value at $`u=\frac{1}{2}\beta \mathrm{}`$ and decreases for $`u>\frac{\beta \mathrm{}}{2}`$ reaching zero at $`u=\beta \mathrm{}`$. $`g(r)`$ gives information about the average packing of solvent particles around the electron. The variational parameter $`\gamma _n`$ measures the strength of the electron fluid coupling. Quantities such as average kinetic energy, potential energy and effective mass etc., can be expressed in terms of $`\gamma _n`$ ref9 ; ref24 ; ref27  as

$`K.E.={\displaystyle \frac{D}{2}}k_BT[1+{\displaystyle \frac{\gamma _n}{\beta m\mathrm{\Omega }_n^2+\gamma _n}}].`$ (25)
$`P.E.=\rho _s{\displaystyle }d𝐫u_{es}(|𝐫|)g(|𝐫|).`$ (26)
$`{\displaystyle \frac{m}{m^{}}}=24{\displaystyle \underset{n0}{\overset{\mathrm{}}{}}}(4\pi ^2n^2+\gamma _n\lambda _e^2)^1.`$ (27)
In this equation $`m^{}`$ is the effective mass of the solvated electron. As mentioned earlier, we need two input for this theory. One is electron-solvent and another is the density-density correlation function of solvent which is related to the structure factor of the solvent. For $`D`$-dimensional hard sphere solvent under cocsideration the Percus-Yevic (PY) equation can be solved analytically for $`D`$ = 1 ref31  and for $`D`$ = 3 ref32 . For $`D=2`$ excellent results of thermodynamic and structural properties have been obtained by Baus and Colot ref33 . In the above $`\widehat{\omega }`$, $`\widehat{c}_{es}`$, and $`\widehat{\chi }`$, are the spatial Fourier transform of $`\omega `$, $`c_{es}`$, and $`\chi `$.

## III Results and Discussions

In presenting our results we mainly focus on the imaginary time correlation function

$`R^2({\displaystyle \frac{1}{2}}\beta \mathrm{})=|𝐫({\displaystyle \frac{1}{2}}\beta \mathrm{})𝐫(0)|^2=4D{\displaystyle \underset{n1}{}}A_n(1\mathrm{cos}\pi n)`$ (28)
and the electron-solvent radial distribution function. Note that $`R(\frac{1}{2}\beta \mathrm{})`$ gives a measure of the physical size of the electron chain or the spread of the wave packet associated with the electron. For a free particle,

$`R({\displaystyle \frac{1}{2}}\beta \mathrm{})=\sqrt{D/4}\lambda _e.`$ (29)
In Fig.1 we plot the reduced correlation length, S $`R(\frac{1}{2}\beta \mathrm{})/\sqrt{D/4}\lambda _e`$ which is the dispersion of the wavepacket associated with the solvated electron relative to the free particle in one, two and three dimension as a function of density for $`\lambda _e=15\sigma `$, $`\alpha `$ = $`\sigma ^1`$ and $`d/\sigma `$ = 0.29 for several values of the attractive interaction, $`\beta ϵ`$. From these figures we find that when the electron-solvent interaction is solely repulsive ($`\beta ϵ=0.0`$), the electron is always gets trapped inside a solvent cage as solvent density is increased in one and two dimensions. In three dimension \[see Fig. (1C)\] the repulsive interaction is not strong enough (because $`d/\sigma =0.29`$) to localize the electron due to ordered structure formed when the solvent density is high. When the attractive interaction between electron and solvent atom is large ($`\beta ϵ100`$), the reduced correlation length of the electron is very small at very low solvent density and it stays almost constant upto $`\rho ^{}0.5`$. This behavior strongly indicates the electron is localized to a single atom irrespective of space dimensionality. As the space dimensionality increases, the range of the constant value of the solvent density decreases for example for $`D`$=1, S is constant upto $`\rho ^{}0.6`$. One interesting feature we have noted from Figs. 1A , 1B, and 1C is that as the attractive interaction increases, the reduced correlation length of the electron decreases as we increase the solvent density in the low-density regime, but the correlation turns upwardto have a maximum in the middle density regime, and finally the electron is trapped in the high density regime. If we carefully look at Figs.1, we observe that as we increase the space dimensionality, we find these effect are pronounsed for large attractive interaction. The similar type of behavior has been observed for the other values of $`d/\sigma `$. When $`d/\sigma `$ is less than $`0.29`$ we found the reduced correlation length of the electron decreases as we increase the attractive interaction i.e. $`\beta ϵ`$ and is self trapped at lower density. Opposite trend is found for increasing $`d/\sigma `$.

From Figs. 1\[A\], 1\[B\] and 1\[C\] it is clear that at the begining when we start increasing the attractive interaction ($`\beta ϵ20`$) the wave packet associated with the electron increases then further increase of electron solvent attractive interaction i.e. when $`\beta ϵ>50.0`$ the electron is localized on a single solvent atom. The reason behind the delocalized state at small vales of attractive interaction ($`\beta ϵ`$) is the cancellation between the repulsive and attractive electron-solvent interaction. In Fig. 2\[A\] to Fig.2\[F\] we plot the electron-solvent radial distribution function for one, two and three dimension at various solvent density and attractive iteraction ($`\beta ϵ`$). Figs. 2 can be explained in consistent manner with the same physical picture given for Figs. . When the attractive interaction is weak, the electron pushes the solvent atoms makes space to self-trap. When the attractive interaction is very large ($`\beta ϵ100.0`$) the electron is trapped on a single atom irrespective of space dimensionality, as evidenced by a large peak in $`g(r)`$ \[see Figs. 2\[F\].

In Fig. 3\[A\], 3\[B\], and 3\[C\] we plot the reduced correlation length as a function of $`\lambda _e/\sigma `$ at $`\beta ϵ=0`$ for $`D`$=1, 2, and 3 respectively. In one and two dimension ( Fig.3\[A\] and Fig.3\[B\] ) we observed as we increase $`\lambda _e`$ the electron is trapped in the solvent cage. If we compare Fig. 3A with Fig. 3B as the space dimensionality increases, the electron is strongly trapped in the solvent cage. One interesting feature we found in Fig. 3\[C\] at $`\rho ^{}0.7`$ there is sharp transition from delocalized to self-trapped state. In less than two dimension when there is no attractive interaction, the electron will always be in a localized stateref16  but in three dimension there is always a tendendency the electron will be in delocalized state. On the other hand, as the temperature is decreased, the electron has tendency to be in the localized state. This competition between vaious length scales probably make the sharp transition around $`\lambda _e16\sigma `$.

In Fig. \[4A\], \[4B\] and \[4C\] we plot the reduced correlation length as a function of density ($`\rho ^{}`$) for various values of $`d/\sigma `$. We found for low $`d/\sigma `$ (for example $`d/\sigma `$ = 0.15) the attractive interaction dominates and the electron is self-trapped in the low density regime. On the other hand, when $`d/\sigma `$ is large the self-trapping of the electron is dominated by repulsive interaction or self-trapping is by cage effect. In the case of electron in helium and xenon, the self-trapping of electron in helium and delocalized state of the electron in xenon has been explained on the basis of $`d/\sigma `$ value ref25  when elecron-solvent interaction is repulsive i.e. $`\beta ϵ=0`$.

## IV Concluding Remarks

The CSR theory for the excess electrons in simple fluid ( consisting of spherical atoms ) has been studied in one, two, and three dimension. The detailed study led us to the following conclusions :

1. The reduced correlation length, S, is very sensitive to the nature of the electron-fluid atom interaction, thermal wavelength of the excess electron ( $`\lambda _e`$ ) and a length associated with the mean volume occupied by each fluid atom which is related to $`\rho ^{^{\frac{1}{D}}}`$, where $`\rho ^{}=\rho \sigma ^D`$, $`\rho `$ being the number density of the fluid atoms and $`\sigma `$ is the effective diameter of a fluid atom. The behavior is interpreted in terms of an interplay among the length scales noted above. When the electron-solvent attractive interaction is very large, the density dependence of the size of the electron polymer relative to the free particle is dominated essentially by the pair attractive interaction, and the electron is trapped in a single solvent atom irrespective of the dimensionality ($`D`$). On the contrary, if the attractive interaction is absent, the electron is trapped in a cage formed by the solvent atoms at higher density. When the attractive interaction is low ($`\beta ϵ10.0`$) due to cancellation between repulsive and attractive interaction, the electron is delocalized irrespective of the dimensionality.
2. In the one dimensional case the reduced correlation length is almost independent of $`d/\sigma `$ when electron-solvent attractive interaction is absent. On the other hand, in two and three dimensions, $`S`$, the reduced correlation length of the electron is sensitive on the value of the $`d/\sigma `$.
3. In three dimension we found when $`d/\sigma =0.29`$ and $`\beta ϵ=0.0`$, the temperature dependence of the reduced correlation length, $`S`$, shows a sharp transition (metal-insulator type) around $`\lambda _e/\sigma 15`$ at $`\rho ^{}0.7`$.
4. The long range electron-solvent attraction dominates at low density, while hard core repulsion dominates at high densities. When both interactions are present, they can counterbalance each other. Consequently, at some intermediate density, the effective electron-solvent interaction can be quite small, resulting in the electron delocalization.

In the present work we have explored the electron self-trapping in simple liquid for various space dimension ($`D`$= 1,2, and 3) which have same model potential (solvent-solvent interaction potential and electron-solvent interation potential) in every space dimension. In every space dimension (i.e. $`D`$=1, 2, and 3) we have in the model system considered here there is repulsive as well as attractive electron solvent interation potential which is not considered in Refs. 25, 27, and 28. We have presented the dimensionality dependent result for self trapping behavior of electron in simple liquid.

###### Acknowledgements.

This work was supported in part by National Science Foundation, the Robert Welch Foundation, and the Institute for Molecular Science, (Japan).

## Appendix A  $`D`$-space dimensional integration

In $`D`$-space we have for any function f(x) depending only on the distance $`x=|𝐱|`$,

$$𝑑𝐱f(x)=S_D_0^{\mathrm{}}𝑑xx^{D1}f(x),$$
$`(A.1)`$
or on $`x`$ and one integration angle $`\theta `$,

$$𝑑𝐱f(x,\theta )=S_{D1}_0^{\mathrm{}}𝑑xx^{D1}𝑑\theta \mathrm{sin}^{D2}\theta f(x,\theta ),$$
$`(A.2)`$
where $`S_D=DV_D`$ is the surface area of the unit sphere of volume $`V_D=\pi ^{D/2}/\mathrm{\Gamma }(1+D/2)`$,where $`\mathrm{\Gamma }(1+z)=z!`$ is the $`\mathrm{\Gamma }`$ function.


format: text

text:
# Annihilation catastrophe: new critical phenomenon in diffusion-controlled dynamics

For the last decade the reaction-diffusion system $`A+B0`$, where unlike species $`A`$ and $`B`$ diffuse and irreversibly react in the bulk of a $`d`$-dimensional substrate, has acquired the status of one of the most popular objects in nonequilibrium statistical physics . Two problems have been investigated most intensively: (i) dynamics of the $`A+B0`$ annihilation in an infinite system with initially homogeneously (randomly) and equimolarly distributed reactants, in which below the critical dimension $`d_c=4`$ a dynamical clustering develops (fluctuation-induced like-particle domains formation) and, as a result, an anomalous reaction deceleration arises (Ovchinnikov and Zeldovich, 1978 , Toussaint and Wilczek, 1983 ); (ii) behaviour of the dynamic reaction front in an infinite system with initially spatially separated reactants (Galfi and Racz, 1988 ), and structure of the steady state front in a finite system, at the ends of which are injected equal currents of $`A^{}s`$ and $`B^{}s`$ particles (Ben-Naim and Redner, 1992 , Cornell and Droz, 1993 ).

In 1994 it was discovered (Shipilevsky, 1994 ) that in another wide class of RD systems, where reaction and diffusion are spatially separated (i.e. reaction proceeds on the surface of the medium and diffusion proceeds in its bulk) the interplay between reaction and diffusion acquires qualitatively new features and leads to the threshold self-organizing dynamics of the $`A+B0`$. It has been found that once particles $`A`$ and $`B`$ diffuse at different mobilities from the bulk of finite medium onto the surface and die on it by the reaction $`A+B0`$, there should exist some threshold difference in the initial numbers of $`A`$ and $`B`$ particles, $`\mathrm{\Delta }_c`$, above which the loop of positive feedback is ”switched on” and the process of their death, instead of usual deceleration, starts to accelerate autocatalytically. Recently, I have shown  that deceleration-acceleration transition is only a prelude to much more nontrivial dynamical effects, which take place far beyond the acceleration threshold. It has been demonstrated that in the diffusion-controlled limit $`\mathrm{\Delta }\mathrm{}`$ a new critical phenomenon develops, which arises as a result of self-organizing explosive growth (drop) of the surface concentrations of, respectively, slow and fast particles (concentration explosion) and manifests itself in the form of an abrupt singular jump in the desorption flux relaxation rate. The aim of this report is to give an insight into the key features of this remarkable phenomenon which I call annihilation catastrophe.

Model. I consider a model, in which species $`A`$ and $`B`$ are supposed to be initially uniformly distributed in the bulk of infinitely extended slab of thickness $`2\mathrm{}`$. Both species diffuse to the surface $`(X=\pm \mathrm{})`$ and desorb as a result of surface reaction $`A_{ads}+B_{ads}AB0`$. The rate of reaction is proportional to the product of surface concentrations $`I=k\rho _A\rho _B=\kappa c_{As}c_{Bs}`$, where the rates of surface-subsurface exchange are supposed to be rather high, so that quasiequilibrium $`\rho _i=f_ic_{is}`$ is always sustained. The system is effectively one dimensional and the boundary conditions are determined from the equality of diffusion $`I^D`$ and desorption $`I`$ flux densities at the surface $`I^D_s=I`$, i.e. it is assumed that the surface layer capacity can be neglected. Thus, introducing index $`\mathrm{"}H\mathrm{"}`$(heavy) for slower diffusing species and index $`\mathrm{"}L\mathrm{"}`$(light) for a faster one, the problem reads ( by symmetry I consider the interval \[$`0,\mathrm{}`$\] only)

$`h/\tau =^2h,l/\tau =(1/p)^2l,`$ (1)
$`h_s=(1/p)l_s=h_sl_s,`$ (2)
with $`(h,l)_{x=0}=0`$ and the initial conditions $`h(x,0)=h_0`$ and $`l(x,0)=l_0`$. Here $`h(x,\tau )=c_H/c_{}`$ and $`l(x,\tau )=c_L/c_{}`$ are the reduced concentrations, $`/x,x=X/\mathrm{}[0,1]`$ is the dimensionless coordinate, $`\tau =D_Ht/\mathrm{}^2`$ is the dimensionless time, $`p=D_H/D_L1`$ is the ratio of species diffusivities, and $`c_{}=D_H/\kappa \mathrm{}`$ is the characteristic concentration scale, from which diffusion starts to play an essential role. Boundary conditions (2) reflect the fact that particles disappear in pairs only, i.e. $`J=h_sl_s=<\dot{h}>=<\dot{l}>`$, where $`J=I/I_{}`$ is the reduced desorption flux density and $`I_{}=\kappa c_{}^2`$ is its characteristic scale, therefore $`<h><l>=\mathrm{\Delta }=\mathrm{const}.`$, i.e., the excess amount stays ”inert” in the bulk ( here $`<h>=_0^1h𝑑x=𝒩_H/𝒩_{}`$ and $`<l>=_0^1l𝑑x=𝒩_L/𝒩_{}`$ are the total reduced numbers of particles in the bulk per unit of surface and $`𝒩_{}=c_{}\mathrm{}=D_H/\kappa `$ is the characteristic particle number scale). I consider here the annihilation dynamics for $`\mathrm{\Delta }>0`$ when the pair number $`N`$ is dictated by the $`L`$ particles number $`N(\tau )=𝒩_{pair}/𝒩_{}=N_L(\tau )`$ so that in the final state all $`L`$ particles disappear $`N(\mathrm{})=l(x,\mathrm{})=0`$ and $`H`$ particles are distributed uniformly with concentration $`h(x,\mathrm{})=\mathrm{\Delta }`$. I am interested first and foremost in the diffusion-controlled limit $`N_0\mathrm{}`$ (note that diffusion-controlled limit $`N_0=l_0=c_L(0)/c_{}\mathrm{}`$ means $`c_{}^1=\kappa \mathrm{}/D_H\mathrm{}`$). Moreover, I focus on the limit $`n_0=N_0/\mathrm{\Delta }1`$ when the initial numbers of $`H`$ and $`L`$ particles are almost equal, $`r=h_0/l_0=(1+n_0)/n_01`$. Precisely in this limit the particles interplay is most ”intensive” and gives rise to the explosive dynamics, that is my prime interest here.

Transition point and self-accelerating annihilation. Assuming $`N_0\mathrm{}`$ and $`=r\sqrt{p}<1`$, I show , that for $`n_01`$ and $`\gamma =11/\sqrt{N_0}`$ at times $`\tau \tau _h=1/N_0^20`$ the annihilation is followed by a sharp drop of $`h_s`$, crossing over to the $`H`$-diffusion-controlled regime, where the desorption flux density is governed by the rate of $`H`$ particles diffusion to the system’s surface. In the range $`\tau _h\tau 1`$, when diffusion of $`H`$ particles proceeds in a semi-infinite medium regime, $`Jh_0/\sqrt{\pi \tau }J_0/\sqrt{\pi }(\tau _h/\tau )^{1/2}`$ and for surface concentrations I find

$$h_s=(r/\gamma \sqrt{\pi \tau })(1+\varrho +\mathrm{}),l_s=l_0\gamma (1\varrho +\mathrm{}),$$

where $`\varrho =2r\sqrt{\tau /\pi }`$ at $`p1`$. At $`\tau >1`$ the slow power-law stage crosses over to the exponential relaxation stage $`J=𝒜e^{\omega _0\tau }(1+e^{8\omega _0\tau }+\mathrm{})`$ with amplitude $`𝒜=2h_0`$, and I show that in the limit $`𝒜\mathrm{}`$ the surface concentrations achieve universal ($`\mathrm{\Delta }`$-independent) asymptotics

$`h_s=\mathrm{\Delta }_c(1+\varphi +e^{8\omega _0\tau }+\mathrm{}),l_s=(𝒜/\mathrm{\Delta }_c)e^{\omega _0\tau }(1\varphi +\mathrm{})`$ (3)
where $`|\varphi |`$ decays as $`O(e^{8\omega _0\tau })`$ in the range $`pp_c=4/9`$ and as $`O(e^{(\pi ^2/p\omega _0)\tau })`$ in the range $`p_c<p<1`$ , and $`\mathrm{\Delta }_c`$ obeys the expression

$$\mathrm{\Delta }_c=\sqrt{\omega _0/p}\mathrm{tan}(\sqrt{\omega _0p})$$

($`\omega _0=\pi ^2/4`$ being the main eigenfrequency of the diffusion field relaxation under absorbing boundary conditions). This important result implies that in the diffusion-controlled limit $`𝒜\mathrm{}`$ the critical asymptotics $`h_s=h_s^c=\mathrm{\Delta }_c`$ appears at any $`\mathrm{\Delta }`$ as a precursor of long-time asymptotics $`h_s_\tau \mathrm{}=\mathrm{\Delta }`$, and precisely this intermediate asymptotics $`h_s^c`$ selects the critical value $`\mathrm{\Delta }=\mathrm{\Delta }_c=h_s^c`$ below and above of which the relaxation dynamics of $`h_s`$ is qualitatively different: at $`\mathrm{\Delta }<\mathrm{\Delta }_c`$ the surface concentration, $`h_s`$, escaping from the asymptotics $`h_s^c`$, decreases, reaching its limiting value $`\mathrm{\Delta }`$ from above, whereas at $`\mathrm{\Delta }>\mathrm{\Delta }_c`$ the surface concentration, $`h_s`$, escaping from the asymptotics $`h_s^c`$, increases, reaching its limiting value $`\mathrm{\Delta }`$ from below . With neglect of rapidly decaying terms of $`O(e^{8\omega _0\tau },e^{(\pi ^2/p\omega _0)\tau })`$ the exact asymptotic expansion in power of $`1/𝒜`$ has the form

$$h_s=\mathrm{\Delta }_c(1+\lambda +\mathrm{}),l_s=(𝒜/\mathrm{\Delta }_c)e^{\omega _0\tau }(1\lambda +\mathrm{}),N=(𝒜/\omega _0)e^{\omega _0\tau }(1\mu \lambda +\mathrm{}),$$

$`J=𝒜e^{\omega _0\tau }(1\sqrt{\omega _0}\mathrm{tanh}\sqrt{\omega _0}(\mathrm{\Delta }\mathrm{\Delta }_c)^1\lambda ^2+\mathrm{}),`$ (4)
with $`\mu =\omega _0/\mathrm{\Delta }_c1p`$ and $`\lambda =\mathrm{\Delta }_c(\mathrm{\Delta }\mathrm{\Delta }_c)𝒜^1e^{\omega _0\tau }.`$ These expressions and more detailed analysis in Ref. suggest that at the crossing of the critical point $`\mathrm{\Delta }=\mathrm{\Delta }_c`$ not only the sign of $`\dot{h}_s`$ changes, but a growth of $`h_s`$ (and with it the whole annihilation process) begins to accelerate autocatalytically, $`\dot{h}_s=\omega _0(h_s\mathrm{\Delta }_c)+\mathrm{}`$, as a result of ”switching-on” of the loop of positive feedback: diffusion-induced growth of $`h_s`$ accelerates the drop of $`l_s`$, the drop of $`l_s`$ accelerates the growth of $`h_s`$ and so on (self-accelerating pairs death). This phenomenon, interpreted in  as a new type of self-organization, develops the brighter the stronger $`\mathrm{\Delta }`$ exceeds $`\mathrm{\Delta }_c`$, and in the limit $`\mathrm{\Delta }\mathrm{}`$ it results in the explosive dynamics - annihilation catastrophe .

Annihilation catastrophe. Scaling and universality. According to (4) in the limit $`\mathrm{\Delta }\mathrm{}`$ at the self-acceleration stage the rates of growth $`\mathrm{\Omega }_{Hs}=+d\mathrm{ln}h_s/d\tau `$ and relaxation $`\mathrm{\Omega }_{Ls}=d\mathrm{ln}l_s/d\tau `$ of surface concentrations are unambiguously defined by the current value of the reduced pair number $`n(\tau )=N(\tau )/\mathrm{\Delta }`$ and are mutually ”compensated” so that with an accuracy to vanishingly small terms $`O(1/n^2\mathrm{\Delta })0`$ the desorption flux relaxation rate $`\tau _J^1=d\mathrm{ln}J/d\tau =\mathrm{\Omega }_{Ls}\mathrm{\Omega }_{Hs}=\omega _0+\mathrm{}`$ is sustained constant. The central result of  is my finding, that at $`\mathrm{\Delta }\mathrm{}`$ the relaxation rate $`\tau _J^1`$ is sustained constant up to the time moment $`\tau _{}`$ of reaching some critical pair number $`n_{}(p)`$, in whose vicinity $`\delta \tau \mathrm{\Delta }^{2/5}0`$ the $`\tau _J^1`$ abruptly jumps from $`\tau _J^1=\omega _0`$ to $`\tau _J^1\mathrm{}`$. I consider here the two consecutive stages of this phenomenon.

i) Concentration explosion. Following the work , I first show that condition $`\tau _J^1=\omega _0`$ in the limit $`\mathrm{\Delta }\mathrm{}`$ leads automatically to the arising of a finite-time singularity. Taking $`J=𝒜e^{\omega _0\tau }`$ and satisfying (1),(2) to the leading order in $`\mathrm{\Delta }`$, I find

$`j=J/\mathrm{\Delta }=\dot{n}=\omega _0(1+n),h_s=\mathrm{\Delta }_c(1+n)/(nn_{}),l_s=\mu \mathrm{\Delta }(nn_{}),`$ (5)
where $`n_{}=(1\mu )/\mu p/(1p)`$. It is seen that $`h_s`$ diverges as $`nn_{}`$ and at the point of singularity exactly $`j_{}=\mathrm{\Delta }_c`$. From equality $`J_{}=2h_0e^{\omega _0\tau _{}}=\mathrm{\Delta }\mathrm{\Delta }_c`$ in accord with the numerical data of  we obtain

$$\tau _{}=(4/\pi ^2)\mathrm{ln}[2(1+n_0)/\mathrm{\Delta }_c].$$

Taking the origin of time at the point $`\tau _{}`$, i.e. introducing a relative time $`𝒯=\tau \tau _{}`$, from (5) we come to remarkably simple laws $`n=(1/\mu )e^{\omega _0|𝒯|}1,j=\mathrm{\Delta }_ce^{\omega _0|𝒯|},h_s=\mathrm{\Delta }_c/(1e^{\omega _0|𝒯|}),l_s=\mathrm{\Delta }(e^{\omega _0|𝒯|}1)`$, and, hence, $`\mathrm{\Omega }_{Hs}=\mu (h_s\mathrm{\Delta }_c)=\omega _0/(e^{\omega _0|𝒯|}1)`$, whence it follows that for any $`p<1`$ at the ”distance” from the singularity point $`|𝒯|\omega _0^1`$ the dynamics of growth of $`h_s`$ changes drastically: the initial exponential increase in the growth rate $`\mathrm{\Omega }_{Hs}=\omega _0e^{\omega _0|𝒯|}(\mathrm{\Omega }_{Ls}\omega _0)`$ turns at $`\omega _0|𝒯|1`$ to a singular one

$`\mathrm{\Omega }_{Hs}=\mathrm{\Omega }_{Ls}=\mu h_s=1/|𝒯|,|𝒯|0,`$ (6)
i.e., concentration explosion develops in the system. The striking feature of this phenomenon is that at times $`\omega _0|𝒯|1`$ the desorption flux practically reaches its critical value $`J_{}`$, therefore the concentration explosion develops, in fact, at a ”frozen” flux $`h_sl_s=J_{},h_s=1/\mu |𝒯|,l_s=\mu J_{}|𝒯|,`$ the consequence of which is synchronization of the growth and relaxation rates $`\mathrm{\Omega }_{Hs}=\mathrm{\Omega }_{Ls}=\mathrm{\Omega }_s`$ during the explosion. A singular growth of $`h_s`$ at the background of the diffusion flux $`J_{}`$, ”frozen” on the explosion time scale, gives rise to a singular growth of the diffusion ”antiflux”, the calculation of which yields  $`\delta J/J_{}1/\mathrm{\Delta }|𝒯|^{3/2}`$. As a result, in the desorption flux relaxation rate $`\tau _J^1=\omega _0+[\tau _J^1]`$ the singularly growing term arises

$`[\tau _J^1]1/\mathrm{\Delta }|𝒯|^{5/2},`$ (7)
whence it follows that the scenario described above holds down to $`|𝒯|_{cat}\mathrm{\Delta }^{2/5}0`$, beyond which a catastrophic growth of $`\tau _J^1`$ begins.

ii) Two universality classes of passing through the point of singularity. Remarkably, the dynamics of passing through the point of singularity can be described analytically (which, in itself, is a rather rare event ) and it appears quite beautiful. I show analytically and confirm numerically that at the singularity point $`𝒯=0`$ the rate of explosion $`\mathrm{\Omega }_{Hs}`$ reaches a maximum $`\mathrm{\Omega }_{Hs}^M`$, the dynamics of passing through which proceeds in one of the two characteristic scaling regimes, depending on the value of parameter $`𝒦=p^{3/2}\mathrm{\Delta }/\mathrm{\Delta }_c`$.

1. Annihilation catastrophe at ”frozen” flux ($`𝒦1`$). In the limit $`𝒦1`$ at times $`|𝒯|𝒯_L=p`$ the particle number $`n=n_{}(1+\omega _0|𝒯|/p+\mathrm{})`$ is, in fact, ”frozen” $`nn_{}`$ and so in the process of explosion in the contracting as $`\delta x_L\sqrt{|𝒯|/p}`$ nearsurface layer there arises a sharp gradient of $`L`$ particles. I show that in this limit up to $`𝒯p`$ the flux and the particle number are retained ”frozen”, $`h_sl_s=J_{}`$ and, hence, the explosion develops synchronouosly, $`\mathrm{\Omega }_{Hs}=\mathrm{\Omega }_{Ls}=\mathrm{\Omega }_s`$, both before and after the singularity point. Equalizing the rates of change of $`H`$ and $`L`$ particles diffusion fluxes leads to the scaling

$`\mathrm{\Omega }_s=\mathrm{\Omega }_s^MS(𝒯/𝒯_f),h_s/h_s^M=l_s^M/l_s=w(𝒯/𝒯_f),`$ (8)
where the scaling functions $`S(\zeta )=1/\sqrt{1+\zeta ^2}`$ and $`w(\zeta )=\zeta +\sqrt{1+\zeta ^2}`$, the characteristic time $`𝒯_f=1/\mathrm{\Omega }_s^M`$, and the amplitudes $`\mathrm{\Omega }_s^M,h_s^M`$ and $`l_s^M`$ at the explosion maximum

$$h_s^M=(2/\mu )\mathrm{\Omega }_s^M=p^{1/4}\sqrt{\mathrm{\Delta }\mathrm{\Delta }_c},l_s^M=p^{1/4}\sqrt{\mathrm{\Delta }\mathrm{\Delta }_c}.$$

Two striking features of this result are symmetrical universalization of $`\mathrm{\Omega }_s=1/|𝒯|`$ beyond the scope of interval $`[𝒯_f,𝒯_f],𝒯_f0`$ and remarkable symmetry $`𝒯𝒯,w1/w`$.

Singular ”flash” of $`\mathrm{\Omega }_s`$ leads to a singular jump in $`\tau _J^1`$, the dynamics of which is described by the scaling

$`[\tau _J^1]=[\tau _J^1]_{}W(𝒯/𝒯_f),[\tau _J^1]_{}=c_1(\mathrm{\Omega }_s^M)^{5/2}/\mathrm{\Delta }\mu ^{5/4}p^{5/8}\mathrm{\Delta }^{1/4},`$ (9)
where in accord with the scaling function $`W(\zeta )=c_2_0^{\mathrm{}}𝑑\theta /\sqrt{\theta }[1+(\zeta \theta )^2]^{3/2}`$ the value of $`[\tau _J^1]`$ reaches at $`𝒯𝒯_f`$ a maximum $`\mathrm{max}[\tau _J^1][\tau _J^1]_{}\mathrm{\Delta }^{1/4}\mathrm{}`$, and then it decays as $`1/\sqrt{p𝒯}`$. Essentially, that $`\mathrm{max}[\tau _J^1]p^{5/8}`$, so the less is $`p`$ the more brightly the effect is displayed.

2.Annihilation catastrophe at uniform $`L`$ particles distribution. Flux breaking effect ($`𝒦1`$). If $`\mathrm{\Delta }\mathrm{}`$ but $`p0`$ so that $`𝒦1`$, the distribution of $`L`$ species is retained uniform $`(l_s=N,\mathrm{\Omega }_{Ls}=h_s)`$ and the final stage of explosion is changed radically. I demonstrate that in this case an intermediate scaling takes place

$`\mathrm{\Omega }_{Hs}=\mathrm{\Omega }_{Hs}^MQ(𝒯/𝒯_{}),h_s=h_s^Mv_H(𝒯/𝒯_{}),l_s=l_s^Mv_L(𝒯/𝒯_{})`$ (10)
where the characteristic time $`𝒯_{}=1/\mathrm{\Omega }_{Hs}^M`$ and amplitudes $`\mathrm{\Omega }_{Hs}^M,h_s^M`$ and $`l_s^M`$ at the explosion maximum

$$h_s^M=q\mathrm{\Omega }_{Hs}^M=b_H\mathrm{\Delta }^{2/3},l_s^M=b_L\mathrm{\Delta }^{1/3},$$

with $`q=2.24,b_H=1.63`$ and $`b_L=0.73`$. It is seen that in the limit $`𝒦1`$ the singular ”flash” of $`\mathrm{\Omega }_{Hs}`$ occurs much sharper, leading to more abrupt jump in $`\tau _J^1`$

$`[\tau _J^1]=[\tau _J^1]_{}V(𝒯/𝒯_{}),[\tau _J^1]_{}=(q1)\mathrm{\Omega }_{Hs}^M\mathrm{\Delta }^{2/3},`$ (11)
which, in turn, gives rise to a qualitatively new phenomenon: the finite jump of the flux by a factor of $`2`$ for the time $`|𝒯|_{break}𝒯_{}\mathrm{\Delta }^{2/3}0`$ (flux breaking effect).

As in the limit $`𝒦1`$ at $`𝒯p`$ the relaxation rate $`\tau _J^1`$ reaches its $`L`$-diffusion-controlled limit $`\tau _J^1\omega _0/p`$, and in the limit $`𝒦1`$ it always grows, I point out finally (and demonstrate numerically) that independently of $`𝒦`$ in the limit $`p0,\mathrm{\Delta }\mathrm{}`$ there arises a most spectacular consequence of the annihilation catastrophe:an abrupt, practically instantaneous (on the scale of $`\omega _0`$) disappearance of the flux.

In conclusion, I believe that annihilation catastrophe may pretend to be not only one of the most dramatic manifestations of the reaction-diffusion interplay, but it is a new type of catastrophe wherein two explosive processes ($`\mathrm{\Omega }_{Hs}`$ and $`\mathrm{\Omega }_{Ls}`$) are developing simultaneously, effectively ”compensating” one another so that for an external observer ($`J`$) the explosion dynamics goes unnoticed up to the point of singularity, in the vicinity of which decompensation of explosions is manifested as a sudden singular jump.
This work was supported by the RFBR Grant No. 02-03-33122-a.
References:  For a review see E. Kotomin and V. Kuzovkov, Modern Aspects of Diffusion Controlled Reactions: Cooperative Phenomena in Bimolecular Processes (Elsevier, Amsterdam, 1996);  A.A. Ovchinnikov and Ya.B. Zeldovich, Chem. Phys. 28, 215 (1978);  D. Toussaint and F. Wilczek, J. Chem. Phys. 78, 2642 (1983);  L. Galfi and Z. Racz, Phys. Rev. A 38, 3151 (1988);  E. Ben-Naim and S. Redner, J. Phys. A 25, L575 (1992);  S. Cornell and M. Droz, Phys. Rev. Lett. 70, 3824 (1993);  B.M. Shipilevsky, Phys. Rev. Lett. 73, 201 (1994);  B.M. Shipilevsky, Phys. Rev. Lett. 82, 4348 (1999);  B.M. Shipilevsky, J. Phys. A 30, L471 (1997);  B.M. Shipilevsky (to be published);  L.P. Kadanoff, Singularities and Blowups, Phys. Today, 11 (1997).


text:
# Polarized Gravitational Waves from Gamma-Ray Bursts

## 1 Introduction

The observations of gamma-ray burst (GRB) afterglows at energies ranging from X-rays to radio have lead to an increased understanding of the possible geometry of the inferred relativistic outflow or ejecta (see van Paradijs, Kouveliotou & Wijers 2000 for a review; also Rhoads 1999; Sari, Piran & Halpern 1999). At least for the class of long bursts (durations $``$ 10 s) there is now significant evidence that the ejecta is collimated in a jet.

The evidence for jets is based on a change of the GRB light curve time-dependence power law index, and there appears to be a correlation in the sense that bursts with the largest gamma-ray fluences have the narrowest implied opening angles. Frail et al. (2001), Piran et al. (2001) and Panaitescu & Kumar (2002) reported that the total gamma-ray energy release, after correcting for the collimation and distance determined by the afterglow observations, are narrowly clustered around $`5\times 10^{50}`$erg, and suggested that the broad distribution in fluence and luminosity for GRBs is largely the result of a wide variation of the opening angles. This result can be obtained assuming that the jet emission inside some angle $`\theta _j`$ is approximately uniform and then drops off abruptly beyond that. However, Rossi, Lazzati & Rees (2002) and Zhang & Mészáros (2002) showed that the observations can also be interpreted in terms of an alternative model where the jet, rather than having a uniform profile out to some definite cone angle, has a universal angle-dependent beam pattern with a luminosity per unit solid angle which is maximal along the axis, and drops off gradually away from the axis. In this model, it is the difference of viewing angles $`\theta `$ which causes a wide variation in the apparent luminosity of GRBs, $`L\theta ^2`$.

Gravitational waves (GW) are also expected from some types of GRB (e.g. Kobayashi & Mészáros 2002). The GW emissivity and its polarization are angle-dependent, and may in principle be measurable, depending on the signal to noise ratio and on the detector alignments. In this Letter we discuss the prospects of exploiting measurements of the gravitational polarization degree to obtain information on the nature and orientation of the GRB system, providing additional constraints on the luminosity of the bursts, and we mention possible implications for the interpretation of dark GRB and X-ray flashes.

## 2 Gravitational Waves from GRB Progenitors

Most GRB models involve a stellar system resulting in a rotating black hole and a massive ($`0.11M_{}`$) disk of accreting matter around it. Such massive disks can form from the fall-back of debris during the formation of the central compact object, which ultimately is likely to lead to a black hole. Several scenarios could lead to a black hole-accretion disk system (e.g. Mészáros, 2002). This includes double neutron star binary mergers (Eichler et al 1989; Ruffert et al. 1997), neutron star - black hole mergers (Paczynski 1991; Janka et al. 1999) and massive stellar collapses (collapsars) (Woosley 1993; MacFadyen & Woosley 1999).

Coalescing compact binaries are one of the most promising GW sources detectable by the laser interferometer gravitational wave observatory (LIGO), and other detectors such as VIRGO, GEO600 and TAMA300. The binary coalescence process can be divided into three phases: inspiral, merger and ringdown (e.g., Flanagan & Hughes 1998). (1) During the inspiral phase, the gravitational radiation reaction time scale is much longer than the orbital period. As the binary losses energy by gravitational radiation, the masses gradually spiral in toward each other. (2) The merger begins when the orbital evolution is so rapid that adiabatic evolution is no longer a good approximation. The two masses then violently merge to form a black hole on a timescale comparable to a small multiple of the dynamical time. A bar-mode instability develops in the early stages of the merger, releasing a fraction of their rest mass energy in GWs. However, a significant fraction of the stellar material retains too much angular momentum to cross the black hole horizon promptly. This creates a temporary disk of debris material around the black hole, whose accretion can power a GRB jet. (3) The black hole is initially deformed, and in a ringdown phase radiates away the energy associated with these deformations as GWs, until it settles down into a Kerr geometry.

Massive rotating stellar collapses and their associated GW emission have been calculated numerically (e.g. Fryer et al. 1999; McFadyen & Woosley 1999; Rampp et al 1998). These numerical estimates are not conclusive, as a number of effects (including general relativity, secular evolution, etc, see Rampp et al 1998) have been neglected, but they suggest that GW emission from massive collapses may be much less important than from compact binary mergers. On the other hand, recent semi-analytical estimates (van Putten 2001; Fryer, Holz & Hughes 2002; Davies et al 2002) have indicated that massive collapses might emit significantly stronger GW signals than expected from the previous numerical estimates. Collapsars, i.e. massive stellar collapses leading to a GRB, require a high core rotation rate to form a centrifugally supported disk around a central, possibly spinning black hole. A high rotation rate, however, is conducive to the development of bar or fragmentation instabilities in the collapsing core or/and in the massive disk (Nakamura & Fukugita 1989; Bonnell & Pringle 1995; van Putten 2002). The asymmetrically infalling matter perturbs the black hole geometry, which leads to ringdown gravitational radiation.

## 3 Polarization of Gravitational Waves

We have seen that various mechanisms, including binary, bar and fragmentation instabilities and oscillations of black holes, could occur in GRB progenitor systems, resulting in GW emission. It is known that the GWs generated by a binary, bar or fragmentation instability are dominated by modes with spherical harmonic index $`l=m=2`$ mode (e.g. Misner, Thorne & Wheeler 1978). The most slowly damped quasi-normal mode of a Kerr black hole also has indices $`l=m=2`$ (Detweiler 1980). Since this mode may be preferentially excited in the presence of binary masses or fragmentation of a massive disk, the ringdown gravitational radiation is also dominated by $`l=m=2`$ mode.

The polarization components of $`l=m=2`$ mode depend on the inclination angle $`\alpha `$ (e.g. Thorne 1980) as

$$h_+(1+\mathrm{cos}^2\alpha ),h_\times 2\mathrm{cos}\alpha .$$
(1)
GRB progenitors emit GWs more strongly along the polar axis than in the equatorial plane, the latter being the orbital plane of the binary, the disk fragments or the bar. Since the GRB jets are launched along the polar (angular momentum) axis (i.e. $`\alpha =\theta `$), the GRB progenitors are stronger GW sources than the average non-bursting merger or collapse, as pointed out by Kochanek and Piran (1993) in the context of a neutron star binary model. When we observe a GRB from a viewing angle $`\theta `$ respect to the polar axis, the GW signal is stronger by a factor of $`(5/32)^{1/2}[(1+\mathrm{cos}\theta )^4+(1\mathrm{cos}\theta )^4]^{1/2}`$ than the value averaged over all possible viewing angles. Since for the angle dependent jet model (Rossi et al. 2002; Zhang & Mészáros 2002) the difference of the viewing angles causes a variation of the apparent luminosity of the GRB, we expect a correlation between the apparent luminosity of GRBs and the amplitude of the associated GWs, even after distance corrections. Also in the case when the jet profile is uniform inside an opening angle $`\theta _j`$, the typical viewing angle is $`2\theta _j/3`$, and the correlation is still expected.

The GWs will be circular polarized when viewed along the polar axis ($`\theta =0`$), while the $`+`$ polarization dominates when viewed along the equatorial plane ($`\theta =\pi /2`$). Therefore, the polarization of the GWs is also correlated with the luminosity of the GRBs. We define a polarization tensor in analogy with the electromagnetic approach (e.g. Landau & Lifshitz 1975). A plane wave at the position of a detector can be written in the form $`h_+=\text{Re}\left\{A_+e^{i\omega t}\right\}`$ and $`h_\times =\text{Re}\left\{A_\times e^{i\omega t}\right\}`$. If the plane wave is monochromatic, the complex amplitudes $`A_+`$ and $`A_\times `$ will be constants, whereas if the wave contains frequencies in a small interval $`\mathrm{\Delta }\omega `$ we take $`\omega `$ to be some average frequency inside this range. $`A_+`$ and $`A_\times `$ are slowly varying functions of time, compared to the wave oscillation period. The polarization tensor $`\rho _{ab}(a,b=+,\times )`$ is defined by using the time averaged values of the amplitudes $`A_aA_b^{}`$ and described by the Stokes parameters $`\xi _1`$ , $`\xi _2`$ and $`\xi _3`$:

$$\rho _{ab}\frac{A_aA_b^{}}{|A_+|^2+|A_\times |^2}=\frac{1}{2}\left(\begin{array}{cc}1+\xi _3& \xi _1i\xi _2\\ \xi _1+i\xi _2& 1\xi _3\end{array}\right).$$
(2)
For the $`l=m=2`$ mode, the parameters are

$`\xi _1`$ $`=`$ $`0,\xi _2={\displaystyle \frac{8\mathrm{cos}\theta (1+\mathrm{cos}^2\theta )}{(1+\mathrm{cos}\theta )^4+(1\mathrm{cos}\theta )^4}},`$ (3)
$`\xi _3`$ $``$ $`P={\displaystyle \frac{2(1+\mathrm{cos}\theta )^2(1\mathrm{cos}\theta )^2}{(1+\mathrm{cos}\theta )^4+(1\mathrm{cos}\theta )^4}}.`$ (4)
The degree of circular polarization and of linear polarization are given by $`\xi _2`$ and $`\sqrt{\xi _1^2+\xi _3^2}=\xi _3P`$, respectively. Thus, GRB progenitors whose polar axes are directed towards the earth would produce brighter GRBs, and circularly polarized GWs. As the viewing angle becomes larger relative to the polar axis, the luminosity of the GRB decreases $`\theta ^2`$, while the degree of linear polarization increases as $`P10^2(\theta /30\text{degree})^4`$.

## 4 Polarization Measurement

The response of an interferometer (interferometer $`1`$) to the gravitational radiation is given by a linear combination of two polarization components $`m_1=F_{+,1}h_++F_{\times ,1}h_\times `$ where the antenna patterns $`F_{+,1}`$ and $`F_{\times ,1}`$ depend on the orientation of the interferometer with respect to the GW source (e.g. Finn & Chernoff 1993). We assume that the interferometer arms are of the same length and that they meet at right angles. Define a right-handed coordinate system with one interferometer arm along the $`x`$-axis and the other along the $`y`$-axis. For simplicity, we assume that the source is in the direction of the $`z`$-axis. Since we can determine the position of a source in the sky by using observations of the GRB and afterglow, it is possible to generalize the following discussion to the case of a source with an arbitrary sky position. The angular momentum vector of the source (direction of the GRB jet) may be oriented in an arbitrary direction. We assume that the projection of the angular momentum vector to the $`xy`$ plane makes an angle $`\zeta `$ with the $`x`$-axis. With these conventions, the antenna patterns are given by $`F_{+,1}=\mathrm{cos}2\zeta `$ and $`F_{\times ,1}=\mathrm{sin}2\zeta `$.

To determine the polarization of GWs, one needs a network consisting of at least two interferometers which have different (i.e. non-parallel) arm orientations. Consider an identical interferometer (interferometer 2) at the same location as interferometer 1. (If we set the arrival time of the GRB signal as the origin of time at each interferometer, we can correct for the actual physical separation so that the two interferometers can be always considered to be at the same location). We assume that the interferometer 2 is rotated by an angle $`\pi /4`$ around the $`z`$-axis with respect to the interferometer 1. The response of the interferometer 2 is $`m_2=F_{+,2}h_++F_{\times ,2}h_\times `$, and the antenna patters are given by $`F_{+,2}=\mathrm{sin}2\zeta ,`$ and $`F_{\times ,2}=\mathrm{cos}2\zeta `$.

The detection of the polarization of GWs requires observations with a high signal-to-noise ratio (SNR) $`\rho `$. A detection is likelier in an optimal case where the wave forms of the polarized components, $`f_+`$ and $`f_\times `$, are known. We define the noise-weighted inner product as

$$s_i,f_a=4\text{Re}_0^{\mathrm{}}\frac{\stackrel{~}{s_i}(f)\stackrel{~}{f_a^{}}(f)}{S_h(f)}𝑑f$$
(5)
where $`\stackrel{~}{s_i}(f)`$ denotes the Fourier transform of the outputs of the two interferometers $`s_i(t)=m_i(t)+n_i(t)`$ $`(i=1,2)`$, $`n_i(t)`$ is the noise of the interferometers, $`S_h(f)`$ is the one-sided noise power spectral density and $`\stackrel{~}{f_a}`$ is the Fourier transform of $`f_a`$ $`(a=+,\times )`$. We normalize the functions $`f_a`$ as $`f_a,f_b=\delta _{ab}`$, hence the dispersion of the noise is unity $`\sigma ^2=\overline{n_i,f_a^2}=1`$ where the overline represents ensemble average. The linear polarization degree $`P=\xi _3`$ is

$$P=\frac{m_1,f_+^2+m_2,f_+^2m_1,f_\times ^2m_2,f_\times ^2}{m_1,f_+^2+m_2,f_+^2+m_1,f_\times ^2+m_2,f_\times ^2}.$$
(6)
Since the outputs of the interferometers $`s_i`$ are always superpositions of signals $`m_i`$ and noises $`n_i`$, the inner products $`m_i,f_a`$ themselves are not observable. We estimate $`m_i,f_a`$ by $`m_i,f_as_i,f_a`$. Substituting this into eq. (6), we can get an estimate $`P^{}`$ of the polarization degree. This approximation introduces an error $`\mathrm{\Delta }`$ in the estimate of the polarization degree $`P`$, given by

$`\mathrm{\Delta }`$ $``$ $`{\displaystyle \frac{P^{}P}{P}}={\displaystyle \frac{1+(AB)/P}{1+A+B}}1,`$ (7)
$`A`$ $`=`$ $`2a_+(\lambda _{1+}\mathrm{cos}2\zeta \lambda _{2+}\mathrm{sin}2\zeta )+\lambda _{1+}^2+\lambda _{2+}^2,`$ (8)
$`B`$ $`=`$ $`2a_\times (\lambda _{1\times }\mathrm{sin}2\zeta +\lambda _{2\times }\mathrm{cos}2\zeta )+\lambda _{1\times }^2+\lambda _{2\times }^2,`$ (9)
where $`a_+=\sqrt{2}(1+\mathrm{cos}^2\theta )[(1+\mathrm{cos}\theta )^4+(1\mathrm{cos}\theta )^4]^{1/2}`$, $`a_\times =2\sqrt{2}\mathrm{cos}\theta [(1+\mathrm{cos}\theta )^4+(1\mathrm{cos}\theta )^4]^{1/2}`$, $`\lambda _{ia}=n_i,f_a/\rho `$ and the SNR $`\rho `$ is $`\rho ^2=h_+,f_+^2+h_\times ,f_+^2`$. The fluctuations $`\lambda _{ia}`$ are of order of $`\rho ^1`$. Thus, a SNR of $`\rho >P^1`$ is required, if one wants to determine the polarization degree to an accuracy of order $`P`$.

To evaluate the error $`\mathrm{\Delta }`$ numerically, we assume that the $`\lambda _{ia}`$ are normally distributed, and that the projection of the angular momentum on the orbital plane $`\zeta `$ relative to the $`x`$-axis is uniformly distributed between $`0`$ and $`2\pi `$. Figure 1 shows the distributions of $`\mathrm{\Delta }`$ for $`10^6`$ random realizations in the case of a $`1\%`$ polarization degree ($`P=10^2`$). The SNR was assumed to be $`\rho =100,300,500`$ or 1000, for which the fractions of the realizations with large errors ($`|\mathrm{\Delta }|>0.5`$) are 80$`\%`$, 44$`\%`$, 20$`\%`$ and 1$`\%`$, respectively. Therefore, $`\rho 10P^1=1000`$ is required to determine a $`1\%`$ polarization degree.

The LIGO interferometers are coaligned, and thus cannot by themselves determine the polarization degree of GWs, although a network consisting of LIGO and other interferometers could in principle do it. We assume here an interferometer identical to the LIGO detectors with an optimal orientation (rotated by $`\pi /4`$ with respect to the LIGO interferometers). If some fraction of GRB progenitors are double neutron star mergers, the SNR is given by $`\rho _{ave}16(d/100\text{Mpc})^1`$ (e.g. eq. (16) in Kobayashi & Mészáros 2002) where $`d`$ is the distance to the binary. This estimate was obtained from an average over different possible orientation of the source and interferometer. The SNR $`\rho `$ is larger by factor of $`(5/4)\left[(1+\mathrm{cos}\theta )^4+(1\mathrm{cos}\theta )^4\right]^{1/2}`$. Then, in an optimal case, we can determine a $`1\%`$ polarization degree up to $`d_{max}7`$ Mpc. If we assume that $`1000`$ GRBs happen in a year within $`3000`$Mpc, the closest event in a year is, on average, at $`d300`$ Mpc. Since the position in the sky is random, the SNR is smaller by a factor of $`\sqrt{2/5}`$ than in the direction of the z-axis. The LIGO is most sensitive around $`f_0150`$Hz and the sensitivity is about $`\sqrt{f_0S_h(f_0)}3\times 10^{23}`$. Since $`d_{max}`$ is proportional to $`S_h^{1/2}`$, a future detector with sensitivity $`\sqrt{f_0S_h(f_0)}4\times 10^{25}`$ would enable us to measure a polarization degree of $`1\%`$ in a timescale of one year.

Frail et al. (2001) show that the observed distribution of jet opening angles (or viewing angles for the angle dependent jet model) is $`f_{obs}\theta _j^{3.5}`$ ($`\theta _j>0.05`$). Assuming this result, the probability that GWs from a GRB have a polarization degree larger than $`P=10^2`$ is given by $`3\times 10^2(P/10^2)^{0.4}`$. However, the $`f_{obs}`$ are estimated from observations of typical GRBs at redshifts $`z1`$, and one might fail to detect a significant fraction of the dim GRBs at large viewing angles. Therefore, when we study GWs from nearby GRBs, the probability that one observes GWs with a large polarization degree could be much larger than the above estimate.

## 5 Discussion and Conclusions

We have argued that GRB progenitors are likely to emit $`l=m=2`$ gravitational waves (GWs), which are circularly polarized on the polar axis, while the $`+`$ polarization dominates on the equatorial plane. Recent GRB studies suggest that the wide variation in the apparent luminosity of GRBs are caused by differences in the viewing angle, or possibly also in the jet opening angle. Since GRB jets are launched along the polar axis of GRB progenitors, correlations among the apparent luminosity of GRBs and the amplitude as well as the polarization degree of the GWs are expected.

At a viewing angle larger than the jet opening angle $`\theta _j`$ (which may be defined also in the case of a universal angle-dependent jet profile) the GRB $`\gamma `$-ray emission may not be detected. However, in such cases an “orphan” long-wavelength afterglow could be observed, which would be preceded by a pulse of GWs with a significant linearly polarized component. An expanding jet with an opening angle $`\theta _j`$ behaves, as long as its Lorentz factor $`\gamma >\theta _j^1`$, as if it were part of a spherical shell, but relativistic beaming effect allows only observers at viewing angles $`<\theta _j`$ to observe the emission from the jet. (By contrast, since GW emission is expected from the central engine itself, GWs are not subject to such extreme beaming). As the jet slows down and reaches $`\gamma \theta _j^1`$, the jet begins to expand laterally, and its electromagnetic radiation begins to be observable over increasingly wider viewing angles. Since the opening angle increases as $`\gamma ^1t^{1/2}`$ (Sari et al 1999), at a viewing angle $`\theta >\theta _j`$, the orphan afterglow begins to be observed (or peaks) at a time $`t_p\theta ^2`$ after the detection of the GW burst. The polarization degree and the peak time should be correlated as $`Pt_p^2`$.

A new type of fast transient source, called “X-ray flashes”, have recently been observed with the BeppoSAX satellite (Kippen et al. 2002). Apart from their large fraction of X-rays ($`210`$ keV), the overall properties of these events are similar to those of GRBs. Recently, Yamazaki, Ioka and Nakamura (2002) and Woosley, Zhang & Heger (2002) suggested that these events may be GRBs with large viewing angles. If this is the case, linearly polarized GWs should be observed prior to the X-ray flashes. The degree of polarization should be positively correlated with longer delays and with the softness of the X-ray flashes, which increase with angle.

Since the degrees of linear and circular polarization depend on the viewing angle, a determination of the polarization degree would be a measure of the viewing angle. Such measurements, which are likely to require the advent of a future generation of detectors, could provide a new tool for estimating the absolute luminosity of GRBs, including its photon component. By comparing the estimated absolute photon luminosity with the apparent luminosity, the distance to the source may be estimated independently of any redshift measurement. No optical afterglows have been found for about half of all the GRBs detected by BeppoSAX (the so called “dark GRBs”), and the present method would have the potential to help determine or constrain the distances to such dark GRBs.

We thank L. S. Finn, B. Owen, I. Jones and P. Sutton for useful discussions. We acknowledge support through the Center for Gravitational Wave Physics, which is funded by NSF under cooperative agreement PHY 01-14375, and through NSF AST0098416 and NASA NAG5-9192.

References
Bonnell, I.A. & Pringle, J.E. 1995, MNRAS, 273, L12.
Davies, M.B., King, A.,Rosswog,S. & Wynn,G. 2002, ApJ, 579, L63.
Detweiler,S. 1980, ApJ, 239, 292.
Eichler, D., Livio,M., Piran,T. & Schramm,D.N. 1989, Nature, 340, 126.
Finn, L.S. & Chernoff,D.F. 1993, Phys. Rev. D, 47, 2198.
Flanagan,E.E. & Hughes, S.A. 1998, Phys. Rev. D, 57, 4535.
Frail,D.A. et al. 2001, ApJ, 562, L55.
Fryer, C.L., Woosley,S.E., Herant,M. & Davies,M.B. 1999, ApJ, 520, 650.
Fryer, C.L., Holz, D.E. & Hughes,S.A. 2002, ApJ, 565, 430.
Janka,H.T., Eberl,T., Ruffert,M. & Fryer,C. 1999, ApJ, 527, L39.
Kippen, R.M., Woods,P.M., Heise,J., in’tZand,J.J.M., Briggs,M.S. & Preece,R.D. 2002, preprint (astro-ph/0203114).
Kobayashi,S & Mészáros,P. 2002, submitted to ApJ, astro-ph/0210211.
Kochanek,C. & Piran, T. 1993, ApJ, 417, L17.
Landau,L.D. & Lifshitz,E.M. 1975, The Classical Theory of Fields, (New York, Pergamon).
MacFadyen,A.I. & Woosley, S.E. 1999, ApJ, 524, 262.
Mészáros, P. 2002, Annu. Rev. Astron. Astrophys., 40,137.
Misner,C.W., Thorne,K.S. & Wheeler,J.A., Gravitation (Freeman, San Francisco 1973).
Nakamura,T. & Fukugita,M. 1989, ApJ, 337, 466.
Rampp,M., Muller,E. & Ruffert,M. 1998, A&A, 332, 969.
Rossi,E., Lazzati,D. & Rees,M.J. 2002, MNRAS, 332, 945.
Paczynski,B. 1991, Acta Astronomica, 41, 257.
Panaitescu,P. & Kumar,P. 2002, ApJ, 571, 779.
Piran,T., Kumar,P., Panaitescu,A. & Piro,L. 2001, ApJ, 560, L167.
Rhoads,J.E. 1999, ApJ, 525, 737.
Ruffert,M., Janka,H.T., Takahashi,K. & Schaefer,G 1997, A&A, 319, 122.
Sari,R., Piran,T.& Halpern,J.P. 1999, ApJ, 519, L17.
Thorne,K.S. 1980, Rev. Mod. Phys., 52, 299.
Yamazaki, R., Ioka,K. & Nakamura,T. 2002, ApJ, 571, L31.
van Paradijs,J., Kouveliotou,C. & Wijers,R.A.M.J. 2000, Annu. Rev. Astron. Astrophys., 38,379.
van Putten, M.H.P.M. 2001 , ApJ Lett., 562, L51.
van Putten, M.H.P.M. 2002 , ApJ Lett., 575, L71.
Woosley,S. 1993, ApJ, 405, 273.
Woosley,S. Zhang,W & Heger,A. 2002 , preprint (astro-ph/0206004).
Zhang,B. & Mészáros,P. 2002, ApJ, 571, 876.