text:
# The Type Ia Supernova Width-Luminosity Relation

## 1 Introduction

Because they are among the brightest optical explosions in the universe, the promise of employing type Ia supernovæ (SNe Ia) as cosmological probes has long been recognized (Kowal, 1968). In the past few years, this promise has been realized by two groups (Schmidt et al., 1998; Perlmutter et al., 1999) who have used observations of SNe Ia to determine luminosity distances to supernovæ at redshifts beyond $`z=1`$. The results have been surprising: ours is a universe which is expanding at an accelerating rate. In several ways this result is in accord with our expectations. The observations imply that the age of the universe is near 14 Gy, comfortably older than the age of its contents. They also imply that the density of ordinary matter is well below the critical density, near $`\mathrm{\Omega }_{matter}=0.3`$, and this too is in agreement with most dynamical measurements Perlmutter et al. (1999). Perhaps the most exciting result, and certainly the most unnerving and controversial, is that the geometry of the universe is today dominated by an energy density which does not arise from ordinary matter but which is best fit by a vacuum energy density for which there are no satisfying explanations (Riess et al., 1998; Garnavich et al., 1998).

To be really useful in the classic cosmological test of the redshift/luminosity-distance relation, the luminosity of a standard candle must have a small dispersion. The set of “typical” SNe Ia used for cosmology exhibits a range of almost a factor of three in peak brightness. This would not be sufficient to enable the current sample, which extends to a redshift of $`z1`$, to place useful limits on cosmological acceleration; the observed departure from a flat, matter-dominated universe amounts to three tenths of a magnitude ($`30`$% in brightness). SNe Ia are useful as cosmological probes only because what variations are observed can, for the most part, be calibrated out, leading to a usefully small dispersion in corrected brightness of about 0.15 mag.

Phillips (1993) demonstrated from a set of nearby SNe Ia that the rate of decline from peak in the optical lightcurve is correlated with peak brightness. Using the larger and more uniform data set provided by the Calán-Tololo supernova survey, Hamuy et al. (1996d) showed that this correlation exists in the B-, V-, and I-band lightcurves (though, interestingly, with different slopes than those determined by Phillips 1993 from the nearest supernovae). They found that not only the decline rate but also the overall width of the lightcurve centered on maximum light is correlated with peak luminosity. Using this same sample, Riess et al. (1995) also showed that the lightcurve shape could be used as a predictor of luminosity. Perhaps most suggestive from a physical point of view, using a sample of high-redshift supernovæ (for which very early time observations are more easily obtained), Perlmutter (1997) have shown that not only is there a correlation of the lightcurve shape with luminosity but that, near peak, the lightcurve is the same in all supernovæ up to a correlated scaling in time and luminosity, although quantitative errors for this relation have yet to be published.

While this width-luminosity relation (henceforward WLR) is very well substantiated observationally, it is an empirical relation which has remained largely unexplained. Höflich et al. (1996) showed that many current explosion models come close to this relation and suggested that the relation was determined primarily by the precise nature of the explosion itself. In this paper we provide the rudiments of a physical explanation for the WLR which is based upon the physics of radiation flow in SNe Ia, one in which the details of the explosion itself are largely masked. If true, this argues strongly against any significant bias arising in the WLR due to evolutionary effects in the supernova progenitors.

SNe Ia are widely believed to result from the thermonuclear incineration of an accreting carbon/oxygen white dwarf. Beyond this, however, little is known with any degree of precision. Present uncertainties include the nature and evolution of the progenitor system, the mass of the dwarf at ignition, and the physics of the subsequent nuclear burning. Given this state of our knowledge, it would at first blush appear hopeless to explain the WLR from physical principles. It may be, however, that the calibrating relation arises naturally from conditions present in many or all such explosions of accreting dwarfs and that the details of progenitor evolution and burning physics are of little consequence to the WLR, though fascinating in their own right.

For want of any more certain starting point, we adopt this point of view and examine how far it may take us toward an understanding of the WLR for SNe Ia. In the first section, we reprise some of the physics of radiation transport in SNe Ia explored in two previous papers (Pinto & Eastman, 2000b, a, hereafter PE I & II). We then develop a schematic description of the lightcurve physics which illustrates the two main properties of the explosion which determine the WLR: the mass of radioactive 56Ni produced in the explosion and the rate of change of $`\gamma `$-ray escape. In the following sections we then present a sample set of lightcurve calculations which demonstrate the essential correctness of this view, and demonstrate that the observed WLR can be reproduced with this physics. We conclude with speculations on the probability that evolutionary effects can materially alter the relation observed in the local universe and lead to false cosmological conclusions.

## 2 Theory

The WLR is the observation that the greater the luminosity at maximum light, the greater the width of the light curve, or, equivalently, the slower the rate of decline from maximum. As discussed in PE I, the lightcurve of a SN Ia is determined by a competition between three effects: the deposition of energy from radioactive decay, the adiabatic conversion of internal energy to kinetic energy of expansion, and the escape of internal energy as the observed lightcurve.

Because the initial configuration is compact and at high density, the explosion is at first tremendously optically thick. The timescale associated with escape of internal energy is very long, and virtually all of the energy liberated in the explosion is converted to kinetic energy on a timescale approximately equal to the elapsed time. Were it not for the continuous energy input from decay of 56Ni and 56Co produced in the explosion, SNe Ia would be both faint and short-lived. Fortunately, 56Ni decay yields an energy of $`3.0\times 10^{16}\text{erg}\text{g}^1`$ and 56Co decay yields a further $`6.4\times 10^{16}\text{erg}\text{g}^1`$.

As the star expands, the optical depth decreases and the timescale over which this energy can escape declines. While trapped, radiant energy is constantly being converted into kinetic energy; the longer the radiation remains in the supernova, the less energy will escape to power the lightcurve. With a constant energy source, the luminosity would increase monotonically, asymptotically to balance the deposition at zero optical depth. In a real supernova, however, the deposition declines as radioactive material decays away and as the $`\gamma `$rays, in which the decay energy is emitted, escape. Peak luminosity thus occurs roughly when the product of the rising escape fraction and the decreasing deposition reaches a maximum. This is complicated somewhat by the fact that the balance is not instantaneous; the energy emerging as UV to IR photons at any time near peak light has been deposited over a significant fraction of the elapsed time. Arnett (1982) first demonstrated that at maximum light the instantaneous bolometric luminosity is approximately equal to the instantaneous rate of energy deposition by radiative decay. Earlier, the expansion timescale is short, the escape time is long, and the deposition rate is high. More internal energy goes into expansion than escapes, and a store of trapped energy is built up. Subsequent to maximum light, more energy escapes than is converted to kinetic energy, and the surviving store of trapped energy is released faster than the deposition.

The WLR is thus seen to be a relation between the maximum-light luminosity and the radiation escape time; the brighter the supernova, the longer its escape time. Radiation escapes from an optically thick medium by diffusion, and we can tentatively identify the escape time with the diffusion time – the average time it takes a photon to random walk its way to the surface, given by the number of scatterings times the time of flight between scatterings.

Following the usual random walk argument, the diffusion time can be written as

$$t_d=\frac{\alpha R\tau _{eff}}{c}$$
(1)
where $`R`$ is a characteristic size of the system (in this case, the radius of the ejecta), $`\tau _{eff}`$ is the effective optical depth, and $`\alpha `$ is a constant of order unity. By the argument above, if the expansion (elapsed) time and escape (diffusion) times are about equal at peak, we can write

$$t_d=t=\frac{\alpha R\tau _{eff}}{c}=\frac{\alpha Vt\tau _{eff}}{c}$$
(2)
which gives a diffusion optical depth at peak of

$$\tau _{eff}30\left(\frac{10^9\text{cm}\mathrm{s}^1}{V}\right)$$
(3)
(taking $`\alpha 1`$). This value is far smaller than the monochromatic optical depth to most trapped photons, which is typically nearer to 3000 for peak conditions in a SN Ia. The resolution to this seeming paradox, of course, is that the typical photons in an atmosphere do not carry the flux – radiation escapes where it can, at energies where the monochromatic optical depth is smaller.

At stellar densities, there is tight collisional coupling between the radiation field and the plasma, keeping the local photon distribution function at its blackbody limit. While the bulk of the energy density exists as UV photons, these photons are trapped by the high UV opacity. The bulk of the flux is transported at longer wavelengths on the Rayleigh-Jeans tail of the distribution, where the optical depth is typically less, “side-stepping” the higher UV opacity. The lower density of these photons is more than overcome by the decreased diffusion time they suffer, and the appropriate opacity to use in describing the flux is the familiar Rosseland mean.

For this picture to work, the relaxation rate of the trapped photon distribution toward a blackbody must be short compared with the rate at which energy escapes from the long-wavelength tail of the distribution. This relaxation is usually mediated by collisions, and its rate must therefore scale as the square of the density. If the optical depth to the majority of photons remains large (so that radiation remains trapped) while the density decreases, eventually collisions will not be able to maintain sufficient energy density in the long-wavelength tail of the photon distribution where the optical depth is less; the system will depart from local thermodynamic equilibrium (LTE).

PE II demonstrated that this has indeed become the case in SNe Ia by maximum light. The escape time is still mediated by the rate for converting energy from short-wavelength, UV photons into long wavelength, OIR photons. The difference from the more typical LTE situation arises from the nature of opacity in rapidly-expanding media and the character of atomic physics in iron-group ions.

In the homologous expansion which develops soon after explosion, the velocity gradient is isotropic. Much as in the expanding universe, there is a one-on-one mapping between path-length and redshift, and a photon traversing the ejecta will accumulate a red shift in the local gas rest frame which will bring it into resonance with line transitions. One can show (Eastman & Pinto, 1993) that under the conditions in SNe Ia, the time a photon spends scattering in a line is very short compared with the time it spends moving between lines. For the purpose of estimating the escape time, then, each line traversed thus acts as a single scattering, and the distance measured in units of a photon mean free path (the optical depth) will be just the number of optically thick lines traversed in accumulating a given Doppler shift; the corresponding effective opacity is known as an “expansion opacity” (cf. Eastman & Pinto, 1993; Wagoner et al., 1991; Karp et al., 1977; Castor et al., 1975).

Using simple analytic estimates, PE I showed that the central temperature in a $`\mathrm{M}_{\mathrm{ch}}`$ explosion at $`t<20`$ days is $`T13,000`$ K and the peak of the blackbody spectrum is in the UV near $`\lambda 2200`$ Å. For an iron-peak composition in the relevant range of temperature and density, the spectral density of optically thick lines in the UV is very large and is a strongly decreasing function of wavelength from the UV through the optical (cf. Figure 2 of PE II). As long as the Sobolev optical depth in these lines remains large, the diffusion optical depth will not change as the density decreases, allowing a departure from LTE not seen at depth in static atmospheres where the optical depth is linearly dependent upon density. The large spectral density of lines at UV wavelengths is far too great for there to be significant transport at UV energies; the time it takes a UV photon to random walk to the surface is long compared to the expansion time. This is substantiated by the observation of very little UV flux in maximum light SNe Ia in spite of their high temperatures.

If the Rayleigh-Jeans tail of the photon distribution (at wavelengths where significant transport can take place) could be re-populated by a mechanism more efficient at low densities than collisional excitation, significant flux could still be developed, leading to an escape time characterized by the required effective optical depth near 30.

When a photon is absorbed in a line, the possible outcomes of this interaction include: 1) re-emission in the same transition, followed possibly by re-absorptions and emissions – a cycle repeated many times until escape from the red wing of the line; 2) collisional destruction, transferring the photon energy by Coulomb interactions to a free electron; and 3) fluorescence, where the excited state populated by absorption de-populates by radiative decay to a state different from the original. In general, the electron density is too low in maximum light SNe Ia for collisional destruction to have much effect. Thus, a UV photon absorbed on the high energy (blue) side of a line will either be re-emitted on the red side, or fluorescence will occur, with the photon replaced by several longer-wavelength photons (usually), but occasionally by a shorter wavelength photon (less often). The efficiency for this conversion process obviously depends on the number and strength of decay channels which would result in the emission of an OIR photon. For the UV radiation trapped in SNe Ia, fluorescence is the most likely outcome. Fluorescence is much more important in SNe Ia than in other astrophysical objects because they are composed in large part of iron-group elements with enormously complex atomic structures and very high spectral densities of line transitions in which fluorescence can occur.

A key point is this: the number of OIR transitions is greatest in the neutral and low ionization states, and decreases with increasing ionization. Thus the UV$``$OIR fluorescence process is more efficient in Co II than in Co III, and more efficient in Co III than in Co IV. Consequently, the conversion process is sensitive to the ionization. In the approximate calculations reported below, we assume LTE ionization and excitation despite densities which are too low for this to be really true. We shall thus use temperature as a proxy for ionization, but the reader should keep in mind that this is an approximation, and it is actually ionization which is the controlling variable.

In light of the importance of fluorescence, we can sketch the escape of radiation from SNe Ia in the following fashion. Imagine a store of deposited energy, trapped in the form of UV photons. These UV photons are created by the absorption and thermalization of decay $`\gamma `$-rays and removed by conversion to kinetic energy and to OIR photons. For the energy density of UV photons, $`E_{UV}`$, we can write the following schematic rate equation:

$$\frac{E_{UV}}{t}=\rho \dot{S}_{dep}(t)\frac{E_{UV}}{t_e}r_fE_{UV}$$
(4)
where $`\rho \dot{S}_{dep}(t)`$ is the time-dependent deposition rate, $`E_{UV}/t_e`$ is the rate of conversion of radiation to kinetic energy, and $`r_fE_{UV}`$ is the rate for converting UV to optical and infrared energy. We have assumed here that the diffusion time suffered by the UV energy is infinitely long; transport in the UV is unimportant. Similarly, for the OIR photons one can write

$$\frac{E_{OIR}}{t}=\frac{E_{OIR}}{t_e}\frac{E_{OIR}}{t_d}+r_fE_{UV}$$
(5)
where $`t_d`$ is the typical OIR diffusion time. Here we have ignored the possibility that deposition results in the direct emission of OIR photons. If optical photons are processed on a short enough time scale – for instance, if $`t_dt_e`$ – then we can set $`E_{OIR}/t0`$ and solve for $`E_{OIR}`$:

$$E_{OIR}=\frac{r_f}{1/t_e+1/t_d+r_f}E_{UV}.$$
(6)
Finally, by adding equations (4) and (5) together, and making use of equation (6), we obtain an equation for the total energy density, $`E_{tot}`$:

$$\frac{E_{tot}}{t}=\rho \dot{S}_{dep}(t)\frac{E_{tot}}{t_e}\frac{E_{tot}}{\stackrel{~}{t}_d}$$
(7)
where the effective diffusion time, $`\stackrel{~}{t}_d`$, is defined as

$$\stackrel{~}{t}_d\frac{t_d+t_e+r_ft_dt_e}{r_ft_e}.$$
(8)
Written this way, it is easy to see that as $`r_f`$ – the UV to OIR fluorescence rate – decreases, the effective diffusion time goes up. In terms of the opacity, $`r_f`$ can be written in LTE as

$$r_f=c\rho \left(\frac{_{\nu _{OIR}}^{\mathrm{}}\kappa _\nu ^{}B_\nu ^{}𝑑\nu ^{}}{_{\nu _{OIR}}^{\mathrm{}}B_\nu ^{}𝑑\nu ^{}}\right)\times \left(\frac{_0^{\nu _{OIR}}\kappa _\nu ^{}B_\nu ^{}𝑑\nu ^{}}{_0^{\mathrm{}}\kappa _\nu ^{}B_\nu ^{}𝑑\nu ^{}}\right),$$
(9)
where $`B_\nu `$ is the Planck function at the local gas temperature, and $`\nu _{OIR}`$ is a frequency which approximately divides the UV from the optical. The first fraction is just the Planck mean opacity in the UV, and the second fraction is the fraction of absorbed energy which is radiated at optical wavelengths. Since, for instance, Co II has a larger opacity in the optical than Co III, increasing the ionization from once- to twice-ionized will decrease the fraction of absorbed energy which is re-radiated in the optical, making $`r_f`$ smaller, and therefore increasing the effective diffusion time.

We can now understand the underlying physics of the WLR: as the amount of 56Ni is increased, the interior temperature and ionization go up, which decreases $`r_f`$ and increases $`t_d`$. Stated another way, brighter SNe Ia (i.e. ones with more 56Ni) are hotter inside, and therefore have longer diffusion times (i.e. evolve more slowly through peak).

This cannot be the whole story, however. As long as most of the 56Ni is unmixed and its mass fraction $`X_{56}1`$ in the 56Ni-dominated layers, increasing the total amount of 56Ni in the star does not increase the local heating rate per gram, and therefore would not, by itself, result in a temperature increase. Overlooked so far is the fact that, if the additional 56Ni is added at higher velocity, then there will be an increase in the effective optical depth. That is, in stars with more 56Ni (in a layered composition), the velocity of the outer edge of the 56Ni is moving faster than is the outer edge of the 56Ni in a star with less 56Ni. The higher the velocity of the outer edge of the 56Ni, the more lines through which a photon must Doppler shift before escaping. This represents an optical depth increase which would increase the escape time and therefore make the gas hotter. Putting iron-peak elements at higher velocity is, by itself, not sufficient to increase $`\tau _{eff}`$ enough to fully explain the observed correlation between light curve width and luminosity. Rather, the slight increase in optical depth provided by putting 56Ni at higher velocities is amplified by its effect on temperature and ionization, and thus the escape time. It is difficult to cleanly separate these two effects, but in the next section we present results which strongly suggest that ionization is the dominant mechanism determining the escape time.

It is worth stressing here that fluorescence is not a one-way process. Energy can be transported to shorter wavelengths by “combining” photons. A photon can be absorbed, raising an ion to an excited state. The ion can then absorb another photon into this state, raising it to yet a higher-energy state from which a UV photon can be emitted. From obvious thermodynamic considerations (embodied in Rosseland’s Theorem of Cycles) this “combining” of photons must in general be less frequent than “splitting” of UV photons. Near the peak in the lightcurve, however, over much of the volume of the supernova these processes are nearly in equilibrium; it is the net rate of transport of energy to longer wavelengths which determines the escape time, not merely a one-way fluorescent cascade. Indeed, (Mazzali2000) has suggested that “combining”, in lower-optical depth regions above the photosphere, is a significant contributor to the UV emission at maximum light.

Another avenue for time-dependence to enter the equations above is through the energy deposition term $`\dot{S}_{dep}(t)`$. The deposition is determined both by how much 56Ni is present and by what fraction of its decay energy is deposited in the gas. Because the Compton opacity to MeV $`\gamma `$-rays is mainly absorptive, the deposition fraction is roughly proportional to $`1exp(\tau _\gamma )`$, where $`\tau _\gamma `$ is the mean gamma-ray optical depth and is in turn roughly proportional to the mean column depth to 56Ni. More-rapid expansion will lead to a lower mean optical depth to $`\gamma `$-rays at a given time, as will mixing radioactive material to higher velocities. Lower column depths will in turn lead to a more-rapid decline in the deposition fraction. If $`\dot{S}_{dep}(t)`$ declines more rapidly, the lightcurve will peak earlier. An earlier peak will have more 56Ni still present, and this will tend to offset the larger escape fraction. The lightcurve will thus evolve more rapidly, leading to a narrower lightcurve for it luminosity.

We have thus identified three properties of the explosion which may have a strong effect upon the WLR. The mass of 56Ni determines the ionization and through it the mean escape time. The expansion velocity (or specific kinetic energy) determines the rate of decline of the energy deposition. The radial distribution of 56Ni can have an effect upon both the escape time of the thermalized radiation and the time evolution of the deposition. In the following sections we present a small set of synthetic lightcurves which illustrate these effects and their relative importance.

## 3 Computations

A complete model for supernova lightcurves would begin with a model for the evolution of the progenitor. This would be evolved to ignition, the subsequent explosive burning followed to completion, and a transport calculation performed to determine the lightcurve. Unfortunately, a fully predictive model is not yet available from such a procedure. While there are schematic models for all of these stages, it is not yet possible to simulate SN Ia explosions without introducing very simple parameterizations to stand in for very complex physics we do not yet understand.

The evolution of the progenitor is undoubtedly the most uncertain part of this process, yet there is the possibility that “convergent evolution” of some sort, perhaps ignition near the Chandrasekhar mass in single degenerate systems, reduces the sensitivity to initial conditions. On the other hand, the initial conditions in the various double-degenerate scenarios might well be the single most important piece of physics determining the subsequent explosion. In the burning phase, the physics of turbulent combustion and the possible spontaneous transition to detonation are probably the most important and least tractable effects. While progress is being made on all of these fronts (Reinecke et al., 1999; Iwamoto et al., 1999), it will be a long time before truly predictive models of progenitor evolution and explosion become available.

Our aim in this work is to explore the radiation transport effects which influence the WLR. We have thus chosen to construct a very simple model for the structure of the explosion. We start with model DD4 of Woosley & Weaver (1991), an $`\mathrm{M}_{\mathrm{ch}}`$ explosion which burned 1.26M to Si or heavier elements to yield $`1.2\times 10^{51}`$ ergs of kinetic energy, producing 0.63 M of 56Ni. We have modified the mass and distribution of 56Ni in this model to explore the effects such modifications have on the lightcurve. The density as a function of velocity in the original model DD4 is preserved in all models. The structure of DD4 is presented by Woosley & Weaver (1991) and by PE II.

An initial grid of models with different 56Ni masses was obtained by mapping between 56Ni and 28Si according to

$`X_{56}(M)`$ $`=`$ $`X_{56}^0(M)\mathrm{exp}\left(\varphi (M)\right)`$
$`X_{28}(M)`$ $`=`$ $`X_{28}^0(M)+X_{56}^0(M),\left(1\mathrm{exp}\left(\varphi (M)\right)\right)`$
where

$$\varphi (M)\mathrm{max}(\alpha (MM_{cut}),0).$$
(10)
To increase the 56Ni mass, all isotopes below $`M=M_{cut}`$ were mapped to 56Ni:

$$X_{56}(M)=1H(MM_{cut})+X_{56}^0(M)H(MM_{cut}),$$
(11)
where $`H(x)`$ is the Heaviside function, and for $`A56`$,

$$X_A(M)=X_A^0(M)H(MM_{cut}).$$
(12)
This simple scheme is actually well-motivated in a physical sense. The total energy released from burning a C/O mixture to intermediate elements in the Si group is very nearly the same as for burning to iron-peak elements. The 56Ni mass thus need have little necessary connection to the explosion kinetic energy. The most stringent constraints on the distribution of elements in velocity come from peak-phase spectroscopy, though this analysis is hampered by our inability to compute time-dependent NLTE models.

As discussed in PE I, the initial thermal energy in the explosion goes entirely into kinetic energy via the second term on the r.h.s. of equation (4) and hence has no effect upon the observable part of the lightcurve. All the models presented here were started with the same initial internal energy.

The lightcurves were computed as described in PE II, using time-dependent, multi-group radiation transport with an LTE equation of state and including the elements He, C, O, Si, S, Ca, Fe, Co, and Ni in ionizations up to Ni XIII. 6000 energy groups were employed, providing moderate resolution; doubling the number of groups does not significantly alter the results, though halving it does leads to marked effects. We are forced to employ LTE level populations because of computer-time constraints – a NLTE calculation would take as long per time step as a whole lightcurve in LTE (several days on a fast workstation). All current NLTE calculations (Höflich, 1995; Nugent et al., 1997; Mazzali & Lucy, 1993) employ temperature structures which are not consistent with the history of the ejecta’s evolution and therefore cannot be expected to give correct time-dependent spectra and lightcurves, as discussed in PE II. While these LTE calculations obviously suffer as well, the ability to include time-dependent physics correctly is evidently necessary to study the effects discussed above.

The WLR is an observed relation between peak magnitude and a measure of width. To obtain an accurate result in the face of measurement error and the temporal sparsity of observed data, both the peak magnitude and the width are measured by fitting observed data to templates determined from well-observed supernovae (Hamuy et al., 1996b; Riess et al., 1995). It is important to employ the same methods to measure the peak magnitude and width in our models. We “observe” our synthetic lightcurves by employing the MCLS (Multi-Color Lightcurve Shape) method of Riess et al. (1995), employing the revised vectors used in Riess et al. (1998) (kindly provided by Dr. Riess). We have also determined the $`\mathrm{\Delta }m_{15}(B)`$ decline-rate parameter for our models using the $`\chi ^2`$ template-fitting technique described by Hamuy et al. (1996c). The results of fitting from both of these techniques are shown in table 1.

We note that the observational difficulty of obtaining early lightcurve data has hindered the construction of templates for times earlier than ten days before peak. Thus, by construction, template-fitting techniques are most sensitive to post-maximum decline rates; their ability to discriminate between different pre-maximum lightcurve shapes is seriously hampered.

Observationally, WLR relations with similar scatter exist in the B, V, R, and I bands – color evolution during the first 30 days is well correlated with luminosity and width. At later times, all SNe Ia show a remarkable homogeneity in color. In the absence of reddening, determining the width of the lightcurve in a single band is sufficient to uniquely determine a supernova’s peak magnitude. The importance of fitting a supernova’s color evolution is primarily in determining extinction. As our models do not suffer from this observational uncertainty, it is sufficient to fit to a single band.

We remind the reader as well that individual supernovae show significant departures from the best-fit templates (cf. Hamuy et al., 1996a); SNe Ia do not form a strictly one-parameter family. While the width of lightcurves measured by the templates is the largest variation from object to object, there are other features in the lightcurve which are not correlated with the decline rate. A good example of this is the pair SN1991T and SN1992bc which have very similar decline rates over the first 15 days past peak but very different rates over the subsequent 45 days (Hamuy et al., 1996c).

The details of the atomic physics and NLTE effects determine the color evolution of the supernova. We have assumed in this study that such effects are less important to the bolometric luminosity, merely reapportioning flux between the B, V, and R bands which comprise most of the luminosity (Suntzeff, 1995). We have thus chosen to take fits to the WLR in the V band only as most representative of the overall relation between width and luminosity. Below we present fits to B, V, and R simultaneously as well, but as an accurate determination of the color evolution depends more upon specific details of the atomic models and NLTE effect, we regard these fits as less reliable indicators of the underlying physics than the trends in a single band.

As the ejecta evolve, the effects of trapping and fluorescence become less important to the bolometric lightcurve. The colors at these later times, however, are very strongly influenced by NLTE effects. LTE calculations cannot be expected to provide an adequate representation of the lightcurve physics as the supernova becomes a nebula. Thus, we fit only to the peak phase of the lightcurve and estimate the sensitivity of fitting to various time intervals by fitting from 10 days before peak (the earliest time in the templates) until 20 and 40 days past peak.

## 4 Results

Figure (1) compares the MCLS $`\mathrm{\Delta }`$ parameter (a measure of width) to the peak V absolute magnitude determined from the models described above. The MCLS method parameterizes the shape of a supernova’s lightcurve as

$$M(t)=A(t)+\mathrm{\Delta }\times L(t)+\mathrm{\Delta }^2\times Q(t)$$
(13)
where $`A(t)`$ is, roughly, the mean lightcurve shape of a sample of 27 well-observed supernovae (Riess et al., 1998) and $`L(t)`$ and $`Q(t)`$ are functions describing a quadratic fit to the departure of individual supernovae from this mean shape as parameterized by $`\mathrm{\Delta }`$ which is thus a measure of lightcurve width. The quadratic relation between $`\mathrm{\Delta }`$ and $`V_{peak}`$ is indicated by the diagonal curves in the figure. The absolute luminosity of the relation is calibrated by the Cepheid distance scale of (Gibson et al., 2000) corresponding to a Hubble constant of 68 km s-1 Mpc-1. The square symbols correspond to the models with 0.27, 0.38, 0.48, 0.55, 0.63, and 0.90 Mof 56Ni in the structure of model DD4. The star symbols represent other models to be discussed below. The top panel shows fits over the interval (-10,20) days relative to peak and the lower panel over the interval (-10,40). The quality of the fit to the mean WLR relation is somewhat better for the scaled-56Ni models for the shorter interval and is probably a consequence of the increasingly inappropriate assumption of LTE at later times.

In both cases, the fits show a dispersion comparable to or smaller than the observations ($`0.15`$ mag) and (remarkably and probably fortuitously) virtually zero departure from the absolute calibration. It is clear from this figure that the V WLR is fit remarkably well by the simple prescription of changing only the 56Ni mass.

Figure (2) shows the same lightcurves fit by the template-fitting ($`\mathrm{\Delta }m_{15}(B)`$) method of Hamuy et al. (1996b) on the interval (-5,26) days. In this method, the models are fit to each of six template lightcurves spanning the range of well-observed supernovae. The templates are characterized by their $`\mathrm{\Delta }m_{15}(B)`$ parameters (roughly the decline in B magnitude over 15 days from peak), while the actual epochs fit are determined by the whole peak phase of the lightcurves. The best-fit $`\mathrm{\Delta }m_{15}(B)`$ is then determined from the minimum of a quadratic fit to the reduced $`\chi ^2`$ as a function of the parameter. Rather than extrapolating, models broader than the broadest template are assigned the minimum value of $`\mathrm{\Delta }m_{15}(B)`$, accounting for the slight pile-up of points at the bright end of the relation.

The filled symbols are the same as those of the previous figure, and the open symbols are fits to 59 well-observed supernovae from Phillips et al. (1999). The curve is the mean $`\mathrm{\Delta }m_{15}(B)`$-$`M_V`$ relation, on the same calibration, from (Phillips et al., 1999). Once again, the models fit the relation well within the scatter of observed supernovae. The significant difference in curvature between the mean $`\mathrm{\Delta }m_{15}(B)`$-luminosity and $`\mathrm{\Delta }`$-luminosity relations is notable and appears to be the result of the differing ways in which the lightcurve shapes are parameterized; both techniques achieve similar reductions in the scatter about the Hubble relation.

Figure (3) shows the V lightcurve family from the six scaled DD4 models compared with the best-fitting MCLS templates. The shapes of the lightcurves depart less than 0.05 magnitude from the templates throughout the first 45 days (about 25 days past V peak) with the exception of the 0.90 Mmodel, where the inflection at 32 days is not seen in the templates. This error is well within the dispersion of data from individual supernovae about the template lightcurves (cf. Riess et al., 1998). At later times, the LTE approximation becomes increasingly poor and may account for the too-shallow decline rates, though real supernovae show similar departures from the templates at this time.

Figures (4) and (5) show the B, V, and R lightcurves of the 0.55 Mand 0.90 Mmodels, respectively. The DD4/55 model displays a generic problem with most of the models; while the B lightcurve is the correct width, near peak it is too flat compared to the templates. The departure from the templates before peak is less disturbing, and is probably more a problem with the templates than the models, as there is comparatively little data so far before maximum. The IR secondary maximum is clearly seen in the R-band curves.

We now turn to variations on these scaled-DD4 models to investigate the sensitivity of the WLR to other changes in the explosion. The models presented so far have nearly pure 56Ni from the center out to a mass coordinate nearly equal to the Ni mass – radioactive material is concentrated as much as possible toward the center. This is a natural result of the explosion physics. White dwarfs are densest at the center. When matter burns explosively, the higher the density, the larger the ratio of the hydrodynamic timescale to the the burning timescale, and the more the distribution of burning products is skewed toward heavier elements. In the absence of significant hydrodynamic mixing after material has cooled below the burning temperature, a natural consequence of explosive burning is that the gradient of atomic weight will point inward, with the heaviest burning product, 56Ni, concentrated toward the core.

To explore the effect hydrodynamic mixing might have on the lightcurve, we have altered the distribution of 56Ni in model DD4 without changing its mass. A significant mass fraction of 56Ni extends approximately to 0.85M in the original model (primarily as a consequence of producing 0.16M of non-radioactive iron-group material). The starred points at $`M_V=19.56`$ in Figures (1) and (2) show the result of homogenizing the composition of this model over the mass range 0.7 to 1.1 Min model M1. This decreases the peak bolometric magnitude by 0.1 mag and leads to significantly narrower B and brighter I lightcurves. The resulting V and R magnitudes are virtually identical, 0.04 and 0.02 mag brighter, respectively. The 56Ni nearer the surface cooled more rapidly. The lower ionization gas converts energy from shorter- to longer-wavelength photons more efficiently, circumventing the effect of the higher expansion optical depth from 56Ni at higher velocities. Altering the distribution of 56Ni would appear to affect the color of the emergent flux more than its overall escape properties. An important conclusion is that, within the limits of fitting 0.63M of 56Ni in an $`\mathrm{M}_{\mathrm{ch}}`$ white dwarf, the WLR is insensitive to the distribution of radioactivity in the explosion.

The scaled-DD4 models traded off production of 56Ni with intermediate mass elements. This is the only possibility afforded if the supernovae are to satisfy nucleosynthetic constraints. For example, the production of large quantities of non-radioactive iron-group elements in the majority of supernovæ is firmly ruled out. Nonetheless, in model NiFe we have converted, in place, half of the mass of 56Ni to 56Fe, thus changing the heating rate per gram of the iron-rich material without changing its distribution, cooling, or $`\gamma `$-ray escape properties. The lower temperature in the 56Ni-rich material leads to lower ionization and hence to more efficient escape. In spite of having 0.32 instead of 0.27M of 56Ni, the model has a much narrower lightcurve than model DD4/27, as shown by the stars near $`M_V=19.01`$ in the Figures 1 and 2. In model M1 56Ni was mixed outward with overlying Si-group material. However the amount of material involved in the mixing was a small fraction of the total mass of 56Ni. Consequently the effectiveness of mixing in lowering the heating rate per gram and narrowing the lightcurve was small.

Interchanging Fe- with Si-group elements makes relatively little difference in the amount of energy released by burning a C/O mixture; at most $`6.2\times 10^{17}`$$`\text{erg}\mathrm{g}^1`$ is released burning to Si & S versus $`7.9\times 10^{17}`$$`\text{erg}\mathrm{g}^1`$ for burning to NSE. It is thus quite possible to imagine a series of explosions in which the Si/Fe ratio varies more strongly than the specific kinetic energy. On the other hand, in most current models a greater production of 56Ni is accompanied by a greater total amount of burning, leading to higher velocities. Model DD3 of Woosley & Weaver (1991) produced 0.96 M of 56Ni. It burned 10% more mass to the Si-group or above in total and has 10% greater kinetic energy than model DD4. In spite of its high 56Ni mass, the lightcurve from this model is considerably narrower than that of model DD4/90 and is shown by the stars in Figures (1) and (2) at $`M_V=19.90`$. The higher kinetic energy in this explosion resulted in a lower column depth, allowing decay $`\gamma `$-rays to escape more readily after peak. By 20 days, $`6`$% of the decay energy in model DD4/90 escaped directly as $`\gamma `$-rays, increasing to 35% by day 40. In model DD3, the escape fraction was more than double at 20 days, at $``$14%, rising to 46% by day 40. A shortened time to $`\gamma `$-ray transparency leads to a shorter rise-time (by 3.5 days), a narrower peak, and a greater decline rate in the lightcurve. The increased $`\gamma `$-ray escape was, however, more than compensated for by the shorter rise-time, and the peak magnitude is even brighter than the difference one might predict based solely on consideration of the difference in total 56Ni mass. It is interesting to note that the $`\mathrm{\Delta }m_{15}(B)`$ fitting procedure did not pick up this variation in lightcurve shape as it weights more heavily the decline rate after peak than the MCLS templates; the two models had very similar post-maximum declines.

## 5 Conclusions

While far from a thorough investigation of all the variations one can imagine on explosions of $`\mathrm{M}_{\mathrm{ch}}`$ white dwarfs, the results presented above are suggestive. They show that it is possible to create a series of models which exhibit, to the best of our current calculation ability yet well within observational errors, the observed WLR.

It seems clear from the arguments of §2 and the results of these simple experiments that the 56Ni yield and the $`\gamma `$-ray escape fraction play the determining roles in the relation between the lightcurve’s width and peak luminosity. The models show a remarkable insensitivity to the distribution of radioactivity and to the composition of non-radioactive material, within limits set by achieving the appropriate explosion energies and acceptable nucleosynthesis.

The mass of 56Ni sets the scale of the peak luminosity and the heating rate of the gas. The resulting temperature determines the ionization and hence the cooling rate of the gas. More 56Ni leads to more heating, higher temperatures, less efficient cooling, and hence broader and brighter lightcurves.

The efficiency of $`\gamma `$-ray escape sets the evolution rate of the energy deposition. Models with higher overall velocities and/or with 56Ni distributions extending to higher velocities exhibit larger $`\gamma `$-ray escape fractions at a given time. This increases the evolution rate of the lightcurve, with greater escape leading to narrower lightcurves. Models with higher 56Ni yields tend to have their radioactivity distributed over a broader range in velocity, leading to larger escape fractions. This greater escape fraction acts to oppose the tendency for broader lightcurves resulting from higher temperatures, decreasing the slope of the WLR.

We have reproduced the observed relation in V by varying only the 56Ni yield. One solution to reproducing the observed WLR is thus to produce an evolution scenario which leads to supernovæ with similar explosion energies yet differing 56Ni yields. Because higher 56Ni yields and higher velocities tend to produce opposing effects on the lightcurve, this solution is almost certainly not unique. A weak dependence of energy on 56Ni yield will likely lead to acceptable lightcurves as well.

While to date the models which best fit lightcurves and spectra are all at the Chandrasekhar mass (with the possible exception of models for SN1991bg), the masses of all SN Ia explosions may not be identical. Varying the explosion mass primarily affects the $`\gamma `$-ray optical depth. As long as a given 56Ni mass is accompanied by a column depth similar to those exhibited by the $`\mathrm{M}_{\mathrm{ch}}`$ models considered here, a variable-mass explosion scenario would not appreciably alter our conclusions. For sub-$`\mathrm{M}_{\mathrm{ch}}`$ model such as those of (Livne & Arnett, 1995; Woosley & Weaver, 1994), the 56Ni produced by detonation in surface layers does not contribute appreciably to the optical lightcurve due to the low $`\gamma `$-ray optical depth. In such models, only the 56Ni produced by carbon detonation in the core C/O dwarf would affect the WLR.

The exact values of the slope and calibration of the WLR obtained in this work must be, to some extent, fortuitous. EP II showed that the conditions experienced by matter in SNe Ia are not consistent in detail with the assumption of LTE. The shape and colors of the lightcurve depend sensitively upon the average branching ratio from UV to OIR photons through a broad range in energies. While many strong lines are unlikely to be missing from our atomic models, all such models are seriously incomplete and the statistical properties of the line distributions and branching ratios are unlikely to be correct in detail. The basic physics, however, will not change significantly, and thus we have some confidence in the general trends the models indicate.

By fitting our models using the same procedures used to calibrate observed SNe Ia, we have demonstrated that the calibration is insensitive to fairly large changes in the underlying explosions. This suggests that the small differences in initial conditions which might arise from evolutionary effects between the local universe and that at $`z1`$ are unlikely to appreciably affect the supernova cosmology results. For example, mixing in model M1 leads to much less than 0.1 magnitude difference from the corresponding unmixed model. Even the unphysically large change in abundances in model NiFe, strongly ruled out by chemical evolution considerations, leads to less than 0.2 magnitude departures from the MCLS and barely 0.3 magnitude departures from the $`\mathrm{\Delta }m_{15}(B)`$ relations; allowable changes in relative iron-peak abundances would lead to insignificant effects.

Altering the specific kinetic energy, as in Model DD3, is the only way we see to achieve significant departures from the WLR. Increasing the kinetic energy by 10% was seen to lead to a much narrower lightcurve (at least when fit by the MCLS templates), primarily by more than doubling the $`\gamma `$-ray escape from 6% at peak in DD4/90 to 14% in DD3. This change is in the wrong direction, however. To account for the 0.3 magnitude deficit in brightness seen at large redshifts, a shift to lower explosion energy would be required. Halving the $`\gamma `$-ray escape from 6% at peak to 3% is a smaller fractional change in deposition and would lead to a smaller change in width. A future paper will explore these effects in greater detail.

We wish to thank M. Hamuy for many useful discussions and for performing the $`\mathrm{\Delta }m_{15}(B)`$ fits to the lightcurves, and A. Riess for supplying his lightcurve templates in advance of publication. This work has been supported by the National Science Foundation (CAREER grant AST9501634), by the National Aeronautics and Space Administration (grant NAG 5-2798). This work was also performed under the auspices of the U.S. Department of Energy by University of California Lawrence Livermore National Laboratory under contract No. W-7405-Eng-48. Philip Pinto gratefully acknowledges support from the Research Corporation though a Cottrell Scholarship.


format: text

text:
# Temperature and Density Gradients Across the Nucleus of M82

## 1 Introduction

M82 (NGC 3034, Arp 337) is considered the prototypical starburst galaxy. Its close proximity ($`D`$=3.25 Mpc, Tammann & Sandage (1968)) and strong CO lines (Rickard et al., 1975) make it an excellent location to study the effects of a burst of star formation on molecular gas properties. Its edge-on orientation make essential the use of optically thin molecular gas tracers for probing the physical conditions in the interstellar medium. Interferometric maps (e.g. Shen & Lo (1995); Neininger et al. (1998)) show that M82 contains three bright lobes: two main lobes on the outer edges, which are thought to be a torus of molecular gas seen edge-on (e.g. Tilanus et al. (1991)), and a central peak which may be housing an AGN (Muxlow et al., 1994; Seaquist, Frayer, & Frail, 1997). Early CO studies showed that M82 had an unusually large 12CO $`J`$=2$``$1/$`J`$=1$``$0 line ratio ($``$2, Knapp et al. (1980); Sutton, Masson, & Phillips (1983)), which suggested that the 12CO gas in M82 is optically thin. As calibration techniques and telescope surface quality improved, the value for the 12CO $`J`$=2$``$1/$`J`$=1$``$0 line ratio has decreased to $``$1 (Tilanus et al., 1991; Wild et al., 1992), which suggests that 12CO is likely optically thick. A high optical depth leads to some problems in the interpretation of the double lobes as a molecular torus. Neininger et al. (1998) point out that if the gas were optically thin, a torus seen from the edge would have a double peaked appearance, while if the torus were optically thick, we should see a elongated, bar-like structure.

Mapping a galaxy in many transitions will help our understanding of the molecular gas dynamics since different transitions and molecules do not necessarily trace the same gas (Petitpas & Wilson, 1998b; Kohno, Kawabe, & Vila-Vilarò, 1999). While M82 has been painstakingly mapped and studied by many authors in the lower-$`J`$ transitions of CO and its isotopomers (see references above), improvements in receiver sensitivity and telescope surface quality have made fully sampled maps of the higher-$`J`$ transitions of rarer CO isotopes more tractable. This paper presents first results from 12CO, 13CO, C18O $`J`$=3$``$2 and 13CO $`J`$=2$``$1 maps of M82. In §2 we describe the observations and discuss the appearance of the spectra in comparison to each other and previously published data. In §3 we use integrated intensity line ratios to derive the physical conditions inside the nucleus of M82 and discuss the results and possible causes for the observed gradients.

## 2 Observations and Data Reduction

One of the difficulties in the interpretation of line ratios is that the different lines are often measured with different telescopes. This method introduces calibration uncertainties (see for example, Figure 6 of Wild et al. (1992)) that can be reduced if all the data are taken with the same telescope (which is often impossible, given the large variation in frequency between some transitions).

A variety of CO $`J`$=3$``$2 observations of M82 were taken using the James Clerk Maxwell Telescope (JCMT) over the period of 1996 April 5 – 10. The 12CO $`J`$=3$``$2 map covers an area of 120″ $`\times `$ 100″, the 13CO $`J`$=3$``$2 map covers 56″ $`\times `$ 35″, and the C18O $`J`$=3$``$2 map covers a 34″ strip along the bar major axis. We also include 12CO and 13CO $`J`$=2$``$1 data taken during CANSERV at the JCMT over 1993 April 22 – 25. The data were reduced using the data reduction package SPECX. The raw data had a linear baseline removed and then the $`J`$=3$``$2 data were binned to 12.5 MHz (10.8 km s-1) and the $`J`$=2$``$1 data to 5 MHz (6.7 km s-1). Since the CO $`J`$=2$``$1 data have a different beam size and grid spacing than the CO $`J`$=3$``$2 data, it was necessary to interpolate and convolve the spectra before measuring the integrated intensity line ratios. The data were exported to the data reduction package COMB, and the CO $`J`$=3$``$2 data were convolved to the same beam size and interpolated to the same location as the nearest CO $`J`$=2$``$1 spectra.

The calibration was monitored by frequently observing both planets and spectral line calibrators. The spectral line calibrators for the CO data had intensities that were within $``$15% of the published values. We therefore adopt the published values for $`\eta _{\mathrm{MB}}`$ of 0.58 for the CO $`J`$=3$``$2 data and 0.67 for the 12CO $`J`$=2$``$1 data. Pointing was checked frequently and was determined to be accurate to within 2″ r.m.s.

The individual 12CO $`J`$=3$``$2, 13CO $`J`$=3$``$2, C18O $`J`$=3$``$2, 12CO $`J`$=2$``$1, and 13CO $`J`$=2$``$1 spectra for the three bright regions of M82 are shown in Figure Temperature and Density Gradients Across the Nucleus of M82. The individual spectra for M82 appear asymmetric on either side of the center, indicating rotation. The line profiles appear similar for each transition which indicates that, at least on large scales, the emission is coming from co-moving material. There are indications of double peaked line profiles at the (0,0) position in the $`J`$=2$``$1 spectra. Since we do not see this structure as strongly in the $`J`$=3$``$2 spectra, it is likely the result of the larger beam overlapping both the NE and SW peaks simultaneously.

## 3 Physical Conditions from Line Ratios

It is well known that the filling factor of molecular clouds in a galaxy is often low. Thus, integrated intensities are generally not a good measure of the physical conditions of individual molecular clouds because the signal is diluted by the large radio telescope beams, which also cover much empty space. We must therefore use integrated intensity line ratios to determine the physical conditions of individual clouds. Using the ratio of two integrated intensities will help cancel out the effects of beam dilution, assuming that similar regions of space are responsible for the emission at both frequencies. The physical conditions recovered from the analysis of line ratios are the average conditions for all clouds within the beam.

The differing morphology in the 13CO and 12CO $`J`$=3$``$2 lines results in a gradient of the 12CO/13CO $`J`$=3$``$2 line ratios, increasing from east to west (Figure Temperature and Density Gradients Across the Nucleus of M82). If both lines are optically thin, this observation indicates a 13CO/12CO abundance gradient across the disk of M82. This same ‘abundance’-type gradient is not visible in the lower $`J`$ transitions (see 12CO/13CO $`J`$=1$``$0 and $`J`$=2$``$1 ratios; Neininger et al. (1998); Loiseau et al. (1990)), which suggests that the 12CO emission is not optically thin and the gradient may be the result of a variation in optical depth. Since optical depth is a function of density and temperature, we have performed a Large Velocity Gradient (LVG) analysis using a code written by Lee Mundy and implemented as part of the Miriad data reduction package. A single component model provided a very nice fit to the observed line ratios. We kept the C18O abundance as a free parameter, so a survey of parameter space was performed using $`T_{\mathrm{kin}}`$ = 50 $`\pm `$ 20 K (Wild et al., 1992), \[12CO\]/\[13CO\] = 50 $`\pm `$ 20 (Tilanus et al., 1991), $`n`$(H2) from 10 to 106 cm-3, and $`N`$(CO)/$`dv`$ from 1015 to 1020 cm-2/km s-1(see Petitpas & Wilson (1998a) for more details). The $`\pm `$1$`\sigma `$ error bars for the line ratios shown in Figure Temperature and Density Gradients Across the Nucleus of M82 were entered into the program and the solutions are shown in Table 1. The 12CO $`J`$=2$``$1/$`J`$=1$``$0 ratios from Tilanus et al. (1991) also overlap our solutions.

In the simplest scenario, the observed line ratio gradients can be explained by varying either temperature or density across the galaxy while holding the other variable constant. Our analysis suggests that either the density varies from $`1.6\times 10^4`$ to $`0.6\times 10^4`$ cm-2 (for $`T_{\mathrm{kin}}`$ = 50 K), or the temperature varies from 70 to 30 K (for $`n`$(H2) = 104 cm-2) from east to west. Of course, the gradients can also be produced by a combination of varying $`T`$ and $`n`$, which suggests we need to include observations at other frequencies to constrain the models.

A similar line ratio gradient was seen in HCO+ $`J`$=4$``$3/$`J`$=1$``$0 by Seaquist, Frayer, & Bell (1998), but since their data combined single dish and interferometry data, they dismissed the gradient, attributing it to an abundance of diffuse emission on the eastern side of M82, which would cause the HCO+ emission to be underestimated. However, our model fits show a higher density and/or temperature in the NE lobe, which could produce an increased HCO+ $`J`$=4$``$3/$`J`$=1$``$0 ratio. Previous 12CO and 13CO $`J`$=1$``$0 interferometric studies by Kikumoto et al. (1998) find 12CO/13CO $`J`$=1$``$0 ratios that increase across M82, but in the opposite direction from our $`J`$=3$``$2 ratios. This result helps rule out the possibility that the observed gradients are the result of a 13CO and C18O abundance gradient across M82. The $`J`$=1$``$0 line ratio gradients are likely the result of the higher temperatures and/or densities in the NE lobe of M82, which bump the CO into the higher $`J`$ transitions, leaving the lower $`J`$ levels less populated.

We can use the 850 $`\mu `$m continuum maps of M82 (Alton, Davies, & Bianchi, 1999) to remove the degeneracy between density and temperature in our models. Since the SW lobe is brighter in the 850 $`\mu `$m continuum maps, which trace the product of the temperature and the column density, we can conclude that either the column density and/or temperature is higher in the SW lobe. In this scenario, the HCO+ maps of Seaquist, Frayer, & Bell (1998) would require that density be higher in the NE. The pair of LVG solutions which best satisfy the combined CO, HCO+, and continuum data are flagged with check marks in Table 1. Utilizing these additional data sets, we therefore conclude that the observed line ratios are the result of a temperature gradient that increases from NE to SW in conjunction with a density gradient that increases from SW to NE.

One possible reason for this density and temperature gradient may be interaction with the nearby galaxy M81. Gravitational interaction may be triggering vigorous star formation in the SW lobe, which could be heating the gas (increasing $`T`$) at the same time that it is being consumed (decreasing the average density by turning the densest regions into stars). While previous galaxy interaction simulations do not show any indication of density gradients across the smaller galaxy, it is unlikely that the current simulations would be able to detect a factor of three difference in the mean density of molecular clouds across the galaxy.

There is some dispute over the nature of the double peaked structure seen in M82. Most studies suggest that it is the result of a molecular torus seen edge-on. This would only be the case if the molecular gas in the clouds were optically thin (Neininger et al., 1998), or there were not enough molecular clouds in the telescope beam to ‘shadow’ the more distant clouds. Our LVG analysis indicates that 12CO $`J`$=1$``$0 transition is on the border between optically thin and optically thick ($`\tau _{J=10}0.13.5`$). However, the $`J`$=2$``$1 and $`J`$=3$``$2 maps also indicate a double peaked structure even though these transitions are likely optically thick (Table 1). This result can only be consistent with the torus model if the filling factor of the optically thick clouds within the beam is low enough that shadowing of one cloud by another is not a problem. We can estimate the filling factor by comparing the predicted and observed line temperatures. The LVG models predict temperatures of 19 K for 12CO $`J`$=2$``$1 while our data show temperatures of $``$ 6 K. In addition, the high resolution CO maps of Neininger et al. (1998) suggest that M82 is only $``$ 10″ thick and thus the filling factor of the CO emission is $``$ 60% of the 22″ JCMT beam. This filling factor is high enough that cloud shadowing may be important in this region. If cloud shadowing plays a large role, the double peaked structure in M82 could not be explained as the result of an edge-on torus of molecular gas. Instead, it could be produced by the accumulation of molecular gas at the Inner Lindblad Resonance radius as it flows inward along the bar (e.g. Kenney et al. (1992)). There is evidence for these double peaks (on similar size scales) in nearly face-on barred galaxies (Kenney et al., 1992; Petitpas & Wilson, 1998b), where we would not expect them to be caused by torus of gas, since the torus would have to be out of the plane of the galaxy.

We wish to thank Lorne Avery and the JCMT staff for taking the 12CO and 13CO $`J`$=2$``$1 data during service observing. The JCMT is operated by the Royal Observatories on behalf of the Particle Physics and Astronomy Research Council of the United Kingdom, the Netherlands Organization for Scientific Research, and the National Research Council of Canada. This research has been supported by a research grant to C. D. W. from NSERC (Canada).


text:
# Magnetization profiles and NMR spectra of doped Haldane chains at finite temperatures

## I Introduction

The physics of low-dimensional spin systems is extremely rich and shows many surprising effects. Using a valence bond picture it has been suggested that open segments of $`S=1`$ spin chains should have $`S=1/2`$ excitations localized at the ends of the chain , that are experimentally observable. In a recent experiment by Tedoldi et al  the presence of these $`S=1/2`$ chain end spins was beautifully demonstrated by measuring the NMR profile of Mg doped Y2BaNiO5, a model $`S=1`$ spin chain system, showing well defined satellite peaks indicative of the chain end excitations. These satellite peaks were, perhaps surprisingly, very well defined at temperatures close to and well above the gap, $`\mathrm{\Delta }100K`$. At temperatures below the gap, where one would expect the signal from the chain end spins to be most clearly observable, the satellite peaks were smeared into a single broad line. The inter chain coupling $`J^{}`$ in Y2BaNiO5 is very small $`(|J^{}/J|5\times 10^4,J285K)`$  and it therefore remains a question why no well defined satellite peaks were observed below the gap. Bulk susceptibility and magnetization measurements  on Zn doped Y2BaNiO5 have also been interpreted in terms of $`S=1/2`$ chain end excitations, but in this case at temperatures all the way down to 4-10K, well below the gap. Very interestingly, a sub-Curie power law divergence was observed at temperatures below 4K hinting at the presence of a gapless phase. In light of these two experiments it is then natural to ask, from a theoretical point of view, if chain end excitations are observable at finite fields and temperatures well above the gap and when a simple one-dimensional model will fail to describe the experiments. In the present paper we systematically address these issues by calculating the on-site magnetization, $`S_i^z`$, at finite temperatures and fields using continuous time quantum Monte Carlo techniques and reconstructing the associated NMR spectra for a given concentration of impurities.

It is by now well established that integer spin chains display a gap above a spin liquid ground-state  as opposed to half-integer spin chains which are gapless. The gap in integer spin chain systems has been observed experimentally in numerous compounds  as well as in many numerical studies . For $`S=1`$ chains the gap is known to be $`\mathrm{\Delta }0.4105J`$ . For integer spin chains, described by a simple antiferromagnetic Heisenberg model, it is possible to understand the low temperature physics using the intuitive valence bond picture. From this point of view the ground-state is close to the valence bond solid (VBS) state. Each integer spin $`S`$ is decomposed into $`2S`$ spin $`S=1/2`$ in the symmetric state and each of the $`S=1/2`$ on neighboring sites are paired into singlets. It is possible to determine Hamiltonians where this VBS state is the exact ground-state. For $`S=1`$ spin chains the corresponding Hamiltonian is rather simple and apart from the usual Heisenberg term includes a biquadratic term $`\beta (𝐒_i𝐒_{i+1})^2`$. For $`\beta =1/3`$, the AKLT point, the VBS state is the exact ground-state as was shown by Affleck-Kennedy-Lieb and Tasaki  (AKLT). Under periodic boundary conditions the VBS state is always a singlet. However, as was noted by Kennedy  for open $`S=1`$ chains, the ground-state is four-fold degenerate in the thermodynamic limit. In the VBS picture this is easy to understand since the $`S=1/2`$ spins at the end of the chain are not paired into a singlet state, leading to a four-fold degeneracy for a $`S=1`$ system. In general a spin-$`S`$ chain should then have spin-$`S/2`$ chain-end excitations. Many real compounds are well described by the antiferromagnetic Heisenberg chain:

$$H=J\underset{i}{}𝐒_i𝐒_{i+1}.$$
(1)
Following the above discussion the ground-state should then be close to the VBS state. For a $`S=1`$ system, the low-energy states of an open segment is then a singlet and a triplet separated by a gap decaying exponentially with L, the length of the segment. For L odd the interaction between the two chain end excitations is ferromagnetic and the triplet is the lowest lying state, for L even the interaction is antiferromagnetic and the singlet is lowest.

In the following we will focus on $`S=1`$ spin chains relevant to the experiments on Y2BaNiO5. Numerically the $`S=1/2`$ chain end excitations have been observed by various methods : exact diagonalisation , Quantum Monte Carlo (QMC) , Density Matrix Renormalization Group (DMRG) . They subsist all through the Haldane phase $`1<\beta <1`$  but as $`\beta `$ is increased from $`0`$ to 1 the peak in the onsite magnetization, $`S_i^z`$, moves away from the first site of the chain. At $`\beta =1`$ a dimerized phase is reached where the $`S=1/2`$ chain end excitations cease to exist. At the AKLT point, $`\beta =1/3`$, it can be analytically shown that the decay of the on-site magnetization away from the end of the chain has a pure exponential form with a length scale exactly equal to the bulk correlation length. Experimentally the $`S=1/2`$ chain end excitations have been observed in EPR measurements in $`S=1`$ chains doped with non magnetic ions  and even in microscopically fractured undoped systems . Initially, experiments measuring the Schottky anomaly in Zn doped Y2BaNiO5  were interpreted in terms of $`S=1`$ excitations but later work  showed that an explanation in terms of $`S=1/2`$ chain end excitations is also possible when anisotropy is taken into account. These studies have mainly dealt with low temperature (i.e. $`T\mathrm{\Delta }`$) behavior of these systems where the only relevant states are the low-lying singlet and triplet. In the recent experiments by Tedoldi et al , a 89Y NMR study of doped system Y2BaNi1-xMgxO5 was performed, clearly indicating the $`S=1/2`$ chain end excitations for temperatures well above the gap where many excited states apart from the low-lying triplet and singlet should be populated. In order to understand the NMR spectra in this high temperature range we simulated the isotropic Heisenberg spin 1 chain using a continuous time Quantum Monte Carlo algorithm at finite temperature and finite field. The QMC method samples the partition function so that all (excited) states are correctly treated. For the purpose of understanding the physics at temperatures close to $`T0`$ we have performed Density Matrix Renormalisation Group calculations of the NMR spectra at $`T=0`$. Our results for temperatures above the gap are in very good agreement with the experimental data and clearly show that the $`S=1/2`$ chain end excitations indeed are observable at such high temperatures. At temperatures below the gap down to $`T=0`$ we find very well defined satellite peaks that at $`T=0`$ only are broadened by the distribution of chain lengths. In the experiments of Tedoldi  et al, no clear satellite peaks were observed at $`T=77K`$, well below the gap. This could be due to a limited resolution at this temperature . In the experiments of Payen et al  on doped Y2BaNiO5 the results were interpreted in terms of the chain end excitations at temperatures down to 4K and it therefore remains unclear when and if higher dimensional physics become important in doped Y2BaNiO5.

The paper is organized as follows. In section II we discuss the numerical methods used. Section III presents our results for the on site magnetization for different temperatures and fields. In section IV we turn to a discussion of the NMR spectra obtained from the magnetization profiles and finally in section V we discuss possible explanations why signals from the $`S=1/2`$ excitations are not clearly seen below the gap.

## II Numerical Method

We use a single-cluster continuous time  version of the loop algorithm . Spins 1 are simulated by two spins $`1/2`$ with symmetrized boundary conditions  in the imaginary time direction. In order to take into account the magnetic field, we use a hybrid method  where field-dependent loop flip probability and global flip of worldlines are used. Each Monte Carlo run consists of $`10^7`$ sweeps (a measurement is taken every 10 sweeps to avoid autocorrelations) preceded by $`10^5`$ sweeps of thermalization. The DMRG  calculations were performed keeping $`m=243`$ states for chain lengths between L=1,100 and $`m=81`$ states for chain lengths L=101,501. Only negligible differences were observed between the results for $`m=81`$ and $`m=243`$. The calculations were performed using the total $`z`$-component of the spin, $`S_{tot}^z`$, and parity, $`P`$, with respect to a reflection about the middle of the chain as quantum numbers. All the magnetization profiles were calculated in the $`S_{tot}^z=1,P=1`$ state.

## III Magnetization profiles

### A Influence of temperature

The first question we address is the temperature dependence of the on-site magnetization $`S_i^z`$. Due to symmetry the on-site magnetization is always zero in the absence of a magnetic field. In order to compare to experiments we therefore perform calculations in the sector with $`S_{tot}^z=1`$. In Fig. 1 we show the magnetization profile for an open chain of 29 spins for various temperatures in the sector $`S_{tot}^z=1`$. Filled circles are T=0 results obtained by DMRG, the other symbols are QMC data. We observe the well-known alternating structure due to the antiferromagnetic interactions which tends to be diminished as the temperature is raised. For the lowest temperature $`T=0.1J`$ (well below the gap), the profile is almost identical to the ground-state result. The decay of the on-site magnetization away from the end of the chain is at $`T=0`$ known to be exponential  with a length scale equal to the bulk correlation length $`\xi 6`$. At high temperature the staggered magnetization decays rapidly from its value at the end of the chain corresponding to a shorter correlation length. At the highest temperatures, the middle spins (around $`i=15`$) have become polarized and show an average magnetization of $`<S_i^z>=\frac{1}{L}0.034`$. At temperatures relevant to the experiments of Tedoldi et al , $`T=0.3J0.8J`$, the structure from the $`S=1/2`$ chain end excitations is still clearly visible. All these results confirm the existence of $`S=1/2`$ excitations at temperature above the gap, even if the magnetization at the end of the chain is significantly reduced above the gap.

### B Influence of magnetic field

Lifting the constraint of working in the $`S_{tot}^z=1`$ sector we now consider the on-site magnetization as a function of the applied magnetic field at temperatures and fields relevant to the experiments. We calculated the dependence on the magnetic field of the profiles for two temperatures : $`T=0.1J\mathrm{\Delta }`$ and $`T=0.526J>\mathrm{\Delta }`$.

For $`T=0.1J`$ (figure 2), the zero-field magnetization is zero for all spins (by symmetry). At this temperature, the only magnetic states are in the sectors $`S_{tot}^z=\pm 1`$ of the triplet (which is the lower state for an odd chain) and are the only ones affected by the field. As we increase the field, the $`S_{tot}^z=1`$ state will be energetically favored and hence the magnetization starts to take on the alternating structure showing the presence of the chain end excitations. At sufficiently large fields we should obtain results comparable to the results discussed in the previous section obtained in the absence of a magnetic field but in the $`S_{tot}^z=1`$ sector. This is clearly the case and for $`H=0.3J`$, we nearly recover the zero-field structure of the $`S_{tot}^z=1`$ sector at $`T=0.1J`$ shown in Fig. 1. If the field is too strong (say $`H>0.3J`$), magnetic states from different sectors above the Haldane gap contribute to the thermal expectation of the on-site magnetization and, even if the alternating structure is still present, we don’t recover the results from the zero-field case. This effect is visible on the $`H=0.4J`$ data (left triangles). For this field, $`<S_i^z>`$ for odd spins (which all have a positive expectation value) is nearly equal to the value in the $`S_{tot}^z=1`$ sector (filled circles) but for even spins the value for $`H=0.4J`$ is greater than in the $`S_{tot}^z=1`$ sector (this is particularly clear for the middle even spins).

At higher temperature ($`T=0.526J`$), slightly above the gap $`\mathrm{\Delta }0.350.4J`$ (see figure 3), end effects are smaller, especially for low fields. At $`H=0`$, the on-site magnetization is zero by symmetry and is not shown in the figure. Increasing the field, we notice two effects: the first is simply progressive magnetization of the middle spins, unaffected by end effects. When we later discuss the NMR spectra this will result in a positive shift of the NMR spectra for Y atoms close to the middle of the chain segments, i.e. of the central peak. The second effect is the progressive magnetization of the end spins with increasing field. The largest value of the field is equal to the field used in the experiments by Tedoldi et al . Clearly at this field and temperature there is a signal from the $`S=1/2`$ chain end spins. Since the temperature is large, the correlation length is small and the end effects are only visible for the first few spins away from the chain end.

### C Estimation of correlation length

As previously mentioned, the on-site magnetization should decay exponentially away from the end of the chain with a characteristic length related to the bulk correlation length. In order to compare to the experiments, we estimated the correlation length for different temperatures for the longest chains studied in a field of $`H=0.08J`$. We fit the exponential decay of the on-site magnetization away from the end of the chain using the following form:

$`S_i^z=`$ $`S_1^z\left(()^{i1}e^{\frac{i1}{\xi (T)}}+()^{Li}e^{\frac{Li}{\xi (T)}}\right)`$ (3)
$`+S_h^z.`$
Here, $`L`$ is the length of the chain and the last term $`S_h^z`$ in this equation corresponds to the magnetization induced by the field at high temperature. Our results for the correlation length obtained using this method are shown in table I and are in agreement with previous numerical results   and the experimental results taken from NMR experiments .

## IV NMR spectra

In the recent 89Y NMR study by Tedoldi et al , of doped Y2BaNi1-xMgxO5 the Y is sitting roughly at an equal distance between two $`S=1`$ chains. The non-magnetic Mg atoms will tend to break the chains into segments of different lengths and the NMR signal should then be proportional to the sum of the two nearest Ni spins sitting on different segments on different chains. Due to the rather large variation in chain lengths and the resulting variation in the on-site magnetization the satellite peaks should be broadened. Averaging over a large number of chain lengths it is then possible to reconstruct the observed NMR spectra up to an overall rescaling factor.

In order to perform this calculation we have simulated open spin 1 segments in a magnetic field of $`H=0.08J`$, corresponding to the experimental situation, for four temperatures between $`0.2J`$ and $`1J`$ equal to the experimental values, for chain lengths between L=1,200. At zero temperature we have computed the profiles in a slightly different fashion by selecting the $`S_{tot}^z=1`$ subspace for chain lengths between L=1,501 using DMRG techniques. For very short even length chains, which are quite probable, this approach might lead to a slight bias since the magnetic field applied in the experimental situation might not be sufficiently strong to select the $`S_{tot}^z=1`$ state as the ground-state. However, we estimate this effect to be rather small.

Assuming that the distribution of the Mg impurities is uncorrelated we use a geometric distribution for chain lengths where the probability of having a chain of length L is given by:

$$p(L)=x(1x)^L.$$
(4)
Here $`x`$ is the impurity concentration. The mean chain length is then

$$\overline{L}=\underset{l=0}{\overset{\mathrm{}}{}}lx(1x)^l=\frac{1}{x}1,$$
(5)
which, if $`x`$ is sufficiently small, reduces to $`\overline{L}=\frac{1}{x}`$. The non magnetic impurity concentration we used is $`x=5\%`$, which corresponds to a mean chain length of 19. One should note that the geometric distribution is quite wide and even though large chain lengths are less probable than short ones they contribute significantly to the observed NMR spectrum. The most probable chain length is in fact $`L=0`$ corresponding to an impurity site and the possibility of having one or more impurities next to each other is non-negligible. Hence, in order to obtain reliable results we have computed profiles for chains of long lengths. For $`T0`$, we calculated profiles for chains of lengths 1 up to 200, and for $`T=0`$ we simulated chains of up to 500 sites .

According to the experimental situation, we randomly choose two lengths of chains and compute the sum of the two chains magnetizations site by site as we go along the chains. The two sites whose magnetizations are added at each step stand for the two next nearest neighbors spins of the NMR probe in the experiments . When the end of one chain is reached we insert a site with zero on-site magnetization corresponding to the impurity. We then select a new chain following the geometrical distribution and continue to go along the system until the end of another chain is reached and we repeat the process till the thermodynamic limit is reached. If the length of the chain is 0 we immediately insert an impurity site resulting in two or more impurities sitting next to each other. The NMR spectrum is then computed by plotting the histogram of the effective field seen on the Y site, i.e. the sum of the on-site magnetization coming from the two sites sitting closest to the Y on the two chains. This procedure should be quite close the experimental situation. In order to compare to the experimental data we rescale the y and x axis of the experimental NMR spectrum by a scale factor. This is done by adjusting to the first satellite peak (for the x axis we use a scale factor of roughly $`0.80.9`$ which is in good agreement with the experimental estimate of hyperfine interaction).

In figures 4 and 5, we show the results of our calculations for the five temperatures and compare them to the experimental spectra. The central peak corresponds to the situation where the effective field at the Y site cancels either due to the fact that the two contributions are of roughly equal and opposite sign or because the Y is sitting close to the middle of both of the two chain segments where the on-site magnetization is very close to zero. If the Y is closest to the first site of a chain segment it should experience a positive shift away from the central peak and assuming that the magnetization on the other chain is close to zero we should observe the most pronounced satellite peak corresponding to the average value of the on-site magnetization on the first site. If we have exponentially localized excitations at the end of the chains the most probably value of the magnetization on the “other” chain is close to zero. However, this is of course not strictly the case, resulting in a broadening of the peaks. If the Y is close to the second site of one of the chain we should observe a satellite peak shifted in the negative direction away from the central peak but with a smaller absolute shift. It is then possible to identify the various satellite peaks with the on-site magnetization as one goes away from the end of the chain. As the distance from the end of the chain is increased the position of the satellite peak should move closer and closer to the central peak. Due to the significant broadening of the peaks only the first few satellite peaks are resolvable. With exponentially localized excitations at the end of the chains there is also a small probability that the Y is close to firsts sites on both chains (for example first-first, first-third …, or second-second, second-fourth…). Theoretically, this should result in a number of secondary peaks beginning at twice the absolute value of the main peaks, also contributing to a general broadening of the spectrum. However, the intensity of such peaks should be rather small. In our calculations at $`T=0`$, a secondary peak corresponding to the Y close to the firsts sites on both chains is clearly visible at twice the value of the first satellite peak (See Fig. 5).

The computed spectra are in extremely good agreement with the experimental results for the three higher temperatures. Note that the central peak is always greater in the NMR experiments than in our simulations, this might be due to a different distribution of impurities in the sample than the simple geometrical distribution (we don’t consider correlations between impurities in this model). There is a central peak for all the temperatures; the position of this peak is increased as temperature is raised due to the fact that (middle) spins are more easily magnetized by the magnetic field for high temperatures, see section III B. The spectrum is wider at lower temperatures since the magnetization profiles are much more developed, as discussed in the previous section, resulting in larger shifts. For each figure, apart from the central one, several peaks can be observed. As explained, the shift from the central peak to the others correspond to the intensity of the magnetization on neighboring sites. On the figures, we have labelled the different satellite peaks according to the proximity to a site on the chain. As can be clearly seen, the neighboring magnetic structure of a site is more and more well defined as the temperature is lowered. For example, at the highest temperature $`T=0.993J`$, only the two first satellite peaks are resolved, the more distant satellites cannot be distinguished and contribute to the central peak. For $`T=0.526J`$, up to six satellite peaks are resolved, and even more at the lowest temperature, $`T=0.27J`$. This temperature is well below the gap and corresponds to the $`T=77K`$ experimental results in Ref. . In this case the experiment show a single broad line and the signal from satellite peaks is lost. Tedoldi et al  argues that one possible explanation of this is the large distribution of the chain lengths which could become more important at temperatures below the gap. This effect is taken correctly into account in our calculations and since we observe well-defined satellite peaks well below the gap and at $`T=0`$ the explanation is likely elsewhere. The resolution of the $`T=77K`$ spectrum is less than for the other temperatures  and it is possible that higher resolution experiments will be capable of resolving the satellite peaks.

Finally, in figure 5 we show results calculated at $`T=0`$ using DMRG techniques. The calculation is performed in zero field in the $`S_{tot}^z=1`$ subspace. Many well defined peaks are now visible. In particular we see secondary peaks corresponding to the case where the Y is between sites which both have large on-site magnetizations ($`1^{}`$).

## V Discussion

Our results clearly indicate that the $`S=1/2`$ chain excitations should be visible at temperatures well above the gap in a finite field. It is well known that the gap in $`S=1`$ chains corresponds to a single magnon branch at $`k=\pi `$. Around $`k=0`$ a two-magnon continuum should begin at energies of $`2\mathrm{\Delta }`$ and a three-magnon continuum at $`k=\pi `$ should begin at an energy of $`3\mathrm{\Delta }`$ . It is natural to expect that chain end effects are visible until the two-magnon and three-magnon continuum become significantly populated, i.e. for temperatures above 300K in agreement with the numerical and experimental results.

The discrepancy between our theoretical results and the experimental results for temperatures below the gap could have several explanations. As suggested by Tedoldi et al  anisotropy effects could be important. However, the single ion anisotropy term $`D=8K`$ is rather small and it is not obvious that it would give rise to a large effect. Such a term would split the low-lying triplet into a doublet and a singlet and could likely cause a systematic shift in the spectra rather than a complete smearing of the profile. Still another explanation could come from the small coupling between the chains. However, it would seem unlikely that this effect could be important at temperatures as high as $`T=77K`$. Payen et al  argue that such effects become important at temperatures around 4K, much lower than the gap, $`\mathrm{\Delta }100K`$. However, these experiments are concerned with bulk susceptibility and magnetization which might be a less sensitive probe. Another simple explanation is that below the gap so many secondary peaks appear that the experimental spectrum appears as a single broad line. This is not unreasonable considering the results shown for $`T=77K`$ in figure 4. Lastly, the actual impurity concentration is not known to a very high precision and a slight modification in this concentration greatly modifies the observed spectra. High resolution NMR experiments at low temperatures will be needed to resolve these issues.

###### Acknowledgements.

We would like to acknowledge F. Tedoldi and M. Horvatić for many fruitful discussions and for kindly letting us reproduce the experimental data.


text:
# Free cooling and inelastic collapse of granular gases in high dimensions

## Abstract

The connection between granular gases and sticky gases has recently been considered, leading to the conjecture that inelastic collapse is avoided for space dimensions higher than 4. We report Molecular Dynamics simulations of hard inelastic spheres in dimensions 4, 5 and 6. The evolution of the granular medium is monitored throughout the cooling process. The behaviour is found to be very similar to that of a two-dimensional system, with a shearing-like instability of the velocity field and inelastic collapse when collisions are inelastic enough, showing that the connection with sticky gases needs to be revised.

PACS numbers: 47.50.+d, 05.20.Dd, 51.10.+y, 47.11.+j

An important difference between a molecular fluid and a gas of mesoscopic or macroscopic grains is the possibility for the latter, associated with the inelastic nature of the collisions, to exhibit clustering and collapse . A large and rapidly growing body of theoretical work is devoted to clustering, which consists in a long wavelength low frequency hydrodynamic phenomenon, and refers to the formation of density inhomogeneities. On the other hand, the phenomenon of inelastic collapse, which is a short wavelength and high frequency singularity inherent to the inelastic hard sphere (IHS) model, seems much less understood, except in one dimension .

In the IHS model, grains are modeled as smooth hard spheres undergoing binary, inelastic and momentum-conserving collisions, which dissipate a constant fraction $`(1r)`$ of the component of the relative velocity $`𝐯_{12}`$ along the center-to-center direction $`\widehat{𝝈}`$. Noting with primes the post-collision velocities, this translates into $`𝐯_{12}^{}\widehat{𝝈}=r𝐯_{12}\widehat{𝝈}`$, while the tangential relative velocity (perpendicular to $`\widehat{𝝈}`$) is conserved.

In an interesting article, Ben-Naim et al. proposed that a freely evolving inelastic gas belongs asymptotically to the universality class of the sticky gas , for which $`𝐯_{12}^{}=\mathrm{𝟎}`$ after each collision. Noticing that the temperature $`T`$ of an inelastic gas is a monotonically increasing function of the restitution coefficient $`r`$ and therefore bounded from below by the totally inelastic case ($`r=0`$), these authors invoked a mapping onto Burgers’ equation to conjecture that the inelastic collapse is avoided for space dimensions $`d>d_c=4`$, and that the standard Haff’s cooling law $`T(ϵt)^2`$ where $`ϵ=(1r^2)/(2d)`$, holds indefinitely above this critical dimension $`d_c`$. Velocity fluctuations and scaling exponents in one dimension can indeed be described by the inviscid Burgers equation . However, in higher dimensions, the completely inelastic version of the IHS model, with $`r=0`$, does not strictly correspond to the sticky gas limit because the tangential relative velocity is not dissipated in a binary encounter ($`𝐯_{12}^{}\widehat{𝝈}=\mathrm{𝟎}`$, but a priori $`𝐯_{12}^{}\mathrm{𝟎}`$), so that it is interesting to test the validity of the above-mentioned predictions. In this article, we report Molecular Dynamics (MD) simulations of IHS gases for $`d=4,5`$ and 6. For each space dimension, the relevant parameters $`\varphi `$ (packing fraction) and $`r`$ (normal restitution) are varied for systems consisting typically of $`N=5.10^3`$ to $`5.10^4`$ particles. The code was first successfully tested by comparing for various densities the MD equation of state (or equivalently the pair correlation function at contact) with the analytical approximation of Song et al . For all investigated dimensions and for high enough dissipation ($`r0.2`$ for $`d=5`$ and $`r0.1`$ for $`d=6`$), the system exhibits the finite time singularity characteristic of the inelastic collapse, in contradistinction to the conjecture of , with a situation closely reminiscent to its two dimensional counterpart : the (hyper)spheres collide infinitely often in a finite time along their joint line of centers. Following Refs.  we probed this multi-particle process occurring through the accumulation of an infinite sequence of binary collisions by a contact criterion: after each collision, the relative distance $`d^{}`$ between the next two colliding partners is monitored; if this interparticle spacing normalized by the diameter $`\sigma `$ has decreased and becomes of the order of machine precision, a three body interaction has occurred, corresponding to an inelastic collapse. The results of a typical run ($`d=5`$) are shown in Fig. 1. When a multi-body interaction commences, $`d^{}`$ decreases geometrically with the number of collisions, as in two dimensions . Enforcing a high precision computation allows to follow the decay over more than 26 orders of magnitude (whereas only 8 orders are accessible with a standard double precision algorithm, see the inset of Fig. 1). After a collapse has occurred, the inaccuracy of the computer disperses the collapsing cluster, before another multi-body event involving different particles occurs at a different location. Our analysis indicates that on a hypothetical infinite precision machine, the collapse would continue forever whereas roundoff errors act as an effective regularization. Throughout a collapse, the time between two successive collisions follows a geometrical decrease very similar to that displayed in Fig. 1.

Let us note that the seemingly low value of the packing fraction $`\varphi `$ in Fig. 1 is a misleading effect of “high” dimensionality. It turns that the reduced density

$$n^{}=n\sigma ^d=\frac{d\mathrm{\hspace{0.17em}2}^{d1}}{\pi ^{d/2}}\mathrm{\Gamma }\left(\frac{d}{2}\right)\varphi $$
(1)
where $`n`$ is the number density and $`\mathrm{\Gamma }`$ the Euler function, is a more relevant parameter to discriminate between “low” and “high” densities. In a simple cubic lattice, the highest packing achievable with spheres at contact corresponds to $`n^{}=1`$. The configuration of Fig. 1 with $`\varphi =0.08`$ corresponds to $`n^{}=0.5`$ and is consequently of high density. We considered the possibility that the conjecture of  applies in the opposite limit of low packing (that would correspond to the so called sticky dust in the case of vanishing both normal and tangential restitution coefficients). We lowered $`\varphi `$ by an order of magnitude which requires to consider large systems to avoid a spurious increase of the mean free path; our simulations with $`N=2.10^5`$ particles and $`\varphi =0.008`$ ($`n^{}=0.05`$) in five dimensions nevertheless exhibit inelastic collapse for $`r<0.05`$.

When the restitution coefficient is above a critical threshold, the inelastic singularity is absent and we can follow the time evolution of the system. Initial conditions for all inelastic runs were equilibrated fluid configurations of elastic hard spheres at the corresponding (uniform) density, with a uniform temperature (coarse-grained kinetic energy) and Gaussian distribution of velocities. The short time regime is then given by the homogeneous cooling state, which is essentially an adiabatically changing equilibrium state where the mean kinetic energy $`E`$ (related to the temperature $`T`$ by $`E=dT/2`$ due to the vanishing of the flow field) is given by

$$E(t)=\frac{E_0}{\left(1+ϵt/t_0\right)^2}.$$
(2)
In Eq. (2), $`t_0`$ is the Enskog mean collision time of elastic particles at the same density  and the inelasticity parameter $`ϵ=(1r^2)/(2d)`$ follows from the assumption of a Gaussian velocity distribution. It has been shown within the framework of Enskog-Boltzmann equation that the corrections due to non Gaussian behaviour on the cooling rate are small for all inelasticities . The internal clock of the system can be parameterized by $`\tau `$, the average number of collisions suffered per particle in a time $`t`$, which reads in the homogeneous cooling state

$$\tau =\frac{1}{ϵ}\mathrm{ln}\left(1+ϵ\frac{t}{t_0}\right).$$
(3)
Denoting $`N_c`$ the total number of collisions having occurred in the system over a time $`t`$, we have $`\tau =2N_c/N`$. Figure 2 shows that the law $`E(\tau )=E_0\mathrm{exp}(2ϵ\tau )`$ expected from Eqs. (2) and (3) is only valid at short times. A crossover, indicated in Fig. 2 by an arrow, is generically observed for all runs in dimensions 5 and 6, provided the inelastic collapse is avoided; the crossover time depends both on the inelasticity and on system sizes. After the crossover, the behaviour of $`E`$ with the number of collisions also depends on system size (see Fig. 3), and $`[\mathrm{ln}(E/E_0)]/\tau `$ does not scale with the inelasticity as $`ϵ`$ (inset of Fig. 3). Moreover, Fig. 4 shows that the energy decays like $`t^2`$, as predicted in ; however, Fig. 4 shows that the prefactor is different at short and large timesWe note that from Eq. (2), the slope of the dotted line in Fig. 4 (short time behaviour) corresponds to the ratio $`t_b/t_0`$ of Boltzmann over Enskog mean collision times, and is an indirect measure of the pair correlation function at contact $`\chi `$ . From Fig. 4, we get $`\chi 1.7`$ while the equation of state of Ref.  yields 1.74 for the density considered. This value is corroborated by an independent elastic run where we measure $`\chi 1.7`$., and (at large times) strongly increases with increasing number of particles; it neither scales with inelasticity as $`ϵ^2`$ (not shown). All these results are at variance with the suggestions of : the energy at large times decays as $`A(N,ϵ)ϵ^2t^2`$, and not simply as $`ϵ^2t^2`$ (Haff’s law).

Before reaching the crossover time $`t_c`$, the velocity distribution of the particles is isotropic and very close to a Maxwellian. Once the crossover time is elapsed, we observe an evolution reminiscent of the shearing instability found in lower dimensions . Fig. 5 shows the distribution of the components of the rescaled velocity $`𝐜=𝐯/\sqrt{T(t)}`$ for $`tt_c`$. For the particular system of Fig. 5, two components are found to be Gaussian (the central peaks) and the remaining directions appear as bimodal and concentrate most of the energy. The corresponding shear bands are difficult to visualize, but insight can be gained by a suitable projection onto a two-dimensional surface. Such a projection is displayed in Fig. 6. The coherent motion we observe in Fig. 6 suggests that the large time dynamics of the system is controlled by the periodic boundary conditions used in the simulations, as analyzed in .

Since, in a 5 dimensional system of $`5.10^4`$ particles, with packing fraction 0.08, the simulation box length $`L`$ is only 10 times as large as the spheres’ diameter, the systems considered here can be considered as small (in the sense that there exists no possible separation of length scales). The behaviour reported here for $`E(t)`$ is therefore consistent with the two dimensional simulations of Orza et al. , showing $`EA(N,ϵ)t^2`$ for small systems, with a prefactor much larger than predicted by Haff’s law. Gaining one order of magnitude for $`L/\sigma `$ would require to increase $`N`$ by $`d`$ orders of magnitude, which is unachievable for $`d>4`$. The validity of Brito and Ernst’s theoretical prediction for large systems $`E\tau ^{d/2}`$  (coinciding with the approach of  for $`d4`$) can consequently not be tested in dimensions higher than 4.

In conclusion, our simulations show the existence of an inelastic collapse in high dimensions, not only for dense systems but also for more dilute ones (even if the situation of very low packing fraction cannot be reached, so that the “sticky dust” limit cannot be tested by numerical simulations). They indicate therefore that the relationship between inelastic and sticky gases put forward in  needs to be refined.

It is a pleasure to thank M. Ernst for his continued interest in our work.