format: text

text:
# Fireball Heated by Neutrinos

## 1 INTRODUCTION

The relativistic fireball (Shemi and Piran, 1990; Rees and Mészáros, 1992; Mészáros and Rees, 1993; Sari and Piran, 1995; Sari et al., 1996) is one of the most promising models of gamma-ray bursts (GRBs). The sudden release of a large quantity of gamma-ray photons into a compact region produces a fireball that is an opaque radiation plasma whose energy is significantly greater than its rest mass. Piran et al. (1993) and Mészáros et al. (1993) investigated the hydrodynamics of fireballs and established basic behaviour in their evolution. Most of the matter and energy is concentrated within a narrow shell. In the radiation-dominated phase, the rest mass and radiation energy density in the fireball behave as $`\rho r^3`$ and $`er^4`$ ($`r`$ is the location of the shell), respectively. The bulk Lorentz factor $`\gamma `$ increases in proportion to $`r`$. After the fireball enters the matter-dominated regime, as long as the shell is optically thick to scattering, the densities behave as $`\rho r^2`$ and $`er^{8/3}`$, and the shell eventually coasts with $`\gamma =\eta E/Mc^2`$, where $`M`$ and $`E`$ are the total baryon rest mass and the total energy of the fireball, respectively. A simple analytical estimate and hydrodynamical simulations indicate that the shell width expands as $`Wr/\eta ^2`$ after the comoving energy density drops below the baryon rest mass density, while the shell width remains constant during the radiation-dominated phase. The actual behaviour of the fireball may be more complex depending on the energy injection process of the central engine. A detailed study of the fireball in various situations will be important to determine the character of the central engine in the future.

The rapid temporal variabilities observed require that the GRB itself must arise from internal shocks within the flow. The shortest time scale is about 1 msec (Bhat et al., 1992). The internal shock is generated by the collision of two shells with different speeds arising from inhomogeneity of the fireball. In this paper we assume that the inhomogeneity in the fireball can be represented by multiple fireballs arising intermittently.

Notwithstanding the very high energy phenomenon ($`10^{52}`$-$`10^{54}`$ ergs for spherical symmetry), the baryon density in the fireball must be extremely small. In order to reproduce GRBs the final Lorentz factors of shells should be $`100`$-$`1000`$. This is the famous baryon contamination problem and still remains unsolved. The central engine of GRBs is still beyond deep mist. The source of GRBs may be one super massive (failed) supernovae (Woosley, 1993; Paczyński, 1998) or may be a merger of two neutron stars or of a neutron star and a black hole (Eichler et al., 1989; Narayan et al., 1992; Mészáros and Rees, 1992a; Katz, 1997; Ruffert and Janka, 1998, 1999).

In these systems the hot core or accretion disk can emit neutrinos and antineutrinos. The neutrino-antineutrino annihilation into electrons and positrons (hereafter neutrino pair-annihilation) is a possible and important candidate to explain the energy source of GRBs (Paczyński, 1990; Mészáros and Rees, 1992b; Janka and Ruffert, 1996; Ruffert at al., 1997; Ruffert and Janka, 1998, 1999; Salmonson et al., 2001). The central engine releases energy in a finite period of time at least longer than $`L/c`$ where $`L`$ is the scale of the central engine. Motivated by the delayed explosion of Type II supernovae, the energy deposition rate due to neutrino pair-annihilation above the neutrinosphere has been calculated (Goodman et al., 1987; Cooperstein et al., 1987; Berezinsky and Prilutsky, 1987). The energy deposition rate is proportional to $`r^8`$ for large $`r`$, and almost all deposition occurs near the neutrinosphere. In hydrodynamic simulations of neutron star mergers, Ruffert and Janka (1999) showed that neutrino pair-annihilation deposits energy in the vicinity of the torus at a rate of $`(3`$-$`5)\times 10^{50}`$ ergs $`\mathrm{s}^1`$. In the failed supernovae model MacFadyen and Woosley (1999) predict that the total energy deposited along the rotational axes of Kerr black holes is $`(1`$-$`14)\times 10^{51}`$ ergs.

Most of the neutrinos’ energy ($`99`$ %) will escape to infinity without annihilation. Since GRBs are very high energy phenomena, it is favorable that the energy of neutrinos is transformed into the fireball much more efficiently. Asano and Fukuyama (2000, 2001) showed that the gravitational effect does not change the energy deposition rate dramatically. There is also another possible process to transfer the neutrinos’ energy to the fireball. In the fireball, high-energy photons produce a large number of electron-positron pairs at the initial stage. The large optical depth due to electron-positron pairs leads to thermal equilibrium (Shemi and Piran, 1990). The neutrino scattering with the thermalized electron-positron pairs (hereafter neutrino scattering) can inject energy into the fireball. This possibility was already pointed out by Woosley (1993).

In this paper we discuss the evolution of the fireball heated by neutrinos and the formation of shells using a simple spherical model. We do not deal with a specific model of the central engine. We discuss only the general evolution of the fireball with energy injection in a finite period of time and the effect of neutrino scattering. While baryonic matter has been considered to be injected together with radiation energy in general, we assume here that the matter is provided around the central engine before neutrinos start to be emitted. In §2 we review the formulae for the energy and momentum deposition rates via neutrino pair-annihilation and derive the deposition rates via neutrino scattering. In §3 we numerically simulate the fireball heated by neutrinos to investigate the evolution of the fireball and the deposition rates via neutrino scattering. Finally, §4 is appropriated for the summary.

## 2 ENERGY AND MOMENTUM DEPOSITION RATE

### 2.1 Pair Annihilation

In this subsection we review the formulae for the energy and momentum deposition rates due to neutrino pair-annihilation. For simplicity, we assume that neutrinos are emitted isotropically from the neutrinosphere of radius $`R_\nu `$. The neutrino number density in the phase space, $`f_\nu d^3p_\nu =n_\nu (\epsilon _\nu )\epsilon _\nu ^2d\epsilon _\nu d\mathrm{\Omega }_\nu `$, is assumed to be conserved along the neutrinos’ trajectory. This is a good approximation as long as the energy efficiency of neutrino pair-annihilation is very small. Defining $`X[1(R_\nu /r)^2]^{1/2}`$, Goodman et al. (1987) derived the energy deposition rate due to neutrino pair-annihilation as

$`{\displaystyle \frac{dE}{cdtdV}}|_{\mathrm{pair}}`$ $`=`$ $`{\displaystyle \frac{K_\mathrm{p}G_\mathrm{F}^2}{12\pi ^2c^2R_\nu ^4}}\dot{N_\nu }\dot{N_{\overline{\nu }}}(\epsilon _\nu \epsilon _{\overline{\nu }}^2+\epsilon _\nu ^2\epsilon _{\overline{\nu }})F_E(r),`$ (1)
where

$`F_E(r)=(1X)^4(X^2+4X+5),`$ (2)
and the Fermi constant $`G_\mathrm{F}^2=5.29\times 10^{44}\mathrm{cm}^2\mathrm{MeV}^2`$. The total emmision rate of neutrinos, $`\dot{N_\nu }`$, and the average energy of neutrinos, $`\epsilon _\nu `$, are written as

$`\dot{N_\nu }\epsilon _\nu ^k=(4\pi R_\nu ^2)(\pi c){\displaystyle \epsilon _\nu ^{2+k}n_\nu (\epsilon _\nu )𝑑\epsilon _\nu }.`$ (3)
If we simplify the energy distribution of neutrinos by using the effective temperature, $`T_{\mathrm{eff}}`$, we can express the total emmision rate as

$`\dot{N_\nu }\epsilon _\nu ={\displaystyle \frac{7\pi ^6R_\nu ^2}{30h^3c^2}}T_{\mathrm{eff}}^4,\dot{N_\nu }\epsilon _\nu ^2={\displaystyle \frac{90\pi ^2\zeta (5)R_\nu ^2}{h^3c^2}}T_{\mathrm{eff}}^5,`$ (4)
where $`\zeta (n)`$ is the Riemann zeta function. The constant $`K_\mathrm{p}`$ is written by the Weinberg angle $`\mathrm{sin}^2\theta _\mathrm{W}=0.23`$. To be exact we should sum up contributions of all species of neutrino. Each species may be emitted at different emission rates. When $`\nu _\mathrm{e}`$ and $`\overline{\nu }_\mathrm{e}`$ are dominant components, $`K_\mathrm{p}=0.124`$. If $`\dot{N_\nu }\epsilon _\nu `$ and $`\dot{N_\nu }\epsilon _\nu ^2`$ of all the six species of neutrino are the same, we can include the contributions of $`\nu _\mu `$ and $`\nu _\tau `$ by setting $`K_\mathrm{p}=0.178`$. The fraction and mean energy of each species of neutrino depend on the model of the central engine (see, e.g., Ruffert at al. (1997)). In most cases the effective $`K_\mathrm{p}`$ may be between $`0.124`$ and $`0.178`$ unless the emission rate of $`\nu _\mu `$ or $`\nu _\tau `$ is much larger than the rate of $`\nu _\mathrm{e}`$. For $`rR_\nu `$, the energy deposition rate behaves as $`T_{\mathrm{eff}}^9r^8`$. Most of the energy is deposited near to the surface of the neutrinosphere. We also obtain the momentum deposition rate in a similar way as

$`{\displaystyle \frac{dP}{dtdV}}|_{\mathrm{pair}}`$ $`=`$ $`{\displaystyle \frac{K_\mathrm{p}G_\mathrm{F}^2}{12\pi ^2c^2R_\nu ^4}}\dot{N_\nu }\dot{N_{\overline{\nu }}}(\epsilon _\nu \epsilon _{\overline{\nu }}^2+\epsilon _\nu ^2\epsilon _{\overline{\nu }})F_P(r),`$ (5)
where

$`F_P(r)=(1X)^4(1+X)(3X^2+9X+8)/4.`$ (6)
### 2.2 Neutrino Scattering

In the fireball high-energy photons produce a large number of electron-positron pairs. The large optical depth due to electron-positron pairs leads to thermal equilibrium (Shemi and Piran, 1990). Neutrino scattering with thermalized electron-positron pairs can transfer the neutrinos’ energy to the fireball. We derive the energy and momentum deposition rates through neutrino scattering in this subsection. The electron-positron pairs are assumed to be relativistic $`\epsilon _\mathrm{e}m_\mathrm{e}c^2`$. For the reaction $`\nu _\mathrm{e}+\mathrm{e}^{}\nu _\mathrm{e}+\mathrm{e}^{}`$, the differential cross-section is given by

$`{\displaystyle \frac{d\sigma }{dy}}={\displaystyle \frac{G_\mathrm{F}^2(p_{i,\nu }p_{i,\mathrm{e}})}{2\pi }}\left[(2\mathrm{sin}^2\theta _\mathrm{W}+1)^2+4\mathrm{sin}^4\theta _\mathrm{W}(1y)^2\right],`$ (7)
where $`yp_{i,\mathrm{e}}(p_{i,\nu }p_{f,\mathrm{e}})/(p_{i,\nu }p_{i,\mathrm{e}})`$ is Lorentz invariant, and $`(p_{i,\nu }p_{i,\mathrm{e}})`$ expresses the inner product of the 4-momenta of a neutrino and an electron. The subscripts “$`i`$” and “$`f`$” denote the initial and final states, respectively. We can write the cross-sections for the other five species of neutrino in a similar way (see, e.g., Kim and Pevsner (1993)).

In the comoving frame of the relativistic fireball, the energy of the neutrinos redshifts as $`\epsilon _\nu ^{}=\epsilon _\nu /\gamma (1+\beta \mathrm{cos}\theta _\nu ^{})`$, where $`\gamma =1/(1\beta ^2)^{1/2}`$ is the Lorentz factor of the fireball and $`\theta _\nu ^{}`$ is the incident angle of neutrinos to the $`r`$axis (see Fig. 1). Hereafter values with a prime are measured in the comoving frame. We have two scattering angles for a given value of $`y`$. The energy a neutrino loses per one scattering depends on the final scattering angle. The average energy an electron/positron obtains for a given value of $`y`$ becomes $`\overline{\mathrm{\Delta }\epsilon ^{}}=y(\epsilon _\nu ^{}\epsilon _\mathrm{e}^{})`$. The average momentum component an electron/positron obtains, which is parallel to the neutrino’s momentum, becomes $`c\overline{\mathrm{\Delta }p_{}^{}}=y(\epsilon _\nu ^{}\epsilon _\mathrm{e}^{}\mathrm{cos}\widehat{\theta ^{}})`$, where $`\widehat{\theta ^{}}`$ is the incidence angle between the two particles. Here we assume that the random motion of electron-positron pairs is isotropic in the fireball frame. When a neutrino interacts with electrons distributed isotropically, the mean vertical momentum the neutrino loses per one scattering is zero. As we have assumed in §2.1, if $`n_\nu (\epsilon _\nu )`$ of all the six species of neutrino are the same, the energy a scattered electron (positron) gains and the total cross section are averaged as

$`{\displaystyle \underset{s}{}}n_{\nu s}(\epsilon _\nu ^{})\overline{\sigma _s\mathrm{\Delta }\epsilon ^{}}`$ $`=`$ $`n_\nu (\epsilon _\nu ^{}){\displaystyle _0^1}{\displaystyle \underset{s}{}}{\displaystyle \frac{d\sigma _s}{dy}}y(\epsilon _\nu ^{}\epsilon _\mathrm{e}^{})dy`$ (8)
$`=`$ $`{\displaystyle \frac{7}{12}}n_\nu (\epsilon _\nu ^{})K_\mathrm{s}G_\mathrm{F}^2(p_{i,\nu }^{}p_{i,\mathrm{e}}^{})(\epsilon _\nu ^{}\epsilon _\mathrm{e}^{}),`$ (9)
where the summation is over all species of neutrino, and $`K_\mathrm{s}(24\mathrm{sin}^4\theta _\mathrm{W}4\mathrm{sin}^2\theta _\mathrm{W}+3)/2\pi =0.533`$. When $`\nu _\mathrm{e}`$ and $`\overline{\nu }_\mathrm{e}`$ are dominant components, $`K_\mathrm{s}`$ should be altered as $`(16\mathrm{sin}^4\theta _\mathrm{W}+8\mathrm{sin}^2\theta _\mathrm{W}+2)/4\pi =0.373`$. To be more precise we have to estimate the deposition rate corresponding to the variation of the distributions of the six species of neutrino. We include the variation by adjusting $`K_\mathrm{s}`$ effectively as we have done in §2.1.

The numbers of neutrinos and antineutrinos are assumed to be the same. Then the reaction rates of electrons and positrons are the same because of the symmetric property. Since the energy distribution of the positrons is the same as the distribution of the electrons in thermal equilibrium, we include the contribution of positrons by doubling the energy distribution of electrons $`n_\mathrm{e}(\epsilon _\mathrm{e}^{})`$. The energy deposition rate in the comoving frame is written as

$`{\displaystyle \frac{dE}{cdtdV}}|_{\mathrm{scat}}^{}`$ $`=`$ $`{\displaystyle 𝑑\epsilon _\nu ^{}𝑑\epsilon _\mathrm{e}^{}𝑑\mathrm{\Omega }_\nu ^{}𝑑\mathrm{\Omega }_\mathrm{e}^{}\underset{s}{}\epsilon _\nu ^2\epsilon _\mathrm{e}^2n_{\nu s}(\epsilon _\nu ^{})n_\mathrm{e}(\epsilon _\mathrm{e}^{})\overline{\sigma _s\mathrm{\Delta }\epsilon ^{}}c(1\mathrm{cos}\widehat{\theta ^{}})}`$
$`=`$ $`{\displaystyle \frac{14cK_\mathrm{s}G_\mathrm{F}^2\pi ^2}{3}}{\displaystyle 𝑑\epsilon _\nu ^{}𝑑\epsilon _\mathrm{e}^{}n_\nu (\epsilon _\nu ^{})n_\mathrm{e}(\epsilon _\mathrm{e}^{})\epsilon _\nu ^3\epsilon _\mathrm{e}^3(\epsilon _\nu ^{}\epsilon _\mathrm{e}^{})}`$
$`\times `$ $`{\displaystyle _0^{\theta _{\mathrm{max}}^{}}}𝑑\theta _\nu ^{}\mathrm{sin}\theta _\nu ^{}{\displaystyle _0^\pi }𝑑\widehat{\theta ^{}}\mathrm{sin}\widehat{\theta ^{}}(1\mathrm{cos}\widehat{\theta ^{}})^2,`$ (11)
where $`\mathrm{cos}\theta _{\mathrm{max}}^{}(X\beta )/(1\beta X)`$, including the redshift effect. Assuming the electron-positron pairs are in thermal equilibrium and their temperature $`T_\mathrm{e}m_\mathrm{e}c^2`$, we get

$`{\displaystyle \frac{dE}{cdtdV}}^{}|_{\mathrm{scat}}=f(T_\mathrm{e}(r))G_E^{}(T_\mathrm{e},\beta ;r),`$ (12)
where

$`f(T_\mathrm{e}(r))`$ $`=`$ $`{\displaystyle \frac{7}{R_\nu ^2c}}{\displaystyle \frac{K_\mathrm{s}G_\mathrm{F}^2}{(hc)^3}}T_\mathrm{e}^4(r)\dot{N_\nu },`$ (13)
and

$`G_E^{}(T_\mathrm{e},\beta ;r)`$ $`=`$ $`\gamma ^2(1X)\{\gamma {\displaystyle \frac{7\pi ^4}{135}}ϵ_\nu ^2[1{\displaystyle \frac{3}{2}}\beta (1+X)+`$ (14)
$`\beta ^2(1+X+X^2){\displaystyle \frac{\beta ^3}{4}}(1+X)(1+X^2)]`$
$`20\zeta (5)ϵ_\nu T_\mathrm{e}(r)[1\beta (1+X)+{\displaystyle \frac{\beta ^2}{3}}(1+X+X^2)]\}.`$
The temperature $`T_\mathrm{e}(r)`$ is measured in the comoving frame. The momentum deposition rate can be obtained in a similar fashion, using the parallel component of momentum to the $`r`$-axis, $`\mathrm{cos}\theta _\nu ^{}\overline{\mathrm{\Delta }p_{}^{}}`$. Carrying out the Lorentz transformation, we obtain the energy deposition rate in the observer frame as $`dE/cdtdV|_{\mathrm{scat}}=f(T_\mathrm{e}(r))G_E(T_\mathrm{e},\beta ;r)`$, where

$`G_E(T_\mathrm{e},\beta ;r)=\gamma \{\gamma {\displaystyle \frac{7\pi ^4}{135}}ϵ_\nu ^2[1\beta (1+X)+{\displaystyle \frac{\beta ^2}{3}}(1+X+X^2)]`$
$`20\zeta (5)ϵ_\nu T_\mathrm{e}(r)[{\displaystyle \frac{1}{2}}(3\gamma ^21)+{\displaystyle \frac{\beta }{4}}(16\gamma ^2)(1+X)`$
$`+{\displaystyle \frac{1}{2}}(\gamma ^21)(1+X+X^2)]\}(1X).`$ (15)
In the case $`T_\mathrm{e}(r)m_\mathrm{e}c^2`$ these formulae overestimate the rates. These equations give the upper limit of the energy and momentum deposition rates. For $`\beta 0`$ and $`rR_\nu `$, introducing the effective temperature of the neutrinos, the energy deposition rate behaves as $`(T_{\mathrm{eff}}^5T_\mathrm{e}^4(r)T_\mathrm{e}^5(r)T_{\mathrm{eff}}^4)r^2`$. The deposition rate depends on the temperature of the fireball. When $`r=R_\nu `$ and $`\beta =0`$, the constant factor, except for the temperature dependence, in the deposition rate via neutrino scattering is larger than the factor via pair-annihilation by a factor of $`5`$. If $`T_\mathrm{e}T_{\mathrm{eff}}`$, the energy injection rate for neutrino scattering becomes comparable to the rate for neutrino pair-annihilation. However, if $`T_\mathrm{e}T_{\mathrm{eff}}`$ the scattering effect cools down the fireball inversely. The deposition rate via pair-annihilation behaves as $`r^8`$. If $`T_\mathrm{e}^4(r)`$ has a flatter distribution than $`r^6`$, much more energy could be injected by scattering than energy injected by annihilation at a great distance.

The momentum deposition rate in the observer frame is obtained by $`dP/dtdV|_{\mathrm{scat}}=f(T_\mathrm{e}(r))G_P(T_\mathrm{e},\beta ;r)`$, where

$`G_P(T_\mathrm{e},\beta ;r)=\gamma \{\gamma {\displaystyle \frac{7\pi ^4}{135}}ϵ_\nu ^2[{\displaystyle \frac{1}{2}}(1+X){\displaystyle \frac{2}{3}}\beta (1+X+X^2)+{\displaystyle \frac{\beta ^2}{4}}(1+X)(1+X^2)]`$
$`+10\zeta (5)ϵ_\nu T_\mathrm{e}(r)[(3\gamma ^2{\displaystyle \frac{5}{2}})(1+X)3\beta \gamma ^2{\displaystyle \frac{\beta \gamma ^2}{3}}(1+2\beta ^2)(1+X+X^2)]\}(1X).`$ (16)
For $`\beta 0`$, neutrinos accelerate the fireball regardless of $`T_\mathrm{e}`$. Neutrino scattering can decelerate the flow for large $`\gamma `$ like Compton dragging.

There is a possibility that the temperature of the fireball could run away via neutrino scattering (Woosley, 1993) and reach $`T_\mathrm{e}T_{\mathrm{eff}}`$. In addition the total energy could be enhanced by the energy deposition at large $`r`$ if $`T_\mathrm{e}(r)`$ distributes flatter than $`r^6`$. However, it seems difficult to lead $`T_\mathrm{e}`$ to the runaway, because the dependence $`T_\mathrm{e}^4`$ makes the deposition rate due to neutrino scattering much weaker than the rate due to neutrino pair-annihilation giving a slightly smaller $`T_\mathrm{e}`$ than $`T_{\mathrm{eff}}`$. To reach higher temperatures the expansion speed of the plasma should be slow enough during the energy injection such that the heating effect via neutrinos is dominant rather than adiabatic cooling. In most cases we need to simulate the evolution of the fireball numerically in order to estimate the conclusive energy deposition via neutrinos.

## 3 NUMERICAL CALCULATION

In this section we numerically simulate the fireball heated by neutrinos. Neutrinos are assumed to be emitted isotropically from the neutrinosphere during a finite period of time. As Piran et al. (1993) assumed in their calculation, we assume the fireball is optically thick. The radiation and matter at each radius behave like a single fluid moving with the same velocity. Since the radiation pressure dominates, the pressure $`p`$ and the energy density $`e`$ of radiation and relativistic electron-positron pairs are related by $`p=e/3`$. The relativistic conservation equations of baryon number, energy and momentum are written by

$`{\displaystyle \frac{1}{c}}{\displaystyle \frac{}{t}}\left(\rho \gamma \right)={\displaystyle \frac{1}{r^2}}{\displaystyle \frac{}{r}}\left(r^2\rho u\right),`$ (17)
$`{\displaystyle \frac{1}{c}}{\displaystyle \frac{}{t}}\left({\displaystyle \frac{4}{3}}e\gamma ^2{\displaystyle \frac{e}{3}}+\rho \gamma ^2\right)={\displaystyle \frac{1}{r^2}}{\displaystyle \frac{}{r}}\left[r^2\gamma u\left({\displaystyle \frac{4}{3}}e+\rho \right)\right]+{\displaystyle \frac{dE}{cdtdV}},`$ (18)
$`{\displaystyle \frac{1}{c}}{\displaystyle \frac{}{t}}\left[u\gamma \left({\displaystyle \frac{4}{3}}e+\rho \right)\right]={\displaystyle \frac{1}{r^2}}{\displaystyle \frac{}{r}}\left[r^2u^2\left({\displaystyle \frac{4}{3}}e+\rho \right)\right]{\displaystyle \frac{1}{3}}{\displaystyle \frac{e}{r}}+{\displaystyle \frac{dP}{dtdV}},`$ (19)
where $`u\gamma \beta `$. The rest mass energy density of baryons $`\rho `$ and the non-baryonic energy density $`e`$ are measured in the comoving frame.

We use the relativistic “Roe-solver” method in the computation of one-dimensional relativistic flow. The Riemann problem is approximated numerically by means of a specially chosen linearized form, requiring less computational effort. This Eulerian code was developed for the relativistic case by Mellema et al. (1991). However, our computation is not so complicated as the computation of Mellema et al. (1991), because we assume spherical symmetry and neglect the gravitational effect. Our code is relatively simple and has passed several standard tests. In addition we have checked our numerical code by reproducing the results of the simulation in Piran et al. (1993). We adopt the spatial and time resolution as $`\mathrm{\Delta }r=R_\nu /3000`$ and $`c\mathrm{\Delta }t=\mathrm{\Delta }r/10`$, respectively. In order to discuss the evolution of shell width, however, we need to simulate with a very high spatial resolution. The expansion/contraction rate of shell width may be limited as $`|dW/dr|1/\gamma ^2`$. It will take very long time to follow the change of shell width accurately for large $`\gamma `$ because of the high resolution. In our simulation, judging from test calculations, the expansion/contraction rate is not so reliable for $`\gamma 10`$. However the rate of change is small enough for large $`\gamma `$, in our calculation, such that we can follow the approximate behaviour of the fireball.

First we investigate the evolution of the fireball in an ideal situation to produce a relativistic flow. Cold baryonic matter is provided around the neutrinosphere before the simulation starts. Then we study the effect of neutrino scattering on the energy deposition rate.

### 3.1 Evolution of the Fireball

In order to produce a relativistic flow, we assume an ideal situation. The initial profile of the matter that we present is compact as $`\rho (r,t=0)(R_0^8+r^8)^1`$ for $`rR_\nu `$. The matter is initially cold and static as $`e(r,0)=0`$ and $`\gamma (r,0)=1`$. In this simulation the energy is injected at the constant rate, $`\dot{E}`$, during $`0t(rR_\nu )/c+t_{\mathrm{dur}}`$, while in previous studies they had provided radiation energy before the simulations started (Piran et al., 1993; Mészáros et al., 1993). Spatial distribution of $`\dot{E}`$ is given by equation (2). We neglect the effect of neutrino scattering in this subsection. We set the initial core radius of the matter as $`R_0=ct_{\mathrm{dur}}`$ hereafter. The rest mass energy is normalized as $`4\pi \eta 𝑑rr^2\rho =\dot{E}t_{\mathrm{dur}}`$, where $`\eta `$ is the ratio of the radiation energy to the rest mass energy. The mean final Lorentz factor of the matter shell should be $`\eta +1`$. Thus, the physical conditions in this simulation are determined by the two parameters $`\eta `$ and $`\tau _{\mathrm{dur}}ct_{\mathrm{dur}}/R_\nu `$.

First of all we show results for the case of $`\eta =60`$ and $`\tau _{\mathrm{dur}}=3`$ as an example. If we introduce physical units as $`R_\nu =10`$ km and $`T_{\mathrm{eff}}=10`$ MeV, the initial matter density at $`r=R_\nu `$ corresponds to $`2.4\times 10^4`$ g $`\mathrm{cm}^3`$ in this case. Figure 2-4 are the profiles of the fireball at different times. At $`\tau ct/R_\nu =3`$ (Fig. 2), namely, at the end of the energy injection, the matter is swept out by the radiation pressure. The information about the initial configuration of the matter is almost lost. The approximation $`\rho e`$, which has been assumed conventionally in discussions about the initial evolution of the fireball, already failed around the matter shell at this time. The matter shell is compressed by the radiation pressure. The FWHM of the matter shell in the observer frame (the profile of $`\rho \gamma `$) is about $`0.1R_\nu `$ and remains nearly constant within the numerical error limit hereafter. The behavior of the fluid around the neutrinosphere agrees with the stationary solution of pure radiation flow with energy injection. We have confirmed that $`e`$ and $`\gamma `$ behind the matter shell in Figure 2 are very close to the numerical solutions obtained by setting the left-hand side of equations (17)-(19) to zero and $`\rho =0`$. Figure 3 is the fireball at $`\tau =6`$. The energy injection has already finished. The pure radiation flow, which had been deposited later, is chasing the matter shell. Around the matter shell the forward and reverse shocks are formed. The reverse shock is radiation dominant, while the forward one is matter dominant. In Figure 4 we display the fireball at $`\tau =30`$. The Lorentz factor becomes a maximum behind the shell, where the radiation is dominant. The radiation flow pushes the matter from behind. The matter shell itself accelerates more slowly than the prediction by the simple analysis, $`\gamma =(r/R_\nu )`$. The behaviour of $`u`$ at the point where $`\rho \gamma `$ has the maximum value is approximated as $`u_\mathrm{M}r^{0.5}`$ at this moment. The value $`u_\mathrm{M}`$ does not accelerate obeying a simple power law. The power law index declines as $`u_\mathrm{M}r^{0.4}`$, $`r^{0.3}`$,… with increasing $`r`$. We need a very long time to reach $`u_\mathrm{M}=\eta `$ in our computation. The maximum values of $`\rho `$ and $`e`$ do not obey simple power-law behaviour.

In Figure 5 we plot $`u_\mathrm{M}`$ for different values of $`\tau _{\mathrm{dur}}`$ and $`\eta =60`$. The acceleration time is roughly in proportion to the duration time of the energy release. When the shell reaches $`r=\eta R_\nu `$, $`\gamma `$ becomes a few tens of a percent of $`\eta `$. The maximum Lorentz factor in the whole system increases in proportion to $`r`$. Extrapolating $`u_\mathrm{M}`$ from Figure 5, we guess $`u_\mathrm{M}`$ may reach $`\eta `$ at $`r\eta ^2\tau _{\mathrm{dur}}R_\nu `$. At the end of the energy injection the distance from $`R_\nu `$ to the head of the shell is about $`\tau _{\mathrm{dur}}R_\nu `$. Therefore, the energy injected at the end will catch up with the matter shell of $`\gamma \eta `$ at $`r\eta ^2\tau _{\mathrm{dur}}R_\nu `$. This is the reason $`u_\mathrm{M}`$ reaches $`\eta `$ at $`r\eta ^2\tau _{\mathrm{dur}}R_\nu `$.

To economize the calculation time we demonstrate the evolution of the fireball for $`\eta =10`$ and $`\tau _{\mathrm{dur}}=1`$. In Figure 6 we plot $`u_\mathrm{M}`$ and FWHM of the profile of $`\rho \gamma `$. Since the initial matter distribution is narrow in comparison with the case $`\tau _{\mathrm{dur}}=3`$, the radiation pressure expands the matter component initially. Then the matter component is compressed. Figure 7 shows that the rear part of the matter shell is faster than the head part at $`\tau =12`$ that implies the shell is compressed. The shell begins to expand again from $`\tau 30`$. In this period the rear part of the matter shell is slightly slower than the head part, as is shown in Figure 8. The expansion rate is $`dW/dr10^3<\eta ^2`$, which is close to the spatial resolution $`\mathrm{\Delta }r/R_\nu `$ in our simulation. This value may not be accurate because of the size of the numerical error. However, the difference of the Lorentz factor in the shell is rather small as $`\delta \gamma \gamma `$. At $`\tau =105`$ the maximum and minimum $`u`$ within FWHM are $`10.28`$ and $`9.89`$, respectively. The maximum relative velocity within FWHM in the observer frame is $`3.8\times 10^4c`$ that suggests $`dW/dr4\times 10^4`$, which is smaller than our error limit in the computation. Thus, we can expect a smaller expansion rate than $`\eta ^2`$ by a factor of $`10`$ or $`10^2`$. As we have assumed, $`u_\mathrm{M}`$ reaches $`\eta `$ at $`r/R_\nu \tau \eta ^2=100`$. We can expect that the behaviour of the fireball for large $`\eta `$ would be similar to the above results.

We refer the simulation in Piran et al. (1993) to compare with the case in which radiation energy has been given before the simulation starts. In their simulation the initial profiles are given as $`e(r,0)(R_0^8+r^8)^1`$ and $`\rho (r,0)=e/\eta `$. According to Figure 2 in Piran et al. (1993) the matter shell expands even in the acceleration regime. At $`\tau 100`$ in our notation the shell width grows to a few times the size of $`R_0`$. The Lorentz factor in the matter shell increasingly distributes with radius and their difference in the shell is rather large as $`\delta \gamma \gamma `$.

The results of this subsection can be summarized as follows: In the acceleration regime the shell is pushed from behind by the radiation pressure. The expansion of the shell is suppressed in this regime. The acceleration declines as the shell evolves. The Lorentz factor of the shell reaches $`\eta `$ at $`r\eta ^2\tau _{\mathrm{dur}}R_\nu `$. After the fireball enters the coasting regime, the shell begins to expand. The velocity distribution in the shell is very flat. We need a more accurate numerical simulation and longer calculation time to estimate the expansion rate. There is a possibility that the expansion rate is much less than $`\eta ^2`$. The simulation in this subsection is merely one example of the evolution of the fireball being affected by the energy deposition due to neutrinos. However, when radiation fluid is deposited inside the initial matter distribution during a finite period of time, the fireball may evolve in a similar manner regardless of the details of the spatial distribution of $`\dot{E}`$ and initial matter.

### 3.2 Effect of Neutrino Scattering

In this subsection we simulate the fireball including the effect of neutrino scattering. Our purpose is to investigate the enhancement of the energy deposition by the scattering effect. The energy deposition by neutrinos including the scattering effect is a nonlinear process. The rate of neutrino scattering depends on the temperature of the fireball. As is shown in the former subsection, the temperature around the neutrinosphere can be expressed by the stationary solution of the pure radiation flow as long as the initial matter density is sufficiently low. It is easily found that the maximum temperature, which is obtained from the stationary solution, is much lower than the effective temperature of the neutrino. In an ideal situation for producing the relativistic shell as in §3.1, the effect of neutrino scattering is negligible.

However the environment is generally very dense around candidates of the central engine: the central region of massive stars, merging neutron stars, and so on. In a very dense environment the matter cannot be swept out promptly. The energy injected by neutrinos may stand near to the neutrinosphere, and the fireball temperature becomes higher than the case we considered in the former subsection. We investigate the effect of neutrino scattering for a dense environment. In order to consider the most effective case on neutrino scattering, we make the matter density constant as $`10^9\mathrm{g}/\mathrm{cm}^3`$ and set $`R_\nu =10`$ km and $`T_{\mathrm{eff}}=15`$ MeV. Similarly, as in the former subsection, the energy deposition due to pair-annihilation is given by equation (1) with $`K_\mathrm{p}=0.178`$. We add the contribution of the scattering effect discussed in §2.2 to equations (17)-(19) with $`K_\mathrm{s}=0.533`$. The temperature of the fireball is obtained from $`T_\mathrm{e}^4=60(\mathrm{}c)^3e/(11\pi ^2)`$. Since our purpose is to investigate the effect of scattering between neutrinos and electron-positron pairs in the extreme environment, we ignore the interaction between baryons and neutrinos. As a matter of fact, for $`T_{\mathrm{eff}}=15`$ MeV, the energy deposition rate due to the neutrino-baryon interaction ($`T_{\mathrm{eff}}^6\rho `$ for $`\beta =0`$) becomes dominant rather than the rate due to pair-annihilation at $`r=R_\nu `$, in a case in which the initial matter density is larger than $`7\times 10^9\mathrm{g}/\mathrm{cm}^3`$. As long as $`\rho `$ is constant, the total energy deposition via the neutrino-baryon interaction is larger than the deposition via pair-annihilation.

In Figure 9, we display the fireball at $`\tau =1`$. We have normalized the length scale and energy density by 10 km and $`(10^5`$ g $`\mathrm{cm}^3)c^2`$, respectively. As we have assumed above, the expansion speed is slow and the radiation energy stands inside the shell. The distribution of $`e`$ inside the matter shell is nearly flat. The typical temperature $`T_\mathrm{e}`$ inside the matter shell, where the radiation is dominant, is about 5 MeV, smaller than the effective temperature of neutrinos. At $`r=R_\nu `$ the scattering effect enhances the deposition rate by only $`5`$ %, while the total enhancement is $`20`$ % at $`\tau =10`$ (see Fig. 10). This implies that a considerable amount of energy is deposited in a wide range of $`r`$ by scattering.

The temperature $`T_\mathrm{e}(R_\nu )`$ increases after the start of the energy deposition (see Fig. 11). However, the huge energy injection ($`T_{\mathrm{eff}}=15`$ MeV) sweeps out the matter immediately. Then adiabatic cooling due to the expansion becomes stronger than the neutrino heating. The temperature $`T_\mathrm{e}`$ around the neutrinosphere declines and reaches the value given by the stationary solution.

If we enlarge the dimension of the region as $`R_\nu =100`$ km, the energy deposition rate ($`R_\nu ^3`$) and the energy efficiency ($`R_\nu `$) due to pair-annihilation increase. The dynamical time scale becomes larger by a factor of 10, while the energy injection rate per volume remains constant for the given value of $`r/R_\nu `$. Therefore, the injected energy per baryon during a dynamical time scale increases, which leads to high temperatures. In our case the fireball temperature $`T_\mathrm{e}`$ increases to $`7`$ MeV (Fig. 11). The energy deposition is enhanced by 60% because of scattering at $`\tau =10`$. The qualitative behaviour in this case is almost the same as in the former case.

On the other hand, if the effective temperature decreases to $`T_{\mathrm{eff}}=5`$ MeV, the smaller energy injection cannot promptly sweep out the matter. The adiabatic cooling effect is weaker than the heating effect, initially. The temperature $`T_\mathrm{e}(R_\nu )`$ increases with time, but remains much smaller than $`T_{\mathrm{eff}}`$ (Fig. 11) in comparison with the case $`T_{\mathrm{eff}}=15`$ MeV. Then the increasing rate of $`T_\mathrm{e}(R_\nu )`$ is saturated. To make matters worse, the energy deposition due to the neutrino-baryon interaction is practically dominant for lower $`T_\mathrm{e}`$. Neutrino pair-annihilation and scattering play secondary roles. Such a situation is not adequate for our purpose. Reduction of the initial matter density leads to a prompt sweeping out of matter. Thus, it is very difficult to make $`T_\mathrm{e}`$ large.

The temperature of the fireball does not run away in our cases and the nonlinear effect is not so strong. The energy deposition due to scattering is not so important even under the extreme condition we adopted.

## 4 SUMMARY

We have simulated the fireball in which radiation energy is deposited in a narrow region ($`\dot{E}r^8`$) during a finite period of time. The baryonic matter is initially prepared around the neutrinosphere. In the acceleration regime the matter shell is pushed from behind by the radiation pressure. The Lorentz factor of the shell reaches $`\eta `$ at $`r\eta ^2\tau _{\mathrm{dur}}R_\nu `$. Mészáros and Rees (2000) estimated that wind becomes optically thin at $`10^6R_\nu `$ with typical parameters. If $`\tau _{\mathrm{dur}}1`$, the shell can be optically thin before $`\gamma `$ reaches $`\eta `$ even if $`\eta `$ is smaller than the critical value $`\eta _{}`$ (see Mészáros and Rees (2000)). The maximum values of $`\rho `$ and $`e`$ do not obey a simple power-law behavior. After the fireball enters the coasting regime, the shell may begin to expand. The velocity distribution in the shell is rather flat in comparison with the case in which the radiation energy has been given before the simulation starts. Although we need a more accurate numerical simulation to estimate the expansion rate quantitatively, $`dW/dr`$ may be smaller than $`\eta ^2`$.

The short time scale of the temporal profile of GRBs requires that the shell width should remain thin at the collision radius. If the shell width expands as $`W=r/\eta ^2`$, the width becomes $`10^{10}`$ cm for $`\eta =100`$ at the typical collision radius $`10^{14}`$ cm. Then the observed emission time scale of one collision becomes $`1`$ sec, which is longer than the observed timescale of the pulse rise. In addition, if the shell expands, the magnetic field in the shell may weaken. It is very difficult to reproduce hard spectra of observed GRBs in this case (Guetta et al., 2001). Therefore, it is favorable that the shell width is constant. The actual shell expansion depends on the Lorentz factor distribution across the shell, which is determined by the energy ejection process around the central engine. Although the actual activity of the central engine is more complex than our case, our simulation shows a possibility that the expansion rate $`dW/dr`$ is smaller than $`\eta ^2`$. Our simulation is one simple example of the evolution of the fireball. If baryonic matter is also injected with radiation energy, the behaviour of the fireball may be close to the behaviour in Piran et al. (1993). In our case baryonic matter should be injected during the quiescent period of the central engine. The matter could be injected from the side if the central engine has a jetlike structure. But we do not know if this is realistic or not. In the future, we will have to examine whether the shell width remains narrow in some kind of model of the central engine. Future observations may constrain shell’s structure and give a clue to determining the central engine.

As for the effect of neutrino scattering, the temperature of the fireball does not run away in most cases that we assume, and the nonlinear effect is not so strong. Although the scattering effect can enhance the energy deposition rate by a few tens of a percent, the effect on the deposition rate is not so important even for the extreme conditions we assumed. When other processes of energy injection, like magnetic reconnection, etc., provide a large number of pairs, the effect of neutrino scattering can be larger. In any case a very dense environment is required for raising the fireball temperature high enough in a wide range of $`r`$. Under the dense environment the baryon contamination is crucial, while the higher temperature increases the energy deposition due to the scattering effect. Since the energy enhancement due to the scattering effect is not significant, it is difficult to make $`\gamma `$ large enough in the situation where the effect of neutrino scattering becomes important.

We thank Dr. A. Iida for his useful advice concerning the hydrodynamical simulation. The authors are supported by a Research Fellowship of the Japan Society for the Promotion of Science.


text:
# Oceanic coastline and super-universality of percolation clusters.

## Abstract

New fractal subset of a rough surface, the “oceanic coastline”, is defined. For random Gaussian surfaces with negative Hurst exponent $`H<0`$, “oceanic coastlines” are mapped to the percolation clusters of the (correlated) percolation problem. In the case of rough self-affine surfaces ($`H0`$), the fractal dimension of the “oceanic coastline” $`d_c`$ is calculated numerically as a function of the roughness exponent $`H`$ (using a novel technique of minimizing finite-size effects). For $`H=0`$, the result $`d_c1.896`$ coincides with the analytic value for the percolation problem (91/48), suggesting a super-universality of $`d_c`$ for correlated percolation problem.

Random surfaces are common objects which can be found in various situations; good examples are deposited metal films Kondev2 , fractured surfaces Sahimi ; Bouchaud , streamfunctions of two-dimensional (2D) turbulent flows MBI  etc. For various applications, it is important to know the statistical geometrical properties (statistical topography) of such surfaces. This problem has been studied with a considerable success, both numerically and analytically Kondev2 ; Sahimi ; Bouchaud ; MBI ; IK ; Kondev1 ; Olami .

Here we introduce a new scale-invariant measure of such surfaces, the fractal dimension of the “oceanic coastlines”. An “oceanic coastline” includes the outer perimeter of the ocean (the “continental coastline”), and the coastline of “oceanic islands”, i.e. of such islands, which are adjacent to the largest (possibly infinite) connected flooded region, when water floods the surface up to a certain level; “oceanic islands” are to be opposed to such islands, which are inside the lakes. Thus, each flooding level defines its own “oceanic coastline”. Here we show that for the so called “monoscale” random surfaces (with bounded height, cf. IK ), the “oceanic coastlines” can be mapped to the percolation clusters of the correlated percolation problem. For rough surfaces (with unbounded height), the fractal dimension $`d_c`$ is calculated numerically, using a novel technique of minimizing the finite-size effects. The results suggest that for the correlated percolation (and hence for the “monoscale” surfaces), there is a super-universality, $`d_c91/48`$.

An important characteristic of random surfaces \[given by the surface height $`\psi (x,y)\psi (𝒓)`$ over the two-dimensional plane\] is the Hurst exponent $`H`$; in the case of Gaussian surfaces, all the statistical topography scaling exponents are functions of $`H`$. The most intuitive definition of the Hurst exponent is given by the scaling law of the surface height drop at a distance $`a|𝒂|`$ (c.f. Mandelbrot ):

$$(\psi (𝒓)\psi (𝒓+𝒂))^2|𝒂|^{2H}.$$
(1)
Here angular braces denote averaging over an ensemble (or over the radius vector $`𝒓`$); we assume that $`a`$ is larger than the lower cut-off scale $`a_0`$, and $`0H1`$. More general definition (not bounded to $`0H1`$) of the exponent $`H`$ is given via the power spectrum $`P_𝒌`$ \[with $`\psi _𝒌\psi _𝒌^{}=P_𝒌\delta (𝒌+𝒌^{})`$, where $`\psi _𝒌`$ is the Fourier transform of the surface height; c.f. MBI \]:

$$P_𝒌|𝒌|^{2H2}\text{for}|𝒌|a_0^1.$$
(2)
According to the value of $`H`$, surfaces can be divided into three classes IK :
(i) $`H<0`$: “monoscale”, surface looks like a landscape with hills of approximately same size, the surface height is bounded;
(ii) $`0H<1`$: “rough”, looks like rocky mountains with smaller and higher peaks and unbounded height;
(iii) $`H1`$: “smooth”, the largest scale defines the slope of the surface.

For “monoscale” surfaces, there is a simple mapping between statistical topography and percolation problem, developed by Ziman and Weinrib Ziman ; Weinrib1  \[the minima of the surface define the sites of an irregular percolation lattice; the saddle-points $`S_{ij}`$ (between $`i`$-th and $`j`$-th minima) correspond to the bonds of the lattice (which can be drawn as the steepest descend paths from $`S_{ij}`$ to the $`i`$-th and $`j`$-th minima); when the surface is flooded by water, a bond is declared broken if the corresponding saddle is (“dry”), i.e. above the water level; the water level itself is the counterpart of the bond-keeping probability $`p`$\]. This mapping can be also applied to the case of rough surfaces. Then, however, the lattice problem is no more the (correlated) percolation problem. In particular, there is no critical probability $`p=p_c`$ at which the global connectivity disappears. Indeed if we “look far enough”, we can find as deep (hence also as wide) connected wet region as we want, because the height is unbounded. Therefore, the mapping between the percolation problem and statistical topography allows us to extend various concepts of the percolation problem outside the range of applicability of the very percolation problem.

A straightforward example here is the hull of a percolation cluster: this is defined as the outer perimeter of a percolation cluster. Since the percolation cluster corresponds to a connected wet region, its outer perimeter corresponds to an isoline of the surface. It has been shown that for uncorrelated percolation, its fractal dimension $`d_h=7/4`$ Saleur . Note that percolation with weak correlation (with $`H<3/4`$) belongs to the same universality class as the uncorrelated percolation Weinrib2 . Meanwhile, for rough surfaces, $`d_h1.50.5H`$ JK , with an exact result $`d_h=1.5`$ for $`H=0`$ Kondev1 . Therefore, for $`3/4<H<0`$, one can expect a non-trivial decreasing function $`d_h(H)`$, with $`d_h(3/4)=1.75`$ and $`d_h(0)=1.5`$. Further, for the percolation problem, one can distinguish between the hull and the unscreened perimeter. This is because the hull is very twisted, and typically, there are narrow “necks” of peninsulas separating “gulfs”. We can postulate that these “necks” are too narrow, so that the coastline of the respective peninsulas is “screened”; the rest of the coastline is called “unscreened perimeter”. For percolation clusters, the fractal dimension of the unscreened perimeter is $`d_u=4/3`$ Saleur . Since for rough surfaces, narrow “necks” are untypical, $`d_u(H)=d_h(H)`$ for $`H0`$. Comparing $`d_u(3/4)=4/3`$ and $`d_u(0)=1.5`$, one can expect an increasing function $`d_u(H)`$ for $`3/4<H<0`$.

Now we are ready to turn to the percolation cluster. Each bond of the cluster is touched by the perimeter (hull) from both sides; thus, the length of the net perimeter scales as the cluster size. However, the perimeter touching it can be either external or internal, see Fig. 1. The sum of external and internal perimeters is the sum of (a) the coastline between the “ocean” (=cluster) and the “continent” (=infinite cluster of broken bonds), and (b) the coastline between the “ocean” and the “islands” (=broken bond clusters inside the cluster of undamaged bonds); we call this the “oceanic coastline”. Note that the coastlines of the “lakes” inside the “islands”, “islands” inside the “lakes” inside the “islands”, etc., do not touch the “ocean”, i.e. the bonds of the cluster, and are not relevant here (including these coastlines, too, would result in the set of all the isolines at the given height, the fractal dimension of which $`d_f=2H`$ for $`0H1`$; for $`H0`$, there is a super-universality, $`d_f2`$ ).

In order to calculate numerically the fractal dimension of the “oceanic coastlines”, we have opted for the following procedure. First, 1+1-dimensional \[(1+1)D\] random surfaces (which are assumed to belong to the same universality class as the statistically isotropic 2D surfaces) were generated on the lattice of the four-vertex model JK . The one-dimensional discretized fractional Brownian functions $`f_H(i)`$ were generated as fractional integrals of uncorrelated random sequence of “spins” $`s_i=\pm 1`$ as follows. The aim function $`g_H(i)`$ is found as

$$g_H(i)=\underset{j=0}{\overset{j}{}}s_j|ij|^{H0.5}.$$
(3)
Further, the discretized function $`f_H(i)`$ is found incrementally, starting from $`H=0`$, by rounding $`g_H(i)`$ to such a nearest integer that

$$f_H(i)f_H(i1)=\pm 1.$$
(4)
Finally, the (1+1)D function is given by the superposition of two 1D-functions,

$$F_H(i,j)=f_{H1}(i)f_{H2}(j).$$
(5)
After generating a 2D surface on a square polygon of side length $`N`$, we started “flooding” the polygon. To this end, at the perimeter of the polygon, the lowest point $`P_0`$ was found. Further, “water” was incrementally “injected” into that point, until the flooded region connected opposite edges of the polygon (either left and right or top and bottom), see Fig. 2. Such a procedure guarantees that only the “ocean” (and not the “lakes”) becomes wet. At that stage, three quantities were recorded: $`L_1`$ — the coastline length, $`L_2`$ — the number of cells touching the coastline, and $`L_3`$ — the coastline length immediately before achieving the critical flood level. The procedure was repeated (for each polygon size and Hurst exponent value) $`10^8`$ or more times. The covered polygon sizes are given by formulae $`N=1282^k`$ and N=$`1962^k`$, with $`k=1,2,3,4`$, and $`N=4096`$ (the largest polygon size was used only for $`H=0`$). The simulations were performed on the cluster of ten AthlonXP-1700 workstations during three months, using an assembly-optimized code.

It was expected that the finite-size scaling of the average values $`l_i(N)=L_i(N)`$ allows us to estimate the fractal dimension $`d_c`$. Such an asymptotic behavior was observed, indeed, for larger Hurst exponents ($`H0.375`$). However, for smaller values of $`H`$ (and in particular, for $`H=0`$), the convergence $`l_i(N)AN^{d_c}`$ was very slow, see Fig. 3 (here $`A`$ denotes a constant). Consequently, if we assume that $`l_i`$ can be expanded asymptotically,

$$l_\nu =\underset{\nu =0}{\overset{\mathrm{}}{}}A_{\nu i}N^{\alpha _{\nu i}},\alpha _{(\nu +1)i}<\alpha _{\nu i},$$
(6)
the non-leading exponents ($`\nu 0`$) have to be very close to the leading term. In order to resolve this difficulty, we made an assumption that for the three recorded quantities $`l_i`$ ($`i=1,2,3`$), the three leading exponents ($`\nu =0,1,2`$) are the same, $`\alpha _{\nu i}=\alpha _{\mu i}`$. The validity of this assumption is evident in the case of the leading exponent $`\alpha _{\nu 0}d_c`$, but not so clear for the other exponents. However, if we accept the assumption, it can be verified later using the numerical results. Thus, if we assume that $`\alpha _{\nu i}=\alpha _{\mu i}`$, one can construct such a linear combination of the three quantities $`l_i`$, for which two terms (out of the three leading ones) cancel out. It suffices to find a root-mean-square fit by minimizing the expression

$$\begin{array}{c}S(d)=\hfill \\ \mathrm{min}_{A_1,A_2,A_3}_{i=i_0}^{i_{\mathrm{max}}}\left[\frac{N^d(i)A_1l_1(i)A_2l_2(i)A_3l_3(i)}{\sigma _i}\right]^2.\hfill \end{array}$$
(7)
Here, $`\sigma _i`$ is the variance of the expression $`A_1l_1(i)A_2l_2(i)A_3l_3(i)`$, and can be expressed via the individual variances of the quantities $`l(i)`$, and the pair-wise correlation coefficients; all this is easily calculated, using the Monte-Carlo simulation results. The resulting curves $`\mathrm{log}_{10}[S(d)/(n4)]`$ are given in Fig. 4; here $`n`$ is the number of terms in the sum $`S(d)`$ (so that $`n4`$ is the number of degrees of freedom). Note that the two smallest data-points (with $`N_0=128`$ and $`N_1=196`$) have been excluded from the sum, otherwise, the residual squares would have been unacceptably large \[resulting in $`S/(n4)1`$\]. Evidently, for $`N196`$, the higher ($`\nu 3`$) terms of the expansion (6) become important. The existence of three sharp minima for all the curves in Fig. 4 confirms the assumption which was made about the equality of the leading exponents. Indeed, if the second and third asymptotic expansion terms of $`l_1`$, $`l_2`$, and $`l_3`$ were linearly independent, they would not cancel out in Eq. (7), and there would be no minima with $`S/(n4)1`$.

This novel method of interpreting the finite-size simulation results allowed us to calculate not only the fractal dimension of the “oceanic coastlines” (see Fig. 5), but also the second and third exponents in the asymptotic expansion of the length of “oceanic coastlines” (Fig. 6). The values $`d_c=1.8965\pm 0.0025`$ and $`\alpha _1=1.58\pm 0.01`$ for $`H=0`$ are relatively close to each other, exactly as it was expected, judging by the finite-size scaling in Fig. 3.

Note that for $`H3/4`$ (uncorrelated percolation), there is analytic result $`d_c=91/481.8958`$ Kapitulnik , which is indistinguishable (at our numerical uncertainty) from our result for $`H=0`$. While the author is not able to give a rigorous proof, it seems natural to assume that the function $`d_f`$ is a monotonically decreasing function of $`H`$. Indeed, a larger $`H`$ means a stronger influence of long scales, which results in less twisted isolines, and a lesser number of small islands (note that these arguments do not apply for $`d_u`$: excessive twists of the hull can efficiently screen and smooth the perimeter). This leads us to the conjecture $`d_c91/48`$ for $`H0`$.

In conclusion, we have introduced the concept of “oceanic coastlines”, and shown that they are the extension of the concept of percolation clusters. The finding that their fractal dimension $`d_c`$ is clearly distinct from both the dimension of single isolines $`d_h`$ and the dimension of the full set of isolines at a given height $`d_f`$ is particularly important, because in the context of Earth’s landscapes, the coastlines of ”oceanic islands” have been often mistakenly treated as the full set of isolines Mandelbrot . Based on the numerical results, we have conjectured that for $`H0`$, the fractal dimension of percolation clusters is super-universal, $`d_c91/48`$. Last but not least, a novel and universal approach to the problem of interpreting the finite-size simulation results has been devised.


text:
# Fermi surface of the heavy fermion superconductor PrOs4Sb12

## Abstract

We have investigated the de Haas-van Alphen effect in the Pr-based heavy fermion superconductor PrOs4Sb12. The topology of Fermi surface is close to the reference compound LaOs4Sb12 and well explained by the band structure calculation based on the FLAPW-LDA+U method, where the 4f electrons to be localized. We have confirmed a highly enhanced cyclotron effective mass 2.4$``$7.6$`m_0`$ which is apparently large compared to the usual Pr-based compounds.

The Pr-based compounds in the filled skutterudites have attracted much attention because of their exotic properties, such as metal-insulator transition at $`T_{\mathrm{MI}}=60`$ K in PrRu4P12Sekine  and unusual heavy fermion (HF) behavior in PrFe4P12Sato ; Sugawara ; Aoki  The quadrupolar interaction is thought to be a key mechanism to explain such anomalous behaviors. Especially for PrFe4P12, both the extraordinarily enhanced effective mass ($`m_\mathrm{c}^{}=81m_0`$) and the non-magnetic low field ordered state below 6.5 K are inferred to originate from the quadrupolar interaction. Sugawara ; Aoki ; Hao

Recently, PrOs4Sb12 was reported to exhibit superconductivity below $`T_\mathrm{C}=1.85\mathrm{K}`$Bauer ; Maple  From the large specific heat jump at $`T_\mathrm{C}`$, $`\mathrm{\Delta }C/T_\mathrm{C}`$$``$ 500 mJ/K$`{}_{}{}^{2}`$ mol, Bauer et al. claimed that PrOs4Sb12 is the first example of the Pr-based HF-superconductor. Bauer ; Maple  They also inferred that the quadrupolar interaction plays an important role in the HF-superconductivity, since the magnetic susceptibility, specific heat $`C(T)`$, and inelastic neutron scattering measurements suggest the crystal electric field (CEF) ground state to be a non-Kramers doublet carrying quadrupole moments. The recent Sb-NQR measurement on PrOs4Sb12 have revealed that the temperature $`T`$ dependence of nuclear-spin-lattice-relaxation rate $`1/T_1`$ shows neither a coherence peak nor $`T^3`$ like power-law dependence below $`T_\mathrm{C}`$Kotegawa  The muon-spin relaxation measurements suggest an isotropic superconducting (SC) energy gap, MacLaughlin  while the anisotropy of thermal conductivity against the magnetic field directions indicates that the SC energy gap has point nodes. Izawa  Such unusual properties cannot be consistently explained neither by a conventional isotropic nor any anisotropic gap superconductivity reported until now, suggesting some novel type of HF-superconductivity of this compound. Bauer ; Maple ; Kotegawa ; MacLaughlin ; Izawa ; Miyake  In the normal state above $`4.5`$ T, Aoki et al. reported an anomalous field-induced ordered phase (FIOP) based on the specific heat measurements; possibly a quadrupolar ordering. Aoki\_JPSJ

In order to understand these unusual properties, it is essential to determine the electronic structure of this compound and to directly confirm the highly enhanced effective mass. In this paper, we report the de Haas-van Alphen (dHvA) experiment on PrOs4Sb12Sugawara\_SCES02  which is the most powerful tool to clarify the Fermi surface (FS) topology and to determine the cyclotron effective mass. The experimental result is compared with the band structure calculation based on the FLAPW-LDA+U method.

Single crystals of PrOs4Sb12 and the reference LaOs4Sb12 were grown by Sb-self-flux method basically same as in ref. Takeda ; Bauer\_CeOs4Sb12 , using high-purity elements, 4N (99.99% pure)-Pr, 4N-La, 3N-Os and 6N-Sb. The typical forms of the single crystals were cubic or rectangular shape with a largest dimension of about 3 mm. The crystal structure of filled skutterudite, belonging to the space group $`\mathrm{𝐼𝑚}`$3̄ ($`T_h^5`$, #204), was verified by the powder neutron diffraction experiments performed by VEGA at KEK Tsukuba, Japan. The lattice constant of PrOs4Sb12 obtained from the Rietveld refinement is $`a=9.30311`$ Å , which is close to the reported value. Braun  Fractional coordinates of Sb at the 24g site are determined as (0, 0.3405, 0.1561) which is used in the band structure calculation. The residual resistivity $`\rho _0`$ and the residual resistivity ratio (RRR) of the present samples are $`\rho _0=8\mu \mathrm{\Omega }`$cm and $`\mathrm{RRR}=55`$ for PrOs4Sb12, and $`\rho _0=2.8\mu \mathrm{\Omega }`$cm and $`\mathrm{RRR}=100`$ for LaOs4Sb12. The dHvA experiments were performed in a top loading dilution refrigerator cooled down to 30 mK with a 17 T superconducting magnet. The dHvA signals were detected by means of the conventional field modulation method with a low frequency ($`f10`$ Hz).

Figure 1 shows (a) the typical dHvA oscillations and (b) the fast Fourier transformation (FFT) spectra in PrOs4Sb12 for the field $`H`$ along $`100`$ direction.

The oscillations become detectable just above $`H_{\mathrm{C2}}=2.2`$ T and a faint peak anomaly appears at $`H_\mathrm{A}4.4`$ T, indicating a slight slope change in the field dependence of magnetization $`M(H)`$Tenya  This field $`H_\mathrm{A}`$ almost agrees with the phase boundary of FIOP determined by the $`C(T)`$ measurement in the magnetic fields. Aoki\_JPSJ  Above $`H_\mathrm{A}`$ a small spin-splitting is observed in the dHvA frequencies reflected as beats in the oscillations, which is frequently observed in the magnetic materials. Onuki  The three fundamental dHvA branches $`\alpha `$, $`\beta `$ and $`\gamma `$ are identified.

Figure 2 shows the angular dependence of dHvA frequency in PrOs4Sb12 along with that in LaOs4Sb12.

$`\alpha `$\- and $`\beta `$-branches have been observed over the whole field directions in the cubic symmetry, indicating closed FSs. The branch $`\gamma `$ observed only in the limited angular ranges centered at $`100`$ suggests a part of the multiply connected one. The striking resemblance of the angular dependence of dHvA frequencies between PrOs4Sb12 and LaOs4Sb12 indicates the closeness of FS topology between the two compounds, and evidences the well localized nature of 4$`f`$-electrons in PrOs4Sb12.

In order to assign the origin of dHvA branches, the band structure calculation is carried out using an FLAPW and LDA+U method, in which the 4$`f`$-electrons in PrOs4Sb12 is treated as localized. Here we assume the $`\mathrm{\Gamma }_1`$ singlet to treat the localized $`4f^2`$ electrons with the cubic symmetry. The details of the calculation are described in ref. Harima\_JMMM\_01 ; Harima\_PhysB\_02 . Figures 3 and 4 show the calculated energy band structure and the FS in PrOs4Sb12, respectively.

The FS is composed of the 48th and 49th band hole sheets. The 48th band forms an almost spherical sheet which slightly projects along the $`100`$ direction centered at the $`\mathrm{\Gamma }`$ point. The 49th band gives a round cubic sheet centered at $`\mathrm{\Gamma }`$ point and a multiply connected one whose main parts are centered at the N points. Figure 5 shows the angular dependence of the calculated dHvA frequencies, which reasonably well explains all the observed dHvA branches. A slight disagreement in the absolute values is probably ascribed to the assumption of the $`\mathrm{\Gamma }_1`$ singlet ground state and/or the effect beyond the LDA+U treatment.

Branch $`\beta `$ originates from the 48th-band spherical FS, branch $`\alpha `$ is due to the 49th-band round cube FS and branch $`\gamma `$ originates from the extremal orbit surrounding the N point of the part of multiply connected one. Many other dHvA branches predicted in the band structure calculation have not been observed in the present experiments. That can be understood as a large reduction of the dHvA amplitude due to the large curvature factors $`A^{\prime \prime }`$ of these branches; i.e., for $`H100`$, $`|A^{\prime \prime }|=0.08`$ for the $`\gamma `$-branch whereas $`|A^{\prime \prime }|=2.72`$ for the dHvA branch exhibiting the frequency of $`F=0.71\times 10^3`$ T.

As the most attractive information, we have estimated the cyclotron effective mass $`m_\mathrm{c}^{}`$ from the temperature dependence of the dHvA amplitude. The comparison of dHvA frequencies and $`m_\mathrm{c}^{}`$ between these compounds is given in Table 1Table1

In contrast with the closeness of the FS topology, the cyclotron effective masses are up to $`6`$ times enhanced compared to those in LaOs4Sb12. The mass enhancement is apparently larger than those in the ordinary Pr-based compounds such as PrIn3 and PrSb, Onuki  although it is not so large as that in PrFe4P12Sugawara  For the comparison with the mass enhancement estimated from the specific heat measurement, we simply estimate the Sommerfeld coefficient from the FS volume and $`m_\mathrm{c}^{}`$ in the present experiments assuming spherical FSs. The estimated value of $`150`$mJ/K$`{}_{}{}^{2}`$mol is still a factor of 2 or more smaller than that in the specific heat measurement; Bauer et al. estimated the Sommerfeld coefficient of PrOs4Sb12 to be $`350750`$ mJ/K$`{}_{}{}^{2}`$mol. Bauer

We have no definitive explanation for this discrepancy only from the present information. Nevertheless, a possible origin may be hidden in the difference of the measuring conditions, such as magnetic field and temperature; the Sommerfeld coefficient in the specific heat was estimated at 0 T, and around $`T_\mathrm{C}`$ close to the Schottky-like peak at 3 K, while the dHvA measurements were made at high fields of $`317`$ T and lower temperatures of $`30400`$ mK. If the peak is associated with the nearby magnetic CEF excited state situated about 8 K above the ground state, Bauer ; Aoki\_JPSJ  an additional enhancement of effective mass could be expected from the magnetic instability. As another origin, the spin-splitting effect should not be discarded. In the present experiments we have estimated $`m_\mathrm{c}^{}`$ only for the larger amplitude spin-direction, since the insufficient experimental resolution prevents to determine the masses for the higher frequencies ($`\alpha ^{}`$ and $`\beta ^{}`$ as shown in Fig. 1(b)). Only for $`\beta `$-branch with the large amplitude, we could roughly estimate $`m_\mathrm{c}^{}`$ of the higher frequency branch $`\beta ^{}`$, which is about 20% enhanced compared with that of the lower one $`\beta `$. The larger mass enhancement for the higher frequency branch might explain the discrepancy of the mass enhancement between the Sommerfeld coefficient and $`m_\mathrm{c}^{}`$.

It might be informative to compare PrOs4Sb12 with another Pr-based filled skutterudite superconductor PrRu4Sb12 as shown in Table 2Sugawara\_ROs4Sb12

The localized character of 4$`f`$-electrons, namely the closeness of the FSs with those in LaRu4Sb12, has been confirmed also in PrRu4Sb12 based on the dHvA experiment. Matsuda  On the contrary, the mass enhancement is quite small in PrRu4Sb12 which is in sharp contrast to that in PrOs4Sb12. For PrOs4Sb12, the CEF ground state was inferred to be a non-Kramers doublet carrying quadrupole moments. Bauer ; Maple  On the other hand, Takeda and Ishikawa inferred PrRu4Sb12 to have the $`\mathrm{\Gamma }_1`$ singlet ground state. Takeda\_JPSJ  On the $`T_\mathrm{C}`$ compared with the La-references, the two compounds have different sense; $`T_\mathrm{C}`$ for PrOs4Sb12 is higher than that for La-reference, that is unusual if we take into account that PrOs4Sb12 contains the magnetic element Pr. It should be also noted that the superconductivity observed in PrRu4Sb12 is the ordinary BCS-type. Takeda\_JPSJ ; Kotegawa  These facts naturally indicate the essential role of the mass enhancement associated with 4$`f`$-electrons in the superconductivity of PrOs4Sb12.

In summary, we have directly confirmed the heavy mass in PrOs4Sb12 by the dHvA experiment, indicating PrOs4Sb12 to be the first member of HF-superconductor as Pr-compound. However, the well localized character of 4$`f`$-electrons, evidenced in the closeness of the FS to those in LaOs4Sb12, indicates the uniqueness of this compounds compared to the typical HF-superconductors such as UPt3 and CeCoIn5Kimura ; Settai  where the itinerant $`f`$-electron model is applicable. From the comparison with another Pr-based filled skutterudite superconductor PrRu4Sb12, the quadrupolar interaction associated with the non-Kramers doublet CEF ground state is thought to be responsible for the HF-superconductivity in PrOs4Sb12.

We thank Dr. H. Kotegawa, Professor Y. Kitaoka, Professor M. B. Maple, and Professor K. Miyake for the helpful discussion. This work was supported by the Grant-in-Aid for Scientific Research from the Ministry of Education, Science, Sports and Culture of Japan.