text:
# A rejection-free Monte Carlo method for the hard-disk system

## I Introduction

Monte Carlo (MC) methods have become more powerful tools with the development of faster and more accessible computers. Many different phenomena have been studied with MC methods Binder ; Landau . The melting behavior of the hard-disk system Alder  is one of such subjects which has been studied with Monte Carlo methods LeeStrandburg ; Sengupta ; Fernandez ; Zollweg ; Jaster ; Weber . Consider the hard-disk system with $`N`$ particles. The Monte Carlo method for the hard-disk system has the following steps. First, choose one particle randomly, i.e. choose it with a probability of $`1/N`$. Then choose a new position for the center of the chosen particle, the new position is chosen uniformly in the circle with radius $`\sigma _s`$ centered on the original position of the particle. The trial move is accepted when the new position has no overlap with other particles, otherwise, the trial move is rejected.

The MC methods introduced above have been used to obtain the equilibrium state of the system. To study the equilibrium state of the system, the value of $`\sigma _s`$, the trial step length, is often chosen in order that the rejection rate is near 50%. Recently, Watanabe et al. Watanabe ; Watanabe2  studied the hard-disk system by observing nonequilibrium behavior using a molecular dynamics (MD) method. In order to compare the results of MD and MC, we have to consider the dynamics of the MC method. The dynamics of the MC method for the hard-disk system can be understood as a Brownian motion, like pollen particles in a fluid. The value of $`\sigma _s`$ corresponds to the amplitude of the external random force. Therefore, it determines a time scale of this system. If values of $`\sigma _s`$ are not identical for all densities, we cannot compare the dynamics, like relaxation phenomena and fluctuations, between different densities. However, keeping $`\sigma _s`$ fixed leads to high rejection rates in the Monte Carlo at high densities. A similar problem occurs for spin systems both at low temperatures and in a strong external field: a rejection rate becomes very high, so a huge number of trials is required to make a change. This problem can be avoided using a different technique called the rejection-free Monte Carlo (RFMC) method. The RFMC method was first constructed for discrete spin systems Bortz , for a review see Novotny . Miyashita and Takano MiyashitaTakano  applied the RFMC method to the kinetic Ising model in order to study dynamical critical behavior. Recently, Muñoz et al. Munoz  proposed a new RFMC algorithm which can treat models with continuous degrees of freedom. In this paper, we construct and utilize a RFMC algorithm for the hard-disk system based on the method of Ref. Munoz .

## II Methodology

### II.1 A Rejection-free Monte Carlo method

A Monte Carlo method is an implementation of a Markov process on a computer, and hence is sometimes called a Markov Chain Monte Carlo. The Monte Carlo method calculates various physical quantities by updating states of a system using random variables. These updating processes can be illustrated by a Markov chain (see the schematic in Fig. 1). Let the current state be at $`S_0`$ and the states possible to move from $`S_0`$ are denoted by $`S_i(i=1,2,\mathrm{},N_s)`$. Define $`E_i`$ as the energy of the state $`S_i`$. The new state $`S_i(i=0,1,\mathrm{},N_s)`$ will be chosen with probability $`P(i|0)`$. One way to ensure that the system will relax to the equilibrium state is to insist that the probability $`P(i|0)`$ satisfy the detailed balance condition Landau .

One of the well-known ways Landau  to satisfy the detailed balance condition is to use a heat-bath transition probability. In the heat-bath method, the probability $`P(i|0)`$ is defined to be

$$P(i|0)=\frac{\mathrm{exp}(\beta E_i)}{_{k=0}^{N_s}\mathrm{exp}(\beta E_k)}.$$
(1)
. When a system has a continuous degree of freedom, the summation of Eq. (1) becomes an integration which is generally difficult to calculate analytically.

Another popular way to satisfy the detailed balance condition is a Metropolis method. In this method, each step contains two parts; selecting a new state and accepting or rejecting the trial to move to the selected state. First, pick up a state $`S_i`$ from all possible states to move to with uniform probability $`1/N_s`$. The probability $`P(i|0)`$ to move from $`S_0`$ to $`S_i`$ is defined to be $`1`$ when $`\mathrm{\Delta }E_i<0`$, otherwise it is $`\mathrm{exp}\left(\beta \mathrm{\Delta }E_i\right)`$ with the energy difference $`\mathrm{\Delta }E_iE_iE_0`$. Therefore, the probability $`P(i|0)`$ in the Metropolis method is

$$P(i|0)=\{\begin{array}{cc}1/N_s& \text{if}\mathrm{\Delta }E_i0,\\ \mathrm{exp}\left(\beta \mathrm{\Delta }E_i\right)/N_s& \text{otherwise}.\end{array}$$
(2)
When a trial is rejected, the configuration of the system is not updated. The probability $`\lambda P(0|0)`$ to stay in the current state after the trial is given by the expression

$$\lambda =1\underset{i0}{}P(i|0).$$
(3)
For some parameters, e.g. under a strong external field and at an extremely low temperature, the value of $`\lambda `$ can be very nearly $`1`$. In such cases, most of computational time is spent on calculating trials which will be rejected. This rejection rate drastically decreases the efficiency of the computation.

In order to overcome this problem, a rejection-free Monte Carlo (RFMC) method is proposed. It is an example of an event driven algorithm Binder ; Landau  and has also been called a waiting time method MiyashitaTakano ; Jesper . Each steps of the RFMC method involves first computing the time to leave the current state (the waiting time $`t_{\text{wait}}`$), and then choosing a new state to move to with the appropriate probability. It does not contain the judgment to accept or reject a trial, and, therefore, it achieves rejectionless updates of the system in each algorithmic step. The waiting time $`t_{\text{wait}}`$ is calculated using a (pseudo) random variable. The probability $`p(t)`$ to keep staying at the current state for $`t`$ steps decays exponentially,

$$p(t)=\lambda ^t=\mathrm{exp}(t\mathrm{ln}\lambda ),$$
(4)
with $`\lambda `$ defined in Eq. (3). Note that $`\mathrm{ln}\lambda <0`$ since $`0<\lambda <1`$. Inversely, the time $`t`$ to stay in the current state is determined with a uniform random number $`r`$ on $`(0,1)`$ to be

$$t_{\text{wait}}=\frac{\mathrm{ln}r}{\mathrm{ln}\lambda }+1,$$
(5)
where $`x`$ denotes the integer part of $`x`$ and the rounding down is introduced to express the discrete time step in MC Novotny ; Munoz .

After the time of the system is advanced by $`t_{\text{wait}}`$, a new state $`S_i`$ is chosen from the all states possible to move to except for the current state with the probability proportional to $`P(i|0)`$ NovotnyPRL ; Novotny ; Munoz . Since all values of $`P(i|0)`$ are required to proceed one algorithm step in the RFMC, the computational cost of one step is higher than the normal MC. However, the waiting time $`t_{\text{wait}}`$, the time which can be advanced in one algorithmic step, becomes large at low temperatures, and consequently the efficiency of the RFMC can become better than that of standard MC.

### II.2 Application to hard-disk systems

Consider a hard disk system with $`N`$ particles with the radius $`\sigma `$. A standard MC method for the system involves choosing a particle, and trying to move the chosen particle within a circle with radius $`\sigma _s`$ centered on the original position of the chosen particle. To apply a RFMC method to the hard-disk system, define $`\lambda _i`$ as the probability that a trial to move particle $`i`$ is rejected (given that particle $`i`$ was chosen as the particle to attempt a move). Using the definition of $`\lambda _i`$, we can construct the algorithm of the RFMC method for the hard-disk system as follows:

1. Calculate the waiting time $`t_{\text{wait}}`$ using Eq. (5) with $`\lambda =\frac{1}{N}_i\lambda _i`$.
2. Advance the time of the system by $`t_{\text{wait}}`$.
3. Choose a particle $`i`$ with the probability proportional to $`1\lambda _i`$, which is the probability that (given that particle $`i`$ was the particle chosen for an attempted move) the trial to move the particle $`i`$ would be accepted.
4. Choose the new position of the chosen particle $`i`$ uniformly from all the points to which the particle $`i`$ is allowed to move.

The steps described above are the same as the RFMC for continuous spin systems Munoz , but the algorithms to calculate $`\lambda _i`$ and to determine a new position of the chosen particle are unique to the hard-disk system. The probability $`\lambda _i`$ can be calculated to be,

$$\lambda _i=1\frac{A_i}{\pi \sigma _{s}^{}{}_{}{}^{2}},$$
(6)
with an area $`A_i`$ into which the particle $`i`$ is allow to move (see Fig. 2). Therefore, if we can calculate $`A_i`$ for all particles, then we can construct the RFMC algorithm for the hard-disk system.

### II.3 Calculation of $`A_i`$

The area $`A_i`$ is the continuous set of positions in which the particle can be placed without any overlaps. Without neighboring particles, the shape of $`A_i`$ would be a filled circle with a radius $`\sigma _s`$. Let’s call it a trial circle. In the general case, the shape of $`A_i`$ is the remaining part of the trial circle after removing the overlap of ‘shadows’ of neighboring particles. The shape of the shadow is a circle with a radius $`2\sigma `$ which is concentric to a neighboring particle. Let’s call this a shadow circle. The area $`A_i`$, thus, consists of areas of arcs of a trial circle and that of shadow circles.

To compute the value of $`A_i`$, we develop a method we call the survival point method. See Fig. 3 (a). The chosen particle is shown as a solid circle, the trial circle is shown as a concentric dashed circle, and the area $`A_i`$ is the shaded region. Each neighboring particle (filled circles) has a shadow circle which is concentric and has radius $`2\sigma `$. An enlargement of the area $`A_i`$ is shown in Fig. 3 (b). It is seen that in this example this area has five vertices which are intersection points of shadow particles, we call them survived vertex points. In Fig. 3 (c), these survived vertex points are shown as small filled circles. Straight lines connect the center of the chosen particle and the intersection points. In this case the area $`A_i`$ is divided into five figures.

Each divided figure is the remaining part of an isosceles triangle with the overlap of a shadow circle removed. It is easy to calculate this area. Thus, all we have to do is to find all survived vertex points which form the area $`A_i`$. First, make a list of all intersection points of shadow circles and the trial circle. Next, remove points which are included in other shadow circles from the list, since these points cannot be vertices forming the area $`A_i`$. After this removal process, we have the vertices which form the area $`A_i`$ (see Fig. 3 (c)). The calculation process of a partial figure which forms $`A_i`$ is shown in Fig. 3 (d). The vertices are denoted by $`P_1`$ and $`P_2`$, and the center of the shadow circle is denoted by $`S`$. The survived vertices $`P_1`$ and $`P_2`$ are on the shadow circle centered at $`S`$, so $`\overline{SP_1}=\overline{SP_2}=2\sigma `$. The area of $`OP_1SP_2`$ can be calculated by summing the two triangles $`OP_1P_2`$ and $`SP_1P_2`$ with Heron’s formula. The area of the chord is $`4\sigma ^2\theta `$. Finally the portion of the area $`OP_1P_2`$ is calculated by subtracting the area of the chord $`SP_1P_2`$ from the area of the quadrilateral $`OP_1SP_2`$. The total area $`A_i`$ is the sum of one such calculation for each survived vertex.

### II.4 Choosing a particle to move

After calculation of $`t_{\text{wait}}`$ and advancing the time of the system by it, we have to choose a particle $`i`$ to move with a probability proportional to $`1\lambda _i`$. With a direct implementation, i.e., with the integration scheme IntegrationSearch , the order of the computation is $`O(N)`$, which is very time consuming. Other approaches are proposed like a three level search for spin systems threelevel . The three-level search improves the efficiency of the search by determining coordinates of a spin to update one by one. However, it is difficult to apply this method for particle systems, since neighbors of particles are not fixed. Here we use a complete binary tree search for the choosing part of the algorithm.

First, calculate the area $`A_i`$ for each of the particles. Since an acceptance probability $`1\lambda _i`$ is proportional to $`A_i`$ as shown in Eq. (6), the particle should be chosen with the probability proportional to $`A_i`$.

Next, construct a complete binary tree as follows,

1. Prepare a complete binary tree with enough height $`h`$, this height $`h`$ should satisfy $`2^{h2}<N2^{h1}`$.
2. Label each node with $`T_n^k`$, which denotes the $`n^{\mathrm{th}}`$ value at level $`k`$. The root node is labeled by $`T_1^h`$. A node labeled $`T_n^{k+1}`$ has branches leading to two nodes $`T_{2n1}^k`$ and $`T_{2n}^k`$.
3. Associate every bottom node $`T_i^1`$ with the value of area $`A_i`$. If the number of bottom nodes $`2^{h1}`$ is larger than $`N`$, the rest of the nodes are associated with zero, namely, $`T_i^1=0(i>N)`$.
4. Associate nodes at higher levels ($`k>1`$) recursively with the sum of the values associated with its two children, namely, $`T_n^{k+1}=T_{2n1}^k+T_{2n}^k`$.

A sample of a compete binary tree is shown in Fig. 4. Each node has the value $`T_n^k`$ and the value of each node at level $`k+1`$ is the sum of the values of its two children nodes at level $`k`$. The root node, which is $`T_1^4`$ in Fig. 4, has the sum of all $`A_i`$, that is,

$$T_1^h=\underset{i}{\overset{N}{}}A_i.$$
(7)
Using this tree, we can choose a particle with the probability proportional to $`A_i`$ in the following way.

1. $`k1`$, $`i1`$.
2. Prepare a random number $`r`$ uniform on $`(0,T_i^k)`$.
3. $`\{\begin{array}{cc}i2i1\hfill & \text{if}r<T_{2i1}^{k1}\\ i2i\hfill & \text{otherwise}\end{array}`$
4. $`ki1`$.
5. if $`k>1`$ then go to 2

Consequently, choosing the bottom node requires $`h1`$ random numbers.

After the above processes, the final value of $`i`$ indicates the index of the particle to move. The order of this search algorithm is $`O(\mathrm{log}N)`$. When the position of particle $`i`$ is moved, the value of $`A_i`$ is also modified. We only have to update part of this tree for the chosen particle and its neighbors. The order of this update is also $`O(\mathrm{log}N)`$, which is much faster than $`O(N)`$ of the direct implementation. Details to implement the complete binary tree search method are described in the appendix.

### II.5 Find a new position of the particle

After choosing a particle, we have to choose the new position for it. It is difficult to choose a point uniformly from the points allowed to move into, since its shape is generally very complicated (see Fig. 2). Therefore, we have chosen to choose the new position using a Monte Carlo rejection method. Namely we generate a random position uniformly over some bounding area that includes all of the area $`A_i`$. \[Such as the dashed circle of radius $`\sigma _s`$ in Fig 3(a) or the rectangle in Fig 3(c).\] If this point is not in $`A_i`$ it is rejected and another uniformly distributed random point over the bounding area is generated. The first point not be rejected is the new position of the particle, since it is in the area $`A_i`$, and this point is used as the new center for the particle. This completes one algorithmic step of the RFMC method.

A value of area $`A_i`$ is very small compared to the trial circle at high density, and hence the Monte Carlo trial to find the new position of the particle to be moved became very inefficient. To improve this, it is effective to limit the trial area for the Monte Carlo by making the bounding area very close to the area $`A_i`$. We outline two different survived point methods, but have only implemented the first.

For the first method, the one actually implemented in this paper see Fig. 3 (c). The solid rectangle is a bounding rectangle which includes the area $`A_i`$. It is easy to obtain the bounding rectangle with the survived vertex points. With the set of survived vertex points $`\{(x_i,y_i)\}`$, a diagonal line of the bounding rectangle is from $`(\mathrm{min}\{x_i\},\mathrm{min}\{y_i\})`$ to $`(\mathrm{max}\{x_i\},\mathrm{max}\{y_i\})`$. Then we can perform Monte Carlo trials for a new position within only this rectangle. The area of the rectangle is on the same order of $`A_i`$, so the probability of success to obtain the new position is drastically improved compared with the direct search over the trial circle.

An alternative method is first use a random number to decide which of the triangles formed with point O and two adjacent survived vertex points the survived point will fall into. This is done analytically since the area of each of these triangles with removed shadow circle chords have been already calculated. Then the shortest side formed with point O and the two survived vertex points (say $`\overline{SP_2}`$ in Fig. 3 (d)) is lengthened to be equal to the longest side ($`\overline{SP_1}`$ in Fig. 3 (d)). The random trial point is then generated within the section of the circle with a radius equal to the longest side ($`\overline{SP_1}`$ in Fig. 3 (d)). Then the point becomes the survived point used for the new location of point O if the trial point is within the shaded area. Otherwise, this procedure repeats in the same extended circular section until a survived point is found.

## III Results

### III.1 Calculation of $`A_i`$

In order to test our method to calculate $`A_i`$ described in Sec. II.3, the values of $`A_i`$ were also evaluated by a Monte Carlo sampling ($`A_{\mathrm{MC}}`$) with trial points uniformly drawn over the trial circle. The density of the system $`\rho `$ is defined to be $`\rho =4N\sigma ^2/L^2`$ with the number of particles $`N`$, the radius of the particles $`\sigma `$ and the linear system size $`L`$, respectively. Throughout this study, the number of particles $`N`$ is set to be $`23288`$ and the periodic boundary condition is taken for both axes. The number of the generated configurations were 3000, and $`10^6`$ MC trial points are taken for each of the configurations to evaluate its area. The density of the system is fixed at $`\rho =0.9`$. The result is shown in Fig. 5. The area $`A_i`$ is normalized by the area of the trial circles (see Fig. 2). The difference between the MC and our survived point method is less than $`0.01\%`$ for all areas, which is the same order as the statistical error of our MC method. The standard deviation of the MC area calculation is determined by dividing the data into 10 groups, each including $`10^5`$ samples. This result shows that the value of $`A_i`$ is properly calculated by our method.

### III.2 Time evolution

To compare the dynamics between the standard MC and the RFMC, we observed the time evolution of the six-fold bond-orientational order parameter $`\varphi _6`$ HalperinNelson . The parameter $`\varphi _6`$ is defined to be

$$\varphi _6=\mathrm{exp}(6i\theta ),$$
(8)
with the bond angle $`\theta `$ which has a definition described in Fig. 6. The average is taken for all pairs of neighboring particles. The parameter $`\varphi _6`$ becomes 1 when all particles are located on the points of a hexagonal grid, and it becomes $`0`$ when the particle location is completely disordered. Therefore $`\varphi _6`$ describes how close the system is to the perfect hexagonal packing. The neighbors in an off-lattice model are strictly defined with the Voronoi construction Jaster , which is a very time-consuming method. In this paper, two particles separated by a distance less than $`2.6\sigma `$ are defined as neighbors. We confirmed that the value of $`\varphi _6`$ is approximately the same value as the value obtained with the Voronoi construction.

At the beginning of the simulation, the particles are set up in a perfect hexagonal order, namely, $`\varphi _6(t=0,\rho )=1`$. The order parameter $`\varphi _6`$ starts to relax to the value of the equilibrium state. With this nonequilibrium relaxation (NER) behavior of order parameters, critical points and critical exponents of various phase transitions can be determined accurately ItoPhysica1 ; ItoPhysica2 ; ItoMastuhisa ; ItoHukushima . This method is called a NER method. Watanabe et al. Watanabe ; Watanabe2  studied two-dimensional melting based on the NER method for the Kosterlitz-Thouless transition NERKT  by observing the relaxation behavior of $`\varphi _6`$. Therefore, the following time evolutions of $`\varphi _6`$ contains information about the two-dimensional melting transition.

Time evolutions of $`\varphi _6`$ are shown in Fig. 7. Solid lines are results of the standard MC simulation and symbols (circles, triangles and squares) are the results of the RFMC. Fig. 7 shows that both behaviors are equivalent for the two methods. This is essentially a check of the program implementation, since the physical dynamic is the same for both the MC and the RFMC methods.

### III.3 Efficiency

The computational times required to achieve $`1000`$ accepted MC steps are shown in Fig. 8(a). Configurations are started from the perfect hexagonal configuration and both measurements are started after $`10^6`$ MC steps. All simulations are performed on an Intel Xeon $`2.4`$ GHz computer. While the computational time of the RFMC (open circles) is almost constant, a longer computational time is required for the standard Monte Carlo (solid circles) at higher density. It shows that the RFMC is more efficient at high densities, in spite of the additional bookkeeping involved in the RFMC method (so one algorithmic step takes much longer than one standard MC step).

The CPU-time ratio of the standard MC to the RFMC is shown in Fig. 8(b). The data are shown as a function of $`1/\epsilon `$, where $`\epsilon (\rho _{\text{cp}}\rho )/\rho _c`$, the closest density is $`\rho _{\text{cp}}`$ and the density of the system is $`\rho `$. The CPU-time ratio, which is the efficiency of the RFMC compared to that of the standard MC, diverges as $`\epsilon ^2`$.

## IV Summary and Discussion

We constructed a rejection-free Monte Carlo algorithm for the hard-disk system. This method conserves the property of the dynamic behavior of the original Monte Carlo method. In other words, the time scales will not depend on the density, but are rather set by some Brownian-motion type of dynamic for all densities. An estimate of the time scales between the MC and physical time can thus be obtained by setting the mean-free path of an isolated particle to be proportional to the value $`\sigma _s`$. Note that strictly this is only true in the limit $`\sigma _s0`$, but it should be a reasonable approximation for a small finite $`\sigma _s`$.

We also find that for a fixed value of $`\sigma _s`$, the RFMC method is more efficient at high density. Therefore, the RFMC method should be useful for studies of two-dimensional solids or studies of high-density glass materials. It may also be possible to make the algorithm even more efficient by utilizing the ideas of absorbing Markov chains (for the MCAMC method for discrete state spaces see Novotny  and references therein). Increased algorithmic efficiencies for the Monte Carlo dynamics of hard disks could be useful to further test physical phenomena using hard disk systems, such as for example the relationship between fluctuations and dissipation of work in a Joule experiment Cleuren .

## Acknowledgement

The authors thank S. Miyashita, P.A. Rikvold, and S. Todo for fruitful discussions. This work was partially supported by the Grant-in-Aid for Scientific Research (C), No. 15607003, of Japan Society for the Promotion of Science, and the Grant-in-Aid for Young Scientists (B), No. 14740229, of the Ministry of Education, Culture, Sports, Science and Technology of Japan, the 21st COE program, “Frontiers of Computational Science”, Nagoya University and by the U.S. NSF grants DMR-0120310, DMR-0426488, and DMR-0444051.

## Appendix

The complete binary tree search can be implemented with an one-dimensional array. To make it simple, let the number of particles $`N`$ be $`2^{h1}`$. The tree with height $`h`$ requires an array $`a(i)`$ with size $`2N1`$. First, associate each bottom node with a corresponding value as

$$a(N+i1)A_i(i=1,2,\mathrm{},N),$$
(9)
which corresponds to $`T_i^1A_i`$. Next, associate parent nodes recursively as

* $`iN`$
* While $`i0`$
* $`a(i)a(2i)+a(2i+1)`$
* $`ii1`$
* next $`i`$

Using this array, we can pick up particle $`i`$ with the probability proportional to $`A_i`$ as,

* $`i1`$
* While $`i<N`$
* Prepare a uniform random number $`r`$ on $`(0,a(i))`$
* $`\{\begin{array}{cc}i2i\hfill & \text{if}r<a(2i)\hfill \\ i2i+1\hfill & \text{otherwise}\hfill \end{array}`$
* next $`i`$
* $`iiN+1`$.

After the above procedure, we obtain the index $`i`$ of the chosen particle. When the value of $`A_i`$ is changed, the tree should be updated. The update process is as follows,

* $`a(N+i1)A_i`$
* $`i{\displaystyle \frac{i+N}{2}}`$
* While $`i1`$
* $`a(i)a(2i)+a(2i+1)`$
* $`ii/2`$
* next $`i`$.

Note that, when the chosen particle is moved, the acceptance probabilities of the neighboring particles of the moved particle are also changed. Therefore, we have to perform the above process for all neighboring particles.


format: text

text:
# GRB Redshift Evolution Within the Unified Jet Model

## 1 Introduction

Most objects at cosmological distances (stars, galaxies, AGN) display evolution of their observable properties with redshift. Since gamma-ray bursts (GRBs) are thought to originate in the core-collapse of massive stars and are observed over a wide range in redshift, it is reasonable to suppose that they might also evolve.

Since the advent of rapid GRB localizations with BeppoSAX and HETE-2, and the consequent follow-up observations, over $`30`$ GRBs have reported redshift measurements. Using 9 BeppoSAX bursts from amati2002 , plus an additional 11 bursts localized with HETE-2, graziani2003  have been able to strengthen earlier indications that GRBs are brighter at higher redshifts.

A uniform-jet model proposed by lamb2003  has been shown to provide a unified picture of GRBs and XRFs. Considering all bursts to have a “standard energy” but a range of jet opening solid-angles that spans five orders of magnitude, this unified jet model can account for the full range of observed burst properties seen by HETE-2 . Here we extend this unified jet model to account for redshift evolution and show that it can also explain the observed properties of BATSE bursts.

## 2 Observations of GRB Evolution

Analysis of the BATSE catalog has revealed evidence that GRBs may evolve strongly with redshift. The use of redshift estimators based on burst variability has provided evidence that bursts are intrinsically brighter at larger redshifts than at smaller lloyd2002 ; reichart2001b . Analyzing a set of 9 GRBs with spectroscopically determined redshifts observed by the BeppoSAX satellite, amati2002  also claim evidence for an increase in the isotropic-equivalent energy ($`E_{\mathrm{iso}}`$) with redshift, but they did not include a discussion of possible threshold selection effects.

Recent results from HETE-2 strengthen the evidence for this relationship. Figure 2 from graziani2003  shows the isotropic-equivalent energies $`E_{\mathrm{iso}}`$ and luminosities $`L_{\mathrm{iso}}`$ as a function of redshift for the HETE-2 and BeppoSAX events. After correcting for threshold effects, graziani2003  find that $`E_{\mathrm{iso}}`$ is correlated with redshift at the $`5.1`$% confidence level, and $`L_{\mathrm{iso}}`$ is correlated with redshift at the $`0.9`$% confidence level. The observed relationship goes roughly as $`E_{\mathrm{iso}}(1+z)^3`$.

## 3 Unified Jet Model Simulations

The unified jet model of GRBs and XRFs lamb2003  provides a natural explanation for redshift evolution in GRBs. Namely, each burst exhibits a “standard energy” frail2001 ; bloom2003 , but the possible range of jet opening angles varies from fairly large values at low redshift, to very small values at high redshift, according to the relationship given above. That is, evolution in $`E_{\mathrm{iso}}`$ is explained by an evolution of the jet opening solid-angle, $`\mathrm{\Omega }_{\mathrm{jet}}`$.

The burst simulations that we have implmented to test the unified jet model also provide a powerful way to explore models including burst evolution with redshift. For each burst, we obtain: (1) A redshift z by drawing from a model of the star-formation rate rr2001 , and (2) a jet-opening solid angle $`\mathrm{\Omega }_{\mathrm{jet}}`$ by drawing from specific distribution range in $`\mathrm{\Omega }_{\mathrm{jet}}`$ that is fixed at $`z=0`$ and shifts to smaller values at higher redshifts. We also introduce three Gaussian smearing functions to generate: (1) A spread in jet energy ($`E_\gamma `$, see Figure 1), (2) a spread in $`E_{\mathrm{peak}}`$ around the $`E_{\mathrm{iso}}`$-$`E_{\mathrm{peak}}`$ relation, and (3) a spread in the timescale T that converts fluence to flux. Using these five quantities, we calculate various rest-frame quantities ($`E_{\mathrm{iso}}`$, $`E_{\mathrm{peak}}`$, etc.), and finally, we construct a Band function for each burst and transform it to the observer frame, which allows us to calculate fluences and peak fluxes and determine if the burst would be detected by various instruments (see Figure 2a).

To obtain the model presented here, we sought to roughly match the observed distribution in the ($`1+z`$, $`E_{\mathrm{iso}}`$)-plane (compare Figure 2b with Figure 2 of graziani2003 ), which displays a range in $`E_{\mathrm{iso}}`$ of $`3\times 10^{49}`$ to $`1.5\times 10^{52}`$ ergs at $`z=0`$ and $`6.4\times 10^{51}`$ to $`3.2\times 10^{54}`$ ergs at $`z=5`$. Assuming the faintest burst at $`z=0`$ corresponds to $`\mathrm{\Omega }_{\mathrm{jet}}=2\pi `$, this translates into a range of $`\mathrm{\Omega }_{\mathrm{jet}}`$ of $`2\pi `$ to $`0.0125`$ steradians at $`z=0`$ and $`0.0291`$ to $`5.79\times 10^5`$ steradians at $`z=5`$.

The observed values of $`E_\gamma `$ are taken from bloom2003 , and their distribution is well-fit by narrow Gaussian frail2001 ; bloom2003 ; lamb2003b  (see Figure 1a). Figure 1b plots the displacement of these values from the central value as a function of redshift and shows no evidence for evolution of $`E_\gamma `$ with redshift lamb2003 . This rules out the possibility that redshift evolution might be explained by an evolution of the “standard energy”.

Since it relies on the random distribution of burst jet axes with respect to the viewing angle, the universal structured-jet model cannot easily accommodate redshift evolution. The only solution would be for the “jet energy”, $`E_\gamma `$, to evolve with redshift, but that would make it difficult to explain Figure 1 and the results of frail2001 ; bloom2003 .

## 4 Results

Figure 4 compares the cumulative distributions of four observable burst quantities against several possible models. We find that the uniform jet model with burst evolution can adequately describe the observed distributions of localized bursts.

In addition, adding evolution to the uniform jet model replicates the observed distribution of peak photon number fluxes as observed by BATSE. Models without redshift evolution (Figure 3a, dotted curve) tend to overpredict the number of high peak flux bursts. However, models with strong evolution (solid curve) provide excellent agreement with the BATSE distribution. Figure 3b shows the differential distribution of redshifts (for bursts detected by WXM) for the models with and without redshift evolution.

This model makes several predictions. Most observed XRFs are predicted to be at $`z<1`$. Bursts at low $`z`$ have $`\mathrm{\Omega }_{\mathrm{jet}}>10^2`$ or $`\theta _{\mathrm{jet}}`$ a few degrees. We require a value of the “standard energy” to be $`E_\gamma 5\times 10^{49}`$ ergs, or about $`50`$ times less than the value reported by bloom2003 . Thus, the fraction of Type Ic supernovae producing GRBs increases from $`0.1`$% at $`z=0`$ to $`10`$% at $`z5`$. Finally, $`70`$% of bursts with $`z>5`$ are detected by the WXM.

## 5 Conclusions

HETE-2 has strengthened the evidence that GRBs evolve with z. The uniform jet model can describe XRFs and GRBs and can accommodate evolution whereas the universal jet model cannot.

Confirmation of this model will require the localization of many more XRFs, the determination of $`E_{\mathrm{peak}}`$ and $`E_{\mathrm{iso}}`$ for many more XRFs and GRBs, and the identification of optical afterglows and the measurement of redshifts for these bursts.

HETE-2 is ideally suited to localize XRFs and study their spectra, but this will be difficult for Swift , which has a nominal threshold of $`E_{\mathrm{min}}15`$ keV and a narrow energy band of $`15\mathrm{keV}<E<150`$ keV. However, Swift is optimized for pinpointing X-ray and optical afterglows, and facilitating spectroscopic redshift measurements. Therefore, it is very important that the HETE-2 mission continue, even after Swift is fully operational. A partnership between HETE-2 and Swift can confirm or rule out GRB evolution with redshift.


text:
# Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context

## 1 Introduction

It is now clear that if stellar core collapse simulations are constrained to be spherically-symmetric the cores do not supernova (Burrows, Hayes, and Fryxell 1995; Mezzacappa et al. 2001; Rampp and Janka 2000,2002; Liebendörfer et al. 2001,2003; Buras et al. 2003; Thompson, Burrows, and Pinto 2003). However, some calculations performed in 2D (Herant et al. 1994; Burrows, Hayes, and Fryxell 1995; Fryer et al. 1999) and in 3D (Fryer and Warren 2002,2003) do explode, while others performed in 2D do not (Mezzacappa et al. 1998; Buras et al. 2003). Even though a consensus is emerging that multi-D effects are crucial to the mechanism of core-collapse supernovae, no censensus has yet emerged concerning which multi-D effects are crucial. The enhancement of the efficiency of neutrino heating behind the stalled shock wave due to neutrino-driven convection seems to play a role, but this conclusion too remains controversial. It is suspected that the approximations to neutrino transport incorporated into the codes which see explosions may have compromised the validity of those results. Furthermore, the degree of neutrino heating in the semi-transparent region behind the stalled shock may depend upon a proper treatment of multi-D transfer in the unstable and convective core. Moreover, the viability of doubly-diffusive instabilites in protoneutron stars hinges upon neutrino transport in non-radial directions (Mayle and Wilson 1988; Bruenn and Dineva 1996). Hence, there is now an emphasis in supernova theory circles on improving the treatment of neutrino radiation transport in two and three dimensions. We have embarked upon such a program of improvement and this paper is the first in a series on our evolving efforts.

In many astrophysical problems matter and radiation are strongly coupled and are described by the system of non-linear radiation hydrodynamics equations. When the optical depth is large everywhere, the diffusion approximation can be used. In the diffusion limit, the radiation field is in local thermodynamic equilibrium (LTE) with the matter and the radiation is described by the non-linear diffusion equation. However, when important subregions of the system are optically thin, a full Boltzmann problem must be solved. A typical example of the latter is the core-collapse supernova problem, where the radiation field consists of several neutrino types, the optical depth varies rapidly with time and position, and the coupling to matter may determine the failure or success of the neutrino-driven explosion mechanism. In that problem, the radiation field depends upon position, direction, frequency (energy), and species, for which any discretization scheme becomes very large and expensive. For this reason, to date most time-dependent numerical radiation/hydrodynamics conducted in the astrophysical context have been limited to one spatial dimension, or if multi-dimensional, has been either flux-limited and gray (“one-group”) (Stone, Mihalas, Norman 1992; Hayes and Norman 2003) or carried out along radial rays (Rampp and Janka 2000,2002; Buras et al. 2003), or both (Burrows, Hayes, and Fryxell 1995). The 2D calculations of Burrows, Hayes, and Fryxell (1995) incorporated flux-limited, gray transport along radial rays, and ignored transverse transport. The 2D calculations of Buras et al. (2003) used a Boltzmann solver, but were conducted in 1D along various radial rays. Their hydrodynamics is, however, done in 2D. These authors correct for transverse radiation pressure and transverse lepton transport using a further update split off from the radial transport. The error introduced by such a procedure is unknown. The 2D calculations of Mezzacappa et al. (1998) actually employed the Multi-Group-Flux-Limited Diffusion (MGFLD) neutrino heating rates from a previous 1D simulation. These deposition rates were incorporated into a 2D hydro simulation, but without the feedback that would have arisen from the large differences between the 1D and 2D density and composition profiles and the shock history, as well as from the angular variation seen in 2D. The ZEUS-2D radhydro code of Stone, Mihalas, and Norman (1992), and its update by Hayes and Norman (2003), solve the zeroth- and first-moment equations of the gray transport equation and use a short-characteristics solution of the static transfer equation to obtain the second moment for closure. These codes are implementations of the variable Eddington factor technique and avoid some of the pitfalls of flux limiters, but since they are gray and not multi-group they are of limited utility for modern supernova simulations. The recent increases in computing power now make possible the development of more sophisticated codes to perform 2D, multi-angle, multi-group, multi-species radiation hydrodynamics simulations.

The vast nuclear engineering literature contains numerous methods for solving the Boltzmann equation (Adams and Larsen 2002). However, these may not be appropriate for astrophysical problems. In astrophysics, we confront stellar configurations in which the optical depth varies by many orders of magnitude from their centers to their peripheries and the density and mean free paths change rapidly with radius. In this paper, we present a test core-collapse simulation using a new algorithm to solve the multi-group Boltzmann equation in one and two spatial dimensions which we incorporated into the broader framework of radiation hydrodynamics. We defer a more complete discussion of the solver, stencils, and difference scheme to a later work (Livne et al. 2004).

We summarize aspects of our Boltzmann transport scheme in one and two spatial dimensions in §2 and §3. For a full radiation/hydrodynamics code, one needs to couple the radiation field to matter through emission, absorption, scattering, and radiation pressure. The resulting time-, space-, energy-group-, and angle-dependent code has both 1D and 2D (axially-symmetric) variants. Section 3 is devoted to a précis of the radiation hydrodynamics scheme in which the new transport module is incorporated. We also touch on various subtleties concerning 2D gravity and momentum conservation and the importance of freeing the inner core to move along the symmetry axis. In §4, we summarize various computational issues concerning problem size and parallelism philosophy. In §5, we first discuss the results of our comparison of 1D stellar collapse simulations using VULCAN/1D with the Feautrier/tangent-ray code “SESAME” 111Spherical, Explicit/Implicit, Supernova, Algorithm for, Multi-Group/Multi-Angle, Explosion Simulations (Burrows et al. 2000; Thompson, Burrows, and Pinto 2003) and then go on to describe our preliminary 2D simulation results in the core-collapse context. This test simulation was carried out for 22 milliseconds of physical time. We focus on the relationship between an anisotropic matter distribution and its corresponding anisotropic radiation field. This is the first 2D time-dependent, multi-group, multi-angle radiation/hydro calculation ever performed in core collapse studies. For these first test simulations, only electron neutrinos are incorporated, though the code is set up to handle the other species. The 2D code can employ unstructured fixed or moving meshes and incorporates the centripetal terms due to rotation about the axis of symmetry. Finally, in §6 we list the accomplishments represented by this paper, as well as the significant work that remains yet to be done.

## 2 The Transport Equations and Approach

The transport equation for the specific intensity ($`I`$)

$$\frac{1}{c}\frac{DI}{Dt}+𝛀I+\sigma I=S$$
(1)
describes for all directions $`𝛀`$ the additions to and subtractions from a beam of massless particles due to the physical processes of absorption, scattering, emission, and streaming. In the notation of Morel et al. (1996), $`I(𝐫,𝛀,\epsilon _\nu ,t)`$ is the angular specific intensity (in energy units), $`\sigma ^a(𝐫,\epsilon _\nu )`$ is the absorption cross section times the number density (inverse absorption mean-free-path), $`\sigma ^s(𝐫,\epsilon _\nu )`$ is the inverse scattering mean-free-path, $`\sigma `$ equals $`\sigma ^a+\sigma ^s`$, $`S`$ (the source term) equals $`S_{em}(𝐫,\epsilon _\nu )+\sigma ^sJ`$ (where $`S_{em}`$ is the emissivity), and $`J(𝐫,\epsilon _\nu ,t)`$ equals $`\frac{1}{4\pi }I𝑑𝛀`$ (the zeroth moment of $`I(𝐫,𝛀,\epsilon _\nu ,t)`$). Scattering is currently assumed to be isotropic, though this is not an inherent feature of our formulation, and we use the transport cross section, $`(1\mathrm{cos}\theta )\sigma _T^s`$, instead of the total scattering cross section ($`\sigma _T^s`$). This approach has a long pedigree and has been shown to work quite well in the 1D core-collapse context (Burrows et al. 2000; Thompson, Burrows, and Pinto 2003). In eq. (1), an energy-group subscript for each energy $`\epsilon _\nu `$ is implied. We solve eq. (1) implicitly and for every group and angle. Note that the term $`\frac{1}{c}\frac{DI}{Dt}`$ is the total, Lagrangian time derivative. In eq. (1), the other velocity-dependent terms, such as the Doppler (radiation compression) and aberration terms (Mihalas and Mihalas 1984), have been dropped. Such velocity effects are not dominant in the core-collapse context, but the compression term in particular would need to be included in a future generation of supernova simulations if neutrino viscosity effects are to be properly treated. Coupling as they do energy groups, angular bins, and spatial zones in the type of global solution required for a fully implicit treatment, these two velocity-dependent terms represent one of the most complicated numerical challenges of multi-D transport. One component of this challenge, as yet unsolved, is to incorporate non-monotonic velocity fields, with the resulting complicated characteristic structures, into an implicit solver. This fact is why to date no time-dependent multi-group, multi-angle astrophysics code has properly included them. We plan to do so in a later version using an operator-split, but explicit, approach that we have developed for the spherically-symmetric code SESAME (Thompson, Burrows, and Pinto 2003).

In spherical symmetry, the one-dimensional transport equation (eq. 1) is

$$\frac{1}{c}\frac{DI}{Dt}+\mu \frac{I}{r}+\frac{(1\mu ^2)}{r}\frac{I}{\mu }+\sigma I=S.$$
(2)
The discretization of the radial transport equation is easiest when eq. (2) is written in conservative form:

$$\frac{1}{c}\frac{DI}{Dt}+\frac{1}{r^2}\frac{(r^2\mu I)}{r}+\frac{1}{r}\frac{((1\mu ^2)I)}{\mu }+\sigma I=S.$$
(3)
This form expresses the fact that the transport operator $`𝛀`$ is the divergence of the flow in direction $`𝛀`$. Straightforward integration of this differential equation over a (small) finite volume ($`2\pi r^2drd\mu `$) provides a conservative scheme in which appears creation and extinction terms, along with the boundary fluxes.

In 2D axisymmetric geometry, the transport equation we solve becomes

$$\frac{1}{c}\frac{DI}{Dt}+\mu \frac{I}{r}+\eta \frac{I}{z}\xi \frac{I}{r\varphi }+\sigma I=S$$
(4)
or in conservative form :

$$\frac{1}{c}\frac{DI}{Dt}+\frac{1}{r}\frac{(r\mu I)}{r}+\eta \frac{I}{z}\frac{1}{r}\frac{(\xi I)}{\varphi }+\sigma I=S.$$
(5)
Here, the directional cosines are just the projections of the unit vector $`𝛀`$ (in the neutrino/photon velocity direction) on the Cartesian axes (Figure Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context):

$$\eta =\mathrm{cos}(\theta )$$

$$\mu =\mathrm{sin}(\theta )\mathrm{cos}(\varphi )$$

$$\xi =\mathrm{sin}(\theta )\mathrm{sin}(\varphi ),$$

where $`\theta `$ is the polar angle and $`\varphi `$ is the azimuthal angle. Generally, we distribute $`\eta `$ evenly from -1 to 1 and $`\varphi `$ evenly from 0 to $`\pi `$, but make the number of $`\varphi `$s a function of $`\eta `$ in order to tile the hemisphere more or less uniformly. We need to calculate only on a hemisphere, and not the entire sphere, because of the imposition of axial symmetry. This is also the reason $`\varphi `$ is bounded by 0 and $`\pi `$, and by not $`2\pi `$.

Our numerical scheme uses two features of the transport operator $`𝛀`$. First, being a first-order differential operator along the direction $`𝛀`$, the transport equation can be integrated along a characteristic that traverses the computational domain. Second, this operator is the divergence of the velocity $`𝐯_\nu =c𝛀`$, so that eq. (1) can be regarded as a continuity equation that reflects a conservation law for the number of particles. Thus, our aim was to construct a conservative difference scheme, unlike the method of short characteristics (e.g., van Noort, Hubeny, and Lanz 2002), that is motivated in part by the prevalence in astrophysics of radiation hydrodynamics problems for which the coupling between the radiation field and matter is stiff. In these cases, small relative errors in the energy budget can cause big errors in the dynamics. The core-collapse supernova problem, in which the neutrino radiation field contains a non-trivial fraction of the energy of the system, is a good example of this type of physical situation.

In VULCAN/2D, we discretize the angular variables using the discrete-ordinates (Sn) method, where the specific intensity depends on a set of representive directions (points on the unit sphere) (Mezzacappa and Bruenn 1993). Multi-group Sn methods can be constructed in an efficient way for general grids in one, two, or three spatial dimensions. The salient characteristic of Sn methods is the representation of the angular variable by a set of direction cosines on the unit sphere. Unlike other methods, such as the tangent-ray method, the parametrization of the angular variable is independent of position and the same set of directions applies for the entire domain. This feature makes it easy to generalize the Sn method to any dimension and to any structured or unstructured grid. A drawback of the Sn method is its tendency at low optical depths to concentrate flux along the rays. However, in the core-collapse context it is the radiation field interior to the shock wave that is most germane, and this region is at low to modest optical depths where the anisotropy of the radiation field is not severe. In this work, we limit the discussion to the 1D and 2D cases only, since the computing resources needed for an accurate 3D solution to the time-dependent transport equation are currently unavailable.

## 3 The Radiation Hydrodynamics Scheme: Matter/Radiation Coupling, Parallelization, and an Aside on Gravity

The numerical scheme used in VULCAN consists of a Lagrangian step, followed by a remapping step to the Eulerian grid (Livne 1993). This makes the code similar in this regard to traditional ALE (Arbitrary-Eulerian-Lagrangian) codes. The hydrodynamical variables are all cell-centered, except for the position and the velocity, which are node-centered. The variables of the radiation field are also cell-centered, so that the interaction between the radiation field and matter is properly centered. For radiation hydrodynamics, it is convenient to define the specific intensity as the energy specific intensity, rather than the particle specific intensity. Thus, $`I_g`$ is the energy intensity in group $`g`$ (averaged over the group) and $`I_g/\epsilon _\nu ^g`$ is the particle intensity in group $`g`$ (averaged over the group), $`\epsilon _\nu ^g`$ being the average particle (e.g., neutrino) energy of group $`g`$.

The time advancement in both the radiative and the hydrodynamical sectors is computed in the Lagrangian step, whereas the remapping step changes only the spatial discretization of the variables over the numerical grid. Therefore, we describe here only the Lagrangian scheme. The transport equation itself, for a given source, is always computed in a fully implicit manner.

We first advance the velocity by half a timestep:

$$𝐯^{n+1/2}=𝐯^n+0.5\mathrm{\Delta }t(\frac{p^n}{\rho ^n}U_G),$$
(6)
where $`U_G`$ is the gravitational potential (but see the end of this section). The position vector is then advanced using

$$𝐫^{n+1}=𝐫^n+\mathrm{\Delta }t𝐯^{n+1/2}.$$
(7)
Denoting by $`V`$ the volume of a cell, Lagrangian mass conservation takes the form:

$$\rho ^{n+1}=\rho ^n\frac{V^n}{V^{n+1}}.$$
(8)
We then solve the adiabatic energy equation for the specific internal energy:

$$e^{}=e^n1/2(p^{}+p^n)(\frac{1}{\rho ^{n+1}}\frac{1}{\rho ^n}).$$
(9)
Equation (9) is iterated to convergence. At this stage, we compute new cross sections and emission sources

$$S_{em}=\underset{g}{}\sigma _g^aJ_g^{eq},$$
(10)
where $`J_g^{eq}`$ is the LTE intensity (a function of density, temperature, and composition). Using those cross sections and sources, we solve the transport equation for $`I_g^{n+1}`$ according to the scheme described in Livne et al. (2004). The net change in the radiation energy density is given by:

$$\mathrm{\Delta }E_r=\mathrm{\Delta }t\underset{g}{}\sigma _g^a(J_g^{eq}J_g^{n+1})$$
(11)
and this is also minus the net change in the matter energy density. Consequently, we can compute the final energy density, pressure, and temperature using

$$e^{n+1}=e^{}\mathrm{\Delta }E_r/\rho ^{n+1}$$
(12)
and the equation of state.

For the supernova problem, we also need to compute the degree of neutronization of matter due to electron capture and other charged-current processes (Burrows and Thompson 2003). We obtain

$$Y_e^{n+1}=Y_e^n\mathrm{\Delta }t[\underset{g}{}\sigma _g^a(J_g^{eq}J_g^{n+1})/\epsilon _\nu ^g]\frac{1}{N_a\rho },$$
(13)
where $`N_a`$ is Avogadro’s number and $`Y_e`$ is the elctron fraction. Finally, we advance the velocity due to the new matter pressure by a further half timestep and due to the radiation pressure ($`𝐅_{rad}^{node}`$) by a full timestep:

$$𝐯^{n+1}=𝐯^{n+1/2}+0.5\mathrm{\Delta }t[\frac{1}{\rho ^{n+1}}(p^{n+1}+2𝐅_{rad}^{node})U_G].$$
(14)
The radiation force at grid nodes is evaluated by a simple averaging process using the radiation force at cells and the definition of the radiation flux:

$$𝐅_{rad}^{cell}=\frac{1}{c}\underset{g}{}\sigma _gI_g𝛀𝑑\mathrm{\Omega }.$$
(15)
Since we are focussed here on doing radiation/hydrodynamics in two spatial dimensions, not three, for reasonable sized grids and EOS and opacity tables we can fit an entire hydrodynamic calculation on a single processor. Because the transport solve consumes far more CPU than the hydro updating, this suggests a path to parallelism that is both simple to implement and almost perfectly scalable. Given that we will handle energy redistribution explicitly, the energy groups are not coupled during a timestep. Hence, we can parallelize in energy groups (and neutrino species) and calculate the hydro on ALL processors concurrently. Each processor handles one energy group for a given neutrino species and redundantly calculates the hydro. Interprocessor communication merely collects the transport results for the neutrino energy deposition and radiation pressure calculations that require integrals over the energy groups and angles. The result is then redistributed over the processors at the next timestep and the process is continued. No domain decomposition is required. For example, if we follow 20 energy groups and three neutrino species, only 60 processors are required. We use MPI to handle the parallelization. In practice, the communication overhead is only 2% to 8% of the total run time.

Gravity is a key force in multi-dimensional astrophysical hydrodynamics. However, many calculations in the past have employed only the monopole term and/or have complemented the gravitational force term ($`F_G`$), written as a gradient of a potential in the momentum equation (eq. 6) with a corresponding $`\stackrel{}{v}F_G`$ term in the energy equation. The latter approach is perfectly reasonable, but, given the inherently approximate nature of finite-difference realizations of the partial differential equations, does not guarantee momentum conservation nor consistency between the momentum and energy equations when written in Eulerian form. To address this, we have implemented a version of the code in which the $`z`$-component of gravity appears as the divergence of a stress tensor (Shu 1992; Xulu 2003). This ensures, in principle, the conservation of momentum in that direction, or at least guarantees that in fact the gravitational force of mass parcel “A” on “B” is equal and opposite to the gravitational force of mass parcel “B” on “A.” This has not been the case in past grid-based astrophysical codes that have addressed the supernova problem. This approach works with any algorithm employed to calculate the gravitational potential (we use a multipole expansion with $``$20 terms) and it, along with a gridding scheme such as is described in Ott et al. (2004), allows the inner core to move on the grid. All previously published grid-based results in supernova theory treated the inner core either in 1D (radially) or as a boundary of finite radius. Surprisingly, no previous such publication of grid-based results allowed the core to move, something that is necessary if we are to address the pulsar kick problem and to treat multi-dimensional hydrodynamic effects with realistic gravitational feedbacks.

## 4 Computational Issues of Multi-dimensional Problems

The central computational difficulty in solving the equations of radiation/hydrodynamics in many dimensions stems from the fact that, at the very best, the number of operations and the number of variables scales with the product of the number of generalized zones in the various dimensions involved. VULCAN/2D treats two spatial dimensions (axially-symmetric), one time dimension, one energy-group dimension, and two angular dimensions ($`\theta `$ and $`\varphi `$). This amounts to a total of 6 dimensions. In fact, the number of angular dimensions is the same as would be required in three spatial dimensions, so a full 3D simulation would involve only one more dimension and one more hemisphere in momentum space (§2).

In VULCAN/2D, the number of spatial cells ($`N_e`$) for a low-resolution simulation might be $`100\times 10010,000`$ and for a higher-resolution simulation might be $`300\times 300100,000`$. For an Sn specific intensity calculation, the number of $`\theta `$ directions might be $`n=1020`$ and the number of $`\varphi `$ directions might average 10, for a total of $`N_d15\times 10=150`$ angular bins. The minimum number of energy groups required for a good core-collapse simulation is $`N_g=1020`$ for each neutrino species (Thompson, Burrows, and Pinto 2003). This results in $`N_e\times N_g\times N_d=`$ 15,000,000 to 300,000,000 specific intensity values at each timestep and for each neutrino type. The radiation solve is always much more costly than the hydro solve and completely dominates any run. Doing the radiation solve for only an inner fraction of the spatial grid has its merits and can increase speed, but we have not yet implemented this. Note that the cost of a corresponding 3D simulation would be larger by $``$100 zones$`\times 2`$ hemispheres.

Currently, on a standard fast processor and for the time-dependent test run described in §5.2 we required $``$0.5 - 2.0 minutes per energy group per timestep. The time-dependent run employed the equivalent of $`16\times 10`$ angular bins and $`102\times 51`$ spatial zones. Since we parallelize in energy groups and neutrino species (§3) and can expect very good scalability on most clusters, this is a benchmark number for the current scheme. Accuracy and stability requirements currently impose a timestep after core bounce of $``$0.1-0.3 microseconds. Decreasing the number of angular bins ($`N_d`$) can increase the speed dramatically, but with a corresponding loss of accuracy. It is the scattering angular redistribution term in eq. (1) that poses the greatest numerical challenges. A problem that involves only absorption and emission, or that does not involve large scattering optical depths, is very much simpler and faster.

## 5 VULCAN Simulations

### 5.1 1D Comparison Test: VULCAN/1D versus SESAME

Though not the focus of this paper, but as a run up to the VULCAN/2D runs, we compared collapses using VULCAN/1D and the 1D Lagrangian radiation/hydro code SESAME, based on Feautrier variables and the tangent-ray method, as described in Burrows et al. (2000) and Thompson, Burrows, and Pinto (2003). The progenitor model used for this comparison was the core of an 11 M model from Woosley and Weaver (1995). The cross sections were obtained by interpolating in a three-dimensional table in $`T`$, $`\rho `$, and $`Y_e`$, generated using the formulations of Burrows and Thompson (2003). The effects of Pauli blocking by final-state electrons and of stimulated absorption, where appropriate, are included. However, only electron-type neutrinos and processes were incorporated for these tests. The equation of state we use is the finite-temperature liquid-drop model of Lattimer and Swesty (1991) which contains nucleons, NSE nuclei, alpha particles, photons, and electrons/positrons, the latter corrected for arbitrary degeneracy and relativity. The spherical (1D) calculations were done with the same number of radial zones (300) and energy groups (20). The latter span an energy range up to 320 MeV. In the VULCAN/1D calculations, 15 Sn angles were used, equally-spaced in $`\mathrm{cos}(\theta )`$.

Small differences in zoning and in time stepping cause slight time shifts during the collapse phase and there is a 20-millisecond shift between the two calculations in the time of bounce (SESAME: 213 msec; VULCAN: 195 msec). The entropies at bounce in the inner 0.5 M differ by 0.1 to 0.2 units. The $`Y_e`$ profiles show a corresponding difference, with the values of $`Y_e`$ differing by as much as 0.01 units (inside) to 0.05 units (at $``$0.6 M) in the inner 1.1 M. The outgoing deleptonization wave and the $`Y_e`$ trough look quite similar and the overall lepton losses to infinity in the two calculations are within 5%.

From the center to an interior mass of 1.0 M, the luminosities for the VULCAN/1D and SESAME calculations were within $``$1%, but exterior to 1.1 M they differed by up to $``$10%. This slight discrepancy in the outer regions can be traced to the inherent differences/strengths of Eulerian and Lagrangian codes, the different stencils, the inherent differences between tangent-ray and Sn methods, and differences in the distribution of radial zones after bounce. Interior to $``$1.25 M, the Eddington factors for the two runs were the same to within better than a percent, but since VULCAN/1D uses the Sn method and Sn can’t resolve a strongly forward-peaked specific intensity distribution without many, many more angular bins, in the outer transparent zones the Eddington factor in the VULCAN/1D run saturated at $``$0.9. The Eddington factors for the SESAME run achieved the correct value of 1.0. However, given that the neutrino-matter energy deposition in the so-called “gain” region occurs only in the inner semi-transparent regions, VULCAN/1D, despite its Sn method, provides in that critical region quite reasonable results. Nevertheless, updates to VULCAN that are still necessary are the energy redistribution due to inelastic neutrino-electron scattering and the Doppler compression term.

The differences between the quantities in the 1D core-collapse simulations recently compared in Liebendörfer et al. (2003) provide a context for our 1D comparisons here. The ORNL(Agile-Boltztran)/Garching(Vertex) codes correspond well qualitatively, and in the main quantitatively. No evolutionary phases are importantly discrepant. However, at a given time entropies differ by a few $`\times `$0.1, flux factors by 0.05 to 0.1, luminosities by 10-20%, shock radii by 10’s of kilometers, and velocities by 10’s of percent. In addition, there seem to be small numerical issues with all codes arising from energy conservation, shock capturing, and the handling of neutrino-matter coupling, to name but a few. Nevertheless, despite the different approaches, be they represented by SESAME, VULCAN/1D, Agile-Boltztran, or Vertex, 1D core-collapse simulations give very similar evolutions and results. Note that all these codes indicate that the naively 1D (spherical) models do not explode.

### 5.2 Two-dimensional, Multi-Group, Multi-Angle Simulation Results with VULCAN/2D

There is little in the supernova literature that addresses the issue of aspherical neutrino radiation fields in 2D (or 3D) matter distributions. Most of it is concerned with the effect of rotation on the flux contrast between the pole and the equator and none of it involves the solution of the transport or Boltzmann equation. None of it is multi-group or multi-angle. LeBlanc and Wilson (1970) performed 2D flux-limited MHD calculations to investigate the possible role of magnetically-driven axial jets during core collapse, but did not address the question of the associated anisotropic radiation field. They found that with very rapid rotation an under-energetic MHD jet punched along the poles. Symbalisty (1984) too was concerned with the role of magnetic fields during 2D collapse, but handled the neutrinos with a crude leakage scheme. Janka and Mönchmeyer (1989a,b) studied the possible ambiguity in the derivation of the total neutrino burst energy of SN1987a due to the anisotropy in the neutrino flux distribution that would be caused by rapid rotation of the protoneutron star left behind. They derived a $`\theta `$ dependence to the neutrino flux from an oblate spheroid by associating it predominating with the solid angle subtended by the spheroid at the observer. Pole-to-equator contrasts of as much as three were obtained, but no multi-D transport calculations was performed.

Recently, Kotake, Yamada, and Sato (2003) and Shimizu et al. (2001) used a similar approximate approach to calculate the flux asphericity for 2D pure-hydro rotating collapse models. They used the code ZEUS-2D, but without including the gray, flux-limited transport it contains. They estimated the position of the neutrinosphere, derived a neutrino flux enhancement at the pole, and concluded that such an asymmetry may play a role in the neutrino-driven mechanism of supernova explosions and in the distribution of the debris. We think that there may be merit in this possibility and are planning to calculate the anisotropy of the neutrino fluxes for rotating cores using the machinery of VULCAN/2D 222Rotation itself, and its effect on the hydrodynamic flow along the poles (not just the consequent enhancement of the driving neutrino fluxes), is also thought to be a potentially important feature in the mechanism of core-collapse supernova explosions (Burrows, Ott, and Meakin 2003)..

The calculations of Fryer and Heger (2000, 2D) and Fryer and Warren (2003, 3D) are the most complete explorations to date of the effects of rotation on the supernova mechanism. However, these authors used a flux-limited, gray diffusion approach and did not publish anything concerning the derived anisotropies of the radiation field, nor on the potential role of such anisotropies in the explosion. Their focus was mainly on the evolution of the angular momentum and on its interaction with the convective motions driven by neutrino heating. Finally, as mentioned in §1, Buras et al. (2003) conducted state-of-the-art 2D supernova simulations, but used multiple 1D radial Boltzmann calculations in lieu of real 2D transport. Our VULCAN/2D developments are being undertaken in part to improve upon such calculations by including the second angular dimension into full multi-group simulations performed in the 2D radiation/hydrodynamic context.

Our first 2D time-dependent radiation/hydro runs with VULCAN/2D, those that we publish here, are to explore aspects of the character of 2D transport in the core collapse context. In a sense, they are intended as a “shakedown cruise” of the code. They are not meant to cover the full evolution from collapse, through bounce, shock formation, and explosion (?). Here, we study the asymmetries of the radiation field using the multi-group, multi-angle code and explore the relationship between 2D matter distributions and the corresponding radiation fields. We also provide a first-of-its-kind glimpse at time-dependent, multi-group, multi-angle “supernova” calculations.

To this end, we mapped a 1D SESAME evolutionary calculation of a Woosley and Weaver (1995) 11 M progenitor core $``$5 milliseconds after bounce into VULCAN/2D and continued the post-bounce calculation in full 2D for a further $``$22 milliseconds. For this study, we employed only 5 $`\nu _e`$ energy groups, with group centers at 2.5, 10.8, 30, 74.6, and 177.9 MeV. In addition, we used the equivalent of $`16\times 10`$ angular bins and $`102\times 51`$ spatial zones, the latter distributed as described in Ott et al. (2004). The calculation covered $``$69000 timesteps, with the timesteps themselves ranging between $``$0.1 to $``$0.3 microseconds.

Figures Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context and Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context are snapshots of the $`Y_e`$ distribution with velocity vectors at 10.2 and 22 milliseconds into the run. The initial perturbations imposed have axial symmetry and are completely artificial. These figures are 240 km$`\times `$240 km and 420 km$`\times `$420 km, respectively, on a side. There is a slight artifact along the poles due in part to the modest angular resolution of the run and to the generic singularity of polar coordinates. The positions of the stalled shock wave (roughly spherical) are clearly seen on Figures Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context and Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context by the jumps in the behavior of the corresponding velocity fields. Inward pointing arrows in the unshocked region indicate infall and accretion. The slight moiré pattern seen in the vectors is a graphical/visual artifact and is not in the numerical data. The red and off-white colors depict regions of higher $`Y_e`$ and proton fraction, while the lowest $`Y_e`$ regions are depicted in green and blue. The characteristic trough in $`Y_e`$ between the shock and the inner core is thereby clearly indicated. The most salient and interesting feature of the flow pattern is the pronounced vortical motion traced by the velocity vectors between 15 and 100 kilometers radius. Such prompt convection is expected right after core bounce (Burrows and Fryxell 1992; Burrows, Hayes, and Fryxell 1995). The motion is clockwise in the positive quadrant and counterclockwise in the negative quadrant, with a reflected counterpart in the bottom hemisphere. The extensions of the high-$`Y_e`$ regions follow the velocity vectors. A slight top/bottom asymmetry imposed in the initial model has translated into a slight top/bottom flow asymmetry which persists, quite naturally, during the run. In Fig. Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context, weaker vortices are seen to be developing further out. Due to the low spatial resolution of the grid employed, it is unlikely that the flow field has converged and only the grossest overall convective/overturning pattern should be considered robust.

However, our purpose here is to determine the 2D radiation fields that result from aspherical matter and composition fields in time-dependent hydrodynamic supernova flows. To depict these we first show the electron neutrino flux spectrum vectors at 10.8 MeV. These are provided at the same times as those chosen for Figs. Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context and Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context. Note that 10.8 MeV is roughly at the peak of the emergent spectrum. Figures Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context and Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context are color maps of $`Y_e`$ at the fiducial times of 10.2 and 22 milliseconds, but with the flux spectrum vectors at 10.8 MeV superposed. The shock wave is in the whitish region (high-$`Y_e`$) just exterior to the outer red band. Its positions are more easily seen in Figs. Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context and Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context. Figure Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context depicts the 10.8-MeV flux vector distribution at 22 milliseconds, but superposed on the corresponding entropy color map.

As expected, at the high opacities in the inner core the fluxes are very small. Elsewhere, the flux vectors generally point outwards, and at large distances they are more or less radially-oriented with magnitudes that diminish as $`\frac{1}{r^2}`$. However, there are many intriguing exceptions to this pattern. In particular, as seen in Fig. Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context and Fig. Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context at 22 milliseconds and between 15 and 80 kilometers (a region that brackets the corresponding neutrinosphere) the 10.8 MeV fluxes generally point away from the region (red) of high $`Y_e`$/proton fraction and higher entropies. In fact, at 30-40 kilometers (at modest to high optical depth), the flux vectors can be 45 off the radius vector. This behavior is a consequence of the higher emissivities due to electron capture in the $`Y_e`$/proton-rich, higher-entropy regime. As Fig. Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context indicates, due to the general vortical circulation pattern set up in this simulation, the entropy distribution roughly follows the $`Y_e`$ distribution. In the thick regions, the fluxes tend (at least in part) to track the local emissivity, while at low optical depths in the semi-transparent region the specific intensities from many angles “vector” add at a given point and result in a smoothed partially radially streaming radiation pattern. The tendency for the radiation field at large distances to be more smooth and spherically-distributed than the matter in the inner regions is a generic feature of multi-D transport. Such an effect is not easily nor accurately captured with simpler transport schemes (e.g., Burrows, Hayes, and Fryxell 1995; Rampp and Janka 2002; Fryer and Warren 2002; Buras et al. 2003). Even so, even at 90-100 kilometers and exterior to the 10.8 MeV neutrinosphere, the radiation field is not radial. As Fig. Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context shows, it is enhanced at angles nearest the inner proton-rich environments near 50-60 kilometers.

One measure of anisotropy is the ratio of the radial component of the flux to its magnitude ($`\mathrm{cos}\theta _r=`$$`F_r/\left|\stackrel{}{F}\right|`$). Figure Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context plots $`\mathrm{cos}\theta _r`$ versus polar angle at 10.8 MeV and 22 milliseconds. As Fig. Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context shows, this ratio at 60 km can vary by more than 40%, while at 96 km it varies by $``$5%. The angular scale of the variation follows the angular scale of the matter asphericity; in this calculation at this physical time this is $``$20-30. As Fig. Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context shows, the degree of asphericity, as measured by $`\mathrm{cos}\theta _r`$, decreases with increasing radius, demonstrating the “smoothing” effect at low optical depths.

The dependence of the flux vector and magnitude distributions on neutrino energy is revealed in Figs. Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context and Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context. These figures are at 22 milliseconds and are similar to Fig. Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context, but are for neutrino energies of 2.5 MeV and 30 MeV, respectively. The low-energy neutrinos are seen to assume a more radial pattern starting deeper in the core, while the fluxes at the higher energies (and those on the tail of the emergent spectrum) are strongly concentrated around the $`Y_e`$ contours and in the optically thick regions. The latter pattern is a consequence of the higher opacities at these higher energies. The lower opacities at 2.5 MeV cause this component to decouple deeper in and to assume a more radial flux (and energy density) pattern at a given radius. Not only are the first-moment (flux) and zeroth-moment (energy density) anisotropies in the interior a function of neutrino energy and radius, but the anisotropy of the emergent spectrum at infinity is a function of energy as well. This means that the spectrum at a given radius is generically a weak to modest function of angle and will reflect the anisotropy of the matter due to rotation and convective overturn interior to the shock and around the respective neutrinospheres. For instance, at 22 milliseconds and 60 km, $`\mathrm{cos}\theta _r`$ varies by as much as $``$40% at 10.8 MeV (as Fig. Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context indicates), but only by $``$25% at 30 MeV and by $``$5% at 2.5 MeV. At 96 km, the corresponding numbers for 2.5, 10.8, and 30 MeV are 0.5%, 5%, and 15%, respectively, while at 125 km they are all less than 2%. Quite naturally, $`\mathrm{cos}\theta _r`$ goes to 1.0 at large distances. However, another measure of anisotropy is the variation with angle of the magnitude of the flux itself ($`\left|\stackrel{}{F}\right|`$) and this does not flatten at large distances. Figure Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context depicts for radii from 60 to 125 km the dependence with polar angle of the magnitude of the flux at 10.8 MeV and 22 milliseconds. In particular, at 125 km $`\left|\stackrel{}{F}\right|`$ for 10.8 MeV varies by as much as 80% from 0 to 40 polar angle and at 250 km it varies by 10-15% from top to bottom. At very large radii ($``$infinity), $`\left|\stackrel{}{F}\right|`$ at 10.8 MeV varies by $``$15% from top to bottom. While not small, a much more anisotropic matter distribution, such as what rotation or a large initial core anisotropy can generate, is required to produce a larger degree of flux anisotropy at infinity. Nevertheless, the variations with angle and radius we identify here are important new components of supernova theory that previous simulations could not calculate.

Figures Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context and Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context depict the associated maps for the energy-integrated zeroth-moment (total $`\nu _e`$ energy density) at 10.2 and 22 milliseconds, respectively. Superposed on these plots are the velocity vectors. These distributions, and not the flux maps, are rough guides to the distribution of neutrino heating thought to be important in the neutrino-driven mechanism of core-collapse supernova explosions (Bethe and Wilson 1985; Burrows, Hayes, and Fryxell 1995; Liebendörfer et al. 2001; Fryer and Warren 2002; Rampp and Janka 2002). However, a more interesting plot is Fig. Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context, which shows the relationship at 10.8 MeV between the flux vector field and the corresponding energy-density pattern. In this figure, the scale is 200 km$`\times `$200 km. Combined in this way we can see more directly the correlation between a flux and an energy-density pattern. In addition, though not obvious in Figs. Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context through Two-dimensional, Time-dependent, Multi-group, Multi-angle Radiation Hydrodynamics Test Simulation in the Core-Collapse Supernova Context, the ray effects inherent in the Sn method do not manifest themselves for this $`N_\theta =16`$ calculation until radii of $``$350 kilometers, significantly exterior to the radius of the stalled shock.

## 6 Conclusions

We have constructed a time-dependent, multi-group, multi-angle radiation code in two spatial dimensions and coupled it to a hydrodynamics code. We have used the product, VULCAN/2D, to simulate the early post-bounce evolution of a supernova core for 22 milliseconds. The purpose of this exercise was to explore for the first time the mapping between the aspherical, anisotropic matter distributions that obtain in the convective protoneutron star and the corresponding radiation field. Highlighted were the electron neutrino flux and energy-density distributions and their angle dependence. We found that these radiation distributions and their anisotropies are systematic functions of neutrino energy and radius and that, though the radiation field depends on the distributions of the physical parameters, the non-local character of the solution to the Boltzmann equation results in smoother radiation patterns in optically-thin regions. In addition, we found that the emergent spectrum should be angle-dependent. Even though the code was constructed to investigate core-collapse supernovae, it can in principle be used for other astrophysical problems. These might include gamma-ray bursts, accretion disk evolution, and star formation.

This paper has been a preliminary exploration and test of the new VULCAN/2D tool. More details on the solution method and numerical approach can be found in Livne et al. (2004). However, this study is but the beginning of a program of development and calculation which will include the incorporation of further velocity-dependent terms, the upgrade of the code to handle energy redistribution in an explicit manner, full simulations with the MGFLD variant, simulations with rotation (cf. Ott et al. 2004), and simulations with more energy groups and neutrino species. Our philosophy has been to develop in systematic fashion, and with a realistic timetable, an evolvingly sophisticated computational capability in multi-dimensional core-collapse supernova studies, in particular, and in astrophysical radiation hydrodynamics, in general. Importantly, we are gearing up to address the neutrino-driven mechanism of core-collapse supernovae by calculating for much longer times, for different progenitors, and with more energy groups and neutrino species. The code’s efficient parallelization in group and species implies that higher-energy-resolution simulations incur little penalty. Finally, since the code follows transport in the angular, as well as the radial directions, it can explore in numerical detail and for the first time the viabilility of neutron-finger or doubly-diffusive instabilities in the core-collapse context (Mayle and Wilson 1988).

We acknowledge discussions with Christian Ott, Jeremiah Murphy, Casey Meakin, Yeliang Zhang, Salim Hariri, Ron Eastman, Marvin Landis, and Stan Woosley. Support for this work is provided in part by the Scientific Discovery through Advanced Computing (SciDAC) program of the DOE, grant number DE-FC02-01ER41184. T.A.T. is supported by NASA through Hubble Fellowship grant #HST-HF-01157.01-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555. In addition, we thank Jeff Fookson and Neal Lauver of the Steward Computer Support Group for their invaluable help with the local Beowulf cluster and acknowledge the use of the NERSC/LBNL/seaborg and ORNL/CCS/cheetah machines. Movies associated with this work can be obtained from AB, or at http://zenith.as.arizona.edu/~burrows after publication.


text:
# On the uniform generation of random graphs with prescribed degree sequences

## I Introduction

In the rapidly growing literature on the modeling of complex networks one of the most important classes of network models is the random graph Bollobas2001 . One well-studied such model is the model consisting of the ensemble of all graphs that have a given degree sequence Bender ; Molloy 1995 ; Molloy 1998 ; Newman 2001 ; Chung\_diameter , and this model has proved useful in understanding a variety of network properties. Realistic applications often require that we restrict ourselves to graphs with no multiple edges between any vertex pair and no self-edges. Unfortunately, both the analytic and numerical study of such networks is known to present challenges Bender ; Snijders ; Rao ; Roberts ; Kannan ; chen ; Itzkovitz ; MSZ02 ; Park . In this short paper we consider computer algorithms for generating graphs uniformly from this ensemble. We are concerned primarily with directed graphs, since the examples we will consider are directed, but the concepts discussed generalize in a straightforward fashion to the undirected case also.

There are two algorithms in common use for the generation of random graphs with single edges. We will refer to them as the switching algorithm Rao ; Roberts ; Newman 2002 ; Maslov 2002 ; Stone ; Shen or 2002 ; Milo 2002  and the matching algorithm Molloy 1995 ; Newman 2001 ; Milo 2002 . We argue that, under certain circumstances, both of these algorithms can generate a nonuniform sample of possible graphs. We then present a new algorithm based on the Monte Carlo procedure known as go with the winners Aldous ; Grassberger , which generates uniformly sampled graphs. We compare the three methods in the context of a particular network problem—estimation of the density of commonly occurring subgraphs or motifs—and show that, in this context, the difference between them is small. This result is of some practical importance, since the “go with the winners” algorithm, although statistically correct, is slow, while the other two algorithms are substantially faster.

## II Algorithms

In this section we describe the three algorithms under consideration.

### II.1 Switching algorithm

First, we describe the switching algorithm, which uses a Markov chain to generate a random graph with a given degree sequence Rao ; Roberts ; Newman 2002 ; Maslov 2002 ; Stone ; Shen or 2002 ; Milo 2002 . For simplicity, we discuss directed networks with no mutual edges (vertex pairs with edges running in both directions between them). The case with mutual edges is a simple generalization Roberts .

The method starts from a given network and involves carrying out a series of Monte Carlo switching steps whereby a pair of edges $`(AB,CD)`$ is selected at random and the ends are exchanged to give $`(AD,CB)`$. However, the exchange is only performed if it generates no multiple edges or self-edges; otherwise it is not performed. The entire process is repeated some number $`QE`$ times, where $`E`$ is the number of edges in the graph and $`Q`$ is chosen large enough that the Markov chain shows good mixing. (Exchanges that are not performed because they would generate multiple or self-edges are still counted to insure detailed balance footnote .)

This algorithm works well but, as with many Markov chain methods, suffers because in general we have no measure of how long we need to wait for it to mix properly. Theoretical bounds on the mixing time exist only for specific near-regular degree sequences Kannan . We empirically find, however, that for many networks, values of around $`Q=100`$ appear to be more than adequate (see Fig. 1).

### II.2 Matching algorithm

An alternative approach is the matching algorithm Molloy 1995 ; Newman 2001 ; Milo 2002 , in which each vertex is assigned a set of “stubs” or “spokes”—the sawn-off ends of incoming and outgoing edges—according to the desired degree sequence. (One can also assign mutual-edge stubs for networks that include such edges.) Then in-stubs and out-stubs are picked randomly in pairs and joined up to create the network edges. If a multiple or self-edge is created, the entire network is discarded and the process starts over from scratch.

This process will correctly generate random directed graphs with the desired properties. Unfortunately, however, many real-world networks have a heavy-tailed degree distribution that includes a small minority of vertices with high degree. All other things being equal, the expected number of edges between two such vertices will often exceed one, making it unlikely that the procedure above will run to completion, except in the rarest of cases. To obviate this problem a modification of the method can be used in which, following selection of a stub pair that creates a multiple edge, the network is not discarded, and an alternative stub pair is selected at random. In general this method generates a biased sample of possible networks King  but, as we will show, not significantly so for our purposes (see Table 1).

### II.3 Go with the Winners algorithm

The “go with the winners” algorithm is a non-Markov-chain Monte Carlo method for sampling uniformly from a given distribution Aldous ; Grassberger . When applied to the problem of graph generation, the method is as follows. We consider a colony of $`M`$ graphs. As with the matching algorithm, we start with the appropriate number of in-stubs and out-stubs for each vertex and repeatedly choose at random one in-stub and one out-stub from the graph and link them together to create an edge. If a multiple edge or self-edge is generated, the network containing it is removed from the colony and discarded. To compensate for the resulting slow decline in the size of the colony, its size is periodically doubled by cloning each of the surviving graphs; this cloning step is carried out at a predetermined rate chosen to keep the size of the colony roughly constant on average. The process is repeated until all stubs have been linked, then one network is chosen at random from the colony and assigned a weight:

$`W_i=2^c{\displaystyle \frac{m}{M}},`$ (1)
where $`c`$ is the number of cloning steps made and $`m`$ is the number of surviving networks. The mean of any quantity $`X`$ (for example, the number of occurrences of a given subgraph) over a set of such networks is then given by

$`{\displaystyle \frac{_iW_iX_i}{_iW_i}},`$ (2)
where $`X_i`$ is the value of $`X`$ in network $`i`$.

## III Comparison of algorithms

In Fig. 2 we show a comparison of the performance of our three algorithms when applied to a simple toy network. The network consists of an out-hub with ten outgoing edges, an in-hub with ten incoming edges, and ten nodes with one incoming edge and one outgoing edge each. Given this degree sequence, there are just two distinct network topologies with no multiple edges, as shown in Fig. 2a and 2b. There is only a single way to form the network in 2a, but there are 90 different ways to form 2b.

We generated $`\mathrm{100\hspace{0.17em}000}`$ random networks using each of the 3 methods described here and the results are summarized in Fig. 2c. As the figure shows, the matching algorithm introduces a bias, undersampling the configuration of Fig. 2a. This is a result of the dynamics of the algorithm, which favors the creation of edges between hubs. The switching and “go with the winners” algorithms on the other hand sample the configurations uniformly, generating each graph an equal number of times within the measurement error on our calculations. The “go with the winners” algorithm truly samples the ensemble uniformly but is far less efficient than the two other methods. The results given here indicate that the switching algorithm produces essentially identical results while being a good deal faster. The matching algorithm is faster still but samples in a measurably biased way.

Now consider the study of network motifs. We are interested in knowing when particular subgraphs or motifs appear significantly more or less often in a real-world network than would be expected on the basis of chance, and we can answer this question by comparing motif counts to random graphs. Some results for the case of the “feed-forward loop” motif Shen or 2002 ; Milo 2002  are given in Table 1. In this case the densities of motifs in the real-world networks are many standard deviations away from random, which suggests that any of the present algorithms is adequate for generating suitable random graphs to act as a null model, although the “go with the winners” and switching algorithms, while slower, are clearly more satisfactory theoretically. The matching algorithm was measurably nonuniform for our toy example above, but seems to give better results on the real-world problem.

Overall, our results appear to argue in favor of using the switching method, with the “go with the winners” method finding limited use as a check on the accuracy of sampling. Accuracy checks are also supplied by analytical estimates for subgraph numbers Itzkovitz . Numerical results in Milo 2002 ; Milo 2004 ; Itzkovitz  were done using the switching algorithm.

## IV Conclusions

In this paper we have compared three algorithms for generating random graphs with prescribed degree sequences and no multiple edges or self-edges. Two of the three have been used previously, but suffer from nonuniformity in their sampling properties, while the third, a method based on the “go with the winners” Monte Carlo procedure, is new and samples uniformly but is quite slow. Of the two older algorithms, we show that one, which we call the “matching” algorithm, has measurable deviations from uniformity when compared to the “go with the winners” method, although for graphs typical of practical studies these deviations are small enough to make no significant difference to most previously published results. The other older algorithm, which we call the “switching” algorithm and which is based on a Markov chain Monte Carlo method, samples correctly in the limit of long times and in practice is found to give good results when compared with the “go with the winners” method. Overall, therefore, we conclude that the switching algorithm is probably the algorithm of choice, with the “go with the winners” algorithm finding a supporting role as a check on uniformity, although its slowness makes it impractical for large-scale use.
We thank Oliver D. King for discussions and for pointing out and demonstrating that the matching algorithm of the supplementary online material of Milo 2002  does not uniformly generate simple graphs.