text:
# Gamma-Ray Bursts - When Theory Meets Observations

## Introduction

Gamma-Ray bursts - GRBs, short and intense bursts of $`\gamma `$-rays arriving from random directions in the sky were discovered accidentally more than thirty years ago. During the last decade two detectors, BATSE on CGRO and BeppoSAX have revolutionized our understanding of GRBs. BATSE has demonstrated BATSE92  that GRBs originate at cosmological distances in the most energetic explosions in the Universe. BeppoSAX discovered X-ray afterglowCosta97 . This enabled us to pinpoint the positions of some bursts, locate optical vanP97  and radio Frail97  afterglows, identify host galaxies and measure redshifts to some bursts Metzger97 .

Since their discovery GRBs were among the prime topics of the Texas Symposia. The high energy release and the rapid time scales involved suggested immediately association with relativistic compact objects. The discoveries of BATSE and BeppoSAX confirmed these expectations. These observations have established the Fireball model demonstrating that GRBs are the most relativistic objects known so far: GRBs involve macroscopic ultrarelativistic flows with Lorentz factors $`\mathrm{\Gamma }100`$. Furthermore, while the “central engines” that drive the relativistic flow and power the GRBs are hidden we have excellent evidence that they involve accretion onto a newborn black hole. GRBs are the birth cries of these black hole.

I review, here, the recent progress in our understanding of GRBs, emphasizing, as appropriate for this conference, their relativistic nature. I begin, in section I with a brief tribute to the 7th Texas symposium. This was the first Texas meeting after the discovery of GRBs was announced and GRBs were the highlight of the discussion there. I continue in section II with a brief exposition of the Fireball model (see P99 ; P00  for details), confronting its predictions with the observations. In III I summarize the implication of the fireball model to the “inner engines”. Concluding remarks, further predictions and open questions are discussed in section IV.

## I The 7th (New York) Texas Symposium

GRBs were the hightlight of the Seventh (New York) Texas symposium that took place in 1974. Five out of the 57 talks were devoted to GRBs (this record was repeated only in this symposium with 3 out of 29 talks): Two observational reviews, a theoretical review, a theoretical model and even a description of an automated system for searching for optical transients accompanying GRBs!

M. Ruderman Ruderman75  reviewed the theory111This review enumerates more than thirty models proposed during the short time passed since the announcement of the discovery. It is remarkable (Ruderman, 1998, private communication) that today we know that none of these models is even remotely relevant.. emphasizing the compactness problem: If GRBs are cosmological then the energy budget and the time scales seem to be incompatible with the observed non thermal spectrum of the bursts. The argument is simple: the variability time scale, $`\delta t`$, imposes an upper limit on the size ($`Rc\delta t`$). The observed flux and the assumed (cosmological) distance determine the energy. Together these yield an extremely large lower limits on the photons density within the source and on the optical depth for pair creation by the energetic photons. Pairs would be copiously produced and the source would be optically thick. The observed optically thin spectrum is impossible. Ruderman points out, however, that relativistic effects would change this conclusion. If the source is moving relativistically the relations between the time scale and the implied distance are modified (by a factor of $`\mathrm{\Gamma }^2`$). Furthermore, photons that are observed with energy below 500$`\mathrm{\Gamma }`$keV have energy below 500keV in the source rest frame and could not produce pairs. These ideas lay the foundation for the current Fireball model. Recent observations have indeed confirmed ultrarelativistic motion in GRBs, showing at least in one case $`\mathrm{\Gamma }100`$.

## II The Fireball Model, Predictions and Confirmations

One can never prove a scientific theory. However we gain confidence in a theory when its predictions are confirmed by observations. I discuss, here, the predictions of the Fireball model (specifically of the Fireball-Internal-External shocks model) and their confirmation by numerous observations. My goal is to demonstrate the success of this model. While some of the specific observations could certainly be interpreted within other theories I strongly believe that the bulk of those observations tell us that this is the correct model.

The Fireball model asserts that GRBs are produced when the kinetic energy (or Poynting flux) of a relativistic flow is dissipated by shocks222Given the low densities involved these shocks, like SNR shocks, must be collisionless.. These shocks accelerate electrons and generate strong magnetic fields. The relativistic electrons emit the observed $`\gamma `$-rays via synchrotron or SSC. There are two variants of this model: The External Shocks model MR93a  assumes that the shocks are between the relativistic flow and the surrounding circumstellar matter. The Internal Shocks model NPP92 ; RM94  assumes that the flow is irregular and the shocks take place between faster and slower shells within the flow. According to the Internal-External shocks model SP97  both kinds of shocks take place: Internal shocks are responsible for the GRB while external shocks produce the longer lasting afterglow (see Fig. 1) . Both shocks occur at relatively large distances ($`10^{13}10^{14}`$cm for internal shocks and $`10^{14}10^{16}`$cm for external shocks) from the source that generates the relativistic flow. The observed radiation from the GRB or from the afterglow reflects only the conditions within these shocks. We have only indirect information on the nature of the “inner engines”.

Fig. 1 depicts a schematic picture of the Internal-External shocks model. An inner engine produces an irregular wind. The wind varies on a scale $`\delta t`$ and its overall duration is $`T`$. The variability scale $`\delta t`$ corresponds to the variability scale observed in the GRB light curve KPS97  thus, $`\delta t1`$ sec. Internal shocks take place at $`R\delta t\mathrm{\Gamma }^2310^{14}\mathrm{cm}(\delta t/1sec)(\mathrm{\Gamma }/100)^2`$. External shocks shocks become significant and the blast wave that propagates into the circunstellar matter and produces the afterglow begins at $`10^{16}`$cm (see SaP95 ; P99  for details). At this initial stage there is also a short lived reverse shock that propagates into the ejecta. This reverse shock is responsible to the prompt optical emission observed in GRB990123.

I review now the predictions of this model and confront them with observations.

$``$ Relativistic Motion - Predictions: Relativistic motion is the key ingredient of the Fireball model. Relativistic motion arises naturally (a Fireball forms) when a large amount of energy is produced in a compact region with $`E/Mc^21`$ Goodman86 ; Pac86 ; SP90 . This relativistic motion overcome the compactness problem. Various estimates FenimoreEH93 ; WoodsL95 ; P95 ; SaLi00 , of the Lorentz factor $`\mathrm{\Gamma }`$ based on the compactness problem lead to comparable values, $`\mathrm{\Gamma }100`$, (see SaLi00  for a critical review).

$``$ Relativistic Motion - Observations I: The radio afterglow observations of GRB 970508 provided the first verification of relativistic motion. The radio light curve (in 4.86Ghz) variesd strongly during the first month. These variations died out later. Even before this transition Goodman Goodman97  interpreted these variations as scintillations. The observation of a transition after one month enabled Frail et al., Frail97  to estimate the size of the afterglow at this stage as $`10^{17}`$cm. It immediately follows that the afterglow expanded relativistically. Additionally, the source is expected to be optically thick in radio KP97  leading to a $`\nu ^2`$ rising spectrum at these frequencies. The observed flux from the source enables us (using the black body law) to estimate the size of the source. As predicted the radio spectrum increases like $`\nu ^2`$. The size estimated with this method agrees Frail97  with the one derived by the scintillations estimate implying as well a relativistic motion.

$``$ The Afterglow - Predictions: The Afterglow - lower frequency emission that follows GRBs was one of the earlier predictions of the Fireball model PaczynskiR93 ; MR93 ; MR97 ; Vietri97 . Paczynski and Rhoads PaczynskiR93  predicted radio afterglow on the basis of the analogy between external shocks and SNRs. Later Meszaros and Rees MR93 ; MR97  performed detailed calculations of multi-wavelength afterglow. Vietri Vietri97  predicted soft x-ray afterglow as a test for the external shocks model. These predictions were done in the context of the external shocks model. According to this model the GRB arises due to a shock between a relativistic flow and a circumstellar medium. In this case the emission observed at time $`t_{obs}`$ arises when $`t_{obs}R/2\mathrm{\Gamma }^2`$. Later emission is related to lower $`\mathrm{\Gamma }`$ (and larger emission radii) and hence to lower observed frequencies. This afterglow would be a direct extrapolation of the GRB and its basic features should be strongly correlated with the properties of the corresponding GRB.

The theory of the afterglow is well understood. Blandford and McKee BM76  have worked out (already in the seventies!) the theory of an adiabatic relativistic blast wave. They show that (as long as the flow is ultrarelativistic, $`\mathrm{\Gamma }1`$), the blast wave is self similar, the relativistic analog of the well known Sedov-Taylor solution. Electrons are accelerated to relativistic velocities by the shocks and their interaction with the magnetic field leads to synchrotron radiation. This provides an excellent model for the observed emission SPN98 . Overall we have a simple theory characterized by five parameters: the total energy, $`E_0`$, the ambient density, $`n_0`$, the ratio of the electrons and magnetic fields energy density to the total energy density, $`ϵ_e`$, $`ϵ_B`$ and the exponent of the electrons’ energy distribution function $`p`$. An additional sixth parameters, the exponent of the circumstellar density distribution, $`n`$, arises in cases when the external matter density ($`\rho r^n`$). Most notable is $`n=2`$ corresponding to a pre-GRB wind expected in some models Chevalier00 . This rather simple theory predicts a robust relations between $`\alpha `$ and $`\beta `$ the exponents describing the flux as a function of frequency, $`F_\nu t^\alpha \nu ^\beta `$. At the high frequencies, above the cooling frequency, we have (for $`n=0`$), $`\alpha =(3p2)/4`$ and $`\beta =p/2`$.

$``$ The Afterglow - Observations: On Feb 28 1997, in a wonderful anniversary celebration for SN87A, BeppoSAX detected x-ray afterglow from GRB 970228. The exact position given by BeppoSAX led to the discovery of optical afterglow vanP97 . Radio afterglow was detected in GRB 970508 Frail97 . By now more than thirty x-ray afgterglows have been observed. About half of these have optical and radio afterglow as well and in most of those the host galaxy has been discovered.

Most x-ray and optical afterglow decay as power laws with $`\alpha 1.2`$ and $`\beta 1.2`$, in excellent agreement with the predictions of the simplest afterglow model: An adiabatic Blandford-McKee hydrodynamics with Synchrotron emission SPN98 . As $`\alpha `$ and $`\beta `$ are determined by $`p`$ the electron’s distribution power law index, these observations suggest that as predicted SPN98 ; SNP96 , $`p2.5`$ and it is fairly invariant from one burst to another PanaitescuK01 . A simultaneous spectral fit for GRB980508, all the way from the radio to the x-ray also agrees with this picture WijersGalama99 .

$``$ The GRB-Afterglow Transition - Predictions: The rapid time variability seen in most GRBs cannot be produced by external shocks SP97 ; Fenimore96 . This leaves internal shocks as the only viable model! Shortly before the discovery of the afterglow from GRB970228, Sari and Piran SP97  pointed out that afterglow should arise also within the internal shocks scenario. The efficiency of the internal shocks depends on the parameters of the flow, most specifically on the variability of the Lorentz factor between different shells KPS97 ; Daigne98 . Even in the most efficient cases a significant fraction of the energy remain as kinetic energy. Sari and Piran SP97  suggested that this energy would be dissipated later by interaction with the surrounding matter and produce an afterglow. Within this Internal-External shocks model the GRB is produce by internal shocks while the afterglow is produced by external shocks. The predictions of this model for the afterglow are similar to those of the External shocks model. However, a critical difference is that here the afterglow is not an extrapolation of the GRB.

The internal shocks take place at a distance $`R_{IS}\delta t\mathrm{\Gamma }^210^{14}`$cm. These shocks last as long as the inner engine is active. The typical observed time scale for this activity $`50`$sec (for long bursts) and $`0.5`$sec (for short ones). External shocks begin at $`R_{Ex}10^{16}`$cm. If $`R_{Ex}/\mathrm{\Gamma }^2T`$ this happens while internal shocks are still going on and the afterglow overlaps the late part of the GRB. At the early time the afterglow emission peaks in the high x-rays contributing also to the observed $`\gamma `$-ray flux. We expect, therefore, a transition within the GRB from hard (pure GRB) to softer and smoother (GRB and afterglow) signal.

$``$ The GRB - Afterglow Transition - Observations: The extrapolation of the x-ray afterglow fluxes backwards generally does not fit the $`\gamma `$-ray fluxes. Moreover there is no direct correlation between the $`\gamma `$-ray fluxes and the x-ray (or optical) afterglow fluxes. This result is in a nice agreement with the predictions of the Internal - External shocks scenario in which the two phenomena are produced by different effects and one should not expect a simple extrapolation to work.

The expected GRB afterglow transition have been observed in several cases. The first observation took place (but was not reported until much latter) already in 1992 Burenin99 . Recent BeppoSAX data shows a rather sharp transition in the hardness that takes place several dozen seconds after the beginning of the bursts Costa99 . This transition is seen clearly in the different energy bands light curves of GRB990123 and in GRB980923 Giblin99 . Connaughton Connaughton98  have averaged the light curves of many GRBs and discovered long and soft tails: the early x-ray afterglow.

$``$ The Prompt Optical Flash - Predictions: The collision between the ejecta and the surrounding medium produces two shocks. The outer forward shock propagates into the ISM. This shock develops later into the self similar Blandford-McKee blast wave that drives the afterglow. A second shock, the reverse shock, propagates into the flow. This reverse shock is short lived. It dies out when it runs out of matter as it reaches the inner edge of the flow. While it is active, it is a powerful source of energy. Comparable amounts of energy are dissipated by the forward and by the reverse shocks SaP95 . We expect that the system is radiative at this stage, namely most of the energy converted by the shock is radiated away.

Sari and Piran predicted at the First Rome Meeting (Oct 1998) an intense (brighter than 11th magnitude) prompt optical flash from this reverse shock SP99c . Previous work MR97  done prior to the discovery of the afterglow considered various possibilities and estimated the magnitude of the prompt optical flash to be anywhere from 9th to 19th magnitude. The observations of the afterglow constrained severely the relevant models and the relevant parameter space. With the new data the constrained model led to a clear prediction with a narrow range of magnitude. At that time this prediction was almost conflicting with upper limits given by systems like LOTIS and ROTSE.

One prediction that was, unfortunately, missed: prompt radio emission from the reverse shock. This radio emission should be short lived, like the burst and should have initially an optically thick component that becomes optically thin later.

$``$ The Prompt Optical Flash - Observations: In Jan 23 1999 just three month after this prediction ROTSE recorded six snapshots of optical emission from GRB990123ROTSE  . Three of those were taken while the burst was still emitting $`\gamma `$-rays. The other three snapshots spanned a couple of minutes after the burst. The second snapshot, taken 70sec after the onset of the burst corresponds to a 9th magnitude signal. A comparison of these optical observations with the $`\gamma `$-rays and x-rays light curves (see e.g. P00 ) shows that the optical emission does not correlate with the $`\gamma `$-rays pulses. The optical photons and the $`\gamma `$ rays are not emitted by the same photons SP99b ; ROTSE . The optical pulses peak some 70sec after the onset of the burst simultaneously with a late peak in the soft x-ray emission.

Radio observations of GRB990123 revealed a short lived radio pulse. This emission can be explained as coming from the reverse shock SP99b . Using the parameters of the reverse shock derived from the optical flash Sari and Piran SP99b  estimated the magnitude of this radio emission. The theoretical curve and the observations are in excellent agreement (see e.g. P00 ). Note that the theoretical curve was calculated just from the optical flash data and it was not “fitted” in any way to the observed data. While prompt optical flashes were not detected in other bursts, short lived radio flashes have been detected in GRB000926 and in GRB970828.

$``$ Relativistic Motion II - Observations: The radio emission from GRB970508 showed relativistic motion in its afterglow. However, the significant observations were done one month after the burst and at that time the motion was only “mildly” relativistic with a Lorentz factor of order a few. The observations of GRB990123 enabled us to obtain three independent estimates of the ultrarelativistic Lorentz factor at the time that the ejecta hits first the ISM SP99b . First the time delay between the GRB and the optical flash suggests $`\mathrm{\Gamma }200`$. The ratio between the emission of the forwards shock (x-rays) and the reverse shock (optical) gives another estimate of $`\mathrm{\Gamma }70`$. Finally the fact that the maximal synchroton frequency of the reverse shock was below the optical band led to $`\mathrm{\Gamma }200`$. The agreement between these three crude and independent estimates is reassuring.

These observations provide us also with an estimate of the position of the external shocks, $`10^{15}`$cm at 70 seconds after the bursts. It is an impressive measurement considering the fact that the distance to this burst is $`3.5Gpc`$. The corresponding angular resolution is $`10^{13}`$ or $`50\mathrm{n}\mathrm{a}\mathrm{n}\mathrm{o}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{s}\mathrm{e}\mathrm{c}`$.

$``$ Jets - Predictions: With redshift measurements it became possible to obtain exact estimate the total energies involved. While the first burst GRB970508 required a modest value of $`10^{51}`$ergs, the energies required by other bursts were alarming, $`3\times 10^{53}`$ergs for GRB981226 and $`4\times 10^{54}`$ergs for GRB990123, and unreasonable for any simple compact object model. These values suggested that the assumed isotropic emission was wrong and GRBs are beamed. Significant beaming would of course reduce, correspondingly the energy budget.

Beaming was suggested even earlier as it arose naturally in some specific models. For example the binary neutron star merger has a natural funnel along its rotation axis and one could expect that any flow would be emitted preferably along this axis. The Collapsar model also requires beaming, as only a concentrated beamed energy could drill a whole through the stellar envelope that exists in this model.

Consider a relativistic flow with an opening angle $`\theta `$. As long as $`\theta >\mathrm{\Gamma }^1`$ the forwards moving matter doesn’t “notice” the angular structure and the hydrodynamics is “locally” spherical P94 . The radiation from each point is beamed into a cone with an opening angle $`\mathrm{\Gamma }^1`$. It is impossible to distinguish at this stage a jet from a spherical expanding shell. When $`\theta \mathrm{\Gamma }^1`$ the radiation starts to be beamed sideways. At the same time the hydrodynamic behaviour changes and the material starts expanding sideways. Both effects lead to a faster decrease in the observed flux, changing $`\alpha `$, the exponent of the decay rate of the flux to: $`\alpha =p/2`$. Thus we expect a break in the light curve and a new relation between $`\alpha `$ and $`\beta `$ after the break Rhoads99 ; SPH99 ; PM99 . The magnitude of the break and the duration of the transition will change if the jet is expanding into a wind with $`r^2`$ density profile PanaitescuK00 . The break is expected to take place at $`t_{jet}6.2(E_{52}/n_0)^{1/3}(\theta /0.1)^{8/3}`$hr SPH99 . Recently numerical simulation Granot01  have shown that the break appears in a more realistic calculations, even though the numerical results suggest that the analytical model developed so far are probably too simple.

$``$ Jet - Observations: GRB980519 had unusual values for $`\alpha =2.05`$ and $`\beta =1.15`$. These values do not fit the ”standard” spherical afterglow model333A possible alternative fit is to a wind (n=2) model but with a unsual high value of $`p=3.5`$. However, these values are in excellent agreement with a sideway expanding jet SPH99 . The simplest interpretation of this data is that we observe a jet during its sideway expansion phase (with $`p=2.5`$). The jet break transition from the spherical like phase to this phase took place shortly after the GRB and it was not caught in time. The light curves of GRB990123 shows, however, a break at $`t2`$days Kulkarni99 . This break is interpreted as a jet break, corresponding to an opening angle $`\theta 5^o`$. Another clear break was seen in GRB990510 Harrison99 ; Stanek99 .

The brightest bursts, GRB990123 and GRB980519 gave the first indications for jet like behaviour SPH99 . This suggested that their apparent high energy was due to the narrow beaming angles. A compilation of more bursts with jet breaks suggests that all bursts have a comparable energy $`10^{51}`$ergs and the variation in the observed energy is mostly due to the variation in the opening angles $`\theta `$ Frail01 ; Kulkarni01 ; PanaitescuK01

## III The Inner Engines

The Fireball model tells us how GRBs operate. However, it does not answer the most interesting astrophysical question: what produces them? which astrophysical process generates the energetic ultrarelativistic flows needed for the Fireball model? Several observational clues help us answer these questions.
$``$ Energy: The total energy involved is large $`10^{51}`$ergs, a significant fraction of the binding energy of a stellar compact object. The “inner engine” must be able to generate this energy and accelerate $`10^5M_{}`$ to relativistic velocities.
$``$ Beaming: Most GRBs are beamed with typical opening angles $`0.02<\theta <0.2`$. The “inner engine” must be able to collimate the relativistic flow.
$``$ Long and Short Bursts: The bursts are divided to two groups according to their overall duration. Long bursts with $`T>2`$sec and short ones with $`T<2`$sec.
$``$ Rates: GRBs take place once per $`10^7(4/\theta ^2)`$yr per galaxy. GRBs are very rare at about 1/1000 the rate of supernovae.

The Fireball Internal-External shocks model provides us with another key clue:
$``$ Time Scales: The variability time scale, $`\delta t`$, is at times as short as 1ms. The overall duration, $`T`$, is of the order of 50sec. According to the internal shocks model these time scales are determined by the “inner engine”. $`\delta t`$ msec suggests a compact object. $`T50`$sec is much longer than the dynamical time scale, suggesting a prolonged activity.444The ratio $`\delta t/T1`$ for short bursts as well NakarP01 . This rules out any “explosive” model that release the energy in a single explosion.

The internal shocks model requires two (or possibly three Ramirez00 ; NakarP01a ) different time scales operating within the “inner engine”. These clues, most specifically the last one suggest that GRBs arise due to accretion of a massive ($`0.1m_{}`$) disk onto a compact object, most likely a newborn black hole. A compact object is essential because of the short time scales. Accretion is needed to produce the two different time scales, and in particular the prolonged activity. A massive ($`0.1m_{}`$) disk is required because of the energetics. We expect that such a massive disk can form only simultaneously with the formation of the compact object. This leads to the conclusions that GRBs accompany the formation of black holes. This model is supported by the observations of relativistic (but not as relativistic as in GRBs) jets in AGNs, which are powered by accretion onto black holes. This system is capable of generating collimated relativistic flows even though we don’t understand how.

An important alternative to accretion is Usov’s model Usov92  in which the relativistic flow is mostly Poynting flux and it is driven by the magnetic and rotational energies of a newborn rapidly rotating neutron star. However this model seems to fall short by an order of magnitude of the energy required.

Several scenarios could lead to a black hole - massive accretion disk system. This could include mergers (NS-NS binaries Eichler\_LPS89 ; NPP92 , NS-BH binaries Pac91  WD-BH binaries Fryer\_WHD99 , BH-He-star binaries Fryer\_Woosley98 ) and models based on “failed supernovae” or “Collapsars” Woosley93 ; Pac98 ; MacFadyen\_W99 . Narayan et al. NarayanPK01  have recently shown that accretion theory suggests that from all the above scenarios only Collapsars could produce long bursts and only NS-NS (or NS-BH) mergers could produce short bursts.

Additional indications arise from afterglow observations. One has to use these clues with care. Not all GRBs have afterglow (for example, so far afterglow was not detected from any short burst) and it is not clear whether these clues are relevant to the whole GRB populations. These clues seem to suggest a GRB-SN connection:
$``$ SN association: Possible association of GRB980425 with SN98bw GalamaSN98  and possible SN signatures in the afterglows of GRB970228 Reichart99  and GRB980326 Bloom99 .
$``$ Iron lines: have been observed in some x-ray afterglows Piro00 . Any model explaining them requires a significant amounts of iron at rest near those GRBs.
$``$ Association with Star formation: GRBs seem to follow the star formation rate. GRB are located within star forming regions in star forming Galaxies Pac98 ; Kulkarni01 .
$``$ GRB distribution: GRBs are distributed within galaxies. There is no evidence for GRBs kicked out of their host galaxies Bloom01 ; Kulkarni01  as would be expected for NS-NS mergers NPP92 .

All these clues point out towards a SN/GRB association and towards the Collapsar model. However, the situation is not clear cut. The association of GRB980425 with SN98bw is uncertain. There are alternative explanations to the bumps in the afterglows of GRB970228 and GRB980326 Esin00 . Iron is produces in Supernovae. But there is no simple explanation what is iron at rest doing around the GRB (see however, Vietri01 ). The association with star formation and the distribution of GRBs within galaxies is real but all that it indicates is short lived progenitors. One cannot rule out a short lived binary NS population TutukovY94  which would mimic this behaviour. Even worse, there are some indication that seem incompatible with the SN association:
$``$ No Windy Afterglow: No evidence for a wind (n=2) in any of the afterglow light curves? Furthermore, most fits for the afterglow parameters show low ambient density PanaitescuK01 ; Kulkarni01 .
$``$ No Jets: Some GRBs don t show evidence for a jet or have very wide opening angles PanaitescuK01 ; Kulkarni01 , this would be incompatible with the Collapsar model.

## IV Conclusions, Predictions and Open Questions

There is an ample observational support For the Fireball model. It also has several other predictions. The most interesting ones are those concerning the very early afterglow and the GRB-afterglow transition. The early afterglow phase is radiative and a detailed look at the first hour of the afterglow should show the radiative to adiabatic transition. It should also show (mostly in radio) small bumps corresponding to refreshed shocks KP00a  which would enable us to learn more on the nature of the flow produced by the “inner engine”. With an operational HETE II and Swift in the not too distant future we hope that this crucial phase will be explored soon. Another prediction of a ring structure of the afterglowGPS99c  will have to wait, however, to futuristic ultrahigh resolution detectors.

We know how GRBs are produced. We are less certain what produces them. We can trace backwards the evolution at the source from the observations of the emitting regions to an accretion disk - black hole system. The traces from this point backwards are less clear. Theoretical considerations NarayanPK01  suggest that only Collapsars can produce the disk-black hole systems needed for long bursts while only NS-NS (or possibly NS-BH) mergers can produce the systems needed for short bursts. These conclusions are supported by the afterglows observations that suggest SN/GRB association for the long burst population. However, the picture is far from clear yet. While the information on the location of the bursts points out towards the SN connection the physical conditions within the afterglow indicate a low circumstellar density and does not show any indication for the almost inevitable pre explosion wind. The origin of the iron lines is still mysterious and confusing.

The Fireball model has still many open questions. Some are concerned with the physics of the fireball model: How do the collisionless shocks work? How are the electrons accelerated and how are the magnetic fields amplified? How do jets expand sideways? What controls the “typical” emission to be in the soft $`\gamma `$-ray region? Other questions deal with more astrophysical issues like: What happens in all the cases (like WD-BH merger) in which a GRB almost form but the conditions are not exactly right? What distinguishes between the prgenitor of a GRB-“failed supernova” and the progenitor of a successful supernova with no GRB? Finally, we have, of course, the sixty four thousand dollars question: How does the “inner engine” accelerates the ejecta to relativistic velocities?.


format: text

text:
# Basis Optimization Renormalization Group for Quantum Hamiltonian

## 1 INTRODUCTION

A method of Hamiltonian diagonalization is suitable for calculating physical quantities associated with wavefunctions, such as structure functions and form factors. However, we cannot diagonalize Hamiltonian directly in quantum field theories because dimension of Hilbert space is generally infinite. To create effective Hamiltonian, we extend a technique of NRG (numerical renormalization group) proposed by K. Wilson .

In spin chain models, any state can be expanded with a direct product of two sets of basis states, each of which is for one of two spin blocks. For a finite lattice, we can calculate wavefunction $`\mathrm{\Psi }_{ij}`$ of the ground state by diagonalizing Hamiltonian. In Ref. , S.White has found that singular values $`D_k`$ of target-state wavefunction $`\mathrm{\Psi }_{ij}=_kU_{ik}D_kV_{kj}`$ are useful for the purpose of creating effective Hamiltonian.

$$|\mathrm{\Psi }=\underset{i}{}\underset{j}{}\mathrm{\Psi }_{ij}|i|j=\underset{k}{}D_k|u_k|v_k,$$
(1)
where $`|u_k_iU_{ik}|i`$ and $`|v_k_jV_{kj}|j`$. If $`|D_k|`$ is large, the corresponding basis state is relevant. However, if $`|D_k|`$ is small, the corresponding basis state is not relevant and can be thrown away from calculation. We can control calculation accuracy using singular values $`D_k`$ when we truncate Hilbert space. Magnitude of $`|D_k|`$ determines the relevance of basis states in the target state. Basis truncation based on singular values works well in one-dimensional spin models such as Heisenberg and Hubbard models . This is an extension of Wilson’s NRG and called DMRG (density matrix renormalization group). White has shown that there exists orthogonal transformation that largely reduces the number of basis states with high calculation accuracy. However, we cannot assure that a number of singular values take values of nearly zero in general cases. In order to create optimized basis states without relying on density matrix, we will find another useful basis transformation. Our aim is to find basis transformation that make coefficients (i.e. transformed wavefunctions) exactly zero when target states are expanded with new basis states.

## 2 FORMULATION

Let us consider $`S=1/2`$ Heisenberg chain. Since we cannot diagonalize infinite dimensional Hamiltonian directly, we start from a finite lattice. We extend a lattice step by step by adding one site to it in each RG transformation. $`|u_i`$ is a complete set of basis states for a finite lattice $`L`$. Hilbert space of an extended lattice with size ($`L+1`$) is spanned with a set of basis states $`|u_i|j`$, where $`|j`$ is an eigenstate of $`z`$-component of a spin operator for the ($`L+1`$)-th site

$$\widehat{S}^z|j=\pm \frac{1}{2}|j\text{for}j=\pm .$$

We want to find a rotation matrix $`U`$ that decrease the number of relevant basis states.

$$|u_j^{}=\underset{i}{}U_{ij}|u_i.$$
(2)
In general, any orthogonal matrix can be expressed as a product of two-dimensional rotation matrices $`R`$.

$$R_{(k,l)}(\theta )=\begin{array}{cccccccc}& & & k& & l& & \\ & (\mathrm{}& & \mathrm{}& & \mathrm{}& & )\mathrm{}\\ k& \mathrm{}& \mathrm{cos}\theta & \mathrm{}& \mathrm{sin}\theta & \mathrm{}\\ & & \mathrm{}& \mathrm{}& \mathrm{}& \\ l& \mathrm{}& \mathrm{sin}\theta & \mathrm{}& \mathrm{cos}\theta & \mathrm{}\\ & & \mathrm{}& & \mathrm{}& \end{array}$$

This matrix rotates the $`k`$-th and $`l`$-th basis states by an angle $`\theta `$. Let us consider a ground state that is an eigenstate of a Hamiltonian.

$$|\mathrm{\Psi }=\underset{ij}{}\mathrm{\Psi }_{ij}|u_i|j.$$
(3)
We consider a part of this target state.

$`|\mathrm{\Psi }_{(1,2)}=\mathrm{\Psi }_1|u_1|+\mathrm{\Psi }_{1+}|u_1|+`$
$`+\mathrm{\Psi }_2|u_2|+\mathrm{\Psi }_{2+}|u_2|+.`$
If we rotate two basis states $`|u_1`$ and $`|u_2`$

$`|u_1=\mathrm{cos}\theta |u_1^{}\mathrm{sin}\theta |u_2^{},`$
$`|u_2=\mathrm{sin}\theta |u_1^{}+\mathrm{cos}\theta |u_2^{},`$
we have

$`|\mathrm{\Psi }_{(1,2)}=\mathrm{\Psi }_1^{}|u_1^{}|+\mathrm{\Psi }_{1+}^{}|u_1^{}|+`$
$`+\mathrm{\Psi }_2^{}|u_2^{}|+\mathrm{\Psi }_{2+}^{}|u_2^{}|+,`$
where $`\mathrm{\Psi }^{}`$ are transformed wavefunctions

$`\mathrm{\Psi }_1^{}=\mathrm{\Psi }_1\mathrm{cos}\theta +\mathrm{\Psi }_2\mathrm{sin}\theta ,`$
$`\mathrm{\Psi }_{1+}^{}=\mathrm{\Psi }_{1+}\mathrm{cos}\theta +\mathrm{\Psi }_{2+}\mathrm{sin}\theta ,`$
$`\mathrm{\Psi }_2^{}=\mathrm{\Psi }_1\mathrm{sin}\theta +\mathrm{\Psi }_2\mathrm{cos}\theta ,`$
$`\mathrm{\Psi }_{2+}^{}=\mathrm{\Psi }_{1+}\mathrm{sin}\theta +\mathrm{\Psi }_{2+}\mathrm{cos}\theta .`$
We define the following two types of rotations. Rotation I is defined to make the fourth component zero $`\mathrm{\Psi }_{2+}^{}=0`$;

$$\left(\begin{array}{c}\mathrm{\Psi }_1\\ \mathrm{\Psi }_{1+}\\ \mathrm{\Psi }_2\\ \mathrm{\Psi }_{2+}\end{array}\right)\stackrel{\mathrm{I}}{}\left(\begin{array}{c}\mathrm{\Psi }_1^{}\\ \mathrm{\Psi }_{1+}^{}\\ \mathrm{\Psi }_2^{}\\ 0\end{array}\right)$$

A condition $`\mathrm{\Psi }_{2+}^{}=0`$ determines the value of $`\theta `$. Rotation II is defined to make the third component zero $`\mathrm{\Psi }_2=0`$ when components $`\mathrm{\Psi }_{1+}`$ and $`\mathrm{\Psi }_{2+}`$ are zero;

$$\left(\begin{array}{c}\mathrm{\Psi }_1\\ 0\\ \mathrm{\Psi }_2\\ 0\end{array}\right)\stackrel{\mathrm{II}}{}\left(\begin{array}{c}\mathrm{\Psi }_1^{}\\ 0\\ 0\\ 0\end{array}\right)$$

When $`\mathrm{\Psi }_{1+}=\mathrm{\Psi }_{2+}=0`$, the transformed components $`\mathrm{\Psi }_{1+}^{}`$ and $`\mathrm{\Psi }_{2+}^{}`$ are always zero. Therefore, we can make $`\mathrm{\Psi }_2^{}`$ zero by choosing $`\theta `$ appropriately.

In order to understand how these two rotations are used to find relevant basis states, let us consider a simple example where a target state is expanded with six basis states.

$$|\mathrm{\Psi }=\underset{i=1}{\overset{3}{}}\underset{j=\pm }{}\mathrm{\Psi }_{ij}|u_i|j.$$
(4)
The following is transformation of the six-dimensional vector of the target-state wavefunction $`\mathrm{\Psi }_{ij}`$.

$$\left(\begin{array}{c}\mathrm{\Psi }_1\\ \mathrm{\Psi }_{1+}\\ \mathrm{\Psi }_2\\ \mathrm{\Psi }_{2+}\\ \mathrm{\Psi }_3\\ \mathrm{\Psi }_{3+}\end{array}\right)\stackrel{\mathrm{I}}{}\left(\begin{array}{c}\mathrm{\Psi }_1^{}\\ \mathrm{\Psi }_{1+}^{}\\ \mathrm{\Psi }_2^{}\\ 0\\ \mathrm{\Psi }_3\\ \mathrm{\Psi }_{3+}\end{array}\right)\stackrel{\mathrm{I}}{}\left(\begin{array}{c}\mathrm{\Psi }_1^{\prime \prime }\\ \mathrm{\Psi }_{1+}^{\prime \prime }\\ \mathrm{\Psi }_2^{}\\ 0\\ \mathrm{\Psi }_3^{}\\ 0\end{array}\right)\stackrel{\mathrm{II}}{}\left(\begin{array}{c}\mathrm{\Psi }_1^{\prime \prime }\\ \mathrm{\Psi }_{1+}^{\prime \prime }\\ \mathrm{\Psi }_2^{\prime \prime }\\ 0\\ 0\\ 0\end{array}\right)$$

We obtain the final form by transforming the vector with two rotations appropriately. If a transformed wavefunction is zero, the corresponding basis state is irrelevant and can be neglected. That is, only the first three basis states ($`|u_1^{\prime \prime }|\pm `$, $`|u_2^{\prime \prime }|`$) are relevant and others ($`|u_2^{\prime \prime }|+`$, $`|u_3^{\prime \prime }|\pm `$) are redundant. The former three relevant basis states can reproduce the target state exactly without other basis states. There exists orthogonal transformation that reduces the number of relevant basis states into three without causing calculation error. This is also true for arbitrary dimensions other than six of this example.

Let us consider a case where we target ground and first-excited states. The followings are wavefunction vectors of the target states.

$$\left(\begin{array}{cc}\mathrm{\Psi }_1^{(1)}& \mathrm{\Psi }_1^{(2)}\\ \mathrm{\Psi }_{1+}^{(1)}& \mathrm{\Psi }_{1+}^{(2)}\\ \mathrm{\Psi }_2^{(1)}& \mathrm{\Psi }_2^{(2)}\\ \mathrm{\Psi }_{2+}^{(1)}& \mathrm{\Psi }_{2+}^{(2)}\\ \mathrm{\Psi }_3^{(1)}& \mathrm{\Psi }_3^{(2)}\\ \mathrm{\Psi }_{3+}^{(1)}& \mathrm{\Psi }_{3+}^{(2)}\\ \mathrm{\Psi }_4^{(1)}& \mathrm{\Psi }_4^{(2)}\\ \mathrm{\Psi }_{4+}^{(1)}& \mathrm{\Psi }_{4+}^{(2)}\\ \mathrm{}& \mathrm{}\\ \mathrm{\Psi }_N^{(1)}& \mathrm{\Psi }_N^{(2)}\\ \mathrm{\Psi }_{N+}^{(1)}& \mathrm{\Psi }_{N+}^{(2)}\end{array}\right)\left(\begin{array}{cc}\mathrm{\Psi }_1^{(1)^{}}& \mathrm{\Psi }_1^{(2)^{}}\\ \mathrm{\Psi }_{1+}^{(1)^{}}& \mathrm{\Psi }_{1+}^{(2)^{}}\\ \mathrm{\Psi }_2^{(1)^{}}& \mathrm{\Psi }_2^{(2)^{}}\\ 0& \mathrm{\Psi }_{2+}^{(2)^{}}\\ 0& \mathrm{\Psi }_3^{(2)^{}}\\ 0& \mathrm{\Psi }_{3+}^{(2)^{}}\\ 0& \mathrm{\Psi }_4^{(2)^{}}\\ 0& \mathrm{\Psi }_{4+}^{(2)^{}}\\ \mathrm{}& \mathrm{}\\ 0& \mathrm{\Psi }_N^{(2)^{}}\\ 0& \mathrm{\Psi }_{N+}^{(2)^{}}\end{array}\right)\left(\begin{array}{cc}\mathrm{\Psi }_1^{(1)^{}}& \mathrm{\Psi }_1^{(2)^{}}\\ \mathrm{\Psi }_{1+}^{(1)^{}}& \mathrm{\Psi }_{1+}^{(2)^{}}\\ \mathrm{\Psi }_2^{(1)^{}}& \mathrm{\Psi }_2^{(2)^{}}\\ 0& \mathrm{\Psi }_{2+}^{(2)^{}}\\ 0& \mathrm{\Psi }_3^{(2)^{\prime \prime }}\\ 0& \mathrm{\Psi }_{3+}^{(2)^{\prime \prime }}\\ 0& \mathrm{\Psi }_4^{(2)^{\prime \prime }}\\ 0& 0\\ \mathrm{}& \mathrm{}\\ 0& 0\\ 0& 0\end{array}\right)$$

The indices $`(1)`$ and $`(2)`$ indicate the ground and first-excited states, respectively. We optimize basis states in two steps. In the first step, we target the ground state $`\mathrm{\Psi }^{(1)}`$ and do not take care of the first excited state $`\mathrm{\Psi }^{(2)}`$. In the same way as the previous example (4), we can find three relevant basis states that can reproduce the ground state exactly by using a combination of rotations I and II. In the second step, we target the first excited state $`\mathrm{\Psi }^{(2)}`$ to optimize basis states further. When we rotate basis states, we don’t touch the first four basis states ($`|u_1^{}|\pm `$, $`|u_2^{}|\pm `$). We apply a combination of rotations I and II to other ($`2N4`$) pieces of basis states ($`|u_3^{}|\pm `$, $`|u_4^{}|\pm `$, $`\mathrm{}`$, $`|u_N^{}|\pm `$). Finally, we obtain seven basis states ($`|u_1^{\prime \prime }|\pm `$, $`|u_2^{\prime \prime }|\pm `$, $`|u_3^{\prime \prime }|\pm `$, $`|u_4^{\prime \prime }|`$) that are relevant for reproducing the ground and first excited states exactly. In the same way, we can also create relevant basis states when arbitrary number of eigenstates are targeted. The number of target states can be chosen according to the desired calculation accuracy. We can use this technique of basis transformation to truncate Hilbert space and create effective Hamiltonian in NRG. Our NRG algorithm for $`S=1/2`$ Heisenberg chain is as follows:

1. Diagonalize $`H_L`$ ($`H_L𝒖_n=E_n𝒖_n`$), which is Hamiltonian for a finite lattice $`L`$.
2. Transform basis states $`𝒖_n𝒖_n^{}`$ using type I and II rotations to make wavefunctions of target states zero as many as possible.
3. Create effective Hamiltonian $`\overline{H}_L=O_L^{}H_LO_L`$, where $`O_L=(𝒖_1^{},𝒖_2^{},\mathrm{},𝒖_M^{})`$. Dimension of $`\overline{H}_L`$ is $`M`$.
4. Add one site to $`\overline{H}_L`$ to provide next Hamiltonian $`H_{L+1}`$. Dimension of $`H_{L+1}`$ is $`2M`$.
5. Regard $`H_{L+1}`$ as $`H_L`$, and go to step 1.

These steps are iterated till vacuum energy per site converges with the desired accuracy. I have applied this algorithm to $`S=1/2`$ Heisenberg chain and calculated vacuum energy per site. An exact value of vacuum energy has been obtained using Bethe ansatz. The difference between our result and exact value is about $`3.9\times 10^4`$ for parameters $`M=256`$ and $`L=2000`$.

## 3 CONCLUSION

We have found basis transformation that largely reduces the number of relevant basis states. Especially, three states are sufficient for reproducing one target state in each RG step. Relevant basis states reproduce target states exactly in each RG step. Therefore, if dimension of effective Hamiltonian $`M`$ is sufficiently large, calculation errors come from non-targeted states. Future problems are improvement of accuracy, calculation of correlation functions, calculation of Haldane gap, and application to other models.


text:
# Criticality in one dimension with inverse square-law potentials

## Abstract

It is demonstrated that the scaled order parameter for ferromagnetic Ising and three-state Potts chains with inverse-square interactions exhibits a universal critical jump, in analogy with the superfluid density in helium films. Renormalization-group arguments are combined with numerical simulations of systems containing up to one million lattice sites to accurately determine the critical properties of these models. In strong contrast with earlier work, compelling quantitative evidence for the Kosterlitz–Thouless-like character of the phase transition is provided.

One-dimensional scalar models with inverse-square potentials have held a special position in the field of statistical mechanics for several decades. The surprisingly rich behavior emerging from these conceptually simple systems has attracted the attention of a variety of workers in the field. The work of Ruelle  and Dyson  proved that pair interactions decaying as $`1/r^2`$ (where $`r`$ denotes the distance between two local variables) represent the boundary between models with and without a phase transition at nonzero temperature; but it took until the 1980s before Fröhlich and Spencer  proved rigorously the existence of a phase transition for the borderline case itself. Meanwhile, early work by Anderson and Yuval  (who were motivated by the connection to the Kondo problem) had suggested that this phase transition has a special character, with an order parameter that exhibits both a discontinuity and a superposed singularity at the critical point. This behavior later gained widespread attention through the Kosterlitz–Thouless (KT) transition in the two-dimensional (2D) $`XY`$ model . Indeed, the vortex-unbinding scenario in the latter case has a direct analog in terms of logarithmically interacting defects (domain walls) in the “kink-gas” representation of the Ising chain considered in . This treatment in the spirit of the renormalization-group (RG) theory was extended in  and most notably by Cardy , who demonstrated that it can be generalized to large classes of discrete models. In this sense, these one-dimensional models can be considered as the simplest examples of the class of “topological” phase transitions, which comprises phenomena like 2D superfluidity and superconductivity, as well as dislocation-mediated melting of crystals.

An exact lower bound on the critical order-parameter jump has been obtained by Aizenman et al.  for Ising and Potts chains with $`1/r^2`$ interactions, and the existence of an intermediate ordered phase, where the correlation function decays like a temperature-dependent power-law, has been proven by Imbrie and Newman . Despite these rigorous results the literature on these models contains various contradictory claims and unconfirmed conjectures, while few of the RG-based predictions for the critical properties have been verified by alternative methods: Thus (i) the critical temperature $`T_c`$ has been estimated by a great variety of techniques, but little consistency has been obtained until now, with estimates differing by up to 30%. Indeed, there is even a conjecture for $`T_c`$ of the $`q`$-state Potts chain that is independent of $`q`$ , despite the fact that this appears to contradict intuitive energy–entropy arguments. (ii) In contrast, it has been claimed on the basis of numerical simulations that the transition in the $`q`$-state Potts chain is continuous for $`qq_c`$ and of first order for $`q>q_c`$ , although this in turn disagrees with the KT-like transition (with $`q`$-dependent exponents) found by Cardy . (iii) It has been conjectured that the power-law exponent of the correlation function in the intermediate ordered phase is directly related to the value of the order parameter . (iv) Numerical simulations have been unable to conclude whether the specific heat diverges at $`T_c`$ or even whether it displays a maximum at $`T_c`$ . (v) No numerical work for a one-dimensional system has been able to reproduce one of the actual hallmarks of the KT universality class, namely the predicted exponential divergence of the correlation length and the susceptibility for $`TT_c`$.

It is the aim of the present work to settle these open questions in a convincing way, although (ii) will be addressed only partially, by combining finite-size scaling relations obtained from the RG predictions (which until now have been ignored in nearly all numerical studies) with accurate numerical simulations that span an unprecedented range in system sizes. Indeed, the long-range nature of the interactions has posed a major hurdle for the numerical study of these models: Whereas almost two decades ago Monte Carlo (MC) simulations could reach a system of $`L=256`$ lattice sites , the most recent work still could only attain $`L=900`$ . By using appropriate techniques, we have now been able to study systems as large as $`L=10^6`$. This allows us not only to answer questions intrinsic to the one-dimensional systems under study, but also to obtain what we believe are the most accurate numerical results to date for a system exhibiting a KT-like phase transition: Indeed, although the 2D $`XY`$ and Villain models have been examined in many numerical studies, it has proven extremely difficult to obtain full consistency with the RG predictions . Because of the exponential divergence of the correlation length $`\xi `$, it is essential to cover a large range in linear system sizes, which makes a one-dimensional system easier to study than a 2D one. In our analysis, we start from the hypothesis that the structure of the scaling functions is correctly described by RG theory, and examine the consistency between our numerical results and the scaling predictions to test the correctness of this hypothesis.

We have first focused attention on the Ising chain (or equivalently the $`q=2`$ Potts chain), described by the reduced Hamiltonian $`\beta =K_{ij}s_is_j/r_{ij}^2`$, where the sum runs over all spin pairs, $`s_i=\pm 1`$, and $`r_{ij}=|𝐫_i𝐫_j|`$. We have carried out accurate MC simulations for chains of $`L`$ lattice sites, with $`10L10^6`$, over a large temperature range. These results could be obtained owing to the use of a cluster algorithm that suppresses critical slowing down and has an efficiency that is independent of the interaction range .

Although the critical coupling $`K_c`$ is a nonuniversal quantity, the large number of available estimates makes it interesting in its own right. In addition, its precise value proves indispensable if one wants to address the issues mentioned above. Since finite-size corrections decay only logarithmically near a KT transition, the determination of $`K_c`$ is notoriously difficult, and a new approach had been found. To this end, we have relied on the observation that the properly scaled order parameter $`\mathrm{\Psi }Km^2`$, with $`m=L^1_{i=1}^Ls_i`$, is expected to exhibit a universal jump at the critical point. Indeed, as mentioned, it has been shown rigorously  that $`\mathrm{\Psi }(K>K_c)\frac{1}{2}`$ \[whereas $`\mathrm{\Psi }(K<K_c)=0`$\]. Furthermore, it was conjectured in  that this inequality is saturated at $`K=K_c`$ (which seems to be implied by the arguments of , but note the comments in ). While it can be easily shown that $`\mathrm{\Psi }`$ is invariant under a (real-space) RG transformation of the kink-gas representation of the Ising chain , it appears not to have been explicitly realized that the critical jump $`\mathrm{\Psi }_c\mathrm{\Psi }(K_c)`$ must be universal, in the sense that it will have the same magnitude in every system with couplings that decay asymptotically like $`r^2`$. We remark that also for the analogous phenomenon in two dimensions (the celebrated universal jump in the superfluid density of a 4He film) this universality was only realized  after the original work by Kosterlitz and Thouless . In addition to the jump discontinuity, $`\mathrm{\Psi }(K>K_c)`$ is predicted to exhibit a superposed singularity:

$$\mathrm{\Psi }(K)=\mathrm{\Psi }_c+a_1(KK_c)^{\stackrel{~}{\nu }}+a_2(KK_c)^{2\stackrel{~}{\nu }}+\mathrm{},$$
(1)
with $`\stackrel{~}{\nu }=\frac{1}{2}`$ for the Ising chain; the ellipsis also includes nonsingular correction terms. In what follows, we demonstrate how to exploit the scaling properties of $`\mathrm{\Psi }`$ to accurately locate $`K_c`$. Figure 1 shows $`\mathrm{\Psi }(K,L)`$ as a function of $`K`$ for different system sizes. While $`\mathrm{\Psi }`$ is already sharpening up for relatively small $`L`$, it is only for very large systems that one observes the development of a discontinuity and an additional singularity.

To proceed further, we recall the lowest-order RG flow equations  for the coupling and the defect fugacity $`y`$, namely

$`dK/dl`$ $`=`$ $`4Ky^2,`$ (2)
$`dy/dl`$ $`=`$ $`y(2K1),`$ (3)
which (under the additional assumption that the system size scales with the rescaling parameter $`e^l`$) allow us to derive the finite-size scaling behavior of $`\mathrm{\Psi }(K,L)`$ in the low-temperature regime. Since this derivation is somewhat involved , we only quote the final result here, namely

$$\mathrm{\Psi }(K,L)=\mathrm{\Psi }(K,\mathrm{})\left[1+b_1L^{2(\overline{\mathrm{\Psi }}1)}+b_2L^{4(\overline{\mathrm{\Psi }}1)}+\mathrm{}\right],$$
(4)
where $`\overline{\mathrm{\Psi }}\mathrm{\Psi }(K,\mathrm{})/\mathrm{\Psi }_c`$. Nonsingular corrections become important only when $`\mathrm{\Psi }>2\mathrm{\Psi }_c`$. The appearance of $`\mathrm{\Psi }(K,\mathrm{})`$ both as an amplitude and as an exponent is typical for a KT transition. A least-squares fit of this equation to $`\mathrm{\Psi }(K,L)`$ has been used to obtain $`\mathrm{\Psi }(K,\mathrm{})`$ for several values of $`K`$, as shown in the inset in Fig. 1. This also yielded $`\mathrm{\Psi }_c=0.49\pm 0.02`$. Keeping $`\mathrm{\Psi }_c`$ fixed at the predicted value $`\frac{1}{2}`$  yielded $`K_c=0.6552(2)`$. Subsequently, Eq. (1) allowed us to determine $`\stackrel{~}{\nu }=0.52\pm 0.03`$. We have supplemented this low-temperature approach with an analysis of data lying within the high-temperature finite-size regime, where $`L\xi `$ and $`K<K_c`$. The algebraic corrections in (4) now become powers of $`1/\mathrm{ln}L`$ and temperature-dependent terms take the form $`L/\xi `$, with $`\xi =\mathrm{exp}[B/(KK_c)^{\stackrel{~}{\nu }}]`$ . We then find $`\mathrm{\Psi }_c=0.496(3)`$, in excellent agreement with the RG prediction, and $`K_c=0.6548(14)`$. Keeping $`\mathrm{\Psi }_c`$ now fixed at $`\frac{1}{2}`$ yields $`K_c=0.6555(4)`$, very close to the value resulting from the low-temperature approach.

Independent estimates for $`K_c`$ and $`\stackrel{~}{\nu }`$ can be obtained from the high-temperature susceptibility, which should display the same exponential divergence as the correlation length. RG considerations regarding higher-order contributions  suggest $`\chi \mathrm{exp}[B_0(KK_c)^{\stackrel{~}{\nu }}+B_2(KK_c)^{\stackrel{~}{\nu }}+𝒪((KK_c)^{2\stackrel{~}{\nu }})]`$. Numerical results for $`L=410^5`$ and $`L=10^6`$ are shown in Fig. 2, along with a truncated high-temperature expansion of order $`𝒪(K^8)`$  and a fit to the RG expression. The agreement with the series expansion is good, but breaks down for $`K0.3`$. The RG expression provides an excellent description of the MC data for $`0.2K0.61`$, confirming the predicted exponential divergence over a large temperature range. Owing to the covariance of the fit parameters, it is difficult to determine them simultaneously: We obtained $`\stackrel{~}{\nu }=0.54\pm 0.03`$ as our best estimate. However, keeping $`\stackrel{~}{\nu }=\frac{1}{2}`$ leads to $`K_c=0.6551(6)`$, which we view as a strong confirmation of our analysis of $`\mathrm{\Psi }`$. The rounding of $`\chi `$ close to $`K_c`$ (see Fig. 2) is due to the finite system size. We estimate that the correlation length $`\xi `$ is $`10^5`$ here, which may be compared with a maximum of $`\xi 140`$ reached in recent $`d=2`$ simulations .

Subsequently, we carried out a similar analysis of $`\mathrm{\Psi }`$ and $`\chi `$ for the $`q`$-state Potts model with $`q=3`$, which has a Hamiltonian of the form $`\beta ^{\mathrm{Potts}}=K^{\mathrm{Potts}}_{ij}\delta _{s_is_j}/r_{ij}^2`$, with $`s_i=1,\mathrm{},q`$. Cardy  has found that, for positive integer $`q`$, this system can be described by RG equations that are very similar to those for the $`q=2`$ (Ising) case, but yielding $`\stackrel{~}{\nu }=2/(q+2)`$. This implies that the entire analysis of the critical behavior can be carried through as before : $`\mathrm{\Psi }`$ should still exhibit a universal jump (of magnitude $`1`$ in units of $`K^{\mathrm{Potts}}`$) with a superposed singularity as in (1) and $`\chi `$ should again display an exponential singularity. The order parameter used for the determination of $`\mathrm{\Psi }`$ and $`\chi `$ is given by $`m^{\mathrm{Potts}}\mathrm{max}_{\alpha =1,\mathrm{},q}[\frac{q}{L}_{i=1}^L(\delta _{s_i\alpha }1)]/(q1)`$, with $`q>1`$. An analysis of $`\mathrm{\Psi }`$ for system sizes up to $`L=10^6`$ yielded $`K_c=1.4105(2)`$, $`\mathrm{\Psi }_c=1.02(2)`$, and $`\stackrel{~}{\nu }(q=3)=0.41(2)`$. A least-squares fit of the data for $`\chi `$ was consistent with an exponential divergence with $`\stackrel{~}{\nu }(q=3)=0.42(2)`$. Not only is $`\mathrm{\Psi }_c`$ in good agreement with the predicted universal value, but to our knowledge this is also the very first independent corroboration of the generalized RG scenario of Ref. , which predicts that $`\stackrel{~}{\nu }(q=3)=\frac{2}{5}`$. As expected from entropic considerations, $`K_c^{\mathrm{Potts}}(q=3)`$ is somewhat higher than $`K_c^{\mathrm{Potts}}(q=2)=1.3104(4)`$, which in turn refutes two conjectures based on a different real-space RG analysis , namely that (i) $`K_c^{\mathrm{Potts}}`$ would be independent of $`q`$ and (ii) $`K_c^{\mathrm{Potts}}`$ would be exactly given by $`12/\pi ^2=1.21585\mathrm{}`$ We note that our results are in strong contrast with earlier numerical work: Both a transfer-matrix study and a very recent MC study for $`q=2`$  and $`q=3`$  concluded that the correlation length diverges algebraically (rather than exponentially) with a $`q`$-dependent power-law exponent.

We now return to some interesting aspects of the Ising chain. First, it has been concluded in Ref.  that the critical spin-spin correlation function should decay logarithmically, $`g(r)=s_0s_rs_0^21/\mathrm{ln}r`$, where now $`s_00`$, unlike the situation for a typical critical phase transition. As illustrated in Fig. 3, we have been able to corroborate this behavior over a considerable range in $`r`$. In the high-temperature phase, the decay of $`g(r)`$ is bounded by the spin-spin interaction, but for $`K>K_c`$ it is expected that $`g(r)r^\theta `$, with $`\theta =4\sqrt{(KK_c)/K_c}`$ . We note that this behavior would be consistent with the conjecture of Ref. , viz. that $`\theta =\mathrm{min}[4(\mathrm{\Psi }\frac{1}{2}),2]`$ for $`K>K_c`$. We have determined $`\theta `$ for a large number of couplings, see Fig. 4, and have subsequently fitted these to a power series in $`\sqrt{KK_c}`$, finding $`\theta =(1.81\pm 0.06)\sqrt{(KK_c)/K_c}`$. This agrees with $`4a_1\sqrt{K_c}=1.73(4)`$, where $`a_1`$ is defined via Eq. (1), corroborating the conjectured relation with $`\mathrm{\Psi }`$. However, the prefactor is approximately half as large as predicted in Ref. . Farther away from $`K_c`$ it becomes increasingly difficult to obtain accurate results for $`\theta (K)`$, owing to the rapid decay of $`g(r)`$. However, since our results corroborate the presumed relation between $`\theta `$ and $`\mathrm{\Psi }`$, we exploit this relation to determine two special values of $`K`$, namely the coupling strength $`K_1`$ for which $`\theta =1`$ and the smallest coupling $`K_2>K_c`$ for which $`\theta =2`$. $`K_1`$ has the special property that the magnetic susceptibility will diverge over the full coupling range $`K_cKK_1`$, as follows from the fluctuation–dissipation theorem (and as is actually confirmed by the numerical results ). We find $`K_1=0.78(1)`$ and $`K_2=1.02(3)`$. Thus, the RG prediction $`K_1=\frac{17}{16}K_c0.696`$  is too low, due to the difference in the prefactor for $`\theta (K)`$. $`K_2`$ marks the end of the regime where $`g(r)`$ displays a leading temperature-dependent exponent.

Finally, we consider the specific heat $`C`$, which is calculated from the fluctuations in the internal energy. The efficiency of the original MC algorithm  is greatly reduced when this energy has to be calculated (although recently this problem has largely been overcome ) and hence only much smaller systems could be accessed. Nevertheless, this turns out not to be a limiting factor, since $`C`$ has already converged within the range of accessible system sizes: Fig. 5 shows unambiguously that $`C`$ is a nondivergent quantity which takes its maximum at a temperature that lies approximately 14% above the critical temperature.

In summary, we have provided compelling evidence that the inverse-square Ising and $`q=3`$ Potts chain display the one-dimensional analog of a Kosterlitz–Thouless transition, in disagreement with earlier numerical studies, but in accurate quantitative agreement with renormalization-group predictions. The (exponentially diverging) correlation length could be examined over more than five orders of magnitude, making this the most tractable KT-like transition.

It is a great pleasure to acknowledge stimulating remarks by Prof. Kurt Binder and helpful comments by Prof. Michael E. Fisher. E.L. gratefully acknowledges the hospitality of the Condensed Matter Theory Group at the Johannes Gutenberg-Universität Mainz, where most of this work has been carried out, and support from the National Science Foundation (through Grant No. CHE 99-81772 to M.E. Fisher) and from the Department of Energy, Office of Basic Energy Sciences (through Grant No. DE-FG02-98ER14858 to A.Z. Panagiotopoulos).


text:
# 56Ni Mixing in the Outer Layers of SN 1987A

## 1 Introduction

Until the explosion of Supernova 1987A in the LMC, computer models of supernovae were usually calculated with assumptions that allowed only the prediction of bolometric light curves. The matter temperature was usually assumed to be equal to the radiation temperature, and grey opacities were used (Arnett et al., 1989; Blinnikov, 1999). Often, the models did not account for the effects of radiation on the hydrodynamic gas flow, or the effects of sphericity on the spectrum (Schmutz et al., 1990).

Observations of SN 1987A produced extensive broad-band photometry, especially during its earliest phase of expansion. This motivated the development of models to predict its spectral and photometric evolution. In our analysis here we make use of the optical data obtained at CTIO (Phillips et al., 1988) and the extensive IUE data which was obtained, re-reduced and analyzed by Pun et al. (1995).

While SN 1987A led to the growth of observational data, so too did it lead to an improvement in theoretical models. Höflich (1991) modeled the early phase of SN 1987A, using a pure hydrogen atmosphere and accounting for non-LTE (Local Thermodynamic Equilibrium). Eastman & Kirshner (1989) modeled the first ten days of the SN 1987A explosion. H I and He I were modeled in NLTE, while metal lines were treated as pure scattering LTE lines. Takeda (1991) used a pure H/He envelope in NLTE and attempted to understand the lineshapes of the Balmer lines.

Our work here utilized a hydrodynamic model developed by Blinnikov and collaborators using the STELLA software package (Blinnikov et al., 1998; Blinnikov, 1999; Blinnikov et al., 2000). This model calculated the light curve for up to six months after the explosion, and includes allowance for time-dependent, multi-group radiation hydrodynamics, monochromatic scattering effects, and the effects of spectral lines on the opacity (Blinnikov et al., 2000). The hydrodynamical models were based on the stellar evolution models of Saio, Nomoto, & Kato (1988), Saio, Kato, & Nomoto (1988), and Nomoto & Hashimoto (1988). The procedure for synthesizing SN 1987A spectra from Blinnikov et al.’s model and the results are described in § 2. Here we concentrate on the fact that, by day 4.52, significant discrepancies begin to appear in the synthetic spectra of SN 1987A with respect to the observed spectra. In § 3 we review earlier works which noted this discrepancy. In § 4 we show that nickel mixing at early times seems to be the only way to make the synthetic spectra agree with observations. In a subsequent paper, we will discuss NLTE effects with regard to the results of Blinnikov et al. (2000).

## 2 Models

### 2.1 PHOENIX

PHOENIX is a general radiation transfer code that computes temperature, opacity, and level populations for each of typically 50–100 radial zones in a moving stellar envelope. PHOENIX allows the user to calculate level populations for a multitude of different nuclear species in LTE or NLTE (Hauschildt & Baron, 1999; Short et al., 1999). For this study, the following species were calculated in NLTE: H I, He I, He II, O I, Ne I, Na I, Mg II, Si II, S II, Ca II, Fe I, Fe II, and Fe III.

The Blinnikov et al. hydrodynamic model consists of an ejecta envelope divided into 300 zones, with compositions defined for each layer (Blinnikov, 1999; Blinnikov et al., 2000). The outer layers of this model were re-zoned onto the 75-zone grid that is used by a modified version of PHOENIX, preserving the compositions in velocity-space. The inner boundary condition was taken to be diffusive, which is well justified in these calculations. The gamma-ray energy deposition rate, which is the energy in gamma-rays deposited in the matter in the envelope per unit mass per unit time, is calculated for each new zone. The new grid was chosen so as to maintain both composition and density profiles with adequate resolution.

### 2.2 Spectral Synthesis

Figures 13 show the spectra resulting from PHOENIX computations on days 1.36, 2.67, and 3.59, respectively, of SN 1987A. The fits between the IUE spectra and the synthetic UV spectra are very good, an indicator of the importance of NLTE effects in the supernova envelope. \[Blinnikov et al. (2000) noted that the predicted UV flux 5–10 days after shock breakout was much greater than the actual flux when computed in LTE.\] For days 1.36 and 2.67, the agreement between the synthetic and observed optical spectra are also reasonably good. However, by day 3.59 a discrepancy in the strengths of the Balmer lines begins to appear in the synthetic spectra, most notably H$`\alpha `$, even when NLTE effects are included. Figure 4 shows the PHOENIX-generated optical/near-UV spectrum for day 4.52 compared with the observations (Phillips et al., 1988; Pun et al., 1995), and here the Balmer lines are clearly much weaker in the synthetic spectrum than they are in the observed spectrum.

## 3 Weak Lines

Schmutz et al. (1990), using various hydrodynamic models, noted the hydrogen atoms appeared to be sufficiently excited to create strong lines, despite the fact that the photospheric temperature had dropped below the ionization threshold ($`T6000K`$). They could offer no explanation for this discrepancy. The pure H/He non-LTE model of Takeda (1991) also resulted in weak Balmer lines, as did a model by Hauschildt & Ensman (1994), which used an earlier version of PHOENIX with H I, He I, Mg II, and Ca II in NLTE and metal line blanketing in LTE. These two groups suggested that differences between the theoretical and actual density distributions might be responsible.

Figure 5 shows temperature vs. velocity for Blinnikov et al.’s models at days 1.36 to 4.52. The diamonds on the graph show where $`\tau _{\text{std}}`$ is approximately equal to 1 in each model ($`\tau _{\text{std}}`$ total continuum optical depth at 5000 Å). Assuming the line-forming region lies in the vicinity of $`\tau _{\text{std}}=1`$, note that for days 1 and 2, the gas temperature is well above the ionization threshold for hydrogen, thus ensuring the presence of strong Balmer lines. (A good rule of thumb is that one sees the strongest lines from the species that is just below the dominant ionization stage. For example, Si II lines are strongest if Si III is the dominant stage, etc.) By day 3, the temperature is beginning to approach the threshold at $`\tau _{\text{std}}=1`$, and at day 4 the temperature is such that sufficient ionization can no longer be assumed in the line-forming region.

Figure 6 shows the gamma-ray energy deposition predicted by Blinnikov et al.’s hydrodynamic model. The total energy deposited per second is shown as a function of $`\tau _{\text{std}}`$ and the integral is performed from the outside proceeding inward. Until we reach the innermost layers of the envelope model ($`\tau _{\text{std}}>100`$), the deposition is negligible, and therefore hydrogen atoms cannot become sufficiently excited through the absorption of gamma rays.

## 4 Nickel Mixing

### 4.1 Model Results

In order for the synthetic Balmer lines to match the strength of the observed lines, either higher matter temperatures (to thermally populate the Balmer states) or higher radiation temperatures (to radiatively populate the Balmer states in NLTE) would be required. Both of these would significantly degrade the quality of the fit, since they would both lead to model spectra that would be significantly bluer than the observed spectrum. A third way to strengthen the Balmer lines is to populate the states non-thermally, via fast electrons created from the deposition of gamma-rays produced by the radioactive decay of 56Ni. In our models this is controlled by the gamma-ray energy deposition rate. In order to obtain a larger energy deposition rate, more mixing of radioactive 56Ni into the supernova’s outer envelope than exists in Blinnikov et al.’s model (Blinnikov, 1999; Blinnikov et al., 2000) is required. In PHOENIX this can be accomplished by replacing the self-consistent energy deposition rate from the Blinnikov et al. model with an ad hoc function that simulates a gamma-ray population created by excess nickel mixing. To simulate mixing, we assume a uniform mass fraction of 56Ni distributed homogeneously (i.e., $`\text{X}_{\text{Ni}}=`$ constant), and local deposition of gamma-rays. The new energy deposition rate is therefore much larger (and smoother) than the original (Figure 6), allowing for more hydrogen ionization and excitation in the line-forming region. Comparisons between the synthetic spectrum and the actual spectrum provide a sensitive probe for how much 56Ni must be mixed into the supernova envelope.

Figure 7 shows the results of modeling the day-4.52 spectrum with constant gamma-ray deposition, resulting from a nickel mass fraction of $`\text{X}_{\text{Ni}}=1.0\times 10^3`$ in the envelope. The fit between the predicted lines and the actual lines is much better, illustrating the validity of this model. A mass fraction of $`\text{X}_{\text{Ni}}=1.0\times 10^3`$ distributed thoughout our model atmosphere with a mass of about $`10^2`$ $`\text{M}_{}`$ corresponds to a requirement of approximately $`(1.01.2)\times 10^5`$ $`\text{M}_{}`$ of nickel to create enough gamma rays above the line forming region of the envelope. The mass is low (and corresponds to a lower boundary in velocity space of about 5000 km s-1), because at these early times the density is high and the continuum optical depth is $`>100`$ already at this mass/velocity point, so that spectrum formation occurs in the very outer layers of the ejecta. At later times we see deeper due to geometrical dilution. Although the model produced by PHOENIX utilized a constant nickel mass fraction (so that the number density of nickel follows the rather steep density gradient), this by no means precludes the possibility of inhomogeneous mixing. In the case of that the nickel is in clumps, our nickel mass is a lower limit (since the clumps could self-absorb gamma rays) and our results indicate that at least some of the clumps must be optically thin to the gamma rays that are produced.

In preliminary lower resolution studies, we found that the synthetic H$`\alpha `$ line in day 4.52 was slightly redshifted with respect to the actual H$`\alpha `$. This was also observed in synthetic spectra for days 3.71 and 5.76 created by the radiative transfer code SYNOW (Jeffery & Branch, 1990). This indicates that when the line strengths are correct, so is their velocity and demonstrates the diagnostic power of detailed spectral synthesis.

### 4.2 More Evidence for Nickel Mixing

The mechanism for the mixing of 56Ni into the envelope is the subject of intense debate, and won’t be settled until the core-collapse supernova mechanism is better understood. Kifonidis et al. (2000b) and Kifonidis et al. (2000a) suggested that nickel mixing into the hydrogen shell would be suppressed in Type II supernovae due to strong deceleration at the He/H boundary. Nevertheless, evidence exists to support the results of the PHOENIX models. Early detection of X-rays (Dotani et al., 1987) and gamma-rays (Matz et al., 1988), the expansion velocities in infrared line widths (Erickson et al., 1988; Witteborn et al., 1989) and in 56Co gamma-ray lines (Barthelmy et al., 1989), the detection of H velocities as low as 800 km s-1 at day 221 (Höflich, 1988), and the late detection of the He I 10830 line (Fassia & Meikle, 1999) all suggest that significant mixing has occurred in the ejecta \[see Arnett et al. (1989), McCray (1993) and references therein\].

Two-dimensional, axially symmetric models calculated by Fryxell et al. (1991) and by Herant & Benz (1991, 1992) suggested that 56Ni could be mixed into the envelope via finger-like structures that are produced shortly after the explosion by instabilities in the matter behind the shock front. By $`t=90`$ days, these fingers of nickel had broken up and expanded through radioactive heating into low-density bubbles of cobalt or iron, the decay products of 56Ni.

Inhomogeneous mixing has already been suggested as an explanation for the so-called “Bochum event,” an emission satellite on the red side of the H$`\alpha `$ line in days 20–100. Utrobin et al. (1995) suggested that a high-velocity clump of 56Ni caused an asymmetric over-excitation of hydrogen in the outer envelope. Recent studies by Chandra have detected the presence of diffuse high-velocity Fe bubbles at the outer edge of the Cassiopeia A supernova remnant. Hughes et al. (2000) noted the detection of Fe bubbles at the outer edge of the supernova ejecta, outside the Si-rich layers, and suggested that their diffuse nature could be explained by radioactive heating from decaying nickel.

Models of the composition of SN 1987A in the nebular phase by Kozma & Fransson (1998a), based on the results of various explosion models, used a number fraction of iron of $`2\times 10^5`$ in the hydrogen layer, with a filling factor $`f0.170.70`$. This corresponds to a mass fraction on the order of $`10^3`$. The results of the simulations of Kozma & Fransson (1998b) indicate that the iron mass fraction should be within a factor of 2 of the predicted value. Li et al. (1993) suggest that the light curve of SN 1987A during the nebular phase was dominated by $`0.078`$ $`\text{M}_{}`$ of 56Co in the form of 50 to 100 clumps. Their models contained nickel clumps that were initially opaque to gamma rays with $`f<0.01`$, accounting for $`1\%`$ of the total envelope mass. Gamma-ray absorption heated the clumps, causing them to expand until $`f=0.5`$. The nickel mass fraction and total mass values predicted by the PHOENIX runs are well within the limits imposed by these two studies.

Detection of the He I 10830 line at times after 10 days suggests that helium is also being re-ionized by 56Ni decay. Fassia & Meikle (1999) applied spectral synthesis models to SN 1987A and noted that the best fit to the helium line required at least 3% of the total 56Ni mass (0.07 $`\text{M}_{}`$) be mixed to above 3000 km s-1, but that 56Ni concentrations should be negligible above $`4000`$ km s-1. Fassia et al. (1998) carried out a similar analysis of the Type II supernova 1995V, where they noted the best fit to the He line required the 56Ni density to be constant up to 650–1150 km s-1, then fall off by a power law of 8–9 (consistent with our mixing assumption). Only $`10^6`$ $`\text{M}_{}`$ of 56Ni was required above the helium photosphere ($`4250`$ km s-1 at day 69) in the SN 1995V models. Burrows & Van Riper (1995), using Monte Carlo calculations of gamma-ray transport in clumpy debris models of SN 1987A’s envelope, suggest that up to 50% of the total 56Ni mass in SN 1987A should remain below 1000 km s-1.

Herant et al. (1994) suggest that significant iron/nickel inversion may be the result of neutrino-driven convection within the fraction of a second just prior to the supernova explosion. Convection would destroy the spherical symmetry of the elemental distribution, possibly creating the iron/nickel clumps predicted and/or detected in later stages. This theory is supported by 2D hydrodynamic simulations performed by Kifonidis et al. (2000a), although in their models the nickel was strongly decelerated at the He/H interface.

## 5 Further Research

Future work with the Blinnikov et al. 1987A model will involve using PHOENIX to synthesize the entire time series of the supernova up to around day 100. The detailed time series will indicate the location of nickel in the entire envelope, since the spectra probe deeper as time goes on due to the geometrical thinning of the envelope. PHOENIX spectra for this period may shed light on the later-time phenomena discussed in § 4. Once the whole time series of SN 1987A is modeled, the synthetic spectra will be applied to the Spectral-fitting Expanding Atmosphere Method (SEAM) in order to calculate the distance to SN 1987A, and hence to the LMC (Baron et al., 1993, 1994, 1995a, 1995b). Our calculations will also be useful to ascertain whether a 1-D representation of SN 1987A is sufficient. So far both breakout models (Blinnikov et al., 2000) and our preliminary calculations here indicate this to be the case.

We thank Claes Fransson for support and helpful discussions. This work was supported in part by NSF grants AST-9731450, NASA grant NAG5-3505, and an IBM SUR grant to the University of Oklahoma; and by NSF grant AST-9720704, NASA ATP grant NAG 5-8425 and LTSA grant NAG 5-3619 to the University of Georgia. SB is supported in Russia by the grant RFBR 99-02-16205. PHH was supported in part by the Pôle Scientifique de Modélisation Numérique at ENS-Lyon. Some of the calculations presented here were performed at the San Diego Supercomputer Center (SDSC), supported by the NSF, and at the National Energy Research Supercomputer Center (NERSC), supported by the U.S. DOE. We thank both these institutions for a generous allocation of computer time.