text:
# Untitled Document

On the paper by M. Valluri and D. Merritt ”Orbital Instability and Relaxation in Stellar Systems”

V.G.Gurzadyan

Yerevan Physics Institute and University of Rome

The paper by M.Valluri and D.Merritt is discussing an interesting and difficult topic, and obviously many questions remain open. I hope, my comments will be helpful to understand the reasons of the mentioned ”debates”.

1. It is not precise to state that ”exponential instability of N-body problem was shown by Miller , Gurzadyan and Savvidy ..”. The exponents in  and  are different quantities and have different content.

Miller  in his pioneering numerical work has shown the exponential growth of separation of coordinates and velocities for N-body systems. It was indeed a remarkable work of early period of computer astronomy. Due to the development of chaotic dynamics during these decades, our present understanding of that problem is much deeper.

In  we have analytically proved the exponential instability of spherical N-body systems in Riemannian space with metric depending on the interaction potential, and the positivity of KS-entropy determined by the curvature (the curvature was zero in Miller’s case). This, according to the theorems of ergodic theory, means the existence of the property of mixing and hence tending to an equilibrium (local). KS-entropy is thus defining the time scale of the tending to equilibrium for large N limit, i.e. the relaxation time; for discussion of the meaning of relaxation time see .

Mixing therefore is not just ”a useful way of thinking” but a strictly defined property of dynamical systems, with sufficient condition (for typical automorphisms) - the positivity of KS-entropy.

Miller’s exponent exists even for N=2, i.e. for two-body problem which is integrable, with time scale obviously equal to the crossing/orbital time.

While exponent in  at N=2 turns to infinity and hence – no mixing, no relaxation. The systems of 32 particles shown to be exponentially unstable by Miller , are not unstable via criterion . Indeed, what is the meaning of relaxation time for 24 or 32 particles?

The time scale obtained numerically by Miller, Valluri and Merritt, as well as by other authors, must always be of the order of crossing/dynamical time, as actually they had obtained, with no direct relation with relaxation time scale. Via such numerical experiments hence, for example, nobody had derived (and cannot derive), for example, the two-body relaxation time scale, though nobody seems to be alarmed by it.

The exponent in  describes mixing and hence can be associated with relaxation, the exponent in  - no. Thus, not every exponential instability means relaxation (for clear example see ). The ’debates’ quoted by Valluri and Merrit are therefore due to misunderstanding of this crucial point.

Meanwhile, the relaxation time scale obtained in  is supported by:

(a) observational data on globular clusters (see refs in );

(b) alternative statistical mechanical methods ;

(c) numerical experiments .

2. The problem of numerical estimation of Lyapunov numbers for N-body gravitating systems is itself a non-trivial one, even forgetting that they are limiting in time quantities. In  a theorem is proved showing the non-equivalence of those Lyapunov numbers and their computer images.

3. The massive central object is known to increase the instability and mixing properties of N-body systems .

 R.H.Miller, ApJ, 140, 250, 1964.

 V.G.Gurzadyan, G.K.Savvidy, Doklady AN SSSR, 277, 69, 1984; AA, 160, 203, 1986.

 V.G.Gurzadyan, in: Ergodic Concepts in Stellar Dynamics, (Gurzadyan V.G., Pfenniger D., eds.), Springer-Verlag, 1994.

 V.G.Gurzadyan, A.A.Kocharyan, AA, 205, 93, 1988.

 G.Meylan (this volume); astro-ph/9912495.

 M.Antoni, S.Ruffo, A.Torcini (this volume); cond-mat/9908336.

 A.A.El-Zant, AA, 326, 113, 1997.

 V.G.Gurzadyan, A.A.Kocharyan, J.Phys.A, 27, 2879, 1994.

 V.G.Gurzadyan, A.A.Kocharyan, Doklady AN SSSR,301, 323, 1988; V.G.Gurzadyan, S.G.Matinian, A.A.Kocharyan, Doklady AN SSSR, 296, 54, 1987; A.A.El-Zant, V.G.Gurzadyan, Physica, D, 122, 241, 1998.


format: text

text:
# I Introduction

## I Introduction

The scientific and technological relevance of a-Si, pure and hydrogenated, is well known and needs not be emphasized here; early work on the atomic structure of the pure amorphous material started more than four decades ago, both experimentally and theoretically. Experimentally, work on the amorphous phases of pure germanium and silicon evolved in parallel, beginning with the electrolytic approach of Szekely  in 1951 and the pioneering efforts of Richter and Breitling  in 1958. Theoretically, the first atomic models of both a-Si and a-Ge appeared in the literature over thirty years ago. Grigorovici and Manaila  and Coleman and Thomas  in 1968 suggested structures based on arrangements of closely packed simplified Voronoi polyhedra that have the shape of truncated tetrahedra both eclipsed and staggered $`60^{}`$ about their common bond, leading to rings of five atoms and to boatlike rings of six atoms, as in the carbon and silicon amorphous clusters that we have recently studied.  As is well known, fivefold symmetry is non-crystallographic and therefore yields an amorphouslike diffraction pattern with broad maxima. The eclipsed configuration with a fivefold symmetry structure is energetically unfavorable in the crystalline phases but occurs in the amorphous form since atomic arrangements with a large internal energy can appear in such frozen structures.

It was soon realized that such pure structures were incapable of being doped with donors or acceptors, due to the presence of dangling bonds that masked the existence of an energy gap in the density of states (DOS) spectrum, and the search began to find a solution to this problem. Spear and Lecomber  discovered that if amorphous silicon is grown in the presence of hydrogen the naturally existing dangling bonds passivate revealing, in this manner, the presence of the gap. It was then possible to introduce dopants whose influence was clearly seen. This radically transformed the study of amorphous semiconductors from an academic subject to a technologically relevant one that has led to the design of devices that are presently used in many applications. In this manner the study of hydrogenated amorphous silicon flourished  and became one of the most important subjects in the area of amorphous semiconductors.

All these materials have been experimentally produced and characterized without theoretical guidelines that atomistic simulations may provide. In particular, several experimental radial distribution functions (RDFs) have been obtained for pure amorphous silicon, and some for hydrogenated amorphous silicon, using a variety of techniques from X-ray to neutron diffraction, in order to better comprehend their atomic structure. A consistent picture has emerged for the structure of a-Si. However, for a-Si:H the neutron diffraction experiments needed to probe the hydrogen presence are relatively more recent  and the information is scarce; nevertheless some interesting structural features have been identified for this system, which makes it particularly appealing for theoretical simulations.

It is clear that, in principle, a full theoretical description of the properties of an amorphous solid depends on a complete knowledge of the atomic structure. However, for these solids there are an infinite number of possible structures and the best we can do in order to characterize their atomic arrangement is to use the RDF, also known as pair distribution function. Defined as $`4\pi r^2\rho (r)dr`$, where $`\rho (r)`$ is the density of particles at point $`r`$, the RDF gives the average number of pairs of atoms separated by a distance lying between $`r`$ and $`r+dr`$. A given structure generates a well defined RDF but a given RDF corresponds to many atomic structures; that is, there is not a one to one correspondence between a given structure and its corresponding RDF. For this reason, clusters with a particular local structure have been used to simulate the amorphous bulk and this has shed light on the effect of specific atomic arrangements on the electronic properties of these materials.

The experimental work of Refs. 10 to 16 for a-Si and Refs. 17 and 18 for a-Si:H will be considered in this paper. For hydrogenated silicon, Ref. 18 reports the only complete study of the total RDF, together with the partials for Si-Si, Si-H and H-H.

On the theoretical side, a considerable amount of work has been done to simulate the atomic structure of a-Si, while the atomic structure of a-Si:H has been the subject of more recent work. A knowledge of their atomic topology leads to better and more realistic calculations of the electronic structure of these systems. These efforts can be classified via two extreme types: i) calculations that are carried out in samples that are constructed essentially “to order” by switching bonds and adjusting dihedral angles, and that use ad hoc classical, parameter-dependent potentials constructed for the specific purpose of describing them; ii) quantum methods, parameterized and ab initio, that can deal from the outset with the thermalization processes that generate the amorphous structure, and the study of their corresponding electronic properties. In between, one finds a variety of hybrid approaches. In the present work we shall not consider the beautiful, computationally cost-effective, work that has been done with classical techniques in very large supercells of a-Si that contain tens of thousands atoms (see for instance Ref. 19 and work cited therein), since such methods are not transferable to other amorphous materials. Moreover, it is necessary to generate new ad hoc parameters and potentials in order to deal with new materials, either from experimental results or from ab initio calculations. Furthermore, it has been established  that even though some of these potentials describe the RDF, they fail in differing degrees to describe other properties of the amorphous lattice, and this is understandable since no matter how good a classical approach may simulate the atomic interactions, there will be quantum properties beyond the capability of these techniques.

Quantum methods answered some of the unsolved questions left by the classical approaches, and were themselves of several kinds. For example, there has been some interesting work done using tight binding methods for pure silicon,  where a transferable model is found by fitting the energies of silicon in various bulk crystal structures and examining functional parameterizations of the tight binding forms. For hydrogenated silicon a transferable model has also been found by fitting it to silane.  On the other hand, there are ab initio methods that attempt to answer all the questions from first principles and are generally applicable without adjustment of parameters, but are very demanding in computer resources and so are limited to handling relatively small amorphous cells. The question is, how generally can an ab initio method that uses a reasonably sized supercell be applied, and how accurate can one expect the results to be? The present work addresses these issues for amorphous silicon and for hydrogenated amorphous silicon.

## II Antecedents

For our purposes, the theoretical calculations described in Car and Parrinello  through to Lee and Chang for a-Si and in Buda et al. through to Tuttle and Adams for a\- Si:H are appropiate since they approach the structural problem by generating amorphous cells using first-principles quantum methods.

More than a decade ago Car and Parrinello performed the first ab initio molecular dynamics (MD) study in an fcc cell with 54 atoms of silicon using their plane wave MD method. In their approach a non-local pseudopotential was used together with the parameterized local density (LDA) form of Perdew and Zunger for the exchange-correlation effects. They obtained good agreement with the experimental RDF of a-Ge, rescaled to simulate a-Si, up to the second radial peak and argued that because of the size of the cell used distances larger that 6 Å could not be studied. They pointed out that the comparison of the simulated atomic structure to the experimental ones should be done with care since a large number of defects were found in their results. The simulation was started above the melting point, at about 2,200 K, and then the liquid was allowed to evolve for $`0.7\times 10^{12}`$s before it was quenched down to $``$ 300 K at a cooling rate of $`2\times 10^{15}`$ K/s. During the initial quenching the volume of the cell was gradually increased to 1080 Å3, the crystalline value. This technique of quenching from the melt has been used in subsequent work although handling the transition from the liquid to the amorphous phase is not an easy task since a volume change has to be dealt with and because liquid silicon is metallic with an average coordination number of between 6 and 7, the quenching preserves some of this overcoordination.

For example, Drabold et al. use a density functional theory, local density approximation (DFT-LDA), molecular dynamics approach developed by Sankey et al. based on a simplification of the Kohn-Sham equations as developed by Harris, and starting with a 64-atom simple cubic cell in the diamond structure with one vacancy, they generate an “incompletely melted” sample by heating it up to 8,000 K. The free evolution of the cell then results in the system acquiring a highly disordered liquidlike structure before final quenching to a solid. They state that their results for the RDF agree well with experiment, without making a detailed, direct comparison. They find only four coordination defects, two dangling and two floating bonds, for appropiate values of the rate of free evolution of the cell. After annealing at 300 K only two defects survived.

A more complete report of the Car and Parrinello results is given in Stich et al., where a cooling rate from the melt of $`10^{14}`$ K/s is used. This slower cooling rate seems to be sufficient to recover a tetrahedral network that nevertheless contains several coordination defects as well as a large fraction of distorted bonds. Annealing at 900 K reduces the defects and the RDF they obtain has two peaks that seem to agree with the first two peaks of experiments. The study performed by Drabold et al. was extended  to the 216 atom periodic supercell of Wooten, Winer and Weaire  and a more complete analysis of the relationship of structural defects, spectral defects and interatomic distances was carried out. Lee and Chang perform ab initio simulations on a 64-atom silicon cell and quench incompletely melted samples as in Ref. 24, but they find that the third peak of the RDF practically disappears, and that more dangling and floating bonds occured than in previously generated samples from liquid-quenched simulations, and than in the samples generated by Drabold et al.

Theoretical work on a-Si:H has been less abundant because the experimental RDF results are limited and it is more difficult to model interactions of H and Si and cope with the time steps needed for computer simulations of hydrogen diffusion. In addition, the strong dependence of the a-Si:H structure on deposition conditions has to be reflected in the simulations, together with the chemical reactivity of hydrogen and its zero point energy; all this requires the use of quantum mechanical methods. The ab initio work performed for this system is found in Refs. 29 to 31. In some of the models H is incorporated into a previously created amorphous network of pure a-Si by hand,  with the consequent inhibition of hydrogen diffusion, whereas Buda et al.  and Tuttle and Adams  do allow bulk diffusion of H by first creating a liquid sample of the material and then quenching it, a procedure that has known shortcomings, such as overcoordination, no definite gap, etc.

Specifically, the plane wave MD Car-Parrinello method was applied to amorphous hydrogenated silicon by Buda et al. using a cubic cell of 64 silicons and 8 hydrogens (11% concentration). They start out with a liquid material containing both silicon and hydrogen atoms that is rapidly quenched, maintaining a density equal to the value of the crystalline material. They report only partial distribution functions and, as is usually the case, the H-H RDF is poorly reproduced and is not compared to the existing experimental neutron diffraction data. The DFT-LDA approach of Sankey and coworkers was applied to this material by Fedders and Drabold using several cells based on a two-defect 63 silicon atom supercell that was constructed in previous work. The hydrogen was introduced into the amorphous silicon sample by hand so that the H-atoms are located near (1.5 Å) the corresponding dangling bond. To eliminate the strained bond defects they removed the silicon considered to be the center of the strain defect and put H atoms near the 4 dangling bonds. Posterior relaxation allowed the hydrogens to be trapped by the dangling bonds. Fedders and Drabold do not report any RDF, total or partial. Tuttle and Adams  also use the Harris functional in the DFT-LDA code developed by Sankey et al. and apply it to a cell of 242 atoms with 11% hydrogen. They generate the structure from a liquid at $``$ 1800 K and quench it to produce an amorphous structure at 300 K. Since they were not concerned with real time dynamics, they set the mass of hydrogen equal to the mass of Si, thereby allowing the use of a large time step in the annealing process (4 fs). However, this means that the RDFs could not reflect real diffusion processes of the hydrogens in the cell, and correspondingly only the Si-Si and the Si-H RDFs are reported. A significant percentage of defects with only 90% of the silicon being fourfold coordinated is found.

Even though tight binding calculations have been labeled as “highly arbitrary and inadequate”,  some of them deserve special mention. Biswas et al.  using this approach studied the electronic structure of dangling and floating bonds in amorphous silicon. Fedders  looked into the energetics of defects and found that tight binding gives surprisingly good results for the energy eigenvalues and the degree of localization of the defect states if some radial dependence is included in the hopping matrix elements. Colombo and Maric  reported the first tight-binding molecular-dynamics simulation of the defect-induced crystal-to-amorphous transition in crystalline silicon. In particular, two recent ones are specially relevant for the present work: Yang and Singh find that the total average energy per silicon atom is a minimum when the hydrogen concentration is in the range 8-14%, the optimum experimental range found; Klein et al.  report a H-H RDF closer to experimental results than previous work. Both studies start out from liquid samples that are fast quenched to generate the hydrogenated amorphous silicon and therefore, as in all other cases, a large percentage of defects, floating or dangling bonds, are naturally created. For a recent account of the state of the art of tight binding methods see Ref. 39 where both parameterized and ab initio approaches are considered.

## III Method

In this work we report the generation of samples of both a-Si and a-Si:H using a new approach (with four different time steps) that leads to structures with a minimum of coordination defects. We use FastStructure,  a DFT code based on the Harris functional, and optimization techniques based on a fast force generator to allow simulated annealing/molecular dynamics studies with quantum force calculations.  We use the LDA parameterization due to Vosko, Wilk and Nusair (VWN)  in the simulations. The core is taken as full which means that an all electron calculation is carried out, and for our simulations a minimal basis set of atomic orbitals was chosen with a cutoff radius of 5 Å, (compare to values of $`2.6`$ Å used by Sankey et al.), a compromise between cost and accuracy. The physical masses of hydrogen and silicon are used throughout and this allows us to see realistic diffusive processes of the hydrogen atoms in the supercell. The default time step is given by $`\sqrt{m_{min}/5}`$, where $`m_{min}`$ is the value of the smallest mass in the system; however, in order to better simulate the dynamical processes that occur in the amorphisation a time step of 10 fs (compare to 2.44 fs, the default value) was also used for the pure silicon and 2 fs was used when hydrogen was introduced, in addition to 0.46 fs, the default time step. The forces are calculated using rigorous formal derivatives of the expression for the energy in the Harris functional, as discussed by Lin and Harris. The evaluation of the 3-center integrals that contribute to the matrix elements in the one-particle Schrodinger equation is the time-limiting feature of FastStructure and each is performed using the weight-function method of Delley.

Since it is clear that quenching from a melt or from partially melted samples generates undesirable structures, we took a different path. Our process, like the ones mentioned above, is not designed to reproduce the way an amorphous material is grown, but has the objective of generating an amorphous sample that would represent adequately the ones obtained experimentally. We amorphisized the crystalline diamond structure with 64 silicon atoms in the cell (a crystalline density of 2.33 g/cm3) by slowly heating it from room temperature to well above the glass transition temperature and just below the melting point, and then slowly cooling it to 0 K. This then was followed by cycles of annealing and quenching at temperatures suggested by experiment, to let the structure obtained adjust to the local energy minimum. Heating, cooling, annealing and quenching cycles were carried out using either the 2.44 fs or the 10 fs time step throughout.

To treat a-Si:H we first constructed two amorphous silicon structures as follows: One was created as described above, with a time step of 2.44 fs, maintaining the crystalline density, and then expanded to reproduce the experimental density, 2.2 g/cm3 of the hydrogenated structure, once hydrogens are incorporated. The other was created using a time step of 10 fs from an expanded crystalline cell so that the density of the hydrogenated cell (2.2 g/cm3) corresponds to the reported experimental values. Once these expanded pure amorphous silicon cells were constructed, we introduced 12 hydrogen atoms at the following relative cell positions: (1/4, 1/4, 1/4), (3/4, 1/4, 1/4), (3/4, 3/4, 1/4), (1/4, 3/4, 1/4), (1/2, 1/4, 1/2), (3/4, 1/2, 1/2), (1/2, 3/4, 1/2), (1/4, 1/2, 1/2), (1/4, 1/4, 3/4), (3/4, 1/4, 3/4), (3/4, 3/4, 3/4), (1/4, 3/4, 3/4), and annealed them at temperatures suggested by experiment.

Specifically, for a-Si the cell was heated from 300 K to 1,680 K, just below the melting temperature of crystalline silicon, in 100 steps of 10 fs, and immediately cooled down to 0 K in 122 steps of 10 fs; this was done also for a time step of 2.44 fs. The heating/cooling rate was $`5.66\times 10^{15}`$K/s for 2.4 and $`1.38\times 10^{15}`$K/s for 10 fs, approximately. The atoms were allowed to move freely within each cell of volume (10.8614 Å)3 with periodic boundary conditions. Once this first stage was complete, we subjected each cell to annealing cycles at 300 K (below microcrystallization ) with intermediate quenching processes.

For a-Si:H we implemented two different procedures. First, we used the amorphous pure silicon cell generated above with a time step of 2.44 fs and then expanded it to a volume of (11.0620 Å)3. Second, we amorphisized a previously expanded crystalline cell of 64 silicon atoms with the same volume of (11.0620 Å)3 using a 10 fs time step and repeating the cycles described above. We then placed the 12 hydrogens distributed throughout the amorphous cells and subjected the system to annealing cycles using a time step of 0.46 fs for the first cell (the 2.44/0.46 cell), and 2 fs for the second one (the 10/2 cell). The annealing cycles consist of two cycles of 50 steps at 300 K for the large time step sample and one cycle of 200 steps for the small time step sample, with in between quenches down to 0 K, to allow the hydrogens to diffuse and move in the cells. This gives a concentration of hydrogen of practically 16%, adequate to compare with existing experimental results.

Once the atomic structures were constructed and their respective RDFs obtained, we analysed their electronic density of states at the $`\mathrm{\Gamma }`$ point of the Brillouin Zone using both FastStructure and the DFT approach of Hohenberg and Kohn and Kohn and Sham implemented in the ab initio DMol3 comercial code; both codes were used to calculate the energy levels and DOS curves of the amorphous atomic structures generated. DMol3 treats the electronic structure of periodic solids by solving the Kohn-Sham self-consistency equations within local or nonlocal density approximations; in our work we carried out just energy calculations of the structures found with FastStructure and used the LDA approximation. We also used a double numerical basis set that includes d-polarization of the atoms (DNP) and the frozen inner core orbitals approximation; a medium grid was used for the numerical calculations. The SCF density parameter that specifies the degree of convergence for the LDA density was set at $`10^6`$. The HOMOs, LUMOs and DOS curves were calculated since they are necessary to study the behaviour of the forbidden energy gaps, the position of the Fermi levels, and the gap levels introduced by either dangling bonds, floating bonds or hydrogen states.

## IV Results and discussion

Defining a dangling or a floating bond is to some extent arbitrary since one has to choose some interatomic distance; Drabold et al. chose an interatomic distance of 2.7 Å, although no clear justification for this choice is given. Lee and Chang use the distance at which the minimum of the RDF occurs, 2.73 Å, which leads to an average number of nearest neighbors in the amorphous cell of 3.9. Were we to use the minimum value of the RDF of Fedders et al., $`1.2\times 2.35=2.82`$ Å, it is certain that the number of defects reported by Drabold and coworkers would have changed significantly. An analysis of structural defects versus spectral defects was performed by these authors and it was concluded that dangling bonds give rise to defects within the gap, that some strained tetrahedral structures may also generate gap states, but that it seems floating bonds do not create gap states. Recent work, however, indicates that floating bonds may generate gap states, and this conclusion is also borne out in the results reported here.

Generating a clean gap in such small cells is a difficult task since defects always appear in the simulations; so it seems that the most one can hope for is to minimize the number of defects and so obtain a reasonably clean electronic gap. On the other hand, generating RDFs, total and partial, that compare favorably with experiment is very important as is reproducing the hydrogen dynamics and its experimental behaviour of passivating the dangling bonds in amorphous silicon. This ’cleans’ the gap, in spite of the limited size of the cells. Herein we address these issues and report total RDFs for pure and hydrogenated amorphous silicon, and partial RDFs for silicon-silicon, silicon-hydrogen and hydrogen-hydrogen in the hydrogenated samples and also report HOMOs, LUMOs and their DOS curves.

### A The atomic structure

For a-Si the RDF obtained from our simulation with 10 fs agrees very well with experiment from 0 to 10 Å, Fig. 1. In this work we fourier-smoothed all the RDF results to have adequate curves to allow comparison with experiment, since the small number of atoms in the cell leads to statistical fluctuations that are not representative of the bulk. In Fig. 1 we do a direct comparison of our results and the upper and lower experimental bounds from Refs. 10 to 16 and, considering that not all experimental results cover the range 0 to 10 Å, the agreement is excellent, including the existence and position of the third and the fourth amorphous peaks. The evolution of the crystalline peaks can be clearly observed in our simulations, Fig. 2, and we can say that the first crystalline peak remains as the most prominent feature in the amorphous material. The second crystalline peak also remains but is widened and highly diminished since it contributes to filling in the first and second crystalline valleys. The third crystalline peak disappears to contribute to the valleys between the second and third and the third and fourth crystalline peaks, and the fourth crystalline peak disappears to contribute to the filling in the valleys between the third and fifth crystalline peaks. The fifth crystalline peak generates a third amorphous peak slightly displaced to smaller distances. This figure presents a comparison of the position of the crystalline peaks and the amorphous structure where all this can be better appreciated.

The total and partial RDF simulations for a-Si:H also compare very favorably with the experimental results given in Ref. 18, where a more complete and more recent study is reported. In Figs. 3 and 4 we show the comparison of the total RDFs and the partial H-H RDFs with experiment for the two cells, 2.44/0.46 (Fig. 3) and 10/2 (Fig. 4), where the peak due to the presence of molecular hydrogen that appears in the simulation of the 10/2 cell has been removed in order to better illustrate the agreement of our results with experiment. In particular, the H-H correlation function that we obtained for the 2.44/0.46 fs cell is close to the experimental results of Bellissent et al. unlike previous ab initio simulations. As we shall see later, it turns out that this sample contains more dangling and floating bonds than the 10/2 fs cell. Since the number of hydrogen atoms in the cells is smaller than the number of silicon atoms, a larger fourier smoothing for the H-H partial was used in order to compare our results with experiment. In Fig. 5 we are able to see the dynamics of the originally evenly distributed hydrogen (see Fig. 6 for the distribution of hydrogens) since for the 10/2 fs cell a peak appears in the total and H-H RDFs due to the formation of molecular hydrogen at an interatomic distance of 0.875 Å, (see Ref. 49 where a molecular radius of 0.86 Å for hydrogen in a crystalline silicon cluster is reported). No formation of molecular hydrogen shows up in the 2.44/0.46 fs cell since the large number of dangling bonds inhibits it. Figs. 7 to 10 depict four structures that clearly indicate the passivation process of the existing dangling bonds in the pure amorphous silicon samples due to hydrogenation: Fig. 7(9) shows 5(2) dangling bonds that existed in the pure amorphous sample for the 2.44(10) fs cell, for a cutoff radius of 2.815(2.743) Å below which the number of neighbors is 4(4); Fig. 8(10) shows the 2.44/0.46(10/2) fs cells indicating that 3(1) bonds were passivated by hydrogens, 2(1) by silicons, and 3(2) new ones appeared; here the same cutoff radius mentioned above for Si-Si was used, whereas the cutoff radius for Si-H was 1.9(1.9) Å which is the position of the valley minimum between the first and second peaks in the Si-H partial RDF (See Figs. 11 and 12). The Si-Si and Si-H partial RDF also agree well with experiment, as indicted in Figs. 11 and 12 where a direct comparison is shown.

The 2.44/0.46 cell shows a large number of floating bonds, 11, after the hydrogen diffusion process, as opposed to 3 floating bonds (fb) for the 10/2 fs cell. Of the 11 fb, 6 are due to hydrogen atoms that bond to silicons in addition to the existing tetrahedral coordination and it is for this cell that the addition of H decreases the size of the electronic gap, as will be shown next, indicating the important role of this type of bonds in reducing the gap.

### B The electronic density of states

Using ab initio methods in order to study electronic characteristics of amorphous materials invariably forces one to consider cells of some hundred atoms at most, due to the limitations in computing power of present day resources. Therefore it is desirable to have practically no defects in these cells since then they would better represent what occurs in the best extended samples. However, as stated before, simulations in cells of 64 to 216 atoms usually have a large concentration of dangling and floating bonds. That is why we have developed a different approach to generate amorphous cells that seem to give a smaller concentration of these defects and that are therefore more in agreement with good bulk samples.

As mentioned above for a-Si our two amorphous cells 2.44/10 have 5 and 3 dangling bonds and 5 and 1 floating bonds for a cutoff radius of 2.763/2.772 Å. At these distances the number of nearest neighbors (area under the first amorphous peak) is equal to 4. Calculations of the DOS curves, and HOMOs and LUMOs with both FastStructure and DMol3 give clean gaps of 0.414 eV and 0.173 eV respectively for the 2.44 fs cell and 0.744 eV and 0.385 eV respectively for the 10 fs cell . The gaps were obtained just by looking at the HOMO and LUMO values; no attempt was made at sorting out the states near the gap due to dangling and/or floating bonds. Figs. 13 and 14 show the DOS curves where the gaps are clearly indicated. It should be kept in mind that density functional calculations using the local density approximation tend to underestimate band gaps.

Two procedures were carried out in order to generate the amorphous cells of hydrogenated silicon of the correct density; both imply an expansion of the cell either before the amorphisation or after, as described above. When the cell is expanded after the amorphisation a larger number of defects is created, 6 dangling bonds and 2 floating bonds, for the cutoff radius of 2.763 Å. Then 12 hydrogens were placed as described above and indicated in Fig. 6 and the cell was annealed as described above. It should be emphasized that in this method the hydrogen atoms were left free to diffuse, unlike other results reported in the literature where the hydrogens are placed in ad hoc positions to saturate the dangling bonds. Figs. 15 and 16 show the DOS curves obtained with FastStructure and DMol3 for hydrogenated amorphous silicon, the gap values are 0.324 eV and 0.142 eV for the 2.44/0.46 fs cell and 0.787 eV and 0.483 eV for the 10/2 fs cell, respectively. The result obtained with the 10/2 fs cell follows the trend found in experiments that the addition of hydrogen increases the gap value with respect to pure silicon. The result for the 2.44/0.46 fs cell indicate that the size of the gap can be diminished by the presence of the 6 hydrogen-related and 5 silicon-related floating bonds; this supports the idea expressed recently by Fornary and coworkers.

## V Conclusions

Dealing theoretically with a covalent amorphous material in bulk is difficult because there are many possible structures for a given RDF. Continuous random-network models have achieved considerable success in generating a-Si structures but these models are not easily generalizable to other glassy materials or systems with different species of atoms. That is why the simulation of covalent amorphous semiconductors in the bulk, both pure and alloyed, has encountered serious obstacles. In the present work, we have been able to perform ab initio simulations of the atomic structure and electronic features of amorphous pure and hydrogenated silicon and simulation of their total and partial atomic RDFs that reproduce many features shown by the experimental data. In particular, our RDFs agree very well with the corresponding experimental distributions up to the third and fourth peaks, as can be seen in the direct comparisons made in the figures. In addition, the dynamics of atomic hydrogen, including the formation of molecular hydrogen under certain conditions, and the passivation of dangling bonds was observed. The electronic density of states of the atomic distributions generated indicate the presence of a gap for the pure amorphous silicon samples that widens when the sample is annealed after hydrogen is distributed within the cell and the number of defects (dangling and floating bonds) is kept low.

The procedure that we used to generate these structures is different from the ones found in the literature and leads to amorphous samples with fewer dangling bonds and fewer overcoordinated defects. It consists of heating a crystalline sample of 64 atoms of silicon just below its melting temperature and then cooling it down to 0 K with subsequent annealing and quenching cycles at temperatures dictated by experiment.

We believe that our procedure can be generalized to other monatomic or diatomic amorphous semiconducting materials and therefore allows ab initio techniques to be used in other areas to provide more representative structures and a more complete description of their atomic and electronic characteristics than is possible with parameterized methods. Furthermore, ab initio simulations can be used to provide the parameters needed in non ab initio techniques to be able to handle larger amorphous cells.

ACKNOWLEDGMENT

AAV thanks DGAPA-UNAM for financial support to spend a sabbatical year at Molecular Simulations, Inc. where this work was begun and for financial support to carry out this research through project IN101798 (Estructura Electrónica y Topología Atómica de Silicio Amorfo Puro y Contaminado). FA thanks CONACyT for supporting his PhD studies. This work was done on an Origin 2000 computer provided by DGSCA, UNAM.

FIGURE LIST

FIG. 1. RDFs for a-Si. The lighter lines are the experimental upper and lower bounds (see text). The dark line is our simulation where the third and fourth peaks are well reproduced.

FIG. 2. Evolution of the crystalline peaks during amorphisation. The RDF for a-Si (dark line) has four well defined peaks that are generated by several crystalline peaks (lighter vertical lines).

FIG. 3. Direct comparison of the simulated (dark lines) total RDF for a-Si:H and partial RDF for H-H (Inset) (with 15.79% hydrogen concentration) to the experimental (light lines) results for the short time steps cell (2.44/0.46 fs).

FIG. 4. Direct comparison of the simulated (dark lines) total RDF for a-Si:H and partial RDF for H-H (Inset) (with 15.79% hydrogen concentration) to the experimental (light lines) results for the long time steps cell (10/2 fs). The peak due to molecular hydrogen has been removed for clarity.

FIG. 5. Total RDF for the 10/2 fs cell of hydrogenated amorphous silicon (dark line) that shows the presence of molecular hydrogen, compared to the experimental results (light line).

FIG. 6. Starting positions of the 12 hydrogens placed within the amorphous cells of pure silicon. The relative coordinates are given in the text.

FIG. 7. Short time step cell of pure amorphous silicon (2.44 fs) that shows the existence of 5 dangling bonds. A cutoff radius of 2.815 Å was used, below which the total number of nearest neighbors is 4.

FIG. 8. When hydrogen is added to the short time steps cell, and the sample annealed, 3 dangling bonds are passivated by hydrogens and 2 by silicons but 3 new ones appear; however, 11 new floating bonds are formed (not shown).

FIG. 9. Long time step cell of pure amorphous silicon (10 fs) that shows the existence of 2 dangling bonds. A cutoff radius of 2.743 Å was used, below which the total number of nearest neighbors is 4.

FIG. 10. When hydrogen is added to the long time steps cell, and the sample annealed, all dangling bonds (2) are passivated, one by hydrogen and one by silicon, and 2 new ones appear. The number of floating bonds now is only 3 (not shown). The appearance of molecular and atomic hydrogen is indicated.

FIG. 11. Partial Si-Si and Si-H RDFs for the hydrogenated short time steps cell. The simulated results are represented by dark lines whereas the experimental ones are represented with lighter lines.

FIG. 12. Partial Si-Si and Si-H RDFs for the hydrogenated long time steps cell. The simulated results are represented by dark lines whereas the experimental ones are represented with lighter lines.

FIG. 13. DOS curves for the 2.44 fs cell of a-Si calculated using FastStructure, curve (a), and DMol3, curve (b). HOMOs and LUMOs are shown that lead to a gap of 0.414 eV for (a) and 0.173 eV for (b).

FIG. 14. DOS curves for the 10 fs cell of a-Si calculated using FastStructure, curve (a), and DMol3, curve (b). HOMOs and LUMOs are shown that lead to a gap of 0.744 eV for (a) and 0.385 eV for (b).

FIG. 15. DOS curves for the 2.44/0.46 fs cell of a-Si:H calculated using FastStructure, curve (a), and DMol3, curve (b). HOMOs and LUMOs are indicated and the gap values are 0.324 eV and 0.142 eV, respectively.

FIG. 16. DOS curves for the 10/2 fs cell of a-Si:H calculated using FastStructure, curve (a), and DMol3, curve (b). HOMOs and LUMOs are indicated and the gap values are 0.787 eV and 0.483 eV, respectively.


text:
# The Probability Density of the Mass of the Standard Model Higgs Boson

## Abstract

The LEP Collaborations have reported a small excess of events in their combined Higgs boson analysis at center of mass energies $`\sqrt{s}\stackrel{<}{_{}}208`$ GeV. In this communication, I present the result of a calculation of the probability distribution function of the Higgs boson mass which can be rigorously obtained if the validity of the Standard Model is assumed. It arises from the combination of the most recent set of precision electroweak data and the current results of the Higgs searches at LEP 2.

preprint:  UPR–T–XXX

Combining all Higgs decay channels and experiments, the LEP Collaborations report a $`2.5\sigma `$ excess in their data . The probability that this is due to an upward fluctuation of the background is 0.6%. Of course, the Higgs boson has been searched for in many different energy bins, and there is an infinitely large energy range out of reach, so that one expects to observe an upward fluctuation somewhere. It is therefore difficult to interpret these numbers, and it would be imprudent to conclude that the Higgs boson has been found with 99.4% probability.

In this communication, I present the answer to a different but related question, which is Given the data, what is the probability that the Higgs boson is within reach of LEP 2? This question can be answered unambiguously once the probability distribution of the Higgs boson mass, $`M_H`$, has been constructed. This is not possible given the Higgs search results at LEP 2 by themselves, regardless of how strong a signal is observed there. The reason is that there is an infinite domain of a priori possible values of $`M_H`$ beyond the kinematic reach of LEP 2. As a result, the $`M_H`$ distribution is improper, i.e., it is asymptotically non-zero. Including the electroweak precision data, however, renders it sufficiently convergent and a proper integration is possible.

The electroweak precision data by itself, especially $`Z`$ pole asymmetries and the $`W`$ boson mass, is now precise enough to constrain $`M_H`$ significantly, as can be seen from Fig. 1.

Including the latest updates as presented at the 2000 summer conferences, I find from a global fit to all data using the package GAPP ,

$$M_H=86_{32}^{+48}\text{ GeV}.$$
(1)
Note, that by definition the central value in Eq. (1) maximizes the likelihood, $`Ne^{\chi ^2(M_H)/2}`$, and that correlations with other parameters, $`\xi ^i`$, are accounted for, since minimization w.r.t. these is understood, $`\chi ^2\chi _{\mathrm{min}}^2`$. The error is the standard $`1\sigma `$ uncertainty ($`\mathrm{\Delta }\chi ^2=1`$). The 68.27% central confidence interval,

$$53\text{ GeV}M_H131\text{ GeV},$$
(2)
differs slightly from Eq. (1) due to the non-Gaussian (asymmetric) distribution of $`\mathrm{ln}M_H`$. The 90% central confidence interval yields,

$$38\text{ GeV}M_H173\text{ GeV}.$$
(3)
The right hand side of Eq. (3), i.e., the 95% upper limit, does not take into account the direct search results of the LEP Collaborations. Negative search results will increase the upper limit, because the probability distribution function is effectively being renormalized. For example, in a previous analysis  we found that the Higgs exclusion curve presented by the LEP Collaborations increased the 95% upper limit by 30 GeV.

The use of the Higgs exclusion curve, however, is only appropriate if no indication of an excess is observed. In general, it is more appropriate to consider the likelihood ratio for the data,

$$Q_{\mathrm{LEP}2}=\frac{(\mathrm{data}|\mathrm{signal}+\mathrm{background})}{(\mathrm{data}|\mathrm{background})},$$
(4)
where both the numerator and the denominator are functions of $`M_H`$. The quantity,

$$\mathrm{\Delta }\chi ^2(\mathrm{direct})2\mathrm{ln}Q_{\mathrm{LEP}2}$$
(5)
can then be added to the $`\chi ^2`$-function obtained from the precision data. If the signal hypothesis gives a better (worse) description of the data than the background only hypothesis we find a negative (positive) contribution to the total $`\chi ^2`$. Note, that this is a consistent treatment also in the case of a large downward fluctuation of the background or even if no events are observed at all. Use of $`Q_{\mathrm{LEP}2}`$ had been originally advocated in Ref.  (see Eq. (23) in that reference).

This treatment can be rigorously justified within the framework of Bayesian statistics , which is particularly suited for parameter estimation. Bayesian methods are based on Bayes theorem ,

$$p(M_H|\mathrm{data})=\frac{p(\mathrm{data}|M_H)p(M_H)}{p(\mathrm{data})},$$
(6)
which must be satisfied once the likelihood, $`p(\mathrm{data}|M_H)`$, and prior distribution, $`p(M_H)`$, are specified. $`p(data)p(\mathrm{data}|M_H)p(M_H)𝑑M_H`$ in the denominator provides the proper normalization of the posterior distribution on the left hand side. Depending on the case at hand, the prior can

1. contain additional information not included in the likelihood model,
2. contain likelihood functions obtained from previous measurements,
3. or be chosen non-informative.

Of course, the posterior does not depend on how information is separated into the likelihood and the prior. As for the present case, I choose the informative prior,

$$p(M_H)=Q_{\mathrm{LEP}2}p^{\mathrm{non}\mathrm{inf}}(M_H),$$
(7)
where the non-informative part of the prior will be chosen as

$$p^{\mathrm{non}\mathrm{inf}}(M_H)=M_H^1.$$
(8)
This choice corresponds to a flat prior in the variable $`\mathrm{ln}M_H`$, and there are various ways to justify it . One rationale is that a flat distribution is most natural for a variable defined over all the real numbers. This is the case for $`\mathrm{ln}M_H`$ but not $`M_H^2`$. Also, it seems that a priori it is equally likely that $`M_H`$ lies, say, between 30 and 40 GeV, or between 300 and 400 GeV. In any case, the sensitivity of the posterior to the (non-informative) prior diminishes rapidly with the inclusion of more data. As discussed before, $`p(M_H)`$ is an improper prior but the likelihood constructed from the precision measurements will provide a proper posterior.

Occasionally, the Bayesian method is criticized for the need of a prior, which would introduce unnecessary subjectivity into the analysis. Indeed, care and good judgement is needed, but the same is true for the likelihood model, which has to be specified in any statistical model. Moreover, it is appreciated among Bayesian practitioners, that the explicit presence of the prior can be advantageous: it manifests model assumptions and allows for sensitivity checks. From the theorem (6) it is also clear that any other method must correspond, mathematically, to specific choices for the prior. Thus, Bayesian methods are more general and differ rather in attitude: by their strong emphasis on the entire posterior distribution and by their first principles setup.

Including $`Q_{\mathrm{LEP}2}`$ in this way, one obtains the 95% CL upper limit $`M_H201`$ GeV, i.e. notwithstanding the observed excess events, the information provided by the Higgs searches at LEP 2 increase the upper limit by 28 GeV.

Given extra parameters, $`\xi ^i`$, the distribution function of $`M_H`$ is defined as the marginal distribution, $`p(M_H|\mathrm{data})=p(M_H,\xi ^i|\mathrm{data})_ip(\xi ^i)d\xi ^i`$. If the likelihood factorizes, $`p(M_H,\xi ^i)=p(M_H)p(\xi ^i)`$, the $`\xi ^i`$ dependence can be ignored. If not, but $`p(\xi ^i|M_H)`$ is (approximately) multivariate normal, then

$$\chi ^2(M_H,\xi ^i)=\chi _{\mathrm{min}}^2(M_H)+$$

$$\frac{1}{2}\frac{^2\chi ^2(M_H)}{\xi _i\xi _j}(\xi ^i\xi _{\mathrm{min}}^i(M_H))(\xi ^j\xi _{\mathrm{min}}^j(M_H)).$$

The latter applies to our case, where $`\xi ^i=(m_t,\alpha _s,\alpha (M_Z))`$. Integration yields,

$$p(M_H|\mathrm{data})\sqrt{detE}e^{\chi _{\mathrm{min}}^2(M_H)/2},$$
(9)
where the $`\xi ^i`$ error matrix, $`E=(\frac{^2\chi ^2(M_H)}{\xi _i\xi _j})^1`$, introduces a correction factor with a mild $`M_H`$ dependence. It corresponds to a shift relative to the standard likelihood model, $`\chi ^2(M_H)=\chi _{\mathrm{min}}^2(M_H)+\mathrm{\Delta }\chi ^2(M_H)`$, where

$$\mathrm{\Delta }\chi ^2(M_H)\mathrm{ln}\frac{detE(M_H)}{detE(M_Z)}.$$
(10)
This effect tightens the $`M_H`$ upper limit by 1 GeV.

I also include theory uncertainties from uncalculated higher orders. This increases the upper limit by 5 GeV,

$$M_H205\text{ GeV}\text{ (95\% CL)}.$$
(11)
The entire probability distribution is shown in Fig. 2. Taking the data at face value, there is (as expected) a significant peak around $`M_H=115`$ GeV, but more than half of the probability is for Higgs boson masses above the kinematic reach of LEP 2 (the median is at $`M_H=119`$ GeV). However, if one would double the integrated luminosity and assume that the results would be similar to the present ones, one would find most of the probability concentrated around the peak. A similar statement will apply to Run II of the Tevatron at a time when about 3 to 5 $`\mathrm{fb}^1`$ of data have been collected.

The described method is robust within the SM, but it should be cautioned that $`M_H`$ extracted from the precision data is strongly correlated with certain new physics parameters. Likewise, the Higgs searches at LEP 2 depend on the predictions of signal and background expectations which are strictly calculable only within a specified theory. This note focussed on the Standard Model Higgs boson.

Acknowledgements:

This work was supported in part by the US Department of Energy grant EY–76–02–3071.


text:
# WIND INTERACTION MODELS FOR THE AFTERGLOWS OF GRB 991208 and GRB 000301C

## 1 INTRODUCTION

The simplest model of the afterglows of gamma-ray bursts (GRBs) involves a relativistic spherical blast wave running into a constant density, presumably interstellar, medium (Mészáros & Rees 1997; Katz 1994; see Piran 1999 for a review). The afterglows are emitted by non-thermal electrons whose energy distribution is usually assumed to be a power-law above some cutoff. The predicted power-law decay of the afterglow emission with time has been subsequently observed at X-ray (Costa et al. 1997), optical (van Paradijs et al. 1997), and radio (Frail et al. 1997) wavelengths, giving basic confirmation to this now “standard” picture.

There are observed features of afterglows that are difficult to accommodate by the simplest model. One of the most noticeable is the steepening of the optical light curves. Good examples include GRBs 990123 (Kulkarni et al. 1999), 990510 (Harrison et al. 1999; Stanek et al. 1999), and 991216 (Halpern et al. 2000). The steepening is usually attributed to a jet-like, instead of spherical, blast wave (Rhoads 1997; Sari, Piran & Halpern 1999).

Radio data are especially useful for analyzing the afterglows of GRBs because the self-absorption frequency, $`\nu _a`$, and the characteristic frequency of the lowest energy electrons, $`\nu _m`$, both typically lie in this range. GRB 970508 was followed extensively at radio wavelengths (Frail, Waxman, & Kulkarni 2000) and Frail et al. (2000) have suggested that there was an initial spherical relativistic expansion in a constant density medium (up to day 25), followed by lateral jet expansion (days 25–100) and spherical non-relativistic expansion (days 100–400). The same data were modeled by Chevalier & Li (2000) as a spherical blast wave running into an $`r^2`$ density medium characteristic of a constant mass loss rate and velocity circumstellar wind, possibly of Wolf-Rayet star origin. The model missed an early low frequency radio data point and some infrared points, but it captured most of the observed behavior over 400 days. The question of a wind vs. a constant density surrounding medium is crucial for the question of the progenitors of GRBs, because massive stars, one of the leading candidates for GRB progenitors (Paczyński 1998), should be surrounded by a wind.

Fairly extensive radio observations are now available for GRBs 991208 and 000301C (Galama et al. 2000; Berger et al. 2000). At optical wavelengths, the light curves of both sources are somewhat unusual: those of GRB 991208 are steeper than $`t^2`$ (Sagar et al. 2000a), and those of GRB 000301C show pronounced steepening (Rhoads & Fruchter 2000; Masetti et al. 2000; Sagar et al. 2000b; Jensen et al. 2000). These two sources have been modeled as a jet running into a constant density medium by Galama et al. (2000) and Berger et al. (2000), based on both radio and optical data. The authors ruled out the simplest, spherical wind model with a standard, power-law electron energy distribution. Here, we show that the spherical wind model fits the observed radio data as well as the constant density jet model. We also demonstrate that a non-standard distribution of electron energy can account for both the steep decline and the steepening of the optical light curves of, respectively, GRBs 991208 and 000301C. In §2, we review the spherical wind model and describe the non-standard electron energy distribution. The model is then applied to the GRBs 991208 and 000301C in §3. We discuss our results in §4.

## 2 WIND MODEL WITH NON-STANDARD ELECTRON ENERGY DISTRIBUTIONS

Our spherical wind interaction model was originally developed to explain the radio observations of GRB 980425/SN 1998bw (Li & Chevalier 1999), and has been subsequently applied to GRBs 980519 (Chevalier & Li 1999) and 970508 (Chevalier & Li 2000). It involves a (trans-)relativistic blast wave propagating into a circumstellar wind. Previously, we have shown that the inferred wind properties are in the range of those expected from Wolf-Rayet stars. We assume that the post-shock material is distributed uniformly inside a thin shell, and approximately determine the blast wave dynamics using shock jump conditions, and particle and energy conservation. The dynamics agree with those from the more exact self-similar solutions of Blandford & McKee (1976) to within a factor of order unity. As usual, we consider synchrotron radiation as the emission mechanism of afterglows and assume that a constant fraction $`ϵ_e`$ ($`ϵ_B`$) of the total energy goes into the radiating electrons (magnetic field). The standard prescription for the energy distribution of electrons is a power-law above some minimum Lorentz factor, $`\gamma _{\mathrm{min}}`$, with a constant power-law index $`p`$. The model takes into account relativistic effects and self-absorption (important at the radio wavelengths), but not cooling. The cooling effects can be accounted for approximately once the cooling frequency $`\nu _c`$ is estimated (Sari, Piran & Narayan 1998; Chevalier & Li 2000). We adopt a flat universe with Hubble constant $`H_0=65`$ km s-1 Mpc-1 for cosmological corrections.

The spherical wind model with the standard, power-law distribution of electron energy works rather well for GRB 970508 (Chevalier & Li 2000). Applying it to GRBs 991208 and 000301C runs into a difficulty: we cannot fit the radio and optical data simultaneously; the optical light curves do not decline with time fast enough. To overcome this difficulty, we are motivated to seek a non-standard electron energy distribution that steepens at high energies. Through trial and error, we find that a wind model with the following broken power-law distribution can fit all of the data for both sources reasonably well:

$$\frac{dN_e}{d\gamma }=C_1\gamma ^{p_1},\mathrm{if}\gamma _{\mathrm{min}}<\gamma <\gamma _\mathrm{b},$$
(1)
$$=C_2\gamma ^{p_2},\mathrm{if}\gamma >\gamma _\mathrm{b},$$
(2)
where $`\gamma _\mathrm{b}`$ is the break Lorentz factor, $`p_1`$ and $`p_2`$ are the power-law indexes for electrons below and above the break $`\gamma _\mathrm{b}`$, and the coefficients $`C_1`$ and $`C_2`$ are related through $`C_2=C_1\gamma _\mathrm{b}^{p_2p_1}`$ to ensure continuity of the distribution. We further assume that the ratio of the break and minimum energies, $`R_\mathrm{b}\gamma _\mathrm{b}/\gamma _{\mathrm{min}}`$, remains constant in time so that the relative shape of the distribution is time invariant. Given $`p_1`$, $`p_2`$ and $`R_\mathrm{b}`$, the minimum Lorentz factor $`\gamma _{\mathrm{min}}`$ and the coefficient $`C_1`$ are determined from the total number and energy of the radiating electrons, which in turn are determined from the dynamic evolution of the GRB blast wave. More elaborate prescriptions of the electron energy distribution may produce better fits to the observed data, but would introduce more free parameters. The implications of the above distribution on particle acceleration will be discussed in §4.

## 3 OBSERVED SOURCES

### 3.1 GRB 991208

GRB 991208 was first detected with the Interplanetary Network (IPN) by the spacecrafts Ulysses, WIND and NEAR on December 8, 1999, at 04:36:52 UT (Hurley et al. 2000). The main properties of its optical and radio afterglows are discussed in depth by, respectively, Sagar et al. (2000a) and Galama et al. (2000). The optical afterglow has one of the most steeply declining light curves, with the usual temporal decay index $`\alpha 2.2`$ (Sagar et al. 2000a), assuming that the flux density evolves as $`F_\nu t^\alpha \nu ^\beta `$. The spectral index on December 16.68, 1999 is determined to be $`\beta =0.75\pm 0.03`$, based on the observed flux in the K-band (Bloom et al. 1999) and extrapolated fluxes in the R- and I-band (Sagar et al. 2000a), with negligible Galactic extinction. A day earlier on December 15.64, the fully calibrated Keck-II spectrum yields an index $`\beta =0.9\pm 0.15`$ between $`3850`$ Å and $`8850`$ Å (Djorgovski et al. 1999). The radio afterglow was observed extensively at a number of frequencies between 1.43 and 250 GHz over a period of two weeks (Galama et al. 2000). The multi-frequency data set allows an approximate determination by the authors of the evolution of several key synchrotron parameters: the absorption frequency $`\nu _at^{0.15\pm 0.23}`$, the “typical” frequency $`\nu _mt^{1.7\pm 0.7}`$, and the peak flux density $`F_mt^{0.47\pm 0.20}`$. The redshift of the source was determined to be $`z=0.706`$ (Dodonov et al. 1999; Djorgovski et al. 1999).

The steep decline of the optical light curves of GRB 991208 resembles those of GRB 980326 and especially GRB 980519 (Halpern et al. 1999; Jaunsen et al. 2000). For GRB 980519, we have developed a spherical wind model with a standard power-law electron energy distribution and found that it fits the available data, from radio to optical to X-ray, reasonably well (Chevalier & Li 1999). To account for the steep light curve decline of GRB 980519 (with $`\alpha =2.05\pm 0.04`$; Halpern et al. 1999), an electron energy power-law index of about $`p=3.0`$ is required. The value of $`p`$ is higher than that found in other GRB afterglows (see Table 1 of Chevalier & Li 2000) but is within the range found in radio supernovae. To explain the even steeper decline of the optical light curves of GRB 991208 (with $`\alpha =2.2\pm 0.1`$; Sagar et al. 2000a) using the standard spherical wind model, an even higher value of $`p3.3`$ would be required (Galama et al. 2000). Such a large value of $`p`$ would be in conflict with the value $`p=2.52`$ determined on December 15.5, 1999 (7.3 days after the burst) from the spectral distribution at several frequencies ranging from 1.43 GHz to R-band (Galama et al. 2000). Since the decline in the observed optical flux (with $`\alpha 2.2`$) is much steeper than that predicted by $`p=2.52`$ (which corresponds to $`\alpha =1.64`$ below the cooling frequency, assuming slow cooling), the value of $`p`$ should be smaller (larger) than $`2.52`$ before (after) day 7.3 to account for the optical flux relative to that in the radio. The apparent variation required in the value of $`p`$ motivates us to seek a non-standard distribution of electron energy with a non-constant power-law index $`p`$. One such non-standard distribution is the broken power-law distribution specified by equations (1) and (2) in §2.

We apply the spherical wind model with a broken power-law electron energy distribution, as described in §2, to the afterglow of GRB 991208. Through trial and error, we find that a “standard” model with the following combination of parameters fits the radio through R-band data reasonably well: the electron energy fraction $`ϵ_e=0.03`$, magnetic energy fraction $`ϵ_B=0.01`$, total blast wave energy $`E=2.7\times 10^{52}`$ ergs, wind mass loss rate $`\dot{M}=4.0\times 10^6`$ M yr-1 for a wind velocity of $`v_w=1,000`$ km s-1, electron energy distribution power-indexes $`p_1=2.0`$ and $`p_2=3.3`$, and the ratio of the break to minimum Lorentz factors $`R_\mathrm{b}=50`$. The model fits are shown in Fig. 1. The radio data are taken from Table 1 of Galama et al. (2000) and R-band data from Table 1 of Sagar et al. (2000a). Note that the above combination of parameters is not unique. Indeed, there is a family of parameter combinations, with $`ϵ_eϵ_B^{1/5}`$, $`Eϵ_B^{1/5}`$ and $`\dot{M}/v_wϵ_B^{2/5}`$, that have the same fitting curves, as long as $`ϵ_B`$ is small enough to keep the cooling frequency above the R-band at the beginning of the afterglow observations (Li & Chevalier 1999) and the electron energy distribution remains the same. Different electron energy distributions are also possible; the parameter $`p_2`$ that describes the distribution of the high energy electrons above the break is particularly ill-constrained by the light curves shown in Fig. 1 alone. It turns out that a “variant” model with a combination of the energy distribution parameters $`p_1=2`$, $`p_2=\mathrm{}`$ (i.e., no high energy electrons above the break at all), and $`R_\mathrm{b}=125`$ provides an equally good fit to the radio data and a slightly better fit to the R-band data.

Additional constraints come from the spectral index $`\beta `$ at the optical wavelengths. As mentioned earlier, a fully calibrated Keck-II spectrum yields $`\beta =0.9\pm 0.15`$ on December 15.64. At this time, the “standard” model with $`p_1=2`$, $`p_2=3.3`$, and $`R_\mathrm{b}=50`$ gives $`\beta =1.04`$, within 1$`\sigma `$ of the observed value, whereas the “variant” model with $`p_1=2`$, $`p_2=\mathrm{}`$ and $`R_\mathrm{b}=125`$ gives $`\beta =1.27`$, about 2$`\sigma `$ away from the observed value. It therefore appears that $`p_2`$ should be close to $`3.3`$. A potential difficulty is the value of $`\beta =0.75\pm 0.03`$ on December 16.68, inferred from the observed flux in the K-band (Bloom et al. 1999) and extrapolated fluxes in the R- and I-band (Sagar et al. 2000a). However, with only two usable data points in the I-band (Sagar et al. 2000a) and possible deviation from a pure power-law decay in the single K-band measurement (Bloom et al. 1999; see R-band data for reference), we regard the index inferred on December 16.68 as less reliable than that measured simultaneously over a wide range of frequencies a day earlier. Therefore, we believe that the overall fit to both the radio and optical data is acceptable, and suspect that the fit may be improved with more elaborate prescriptions of the electron energy distribution than the simple broken power-law adopted here.

Panel (a) of Fig. 1 should be compared with Fig. 1 of Galama et al. (2000), where they fit the same radio data with a spectral form given by Granot et al. (1999a,b), based on the Blandford-McKee self-similar solution of a spherical, ultra-relativistic blast wave propagating in a constant density medium, and the best-fit power law time evolution of $`\nu _a`$, $`\nu _m`$, and the peak flux $`F_m`$. Inspection by eye reveals that the two model fits are of comparable quality; both models fit the three highest frequency light curves rather well. The fits to the three lowest frequency light curves are less satisfactory in both models, presumably because the lower frequency data are affected more by interstellar scintillation (ISS; Goodman 1997; Galama et al., in preparation). It is partially based on the low frequency data, which fix the self-absorption frequency $`\nu _a`$, that Galama et al. (2000) ruled out the spherical wind model for this source. They inferred that $`\nu _a`$ evolves as $`t^{0.15\pm 0.23}`$, which is incompatible with the evolution, $`\nu _at^{3/5}`$, predicted by a wind model. The fitted slope, $`0.15\pm 0.23`$, hinges, however, to a large extent on the non-detection on December 22.96 at 1.43 GHz (see the lower-right panel of their Fig. 2), the lowest observing frequency that is most affected by ISS. The poor quality of the power-law fit to the inferred values of $`\nu _a`$ is reflected in the relatively large reduced chi-square: $`\chi _\mathrm{r}^2=3.5`$ for 2 degrees of freedom (see Table 3 of Galama et al. 2000). Indeed, if one were to ignore the value of $`\nu _a`$ based on the non-detection (the last data point in the leftmost panel of their Fig. 4), the remaining three points would yield a much steeper slope. Therefore, this piece of evidence against the wind model is not compelling in our opinion. As noted by Galama et al. (2000), the inferred evolution for two other key parameters, $`\nu _mt^{1.7\pm 0.7}`$ and $`F_mt^{0.47\pm 0.20}`$, is compatible with the wind model. A second objection raised by Galama et al. (2000) was that the spherical wind model requires an unusually large electron energy index, $`p3.3`$, to account for the steep decline of the optical light curves. This would indeed be the case if the electron energy distribution were a single power-law, as commonly assumed. The objection motivated us to seek a non-standard, broken power-law distribution of electron energy which, as demonstrated by Fig. 1, explains both the radio and R-band data reasonably well.

Based on the same set of optical and radio data, Galama et al. (2000) reached a different conclusion. They favored a constant density jet model for GRB 991208. As noted by these authors, the afterglow properties inferred from observations differ from the asymptotic model predictions (Sari, Piran & Halpern 1999) on two accounts. First, the inferred decrease of the peak flux density with time, $`F_mt^{0.47\pm 0.20}`$, is substantially slower than the predicted $`F_mt^1`$. Second, the inferred electron energy index $`p=2.52`$ on December 15.5 is larger than the predicted $`p2.2`$, based on the temporal decay index of $`\alpha 2.2`$ in R-band. The authors argued that a slow transition to the fully asymptotic regime (Kumar & Panaitescu 2000; Moderski, Sikora, & Bulik 2000) may resolve these discrepancies, since during the transition the rate of the peak flux decay should be between that of a spherical model, $`F_mt^0`$, and that predicted asymptotically, $`F_mt^1`$, and the decay index $`\alpha `$ of the R-band light curve should be greater than $`p`$ (or –2.52). While this explanation is plausible, it is far from being proven; detailed modeling of the transition is required.

The continued evolution of the radio flux provides a test of the models. Galama et al. (2000) note that in their model $`\nu _m`$ should pass 8.46 GHz at $``$12 days and the flux should then decay rapidly, $`F_\nu t^{2.22.5}`$. In our model, the flux evolution at this frequency should tend to $`F_\nu t^{1.25}`$ (see Fig. 1). Continued monitoring of the source was apparently undertaken (Galama et al. 2000).

### 3.2 GRB 000301C

GRB 000301C was first detected with the RXTE All Sky Monitor and the IPN spacecrafts Ulysses and NEAR on March 01, 2000, at 09:51:37 UT (Smith, Hurley & Cline 2000). The IR/optical properties of its afterglow are discussed in Rhoads & Fruchter (2000), Masetti et al. (2000), Sagar et al. (2000b), and Jensen et al. (2000). Two features of the IR/optical light curves stand out: (a) large amplitude, short time scale variations and (b) a relatively abrupt steepening from a temporal decay index of $`\alpha 1`$ to roughly $`2.7`$. The spectral index $`\beta `$ was determined in the IR/optical spectral region at several times between day 2 and 14, with values ranging from $`0.5`$ to $`1.5`$ (Rhoads & Fruchter 2000; Sagar et al. 2000b; Jensen et al. 2000). The redshift of the source was determined to be $`2.03`$ (Smette et al. 2000; Castro et al. 2000; Jensen et al. 2000). Radio light curves of varying degree of coverage are available at four frequencies (4.86, 8.46, 22.5 and 250 GHz; Berger et al. 2000), which allow for detailed modeling.

As pointed out by Berger et al. (2000), the steepening of the IR/optical light curves rules out the simplest spherical wind model with a standard, power-law distribution of electron energy for this source. A spherical wind model with a non-standard distribution of electron energy is still possible. Indeed, the model with a simple broken power-law distribution, outlined in §2, can fit the IR/optical light curve steepening, as well as the radio data reasonably well. This is demonstrated in Fig. 2, where we compare the observed radio and R-band light curves with those predicted from a wind model with $`ϵ_e=0.04`$, $`ϵ_B=0.01`$, $`E=2.3\times 10^{52}`$ ergs, $`\dot{M}/v_w=4.5\times 10^6`$ M yr-1/$`10^3`$ km s-1, $`p_1=2.2`$, $`p_2=15`$, and $`R_\mathrm{b}=140`$. Inspection of panel (b) reveals that the model fits rather well the overall steepening of the R-band light curve, but not the short time scale variations, most noticeably the bump around day 4 and, to a lesser extent, the bump around day 7. The first bump has been interpreted by Garnavich, Loeb & Stanek (2000) as due to microlensing, although it could also be caused by a local enhancement in the ambient density (Berger et al. 2000) or an impulsive energy input (Panaitescu, Mészáros & Rees 1998; Li & Chevalier 1999; Sari & Mészáros 2000; Dai & Lu 2000). As noted by Berger et al. (2000), the physical process responsible for the first bump on the R-band light curve may also explain the factor-of-two discrepancy between the 250 GHz data taken around the same time and the model predictions (see their Fig. 2 and panel \[a\] of our Fig. 2). Panel (a) of Fig. 2 should be compared with Fig. 2 of Berger et al. (2000), where the same radio data set is fitted with a constant density jet model. Judging by eye, we find the quality of the two model fits comparable: the jet model fits the 250 and 22.5 GHz data slightly better, whereas the wind model fits the 8.46 and 4.86 GHz data slightly better. We note again that the combination of parameters listed above is not unique. In particular, the power-law index $`p_2`$ for the high energy electrons above the break $`\gamma _\mathrm{b}`$ is not well constrained. A relatively large value of $`p_2`$ is required nevertheless, to explain the rapid decline of the R-band flux at late times (with a decay index of $`\alpha 2.7`$) in the context of the spherical wind model.

Berger et al. (2000) ruled out the possibility that the steepening in the light curves of GRB 000301C is due to a time-varying $`p`$, based on the fact that the single value, $`p=2.70`$, that they inferred from a global fit to the radio/IR/optical data set, appears to fit the spectral flux distributions simultaneously at both day 4.26 and 12.17 reasonably well. Since their fitting is based on a model that assumes a single power-law distribution for the electron energy (i.e., a constant value of $`p`$ at any given time), it does not necessarily rule out a curvature in the slope of the energy distribution (i.e., a $`p`$ that varies with energy at any given time) as the cause for the light curve steepening. Indeed, there is some evidence for the curvature in the spectral energy distributions (SEDs) in the near IR to near UV region compiled by Rhoads & Fruchter (2000; the left panel of their Fig. 3). These SEDs are reproduced in panel (a) of our Fig. 3, with the modest amount of Galactic extinction ($`E_{BV}=0.053`$; Schlegel, Finkbeiner & Davis 1998) already corrected for. Note, however, that the data at different frequencies are not taken simultaneously and considerable uncertainties are involved in interpolating their fluxes to a common time for this source with known large amplitude, short time scale variabilities. The relatively shallow distribution at day 7.59 could, for example, be due to the previously mentioned, second (smaller) bump on the R-band light curve around day 7, which may or may not be present in the K-band (where the data points are too sparse to tell). Nevertheless, there appears to be a general trend that the SED steepens with time (see also column 5 of Table 3 of Rhoads & Fruchter). We expect such a steepening to occur in our model, as electrons with higher and higher energy emit into the frequency band. The spectral steepening is a robust signature of the curvature in the electron energy distribution and could in principle be used to distinguish our spherical wind model from the jet model. In panel (b) of Fig. 3, we plot the predicted SEDs of the wind model of GRB 000301C shown in Fig. 2. As expected, the SEDs steepen with time, although they are generally shallower than their observed counterparts in panel (a). A better agreement is reached if we invoke a modest amount of host galaxy extinction of $`A_\mathrm{V}=0.09`$ of the Small Magellanic Cloud (SMC) type (Pei 1992), as recommended by Rhoads & Fruchter (2000). A somewhat different value of $`A_\mathrm{V}=0.14`$, also of the SMC-type, is favored by Jensen et al. (2000) on day 3.

An alternative model for the afterglow of GRB 000301C was proposed by Kumar & Panaitescu (2000b). These authors attribute the steepening of the R-band light curve to a sudden, large drop in the density of the ambient medium into which the GRB blast wave propagates. Initially, the light curve steepens continuously as the blast wave expands freely. It tends to a constant decay rate when the observed flux is dominated by the emission from the high latitude parts of the blast wave away from the line of sight to the explosion center. This model and ours are similar in that both invoke a sudden change in either the ambient density or the electron energy distribution and that the finite light-travel time effect plays an important role in shaping the observed light curves. A potential problem with the model of Kumar & Panaitescu is that it predicts a rather steep flux decay of $`t^{2.8}`$ for frequencies above $`\nu _\mathrm{m}`$ but below $`\nu _\mathrm{c}`$; the late-time decay of the observed radio flux at 8.46 GHz does not appear to be that steep (see Fig. 2), although additional observations at even later times are needed to draw a firmer conclusion.

Finally, we note that in our model even though the change in the power-law index $`p`$ of the electron energy distribution is discontinuous (jumping from $`p_1`$ to $`p_2`$), the variation of the spectral index $`\beta `$ with time is smooth (see Fig. 3b). This makes inferring the instantaneous value of $`p`$ from the value of $`\beta `$ at any given time, as has been done by, e.g., Sagar et al. (2000b), difficult.

## 4 DISCUSSION

Our study of these two GRBs can be separated into the radio evolution and the optical/IR evolution. ISM + jet models have been proposed for the radio emission from both sources, but we find that spherical wind interaction models can produce fits of comparable quality. The same is true for GRB 970508 (Frail et al. 2000; Chevalier & Li 2000). In the case of GRB 991208, the approximate power law behavior observed over a factor 10 in time is intermediate between that expected for a spherical explosion in a constant density medium and the asymptotic lateral expansion of a jet (Galama et al. 2000). Kumar & Panaitescu (2000a) found that the transition to the asymptotic jet evolution should take at least an order of magnitude, perhaps several orders of magnitude, in observer time. Yet the model of Berger et al. (2000) for GRB 000301C assumes a rapid transition to the asymptotic jet evolution. The uncertainties in the jet models give them additional parameters that are not present in the wind interaction models which assume spherical symmetry throughout the evolution. Kumar & Panaitescu (2000a) found that for jet expansion in a wind medium, any break in the observed light curves would be unlikely to be detected, even for a narrow jet with an opening angle of about a few degrees. The spherical assumption, which we make, should be adequate even if a jet is present.

The problem with the simplest wind model is that the optical light curves of both GRB 991208 and GRB 000301C drop more rapidly than predicted. We have solved the problem by invoking a steepening of the electron spectrum. We attribute the steepening to the acceleration process, although the details of how this occurs are not clear. Bednarz & Ostrowski (1998) have studied acceleration in relativistic shock waves and found that at high shock Lorentz factors, the spectrum tends to a power law with $`p2.2`$. However, at low Lorentz factors ($`\gamma 3`$), which are characteristic of afterglow shock waves, there are a number of possibilities, including $`p<2`$, depending on the shock parameters. A particle spectrum that is this flat must steepen at high energy in order to have finite total energy. The spectral steepening with $`p`$ increasing to a relatively large value of $`3.3`$ or more that we propose for GRBs 991208 and 000301C would be clearer if X-ray observations of the sources were available. Unfortunately, these two sources are not constrained by X-ray observations. Spectral steepening with $`p`$ increasing to a smaller value may also be possible. In such a case, the X-ray flux would be less affected.

The difference between an intrinsic spectral break and a laterally expanding jet is that, in the latter case, the break in a light curve should be achromatic. It appears to be difficult to decide this point solely on the basis of optical/IR observations because of the faintness of the sources and the sparse sampling. The difference between the radio and optical evolution might be clearer. The late radio light curve of GRB 000301C is most complete at 8.46 GHz and the wind model appears to give a better representation of these data than the jet model (compare Fig.  to Fig.  of Berger et al. 2000), although the uncertainties are too large to draw any firm conclusions. A similar test may be possible for GRB 991208 (see § 3.1). Another test provided by radio observations is at low frequencies. In the wind interaction model, self-absorption is high at early times so that the low frequency flux should be low. Both GRB 970508 and GRB 000301C have one early 1.4 GHz radio observation that is higher than expected in the wind interaction model (Chevalier & Li 2000; Fig. ); the 1.4 GHz observations of GRB 991208 show a decrease that is not expected in any model. More extensive early, low frequency observations would provide a useful test of models.

We thank J. Rhoads for providing us with the SED data of GRB 000301C shown in Fig. 3a. This work was supported in part by NASA grant NAG5-8232.