Gamma ray bursts (or GRBs) are sudden and intense bursts of radiation, peaking in the gamma ray band. GRBs are among the most luminous events observed from Earth. They produce tightly directed jets of highly relativistic matter, and the total energy released is on the order of 10^51 to 10^53 ergs. The majority of this energy is concentrated in the jet, so GRBs as observed from Earth appear significantly more luminous when viewed from directly on-axis. (If the energy flux within the jet were emitted isotropically, the total energy emitted would exceed 10^54 ergs--on the order of one solar mass!)

The intense energy involved suggests that GRBs result from stellar-scale explosions--supernovae producing a dense stellar end state such as a neutron star or black hole. The precise mechanism that produces a GRB--rather than an ordinary supernova lacking relativistic jets--is still poorly understood. However, numerous methods exist to model their formation and explain the radiation spectrum GRBs emit. Despite their name, gamma ray bursts do not solely involve gamma rays. Some GRBs are followed by a an afterglow at lower frequencies, drifting through the X-ray, optical, and radio bands.

One model of GRBs injects relavistic matter at the center of a star as a "central engine" operating over a specific opening angle and duration. When the injected jet of matter emerges from the star, it forms a highly relativistic shockwave which expands into the medium around the star. The interaction between this shock and magnetic fields produces synchrotron radiation (or, by a related process, jitter radiation). As the shock expands and cools, the peak frequency of the spectrum decreases.

Previous research has modeled the radiation emitted by an isotropic, relativistic explosion. However, GRBs are highly anisotropic, and this must be considered to accurately model the evolution of a GRB's spectrum as viewed from Earth. The goal of this research is thus to apply radiation-modeling techniques to a simulated (and therefore anisotropic) gamma ray burst viewed from an arbitrary angle.

Pre-existing simulations [1] model gamma ray bursts as a "central engine" placed at the core of a massive star. The engine injects jets of relativistic matter into the star, given a specified energy content, duration of the jet, and opening angle. The simulations track the evolution of the system during and after the period when the central engine is active.

The following images are generated from one such simulation, in which 3*10^51 ergs are released over a 15 second duration and an opening angle of 10 degrees. The images are cross-sections of the upper hemisphere of the star; lighter color corresponds to higher density.

At t=0s, the central engine has not had a chance to disrupt the star, so it appears as a basic sphere.

At t=10s, the energetic jet has burrowed through much of the star. The material is not moving at relativstic speeds since it must travel through the relatively unenergetic, dense material of the star.

At t=11s, the energetic jet has just emerged from the star, generating a shockwave in the rarefied medium surrounding the star.

At t=16s, the shock has expanded to roughly 4 times the diameter of the star. Note that the shockwave is expanding relativistically; it took only 5 seconds to expand to four times the star's diameter, compared to the 10 seconds needed for the jet to burrow through the star's radius.

The above density cross-sections are taken from various times during the central engine simulations. The final frames of these simulations represent the star and shockwave system at late times, after the central engine has been deactivated (so no additional energy is being added to the system) and the shockwave has been allowed to expand until it is much larger than the original star.

In addition to density, the simulations also track pressure and velocity. This is enough information to calculate the kinetic energy density as a function of spatial coordinates. The cross-sections from the central engine simulations are divided into a fine grid. We assume the density, pressure, and velocity to be uniform within each individual grid cell. We also assume the system exhibits symmetry when rotated about the vertical axis and when reflected about the equatorial plane. Given this, we can calculate the energy distribution in spherical coordinates by the following formulae.

The adiabatic index (denoted as gamma in much of the literature, and subscripted "ad" here to slightly reduce confusion) is equal to the ratio of specific heats; for a relativistic, monatomic gas, this equals 4/3.

The Lorentz factor at infinity describes the velocity at which the gas will be moving after maximal adiabatic expansion. It is found by multiplying the Lorentz factor gamma at the present time by a correcting factor, equal to one plus the product of 3, the adiabatic index, and the ratio of pressure over density.

The energy content of a cell (in the shape of a hollow cylinder due to rotational symmetry) is equal to the energy density multiplied by the volume of the cell. The kinetic energy density is equal to the mass density rho, multiplied by the Lorentz factor at infinity minus one, and multiplied by the Lorentz factor at current times to correct for relativistic length contraction. For simplicity in these calculations, c is set to be one.

Since each cell in the grid has a corresponding energy content and velocity, it is possible to arrange energy as a function of velocity. It is often useful to use the product of beta (velocity in units of c) and gamma (the Lorentz factor) as the velocity metric. Beta must always be between zero and one, and gamma must always be greater than one--both therefore can only be graphed over a limited domain. Using gamma*beta provides a much broader domain on which to consider the energy distribution. From this, we can construct cumulative energy distributions like those graphed below.

Above are plotted the cumulative energy distributions from three simulations with 3*10^51 ergs, a 10-degree opening angle, and a duration of 4.0 seconds, 7.5 seconds, and 15.0 seconds respectively. For a given gamma*beta on the horizontal axis, the value on the vertical axis represents the fraction of the total energy corresponding to a gamma*beta ABOVE the given gamma*beta. This graph illustrates that the duration of the central engine strongly influences the velocity distribution of the resulting explosion.

A short-duration (4.0 second) engine produces an energy distribution which is largely sub-relativistic (low gamma*beta). All the central engine's energy is input to the system before the jets of material can burrow their way out of the star, so all the material expends its energy escaping the star in a slower, more uniform explosion--in many ways similar to a supernova which does not produce a gamma ray burst. However, a longer-duration (15.0 second) central engine produces an extremely relativistic distribution (gamma*beta above 100). The longer central engine duration allows the earliest injected material to expend its energy emerging from the star, allowing the material injected later to emerge at highly relativistic velocities with minimal energy loss. This directed energy outflow resembles a typical gamma ray burst.

A 1976 paper by R. D. Blandford and C. F. McKee [2] describes a model for the evolution of a uniform, relativistic explosion, given several initial conditions including the energy injected at the explosion's epicenter and the pressure and density of the medium into which the explosion expands. This model can be applied to the explosions in our gamma ray burst simulations, since the star--the epicenter--is significantly smaller than the relativistic shockwave of the explosion.

The Blandford-McKee solution describes an isotropic explosion, but it can still be used for our highly anisotropic gamma-ray bursts. The shockwave expands radially outward at relativistic velocity but exhibits negligible motion in theta or phi. We can therefore divide the explosion into a grid in theta and phi, and treat each subdivision as if it were a segment of an isotropic explosion. For each explosion, we assume all conditions such as pressure and density are the same, and that only the isotropic equivalent energy varies. The isotropic equivalent energy is given by the formula below, where sigma and E subscripted "slice" represent the solid angle covered by and energy contained in a section of the explosion bounded by the particular limits in theta and phi.

To correctly calculate the spectrum emitted from the system, we must consider the shape of the expanding shock as observed by a distant viewer. Since the shock is moving relativistically, the spectrum seen by an observer at a given time is composed of radiation emitted over a broad range of times in the frame of the shock. To account for this, we calculate the apparent shape of the shock from an observer's perspective. A plot illustrating the extreme distortion of the roughly spherical shock is shown below.

Above is a plot of the location of the edge of the shock in Cartesian coordinates, as seen by an observer located far along the positive x-axis looking toward the origin. The shock considered here actually results from an isotropic explosion, resulting from 10^52 ergs injected uniformly. From the observer's perspective, it is distorted into something resembling an elongated football or highly eccentric ellipse. Note that the axes are not equivalent; the horizontal axis covers 10^17 cm, while the vertical axis covers only 10^16 cm. The shock's true shape in the observer's frame is thus far more distorted than it appears here. The distortion is even more severe for anisotropic gamma ray bursts.

Once we have determined the time-evolution of the shock's structure, we can calculate an observed spectrum by modeling the emission of synchrotron radiation, using a similar technique to that employed by Morsony et. al. in 2009 [3]. We begin by dividing the interior of the shock into a grid, then determining the conditions of pressure, density, and velocity within each grid cell. This step is complicated by the fact that the spectrum observed at a given time was emitted over a broad range of times in the shock's frame, as mentioned earlier. Once these conditions have been determined, we implement a model of synchrotron radiation to calculate the absorption and emissivity at each cell of the grid, then integrate along the line of sight to find the total power as a function of frequency for various observation times. We can compare these calculated spectra to observed spectra of actual gamma ray bursts.

The spectra of observed gamma ray bursts often exhibit an afterglow which decays slowly through the electromagnetic spectrum. As the shockwave expands, it slows and cools, and the peak frequency of the emitted radiation consequently decreases. In the X-ray band, the flux at a particular frequency decays with time according to a simple power law:

In some cases, the decay rate of the X-ray flux will increase at late times (that is, the exponent in the above equation becomes more negative), roughly 10^5 seconds after the activation of the central engine. This is referred to as a "jet break." Jet breaks are not always observed; we hypothesize that the time at which the jet break occurs (or whether it occurs at all) may depend on the angle between the axis of the gamma-ray burst and the line of sight to the observer.

It is possible to crudely mimic the structure of a narrowly directed jet by constructing a jet which inputs energy uniformly, but only between zero and ten degrees of the vertical axis. The analog is imprecise, since the "top-hat" jet lacks the fine structure and smoothed edges of the jets from an actual gamma ray burst. (The top-hat jet is named for the shape of its energy versus angle graph, above; note that the units on the vertical axis are incorrect and should be divided by the product of 4 and pi.) However, the top-hat jet can still be used to roughly calculate the spectrum which would be observed from a gamma ray burst. When input to our model for synchrotron radiation, the top-hat jet yields the power-versus-frequency graphs at various times animated below.

View each graph individually: 1s | 10s | 100s | 1000s | 10^4s | 10^5s | 10^6s

The lower-frequency slope to the left of the peak is proportional to frequency squared, and is due to self-absorption of lower-frequency radiation shortly after emission. The higher-frequency drop-off right of the peak is proportional to frequency to the negative 5/4; this is a combination of dimming with time and the obvious movement of the peak frequency to the left.

We are continuing to analyze gamma ray burst simulations, and to optimize the interface between the GRB simulations and the programs which calculate the emitted spectrum. We are also modifying the program used to transform coordinates and view GRB spectra from off-axis. Once both projects are complete, we will be able to generate a much larger set of spectra for various conditions and viewing angles. We can then begin making thorough comparisons between our simulations and observed spectra, as well as looking for jet breaks. As of yet, we do not have sufficient simulations to identify the occurence of jet breaks, and expanding the bank of analyzed spectra will also further this goal.

[1] D. Lazzati et. al. “Unifying the Zoo of Jet-Driven Stellar Explosions.” Retrieved from arXiv: http://arxiv.org/abs/1111.0970. Uploaded 13 March 2012.

[2] R. D. Blandford and C. F. McKee. “Fluid Dynamics of Relativistic Blast Waves.” *Physics of Fluids*, Vol. 19, No. 8, Aug. 1976

[3] B. J. Morsony et. al. “Jitter Radiation Images, Spectra and Light Curves from a Relativistic Spherical Blastwave.” *Mon. Not. R. Astron. Soc.* **392**, 1397-1402 (2009)

This research was supported by the National Science Foundation and the University of Wisconsin-Madison REU Program.

Special thanks to my mentor, Dr. Brian Morsony.

Thanks also to Dr. Morsony's collaborator, Dr. Davide Lazzati, for providing some of the code used in generating spectra from our simulations.