text:
# Near–infrared imaging of the host galaxies of intermediate redshift steep spectrum radio quasarsBased on observations collected at the European Southern Observatory, La Silla, Chile.

## 1 Introduction

Detailed studies of orientation-independent properties of quasars, such as their host galaxies and environments, are a fundamental tool to understand the quasar phenomenon, and AGN in general. Such studies can crucially test models which attempt to unify different observed classes of AGN either through evolution (e.g. Ellingson, Yee & Green 1991) or orientation (e.g. Antonucci 1993; Urry & Padovani 1995). They can shed light on whether interactions (as evidenced by close companions or disturbed morphology) are important for triggering and fueling of the quasar activity (e.g. Smith & Heckman 1990; Hutchings & Neff 1992). These studies allow one to investigate the possible reciprocal effect of the AGN on the properties and evolution of the hosts. Possible cosmological evolution of the AGN population (e.g. Small & Blandford 1992; Silk & Rees 1998) can be constrained by comparing the properties of the hosts and environments of AGN at different redshifts.

The host galaxies of low redshift (z $`\stackrel{<}{}`$0.3) quasars have been thoroughly investigated using both ground-based imaging (e.g. McLeod & Rieke 1995; Taylor et al. 1996; Percival et al. 2000) and the Hubble Space Telescope (HST; e.g. Bahcall et al. 1997; Boyce et al. 1998; McLure et al. 1999). These investigations have shown that most quasars live in galaxies at least as bright as the Shechter function’s characteristic luminosity L (e.g. Mobasher, Sharples & Ellis 1993). While most quasar host galaxies are brighter than L, and in many cases comparable to brightest cluster member galaxies (BCM; Thuan & Puschell 1989), some undetected or marginally detected hosts may be under-luminous. Recent morphological studies of host properties at low $`z`$ (e.g. Taylor et al. 1996; Percival et al. 2000) have concluded that while radio-loud quasars (RLQ) are found exclusively in giant elliptical galaxies, radio-quiet quasars (RQQ) reside in both elliptical and spiral (disc dominated) galaxies. It has been suggested that the morphological type may depend on the power of the quasar, with the most luminous quasars found only in spheroidal host galaxies (Taylor et al. 1996). Close companions have been found around many quasars (e.g. Hutchings 1995; Bahcall et al. 1997; Hutchings et al. 1999), some with signs of interaction, but the physical association has been confirmed only in some cases through spectroscopic measurements (e.g. Heckman et al. 1984; Ellingson et al. 1994; Canalizo & Stockman 1997).

Our knowledge of the host properties is basically limited to low redshift, because of the increasing difficulty of both resolving the quasar and the rapid cosmological dimming ($`\alpha `$ (1+z)4) of the host in contrast to the nucleus. At intermediate redshift (0.5 $`<`$ z $`<`$ 1), quasar hosts have been little investigated (e.g. Carballo et al. 1998; Rönnback et al. 1996). Recently, we studied in the near-infrared (NIR) the host galaxies of a complete sample of 20 flat spectrum radio quasars (FSRQ) (most of which are at 0.5 $`<`$ z $`<`$ 1) (Kotilainen, Falomo & Scarpa 1998; hereafter KFS98). We were able to resolve the host galaxy in a substantial fraction ($``$50 %) of them and found these highly polarized quasars to reside in galaxies of M(H)$``$–27.

While most FSRQs are characterized by rapid variability, high and variable polarization, high brightness temperatures, core–dominated radio emission and confirmed or suspected superluminal motion (e.g. Impey & Tapia 1990; Quirrenbach et al. 1992; Padovani & Urry 1992; Vermeulen & Cohen 1994), very few of these characteristics are shared by the more common steep spectrum radio quasars (SSRQ). This difference is usually explained in the unified model in terms of synchrotron radiation strongly relativistically beamed close to our line–of–sight in FSRQs, while SSRQs are viewed further away from the beaming axis (Blandford & Rees 1978).

We present here deep high spatial resolution ($``$0$`\stackrel{}{.}`$3 px-1, $``$ 1′′ FWHM seeing) NIR $`H`$-band (1.65 $`\mu `$m) imaging study of the host galaxies of a sample of SSRQs in the redshift range 0.5 $`<`$ z $`<`$ 1. We thus explore the host galaxies in the rest-frame 0.8 – 1.1 $`\mu `$m, offering many advantages over the optical wavelengths, including larger host/nucleus contrast, negligible scattered light from the quasar, small K-correction and insignificant extinction (see KFS98 for details). For practically all the quasars in the sample we present the first high quality NIR imaging observations, and the first detection of the host galaxy.

The SSRQ sample was chosen from the Véron-Cetty & Véron (1998) AGN catalogue to match the FSRQ sample (KFS98) in terms of redshift and optical and radio luminosity distribution, while keeping the two samples separated regarding the radio spectral index. The selection criteria applied were: 0.5 $`<`$ z $`<`$ 1.0, $`\delta `$ $`<`$ +20°, MB $`<`$ -24.0, $`\alpha `$(6-11 cm) $`>`$ 0, and L(11 cm) $`>`$ 1044 erg sec-1. This selection yielded a total of 54 SSRQs. Out of these, 33 were observable during the observing runs. A total of 19 SSRQs out of these were observed, and their general properties are given in Table 1. The properties of the SSRQ and FSRQ samples are compared in Fig. 1 and Table 2. The MB distribution is very similar for the two samples. The FSRQ sample has more high redshift quasars, but the average value is in agreement with the SSRQs. The samples are reasonably well separated in $`\alpha `$, although some overlap exists. Finally, there is a marked difference in the L(11 cm) distribution, in the sense that the FSRQ distribution extends to much higher radio luminosities than the SSRQ sample. This is due to the unfortunate circumstance that there were not enough very radio-luminous SSRQs available at the time of the observations. This caveat should be kept in mind in the following discussion.

The outline of the paper is as follows. In section 2, we briefly describe the observations, data reduction and the method of the analysis. Our results concerning the host galaxies and the nuclear component are given in section 3 and compared with the results for FSRQs. In section 4, we analyze the close environment of SSRQs and FSRQs and, finally, in section 5 summary and main conclusions of this study are drawn. In the Appendix, we give comments for individual quasars. Throughout this paper, H0 = 50 km s-1 Mpc-1 and q0 = 0 are used.

## 2 Observations, data reduction and modeling of the luminosity profiles

The observations were carried out during two observing runs. The first was at the 2.2m telescope of European Southern Observatory (ESO), La Silla, Chile, between 30 December 1997 and 3 January 1998, using the 256x256 px IRAC2 NIR camera (Moorwood et al. 1992) with pixel scale 0$`\stackrel{}{.}`$278 px-1, giving a field of view of 1.2 arcmin2. The second run was at the ESO 3.5m New Technology Telescope (NTT) on the nights of 31 January and 1 February 1999, using the 1024x1024 px SOFI NIR camera (Lidman & Cuby 1998), with pixel scale 0$`\stackrel{}{.}`$292 px-1, giving a field of view of 5.0 arcmin2. Details of the observations and NIR photometry are given in Table 1. The observations, data reduction and modeling of the luminosity profiles were performed following the procedure described in KFS98. Briefly, targets were always kept in the field by shifting them in a grid across the array, with typical offsets of 30′′. Individual exposures of 60 sec duration were coadded to achieve the final integration time reported in Table 1.

Data reduction, performed using IRAF111IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation., consisted of correction for bad pixels by interpolating across neighboring pixels, sky subtraction using a median averaged, scaled sky frame from all observations of a given target, flat-fielding using illuminated dome exposures, and combination of all images of the same target, using field stars or the centroid of the light distribution of the quasar as a reference point. Standard stars from Landolt (1992) were used for photometric calibration, for which we estimate an accuracy of $``$0.1 mag, based on comparison of multiple observations at different nights with IRAC2 as listed in Table 1. K–correction was applied to the host galaxy magnitudes following the method of Poggianti (1997). The applied K–correction for each quasar is reported in Table 3. This correction, insignificant at low redshift (m(H) = 0.03 at z = 0.2), has to be taken into account at moderate redshifts (m(H)$``$0.14 at our median redshift of z = 0.65). No K–correction was applied to the nuclear component.

After masking all the regions around the target contaminated by companions, azimuthally averaged radial luminosity profiles were extracted for each quasar and field stars out to the background noise level, corresponding to a surface brightness of $`\mu `$(H) = 23–24 mag arcsec-2. In the small field of view IRAC2 images, bright stars were present in the field of only a few quasars. For most quasars, the core of the PSF was derived from faint field stars, while the wing was derived from separate observation of a nearby bright star or extrapolated using a suitable Moffat (1969) function obtained from fitting bright stars in other frames of similar seeing during the same night. In the much larger field–of–view SOFI images, all quasars had bright stars in the field, enabling a good PSF determination.

The luminosity profiles were fitted into a point source (PSF) and a galaxy (de Vaucouleurs law, convolved with the PSF) component by an iterative least–squares fit to the observed profile. For these moderately high redshift quasars, data quality does not allow one to discriminate between elliptical and disk morphologies for the host galaxy, and the elliptical model was assumed, as demonstrated to be the case for low redshift RLQs (e.g. Taylor et al. 1996; Bahcall et al. 1997). For quasars with no host galaxy detection, we determined an upper limit to the brightness of the host galaxy by adding simulated “host galaxies” of various brightness to the observed profile until the simulated host became detectable within the errors of the luminosity profile. For a detailed discussion of potential error sources in this fitting procedure, see KFS98. We estimate the uncertainty of the derived host galaxy magnitudes to be $``$$`\pm `$0.3 mag for the clearly detected hosts. For the marginally detected hosts, we can only assess a lower limit to the error margin as $``$$`\pm `$0.5.

For the clearly detected hosts, there is expectedly a good agreement between the total magnitude of the quasars derived from aperture photometry (Table 1, column 11) and that resulting from the profile fit (Table 3, columns 5 and 6), with average magnitude difference of –0.03$`\pm `$0.12.

## 3 Results and discussion

In Fig. 2, we show the $`H`$-band contour plots of all the SSRQs, after smoothing the images with a Gaussian filter of $`\sigma `$ = 1 px. In Fig. 3 we show the radial luminosity profiles of each quasar, with the best–fit decomposition overlaid. We detect the host galaxy clearly for 10 (53 %) SSRQs and marginally for 6 (32 %) more. The host remains unresolved for 3 (16 %) SSRQs. We summarize the best–fit model parameters of the profile fitting and the derived properties of the host galaxies in Table 3, where column (1) gives the name of the SSRQ, (2) the redshift, (3) the applied K–correction, (4) the scale-length of the host galaxy, (5) – (6) the apparent magnitude of the nucleus and the host, (7) the nucleus/host luminosity ratio, (8) – (9) the absolute magnitude of the nucleus and the host, and column (10) note about the detection of the host.

Table 4 presents a comparison of the $`H`$-band absolute magnitudes of the SSRQ hosts with relevant samples from previous studies in the literature, for which we report the average values after correcting the published values for colour term and to our cosmology. In the Appendix, we compare our NIR photometry with the few previous studies, and discuss in more detail individual quasars, including comparison with the few previous optical/NIR determinations of the host galaxies.

### 3.1 The host galaxies

In Fig. 4 (upper panel) we investigate the location of the SSRQ and FSRQ hosts in the $`H`$-band absolute magnitude vs. redshift diagram. In Fig. 4 (lower panel) we show the same diagram for the mean value of various samples of RLQ and RQQ hosts from the literature. The average $`H`$-band absolute magnitude of the resolved SSRQ host galaxies is M(H) = –27.2$`\pm `$1.1 and the average bulge scale-length R(e) = 9.0$`\pm `$1.7 kpc, while the values after adding the marginally resolved hosts are M(H) = –27.0$`\pm `$1.2 and R(e) = 8.6$`\pm `$1.9 kpc. The SSRQ hosts are therefore large (all have R(e) $`>`$ 3 kpc, the empirical upper boundary found for normal local ellipticals by Capaccioli, Caon & D’Onofrio 1992), and very luminous, much brighter than the luminosity of an L galaxy, which has M(H) = –25.0$`\pm `$0.3 (Mobasher et al. 1993). It is therefore evident that the SSRQ hosts are preferentially selected from the high–luminosity tail of the galaxy luminosity function (the derived upper limits for the unresolved hosts are also consistent with this; see Table 3). Indeed, we find no case of an SSRQ host with M(H) $`>`$ –25, indicating that for some reason these quasars cannot be hosted by a galaxy with L $`<`$ L, similarly to what was found by Taylor et al. (1996) for low redshift RLQs. Although there is a large spread, the SSRQ hosts are on average of similar luminosity to those of FSRQ hosts (KFS98).

There is a suggestion of a positive correlation of host luminosity with redshift (Fig. 4, lower panel) as the FSRQ and SSRQ hosts fall between the luminosities of lower redshift (M(H) $``$–26) and higher redshift (M(H) $``$–29) RLQs. This is consistent with what is expected from passive stellar evolution models for elliptical galaxies (e.g. Bressan et al. 1994; Fukugita, Shimasaku & Ichikawa 1995), and the evolution of galaxies in clusters (Ellingson et al. 1991), although the scatter is quite large. On the other hand, it is not consistent with models of hierarchical galaxy formation (e.g. Kaufmann & Hähnelt 2000), which predict $``$L hosts at z = 2 – 3 that afterwards undergo mergers to become present-day giant elliptical hosts. The detected correlation, therefore, suggests that major merging has already happened at z = 2 - 3. Also, this is consistent with the discovery of ellipticals with old stellar populations even at z$``$2 (e.g. Spinrad et al. 1997; Stiavelli et al. 1999). However, this scenario seems not to be valid for RQQs, which show no or little evolution with redshift, as evidenced by the non-detection of high redshift RQQ hosts (Lowenthal et al. 1995).

On the other hand, the detected correlation can be interpreted as evidence for a relationship between optical and radio luminosities. Ledlow & Owen (1996) found such a correlation for radio galaxies and in the unified model (e.g. Urry & Padovani 1995), where RLQs and radio galaxies are identical objects seen from different viewing angles, this relationship is expected also to exist between the radio luminosities and host galaxy absolute magnitudes for RLQs at different redshifts.

Most of the available comparison data from literature are for low and intermediate redshift RLQs. These samples span a moderately large range in redshift from z$``$0.1 up to z$``$0.7, which is the average redshift of the SSRQ sample. The average host galaxy magnitudes for the various samples are given in Table 4 in order of increasing average redshift. Considering all these samples together gives average host magnitude of M(H) = –26.2$`\pm `$0.7. As can be seen from Fig. 4, there is no significant difference among the average values of these low z samples. The average M(H) = –27.0$`\pm `$1.2 for the SSRQ hosts is therefore $``$1 magnitude brighter than the luminosity of low z RLQ which suggests evolution in the host brightness with redshift, and/or a relationship of the host luminosity with the nuclear luminosity (see section 3.2).

It is worth to compare these data also with those for samples of low redshift RQQs extracted from the literature. The average host galaxy magnitudes for the various samples are given in Table 4. Considering all these RQQ samples together yields an average host magnitude of M(H) = –25.4$`\pm `$0.4, $``$1 magnitude fainter than comparable redshift RLQ hosts and $``$1.5 mag fainter than SSRQ hosts at intermediate redshift, indicating that RLQs and RQQs inhabit different types of galaxies. However, this difference may also be due to the above mentioned correlation between optical and radio luminosities (Ledlow & Owen 1996), in the sense that RQQ surveys select objects according to the optical luminosity, while RLQ surveys select objects with high radio luminosity and, therefore, a bright host galaxy.

Lehnert et al. (1992, 1999) have reported spatially resolved structures around RLQs at z$``$1.5 and z$``$2.3 that, if interpreted as host galaxies, would correspond to extremely luminous galaxies with average M(H) = –29.1$`\pm `$1.1, $``$2 mag brighter than the FSRQs and SSRQs at z$``$0.7. However, within the scatter involved in these numbers, our results appear to be consistent with those of Lehnert et al., for the trend between the nuclear and host galaxy luminosities (see section 3.2), and is supporting evidence for the existence of a real upturn in the host luminosity occurring between z$``$0.5 and z$``$2, leading from L$``$L\* hosts at low redshift to the host galaxies of high redshift quasars that are several magnitudes brighter than L\* (see Fig. 4). While this type of change is consistent with evolution of the stellar population in the elliptical hosts (as argued for high redshift RGs by Lilly & Longair 1984), or being intrinsic AGN luminosity effect (as argued for high redshift RGs by Eales et al. 1997), there are many caveats in this comparison, most notably differences in the intrinsic quasar luminosity of the various samples. In addition, optical and NIR imaging by Lowenthal et al. (1995) failed to detect extended emission in a sample of six RQQs at z$``$2.3. Their upper limits (average m(H) $`>`$ 18.8$`\pm `$0.4 and M(H) $`<`$ -28.4$`\pm `$0.6) indicate that the RQQ hosts at high redshift must be $``$3 mag brighter than L\* and $``$1 mag fainter than those found in the Lehnert et al. (1992) sample of RLQs at similar redshift and with roughly similar nuclear luminosity, again suggesting that RLQs and RQQs are different types of objects.

It is well established that elliptical galaxies form families of homologous systems with characteristic parameters R(e), $`\mu `$(e) and velocity dispersion $`\sigma `$. These are commonly represented in the Fundamental Plane (e.g. Djorgovski & Davis 1987). These relations have been claimed to be related to the morphological and dynamical structure of the galaxies, and to their formation process. We have investigated the properties of the SSRQ and FSRQ hosts in the projected Fundamental Plane (F-P) concerning the effective surface brightness $`\mu `$(e) and the effective radius R(e). Surface brightness $`\mu `$(e) was corrected for Galactic extinction, K–correction and for the (1+z)4 cosmological dimming. Fig. 5 shows the correlation between $`\mu `$(e) and log R(e) for the SSRQ and FSRQ hosts. It can be seen that the behavior of the SSRQ and FSRQ hosts are similar. Both follow the Kormendy (1977) relation for giant massive ellipticals (e.g. Capaccioli et al. 1992). No host galaxy is in the (scatter) area at log R(e) $`<`$ 0.5 kpc. This confirms that the SSRQ and FSRQ hosts are exclusively drawn from the population of giant ellipticals.

### 3.2 The nuclear component

The average absolute magnitude of the nuclear component for all SSRQs is M(H) = –28.2$`\pm `$1.2, that is $``$1.5 mag brighter than low redshift RLQ nuclei (average M(H) = –26.8$`\pm `$0.9) and $``$2 mag brighter than low redshift RQQ nuclei (average M(H) = –26.4$`\pm `$1.0). This is likely due to a selection effect present in the original samples induced by the different average redshift. On the other hand, the SSRQ nuclei are $``$1.5 mag fainter than the nuclei of FSRQs (KFS98; M(H) = –29.7$`\pm `$0.8) in the similar redshift range. This difference is even more evident when considering the nucleus/galaxy (LN/LG) luminosity ratio, shown in Fig. 6. The average LN/LG ratios are 3.8$`\pm `$3.2 and 21$`\pm `$11 for the SSRQs and FSRQs, respectively. The large range in the luminosity ratio can be due to differences in the intrinsic nuclear or host luminosity, or a difference in the beaming factor from one object to another. Whereas the majority of the FSRQs have LN/LG $`>`$ 10 in the $`H`$-band, only one SSRQ is above this limit, and the distribution of the LN/LG ratio for the SSRQs is similar to that of the low redshift RLQs (e.g. Taylor et al. 1996).

From Fig. 4 it appears that the host galaxies of the various samples are not dramatically different in intrinsic luminosity, especially if stellar evolution in the elliptical host galaxies is taken into account. Therefore, Fig. 6 clearly indicates that SSRQs exhibit a nuclear component which is systematically brighter than that of lower redshift RLQs and RQQs, but fainter than that of the FSRQs. This is consistent with the beaming model with larger Doppler amplification factor for FSRQs than SSRQs, that makes the observed differences understandable.

In Fig. 7, we show the relation between the luminosities of the nucleus and the host galaxy for the SSRQs and FSRQs. We confirm for the SSRQs the tendency noted for FSRQs by KFS98 for the more powerful quasars to reside in more luminous hosts, in the sense that there appears to be a lower limit to the host luminosity, which increases with the quasar luminosity. All the SSRQs and FSRQs appear consistent with this boundary limit, first proposed by McLeod & Rieke 1995).

Recent dynamical investigations of nearby galaxies have reported evidence for supermassive black holes in their nuclei (Kormendy & Richstone 1995; Magorrian et al. 1998). Assuming that quasar activity results from accretion of material onto a supermassive black hole, the nucleus-host relationship is in agreement with that found by Magorrian et al. (1998) from HST kinematic study between the mass (luminosity) of the black hole and the mass (luminosity) of the spheroidal component in nearby galaxies. Our results for SSRQs and FSRQs therefore suggest that the Magorrian et al. relationship extends to host galaxy masses at cosmological distances. On the other hand, the much weaker correlation found for more nearby, lower luminosity AGN (e.g. McLeod & Rieke 1995; McLeod, Rieke & Storrie-Lombardi 1999) may indicate that the onset of the correlation occurs only after a certain level in nuclear and/or galaxy luminosity has been reached. For example, at highest luminosities (highest black hole masses) quasars may accrete near the Eddington limit, triggering the onset of the host - black hole correlation (Lacy et al. 2000).

However, we advise caution about obvious selection effects that may de-populate the upper left hand (faint nuclei are difficult to detect against luminous hosts) and lower right hand (faint hosts are difficult to detect under luminous nuclei) corners of the diagram. Moreover, it is not clear if this boundary limit depends on the radio loudness since a recent investigation of luminous RQQs at z $``$ 0.3 (Percival et al. 2000) does not find any correlation between host and nuclear luminosities.

## 4 The close environment of SSRQs and FSRQs

It has been noted in many previous studies of quasar environments (e.g. Yee & Green 1984; Stockton & MacKenty 1987; Hutchings & Neff 1990; Hutchings 1995; Hutchings et al. 1999) that quasars have often companions in their immediate environments and they are sometimes associated with disturbed morphology. Spectroscopy of these companions (e,g. Stockton 1978; Heckman et al. 1984) has shown that the companions are at the redshift of the quasar and are therefore physically associated. These observations have sustained the potentially important idea that nuclear activity can be triggered and/or fueled by strong tidal interactions and/or galaxy mergers (e.g. Heckman 1990)

Although both RLQs and RQQs seem to occur in dense groups of galaxies (and only rarely are located in rich galaxy clusters), a systematic difference of environment is found between RLQs and RQQs in the sense that the latter are generally found in poorer fields (e.g. Ellingson et al. 1991). There have also been some suggestions (Stockton 1978; Hutchings & Neff 1990) that FSRQs (and in general, core-dominated RLQs) have a smaller frequency of companions with respect to the more common SSRQs. These differences, if confirmed, are important for understanding the connection between the characteristics of the nuclear activity and the environment.

To investigate these differences, we have compared the frequency of companions found around the matched samples of SSRQs and FSRQs. For both samples, we have counted the number of resolved companion objects (i.e. galaxies) within a specified radius and brighter than m(H) + 2, where m(H) is the apparent magnitude corresponding to M(H) = –25 at the redshift of each quasar. This limit is in all frames at least 1 mag brighter than the magnitude limit of the images. To evaluate the number of companion galaxies around the quasars, we have chosen two radii corresponding to projected distances of 50 and 100 kpc. At the average redshift of $``$0.7 of the samples, 1 arcsec corresponds to $``$9.4 kpc. In Table 5 we give for each quasar the number of resolved companions N(r).

The average number of companions around FSRQs is 0.19$`\pm `$0.39 and 1.19$`\pm `$1.01 within 50 and 100 kpc projected distance from the quasar, respectively. The corresponding values for the SSRQs are 0.05$`\pm `$0.22 and 0.53$`\pm `$0.68. The fraction of FSRQs with no close companion is 81% and 25% within 50 and 100 kpc from the quasar, respectively, while the corresponding fractions for the SSRQs are 95% and 58%. Contrary to previous suggestions (Stockton 1978; Hutchings & Neff 1990), we find that the immediate environment of FSRQs is comparably rich (and possibly even richer) than that of SSRQs. This argues against a possible tight connection between the formation/evolution of the radio structure and the quasar galaxy environment.

## 5 Conclusions

We have presented near–infrared $`H`$-band images of 19 SSRQs in the redshift range 0.5 $`<`$ z $`<`$ 1 with the aim of studying the quasar environment (host galaxies and companions) in comparison with our previously studied matched sample of FSRQs with respect to redshift and optical and radio luminosity. We are able to clearly detect the host galaxy in 10 (53 %) SSRQs and marginally in 6 (32 %) others, while the host remains unresolved in 3 (16 %) cases. The galaxies hosting the SSRQs are large (average bulge scale-length R(e) = 9.0$`\pm `$1.7 kpc) and luminous (average M(H) = -27.2$`\pm `$1.1). They are, therefore, $``$2 mag more luminous than the typical galaxy luminosity L\* (M\*(H) = –25.0$`\pm `$0.2), and $``$1 mag more luminous than the brightest cluster galaxies (M(H) = –26.3$`\pm `$0.3). The SSRQ hosts appear to have similar luminosity as those of the FSRQ hosts (M(H)$``$–27). On the other hand, the average nucleus–to–galaxy luminosity ratio of SSRQs (LN/LG = 3.8$`\pm `$3.2) is significantly smaller than that found for the FSRQs (LN/LG = 21$`\pm `$11), in agreement with the idea that the latter are beamed. We find for SSRQs a positive trend, noted earlier for the FSRQs, between the host and nuclear luminosity. Finally, we find that the number of close companion galaxies is not strongly influenced by the characteristics of radio emission. Comparable (and possibly larger in FSRQ) densities of companions are observed for FSRQs and SSRQs,

###### Acknowledgements.

This work was partly supported by the Italian Ministry for University and Research (MURST) under grants Cofin 98-02-32 and Cofin 98-02-15. This research has made use of the NASA/IPAC Extragalactic Database (NED), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.
## Appendix: Notes on individual quasars and comparison with previous NIR photometry.

PHL 923 = PKS 0056-00. The host remains unresolved with upper limit MH $`>`$ -27.0, assuming R(e) = 10 kpc.

3C 57 = PKS 0159-11. Hyland & Allen (1982) measured H = 14.78 in a 7 ′′ aperture, fainter than our H = 14.1 – 14.2. The profile fit is rather poor with PSF only ($`\chi ^2`$ = 3.7). Using a PSF + elliptical model results in a better fit ($`\chi ^2`$ = 1.0), but there remains excess emission with respect to the first 2 pixels, indicating a bad match with the PSF. The host is thus only marginally resolved with R(e) = 8.8 kpc and MH = –25.4.

3C 95 = PKS 0349-14. Neugebauer et al. (1979) measured H = 14.73 in a 7.5′′ aperture, fainter than our H = 12.9.

3C 110 = PKS 0414-06. Hyland & Allen (1982) measured H = 14.71, while Sun & Malkan (1989) measured H = 14.89, slightly fainter than our H = 14.6. This quasar was studied in the $`R`$-band by Rönnback et al. (1996) who found it to be unresolved, with M(R) $`>`$ -20.9, whereas we resolved the host with M(H) = –26.6.

3C 138 = PKS 0518+16. The host galaxy is marginally resolved, with MH = -25.4 and Re = 8 kpc. This quasar was imaged with HST in the optical by De Vries et al. (1997), who claim it to be slightly extended (although they did not attempt any modelling).

3C 175 = PKS 0710+11. Simpson & Rawlings (2000) measured H = 14.72$`\pm `$0.01 in a 3 arcsec aperture, fainter than our H = 13.1. The host is marginally resolved with MH = –28.7 and re = 8 kpc. This quasar has not been resolved by previous ground-based optical studies (Malkan 1984; Hes, Barthel & Fosbury 1996) but is marginally resolved with HST (Lehnert et al. 1999), showing roughly E-W orientation and a plume to S-SE up to 2 arcsec from the nucleus.

PKS 0825-202. The host remains unresolved with M(host) $`>`$ –25.0.

3C 207 = PKS 0838+13. Simpson & Rawlings (2000) measured H = 16.11$`\pm `$0.04 in a 3 arcsec aperture, in agreement with our H = 15.9. This quasar was not resolved by HST (Lehnert et al. 1999).

PKS 0903-57. The profile fit is poor, resulting in a marginally resolved host with MH = -27.3 and Re = 4.3 kpc.

PKS 0959-443. The profile fit is poor with PSF only. The quasar is marginally resolved with a reasonably good PSF + elliptical fit, yielding MH = -27.4 and Re = 11.1 kpc. Note that there is a bright star close to the quasar, and although it was removed prior to the fit, the possible detection of a host galaxy needs to be approached with caution.

PKS 1116-46. This quasar was studied in the $`R`$-band by Rönnback et al. (1996), who found it to be marginally resolved with M(R) = –23.0.

PKS 1136-138. The rather short exposure time for this quasar resulted in limiting surface brightness of only $`\mu _{lim}`$ = 21. The quasar remains unresolved with a poor fit of the PSF in the center.

PKS 1327-21. Glass (1981) measured K = 14.28$`\pm `$0.24 in a 12′′ aperture, brighter than our H = 15.1. Again, the fit with the PSF does not match well in the central pixel. The quasar is marginally resolved with a reasonably good PSF + elliptical fit, with MH = -25.4 and Re = 8.0 kpc. Wyckoff, Wehringer & Gehren (1981) did not resolve the host in the optical with M(R) $`>`$ –21.9, whereas Véron-Cetty & Woltjer (1990) claim the host to be resolved with M(V) = –23.7.

PKS 1335-06. Romanishin & Hintzen (1989) resolved the host galaxy in the optical with M(B) = -23.7 and R(e) = 7.2 kpc. The radio morphology of this quasar is a bent triple with bridges connecting the components (Price et al. 1993, Bogers et al. 1996). It appears that we have detected NIR emission from the NW hot spot and its associated radio lobe, to the NW and N of the quasar (Fig. 2). However, the SE hot spot and lobe are not detected in our image.

References

Antonucci,R.R.J., 1993, ARA&A 31, 473
Bahcall,J.N., Kirhakos,S., Saxe,D.H., Schneider,D.P., 1997, ApJ 479, 642
Blandford.R.D., Rees,M.J., 1978, Pittsburgh Conference on BL Lac Objects (ed. A.M.Wolfe), 328
Bogers,W.J., Hes,R., Barthel,P.D., Zensus,J.A., 1996, A&AS 105, 91
Boyce,P.J., Disney,M.J., Blades,J.C. et al. 1998, MNRAS 298, 121
Bressan,A., Chiosi,C., Fagotto,F., 1994, ApJS 94, 63
Capaccioli,M., Caon,N., D’Onofrio,M., 1992, MNRAS 259, 323
Carballo,R., Sanchez,S.F., Gonzalez-Serrano,J.L., Benn,C.R., Vigotti,M., 1998, AJ 115, 1234
de Vries,W.H., O’Dea,C.P., Baum,S.A. et al., 1997, ApJS 110, 191
Djorgovski,S., Davis,M., 1987, ApJ 313, 59
Eales,S., Rawlings,S., Law–Green,D., Cotter,G., Lacy,M., 1997, MNRAS 291, 593
Ellingson,E., Yee,H.K.C., Green,R.F., 1991, ApJ 371, 49
Fukugita,M., Shimasaku,K., Ichikawa,T., 1995, PASP 107, 945
Glass,I.S., 1981, MNRAS 194, 795
Heckman,T.M., 1990, in IAU Colloquim 124, p. 359
Heckman,T.M., Bothun,G.D., Balick,B., Smith,E.P., 1984 AJ 89, 958
Hes,R., Barthel,P.D., Fosbury,R.A.E., 1996, A&A 313, 423
Hooper,E.J., Impey,C.D., Foltz,C.B., 1997, ApJ 480, L95
Hutchings,J.B., 1995, AJ 110, 994
Hutchings,J.B., Neff,S.G., 1990, AJ 99, 1715
Hutchings,J.B., Neff,S.G., 1992, AJ 104, 1
Hutchings,J.B., Crampton,D., Morris,S.L., Durand,D., Steinbring,E., 1999, AJ 117, 1109
Hyland,A.R., Allen,D.A., 1982, MNRAS 199, 943
Impey,C.D., Tapia,S., 1990, ApJ 354, 124
Kauffmann,G., Hähnelt,M., 2000, MNRAS 311, 576
Kormendy,J., 1977, ApJ 218, 333
Kormendy,J., Richstone,D., 1995, ARA&A 33, 581
Kotilainen,J.K., Falomo,R., Scarpa,R., 1998, A&A 332, 503 (KFS98)
Lacy,M., Bunker,A.C., Ridgway,S.E., 2000, AJ 120, 68
Landolt,A., 1992, AJ 104, 340
Ledlow,M.J., Owen,F.N., 1996, AJ 112, 9
Lehnert,M.D., Heckman,T.M., Chambers,K.C., Miley,G.K., 1992, ApJ 393, 68
Lehnert,M.D., Miley,G.K., Sparks,W.B., et al., 1999, ApJS 123, 351
Lidman,C., Cuby,J.G., 1998, SOFI User Manual, ESO
Lilly,S.J., Longair,M.S., 1984, MNRAS 211, 833
Lowenthal,J.D., Heckman,T.M., Lehnert,M.D., Elias,J.H., 1995, ApJ 439, 588
Magorrian,J., Tremaine,S., Richstone,D. et al., 1998, AJ 115, 2285
Malkan,M.A., 1984, ApJ 287, 555
McLeod,K.K., Rieke,G.H., 1994a, ApJ 420, 58
McLeod,K.K., Rieke,G.H., 1994b, ApJ 431, 137
McLeod,K.K., Rieke,G.H., 1995, ApJ 454, L77
McLeod,K.K., Rieke,G.H., Storrie-Lombardi,L.J., 1999, ApJ 511, L67
McLure,R.J., Dunlop,J.S., 2000, MNRAS, in press (astro-ph/9908214)
McLure,R.J., Kukula,M.J., Dunlop,J.S., et al., 1999, MNRAS 308, 377
Mobasher,B., Sharples,R.M., Ellis,R.S., 1993, MNRAS 263, 560
Moffat,A.F.J., A&A 3, 455
Moorwood,A.F.M. et al., 1992, ESO Messenger 69, 61
Neugebauer,G., Oke,J.B., Becklin,E.E., Matthews,K., 1979, ApJ 230, 79
Padovani,P., Urry,C.M., 1992, ApJ 387, 449
Percival,W.J., Miller,L., McLure,R.J., Dunlop,J.S., 2000, MNRAS, in press (astro-ph/0002199)
Poggianti,B.M., 1997, A&AS 122, 399
Price,R., Gower,A.C., Hutchings,J.B. et al., 1993, ApJS 86, 365
Quirrenbach,A. Witzel,A., Krichbaum,T.P. et al., 1992, A&A 258, 279
Romanishin,W., Hintzen,P., 1989, ApJ 341, 41
Rönnback,J., van Groningen,E., Wanders,I., Örndahl,E., 1996, MNRAS 283, 282
Silk,J., Rees,M.J., 1998, A&A 331, L1
Simpson,C., Rawlings,S., 2000, MNRAS, in press (astro-ph/0005570)
Small,T.A., Blandford,R.D., 1992, MNRAS 259, 725
Smith,E.P., Heckman,T.M., 1989, ApJ 341, 658
Spinrad,H., Dey,A., Stern,D., et al., 1997, ApJ 484, 581
Stiavelli,M., Treu,T., Carollo,C.M. et al., 1999, A&A 343, L25
Stockton,A., 1978, ApJ 223, 747
Stockton,A., MacKenty,J.W., 1987, ApJ 316, 584
Sun,W.H., Malkan,M.A., 1989, ApJ 346, 68
Taylor,G.L., Dunlop,J.S., Hughes,D.H., Robson,E.I., 1996, MNRAS 283, 930
Thuan,T.X., Puschell,J.J., 1989, ApJ 346, 34
Urry,C.M., Padovani,P., 1995, PASP 107, 803
Vermeulen,R.C., Cohen,M.H., 1994, ApJ 430, 467
Véron–Cetty,M–.P., Woltjer,L., 1990, A&A 236, 69
Véron–Cetty,M–.P., Véron,P., 1998, ESO Scientific Report
Wyckoff,S., Wehinger,P.A., Gehren,T., 1981, ApJ 247, 750
Yee,H.K.C., Green,R.F., 1984, ApJ 280, 79


format: text

text:
# A nonlinear hydrodynamical approach to granular materials

## I Introduction

The nature of the theory which describes the macroscopic transport in granular materials overview  remains an unsolved problem. We propose here a candidate theory, inspired by features of nonlinear hydrodynamics (NH). Given its success in treating transport in a variety of complex systems NH , NH seems natural to use in this context. In developing any continuum model of granular materials, however, one is challenged by the need to incorporate some remnant of the discrete nature of the underlying material.

The development of our theory is based on the hypothesis that the states of a system can be specified in terms of a few collective variables (at sufficiently long length scales), and that these states are connected in time by the local conservation laws supplemented with constituitive relations. Typically, the collective variables used are the conserved densities which, in the simplest view of granular materials, are conservation of mass and momentum. Conservation of energy is more complicated in this case than in simple fluids, due to locally inelastic processes collapse  which may transfer energy into internal degrees of freedom. Although we know how to include the energy density into the description, we begin with a more primitive theory which uses only the mass density $`\rho `$ and the momentum density $`\stackrel{}{g}`$.

The key difference between simple fluids and granular materials is that fluids organize over short time scales to be spatially homogeneous, while granular materials can exist over very long times in spatially inhomogeneous metastable states: metastable, because individual grains do not always pack together efficiently; long-lasting, because the large masses of the grains prevent thermal fluctuations from adjusting particles into a tighter configuration. Any complete description of granular materials must be able to explain the nature of these *quenched* inhomogeneous states. Experiment suggests that a sandpile contains “force chains” forcechains  which serve to support the pile in the presence of external stress; these chains are surrounded by regions of sand which are comparatively unstressed, and which allow the pile to flow and to settle under vibrations.

Since these sandpiles (in the absence of forcing) are static, the momentum is negligible, and we are left with the density field as our sole tool in describing the pile. Our primary hypothesis is that the density field alone is sufficient to capture the metastable nature of the system when at rest. Specifically, we describe the force chains as regions of sand in which the packing of the grains is ideal, having the maximum possible density. We will call these regions “close-packed”. The other, “loose-packed” regions have their grains arranged in some non-optimal configuration, with a slightly smaller density. We propose that these small variations in the density field are sufficient to describe the metastable nature of sand.

Using dynamical equations based on local conservation of mass and momentum, this *quenched* pattern in the density can be prepared in a dynamic process driven by an effective free energy with minima characterizing the loose-packed and close-packed regions. This presents the opportunity to grow a sand pile under the influence of gravity starting from an initial homogeneous state, and in the construction of the sand pile one builds up, in a self-consistent manner, an inhomogeneous equation of state characterized by a nontrivial stress tensor.

The basic structure of our theory, based on local conservation laws, is straightforward. More difficult is the inclusion of these competing close and loose-packed states. Clearly these require a detailed description of the shortest length scales treated in the model; there is a competition of length and energy scales of a type not encountered in simple fluids. Such a short-range description will depend on details of the particular granular material; we will be satisfied with a theory that demonstrates some generic properties of sand, as described in the next section. It seems reasonable that one may then be able to refine the model to include variations in the shapes and types of individual grains.

While the short-range length scales give difficulty, the large-scale nature of our model presents us with the advantage of a description close in scale to those phenomena which are most visible to experiment and casual observation. One ultimate goal is to address the existing Knight  macroscopic shaking and rotation experiments. Also, such a theory should eventually be useful in analyzing the surface states produced in large-amplitude shaking experiments surfacestates .

There have been a number of attempts Hydro  to use hydrodynamics to understand the dynamics of sand piles. These efforts differ from the one advanced here in that they do not propose to follow the full evolution of the system. There has also been a significant theoretical effort to describe granular materials through the use of kinetic theory kinetictheory . This approach is organized at a more microscopic level, where the connection to the macroscopic, static, metastable phase of the system is less direct foot-kinetic  and coupling to experiment more difficult. In many of these cases it is necessary to “liquefy” the granular material in question, giving it so much energy that it loses its metastable quality which seems to be crucial in describing certain aspects of a sandpile’s behavior.

In this paper we begin by describing some generic features of sand which we hope to incorporate into our model. We then proceed to a development of our dynamical equations from a generalized Langevin equation. This will lead into a discussion of the choice of a driving free energy, which will be crucial in creating the quenched pattern in the density field.

To test the validity of the theory we turn to simulation on a two-dimensional lattice. We begin with a rectangular box containing a homogeneous distribution of sand under the influence of gravity, and we show that it does indeed form a sand pile with an inhomogeneous distribution of loose and close-packed states. We then turn to simulations in a circular container, first forming a sand pile using the same technique, and then rotating the pile at different speeds. We compare the resulting angles of repose and oscillations about those angles to corresponding experimental data in the literature.

## II Properties of Granular Materials

We want to construct our model to be compatible with the following aspects of granular materials:

1. Granular materials have a clumping property; simulations of simple systems of inelastically colliding balls provide evidence for this phenomenon. One possible explanation for this clumping is the theory of *inelastic collapse* collapse : when two particles collide inelastically, they lose energy from their translational degrees of freedom and so recede more slowly than they approach each other. On average, the particles stay closer together than if they had collided elastically, and regions of higher density build up. This can also be described in terms of a hydrodynamical instability Goldhirsch .
2. Thermal energies in a sand pile are very small compared to, for instance, the average gravitational energy. Thus we should be able to ignore thermal noise.
3. Because of the large masses of the grains, one does not expect a significant vapor pressure above the interface as is found in liquid-gas systems. Indeed we expect a *very* dilute gas of grains above the pile, whose density exhibits a Boltzmann distribution due to the effect of gravity.
4. Sand piles are strongly driven by gravity and, because of the lack of thermal effects, the dense sand pile is separated from the dilute gas above by a sharp interface.
5. Granular materials are strongly disordered, with metastable structures forming upon creation of a pile. This is also due to the lack of thermal effects. As mentioned in the introduction, there is evidence of stress chains running through the bulk of sand piles, which may be associated with the ramified clumps formed in the more dilute systems studied using kinetic theory. The important point for us is that there is competition between somewhat more dense domains (stress chains and arches) and other, more loosely-packed regions in the pile.
6. Granular materials are stiff: when moved or jostled the pile can for a time behave as a solid object. Thus the sand pile should have a relatively uniform density.
7. One specific example of this stiffness is that a sand pile can maintain a non-horizontal surface. For example, when a pile of sand in a container is rotated about a horizontal axis, it does not flow immediately but waits until its surface passes a certain critical angle (called the angle of repose) with the horizontal. According to Reynolds Reynolds , this effect is due to the static interlocking grains that make up a granular material: the pile must dilate first, giving these grains freedom to move about, before flow is possible.
8. It is well known that granular materials can be very sensitive to boundary and finite-size effects; the famous *Brazil-nut-to-the-top* phenomenon Brazil-nut , for instance, is strongly influenced by both.

We will touch on these points while developing our model.

## III Dynamical Equations

Inspired by nonlinear hydrodynamics, we begin with the generalized Langevin equation Ma ; Das

$$\frac{\psi _\alpha (\stackrel{}{x},t)}{t}=V_\alpha [\psi (\stackrel{}{x},t)]\mathrm{\Gamma }_{\alpha \beta }(\stackrel{}{x})\frac{\delta F[\psi (\stackrel{}{x},t)]}{\delta \psi _\beta (\stackrel{}{x},t)}.$$
(1)
(Throughout this section, summation is implied over all repeated indices.) The variables $`\psi _\alpha `$ are the “slow” fields of interest in the problem: as mentioned earlier, these are the density field $`\rho (\stackrel{}{x})`$ and the momentum field $`\stackrel{}{g}(\stackrel{}{x})`$. The first term $`V_\alpha [\psi ]`$ is the *streaming velocity*, and corresponds to the reversible terms in typical hydrodynamical equations. As we will see, it depends on the Poisson brackets of the slow fields, as well as the derivatives of the Landau-Ginzburg-Wilson (LGW) effective free energy $`F`$. The second term is dissipative in nature, and $`\mathrm{\Gamma }_{\alpha \beta }`$ is the symmetric matrix of dissipative coefficients. Because granular materials are essentiallly zero-temperature systems (property 2), there is no thermal noise driving the system.

For our system, the free energy $`F`$ can be broken up into kinetic, potential, and external parts:

$$F=F_K+F_P+F_E.$$
(2)
From thermodynamics, one has the result that the variable conjugate to the momentum density is the velocity field:

$$\frac{\delta F}{\delta g_i(\stackrel{}{x})}v_i(\stackrel{}{x})\frac{g_i(\stackrel{}{x})}{\rho (\stackrel{}{x})},$$
(3)
where $`i`$ is a vector label. The kinetic energy contribution to the free energy thus has all of the momentum dependence:

$$F_K=d^d\stackrel{}{x}\frac{\stackrel{}{g}^2(\stackrel{}{x})}{2\rho (\stackrel{}{x})},$$
(4)
while $`F_P`$ and $`F_E`$ are functions only of the density. We will make further assumptions as to the form of $`F_P`$ below. $`F_E`$ is the free energy due to external forces: for instance, in the presence of a uniform gravitational field

$$F_E=d^dxg\rho (\stackrel{}{x})(zz_0),$$
(5)
where the scalar $`g`$ is the acceleration due to gravity and $`z_0`$ is the bottom of a confining box.

The streaming velocity $`V_\alpha `$ (in Eq. (1)) is given by the equation

$$V_\alpha [\psi ]=\{\psi _\alpha ,\psi _\beta \}\frac{\delta F}{\delta \psi _\beta }$$
(6)
where the indices $`\alpha ,\beta `$ run over the set $`\{\rho ,\stackrel{}{g}\}`$. We calculate the Poisson brackets by identifying the fields $`\psi _\alpha `$ with microscopic variables, evaluating their Poisson brackets, and expressing the results in terms of the $`\psi _\alpha `$, getting the standard results Das

$$\{\rho (\stackrel{}{x}),g_i(\stackrel{}{x}^{})\}=_i[\delta (\stackrel{}{x}\stackrel{}{x}^{})\rho (\stackrel{}{x})]$$
(7)
and

$$\{g_i(\stackrel{}{x}),g_j(\stackrel{}{x}^{})\}=_j[\delta (\stackrel{}{x}\stackrel{}{x}^{})g_i(\stackrel{}{x})]+_i^{}[\delta (\stackrel{}{x}\stackrel{}{x}^{})g_j(\stackrel{}{x})].$$
(8)
With these and the appropriate derivatives of the free energy, we can calculate the streaming velocities:

$$V_\rho =\stackrel{}{g}$$
(9)
and

$$V_{g_i}=\rho (\stackrel{}{x})_i\frac{\delta F_P}{\delta \rho (\stackrel{}{x})}\rho (\stackrel{}{x})g\delta _{iz}_j\left(\frac{g_i(\stackrel{}{x})g_j(\stackrel{}{x})}{\rho (\stackrel{}{x})}\right).$$
(10)
By choosing $`\mathrm{\Gamma }_{\rho \beta }=0`$, one easily obtains the usual continuity equation

$$\frac{\rho }{t}=\stackrel{}{g}.$$
(11)
The components of the damping tensor associated with the momentum density,

$$L_{ij}\mathrm{\Gamma }_{g_ig_j},$$
(12)
can be expressed in terms of the viscosity. If $`L`$ were independent of the fluctuating fields, then we could write it in the general form

$$L_{ij}=\eta _{il,kj}_l_k$$
(13)
where

$$\eta _{il,kj}=\eta _0(\delta _{il}\delta _{jk}+\delta _{ik}\delta _{jl})+(\zeta _0\frac{2}{3}\eta _0)\delta _{ij}\delta _{kl}$$
(14)
is the viscous tensor, and $`\eta _0`$ and $`\zeta _0`$ are the “bare” viscosities. It is more realistic, however, for the dissipation to depend on the density of the sand, and so we introduce a function $`\varphi (\rho )`$ into the damping tensor:

$$L_{ij}=\eta _{il,kj}_l\left(\varphi (\rho )_k\right).$$
(15)
We shall say more about $`\varphi (\rho )`$ below. It should be emphasized that, in a system with strong nonlinearities and spatial inhomogeneity such as ours, there exists no simple relationship between these bare viscosities and their physical counterparts.

We can now write the remaining Langevin equation:

$`{\displaystyle \frac{g_i}{t}}`$ $`=`$ $`V_{g_i}\mathrm{\Gamma }_{g_ig_j}{\displaystyle \frac{\delta F}{\delta g_j}}`$ (16)
$`=`$ $`\rho _i{\displaystyle \frac{\delta F_P}{\delta \rho }}\rho g\delta _{iz}_j(g_ig_j/\rho )L_{ij}{\displaystyle \frac{g_j}{\rho }}`$
At this point, we must choose a form for $`F_P`$. In this simple theory, we make the assumption that $`F_P`$ is of the square-gradient form

$$F_P=d^dx\left[f(\rho )+\frac{1}{2}c(\rho )^2\right],$$
(17)
where $`f(\rho )`$, the free energy density, is a local function of $`\rho `$, and $`c`$ is a positive constant foot-sqgradient . Thus we have

$$\frac{g_i}{t}=\rho _i\frac{f}{\rho }+c\rho _i^2\rho _j(g_ig_j/\rho )+\eta _{ij,kl}_j\left[\varphi (\rho )_k(g_l/\rho )\right]g\rho \delta _{iz}.$$
(18)
It should be noted that we can rewrite this equation as the divergence of a stress tensor plus a term representing external forces:

$$\frac{g_i}{t}=_j\sigma _{ij}+F_i^E.$$
(19)
This is the usual continuity equation for conservation of momentum, just as our expression for the density’s evolution is the continuity equation for conservation of mass.

For calculational reasons, we prefer to work not with the momenta $`\stackrel{}{g}`$, but the velocity fields $`\stackrel{}{v}`$, defined by Eq. (3) above. In terms of the density and velocity fields, the equations of motion become

$$\frac{\rho }{t}=(\rho \stackrel{}{v})$$
(20)
and

$$\frac{v_i}{t}=_i\frac{f}{\rho }+c_i^2\rho v_j_jv_i+\frac{1}{\rho }\eta _{ij,kl}_j\left[\varphi (\rho )_kv_l\right]g\delta _{iz}.$$
(21)
To proceed further, one needs to choose forms for $`\varphi (\rho )`$ and the viscosities. For very low densities we expect that the dissipative part of the equation should vanish; for this to happen, $`\varphi (\rho )`$ must go to zero at least as fast as $`\rho ^2`$. Higher powers of $`\rho `$ seem to make our numerical calculations more stable, so we have used $`\varphi =\rho ^4`$. This issue is not very important in the dense sand pile where $`\rho 1`$. We choose $`\zeta _0=\frac{5}{3}\eta _0`$ in Eq. (14) so that $`\eta _{ij,kl}`$ is isotropic:

$$\eta _{ij,kl}=\eta _0(\delta _{ij}\delta _{kl}+\delta _{ik}\delta _{jl}+\delta _{il}\delta _{jk})$$
(22)
In this theory, the input most crucial in specifying sand-like behavior is the free energy density $`f(\rho )`$. With an appropriate choice of $`f`$ we will be able to create a system that is very different from conventional fluids.

## IV Constructing the Free Energy Density

We will construct the free energy $`f(\rho )`$ using three steps, as outlined in Figure 1. The mathematical details of the description of the free energy can be found in the Appendix.

We can think of $`f(\rho )`$ as a potential where the system picks out values of the density which correspond to minima of $`f`$. A key assumption is that, from Property 1, grains of sand tend to clump together, and so the potential has a minimum which we arbitrarily place at $`\rho =1`$. The grains themselves are incompressible, so the potential rises rapidly for densities higher than $`\rho =1`$. Similarly, the potential becomes large as one approaches zero density, to prevent unrealistic negative values of the density anywhere in the pile.

Calculations show that this potential, so far, is enough to capture some properties of sand: the sharp interface of a sand pile, for instance. However, it is still missing one key element, and that is the metastability found in experiment. Any pile formed using this free energy would be homogeneous in density, and show no evidence of the close-packed and loose-packed regions mentioned in property 5. Without these regions, our material would be little more than a slightly compressible liquid.

In Figure 1b, we introduce this metastability into the system by placing a barrier into the potential just to the left of its minimum. We now have the following picture: the minimum represents the optimal packing for the sand grains. As sand comes together, it increases the density of the pile until it reaches the density of the barrier where it is frustrated: it can no longer proceed to the optimal packing. However, as the sand begins to pile up, pressures in the pile can provide enough of a push so that some regions may pass the barrier and reach the minimum of the potential.

This potential has all of the right features, but in simulations it proves to make a sandpile that is too soft and liquid-like. To stiffen the pile, we introduce two wells, as in Figure 1c: one for the close-packed states and one for the loose-packed. With the wells, this is the free energy we have used in our calculations. Future possible modifications to this free energy will be discussed in our conclusion.

## V Numerical Analysis

Because the proposed model has a strongly nonlinear nature, we are forced to study this system’s evolution numerically, hoping that this will inspire subsequent analytical work. The first order of business is to determine numerically whether this system of equations, running forward in time, can prepare a state that looks like a sand pile. In principle, this should be straightforward to do. In practice, there are a number of difficulties in putting these partial differential equations onto a space-time lattice.

1. Lattice spacing: It is understood that, to solve PDE’s numerically, one must keep the time step much smaller than the lattice size; otherwise, the system develops physical instabilities. We have chosen a typical convention of setting $`\mathrm{\Delta }x=1`$ and $`\mathrm{\Delta }t=0.001`$.
2. We know that there are important finite-size effects in granular systems, and so we must be prepared to examine our model in various sizes and shapes of containers. This paper focuses on two such containers: a rectangle 200 units high and 100 units wide, and a circle with a 100 unit diameter. Other containers have been tried as well, yielding similar results.
3. One of the discouraging aspects of this approach is the large number of parameters governing the system. There are parameters to control the intrinsic free energy, the square-gradient energy, as well as those for external forces like gravity. There is the finite size of the box and the nature of the boundary conditions. Dynamically, there are the viscosities and the initial level of the density and velocity fields. In this paper, the values of these parameters have been chosen to obtain structures that resemble sand piles. Most of these choices are reported in the appendix where the form of the free energy is described in detail; in addition, we take $`c=10`$ and $`\eta _0=12`$. Further work is needed to map out the range of parameters for which we obtain physical behavior, and determine which parameters are ultimately important.
4. This numerical problem would be rather straightforward if not for a set of persistent unphysical instabilities; without care the system can and will explode in a very unnatural fashion. These explosions are of two, closely related types: *runaway velocities* and *negative densities*. To alleviate the first, we have taken the pragmatic step of locally averaging the velocity if it gets larger than some empirically determined cutoff.

Our free energy density has been designed with a barrier at $`\rho =0`$, but negative densities do tend to show up in our calculations, specifically in the dilute region above the sand pile. We take two steps to counter them: when the density of a site is small and positive ($`\rho <0.05`$), we average the site’s velocity with half of the neighbors’ average velocity. In addition, when the density dips below zero, we bring in density from its four neighbors to bring it above zero. In treating these instabilities, we have been very careful not to violate conservation of mass.

## VI Forming a Sand Pile

As a first example of how our model works, we begin with a nearly uniform distribution of sand which, under the influence of gravity, forms a sand pile at the bottom of the container. Our container is a two-dimensional box, overlaid with a square lattice 100 units wide and 200 units high. Boundary conditions on the walls of the container are non-slip. Our initial conditions are $`\rho =0.5+\mathrm{\Delta }\rho `$ and $`\stackrel{}{v}=\mathrm{\Delta }\stackrel{}{v}`$, where $`\mathrm{\Delta }\rho `$ and $`\mathrm{\Delta }\stackrel{}{v}`$ are both small (0.001 rms) random perturbations in the system. Because mass is conserved and because the density inside the formed pile should have a density about $`1`$, we expect the interface to form in the middle of the box, around $`z=100`$. In principle, one should average over ensembles of initial conditions. That is not necessary for our purposes here, but will become important in more quantitative work.

The first question is whether the model captures the overall dynamics of the system: our physical picture of the situation has the sand particles all under the influence of gravity, so at first there will be a net acceleration. As particles begin to hit the bottom of the container, their velocity drops to zero, and the system as a whole decelerates until it ends up at rest. To see if this is true we measure the average kinetic energy

$$KE=\frac{1}{V}\underset{z,x}{}\rho (z,x)v^2(z,x)$$
(23)
where the sum is over the entire box, and $`V`$ is the total number of lattice sites. Figure 2 shows a plot of the kinetic energy for two runs with different gravitational accelerations.

In both cases the behavior is just as expected, with an immediate acceleration followed by a deceleration to rest. Also interesting to note is the difference between the curves: the run with stronger gravity causes the pile to form faster, with a higher maximum kinetic energy resulting from the higher initial potential energy of the system.

To get a clearer picture of the stationary state our system is settling into, we consider the density profile:

$$P(z)=\frac{1}{N_x}\underset{x}{}\rho (z,x)$$
(24)
where $`N_x`$ is the number of sites across the lattice.

Figure 3 superimposes snapshots of this profile at seven instances in time. Initially, the profile is a flat line $`P(z)=0.5`$, representing the initial homogeneous distribution. As time progresses, the density begins to increase for low $`z`$, indicating pile formation; at the same time, the density at higher altitudes is decreasing. By $`t=30`$ the system has settled into a phase-separated system of high-density pile below and low-density “gas” above, with a sharp interface between the two (agreeing with property 4 in section II.)

The height of the pile can be defined as that value of $`z`$ for which the profile reaches $`P(z)=0.5`$, or half the average density of the pile foot-profile . This is plotted as a function of time in Figure 4, along with the vertical coordinate of the system’s center of mass:

$$z_{CM}=\frac{\underset{x,z}{}z\rho (z,x)}{_{x,z}\rho (z,x)}$$
(25)
with the sums over all the points in the lattice.

At first the center of mass moves downward in time parabolically, as it should with the system in free fall. As the sand starts to pile up, less and less of the sand is in motion and the center of mass approaches a constant height of 52.6. The height of the interface moves contrary to this, of course, settling in at 102, about halfway up the container. Since the average density of the dense pile should be about $`1`$, and we began with a uniform distribution of density $`0.5`$, the position of the interface is just a little bit higher than our expectation that the container be half full of sand, suggesting that the pile has an average density slightly less than one. The center of mass is slightly higher than the midpoint between the base and the interface of the pile, suggesting that the pile is somewhat top-heavy, which is surprising.

Since on the large distance scale our system looks like a sand pile, we turn our attention to the small-scale features inside the pile. Specifically, we want to investigate the competing regions of loose-packed and close-packed states, mentioned in property 5 of section II, corresponding to the two wells in the free energy. Recall that these wells are separated by a barrier, which we have placed at $`\rho =0.99`$; for definiteness, we say that all sites with $`\rho >0.99`$ are close-packed and all sites with $`0.85<\rho 0.99`$ are loose-packed.

In Figure 5 we count the number of loose and close-packed sites as a function of time, for two different strengths of gravity. Naturally, the loose sites start to form first; some of these are then pushed across the barrier and become close-packed due to the weight of sand on top of them. Both graphs show the numbers of loose and close-packed sites approaching a constant as the pile settles. Notice that the final ratio of low-density to high-density sites is sensitive to parameters such as gravity: a stronger gravitational field has a greater capacity to compact the sand and create more close-packed sand.

Figure 6 shows how these loose and close-packed regions are distributed through the pile for the same two values of $`g`$. The top halves of each pile share similar features: both have a layer of loosely-packed sand on top, where there is no other sand to weigh it down and compress it. Immediately below the surface in the center there is a region of tightly-packed sand, which dips farther below the surface in the center of the pile than to either side. This region is flanked by columns of loosely-packed sand, which are in turn flanked by dense sand up against the edges of the box. The differences in the two piles are mostly at the bottom: the run with less gravity is mostly loose-packed in the bulk of the pile, while the run with more gravity is tighter. Note however that in neither case is the density uniform: one can find patches of both types in both piles. In particular, the patterns in the center of the $`g=0.5`$ pile may hint at the stress chains that are seen in experiment.

## VII Rotating a Sand Pile

### VII.1 Set-up

One of the standard ways to probe a granular system is to rotate it about some horizontal axis. This method demonstrates several characteristics of a granular material, and the most important of these is the angle of repose. Actually, there are (at least) two such angles associated with rotation. If one begins with sand in a cylindrical container, having a horizontal free surface, and then begins to rotate the cylinder about its axis, the surface of the sand will increase in slope along with the container’s rotation until it reaches some critical slope, at which point an avalanche will restore the surface to a more horizontal state. This we might call a “static” angle of repose, $`\theta _s`$. If one continues to rotate the cylinder, then one of two things will happen:

1. If the rotation is slow enough, then the flow along the surface will come to a stop before a further increase in slope should cause an avalanche again. The result is a periodic, intermittent flow, and the interface oscillates about some average angle.
2. If the rotation is faster, the flow from the first avalanche does not have time to stop before the next avalanche begins. The system settles into a state where there is a continuous flow along the surface, and the interface maintains a constant angle with the horizontal.

In both cases, we may take the mean angle that the surface makes with the horizontal over time and call it a “dynamic” angle of repose, $`\theta _d`$. In the periodic case, a quantity as important as the angle of repose is the size of the fluctuation about that angle, $`\delta _d`$. The actual values of these angles seems to depend experimentally on a number of factors, including particle size and shape, the humidity of the air, how the pile is formed, boundary conditions, and so forth. The experimental results of Jaeger et al. JLN , for instance, find that spherical glass beads with diameters of about half a millimeter show a dynamic angle of repose of 26, with a fluctuation of 2.6; while rough aluminum-oxide particles with the same diameter show a higher angle, 39, with a fluctuation of 5. These seem to be typical values. What will be most important in our qualitative analysis is that these angles exist and are non-negligible.

We set out to find evidence of these two angles in our model. To better match experiments in this field, we move from a rectangular box to a circular one having a diameter of 100 lattice sites and the same no-slip boundary conditions. We grow a sand pile here as we did before: Figure 7 shows the pattern of loose and close-packed states in a pile formed in our circular container.

To best mimic rotation under experimental conditions, one would ideally set up rotating boundary conditions, giving a constant velocity to the sand at the edge of the container. This has turned out to be difficult to do in practice: one reason is that, due to the lattice, our boundary is not a perfect circle, and there is a tendency for mass to “leak” in or out of the boundaries when the boundary conditions are not strictly no-slip. Because of this, we choose to implement rotation by rotating gravity with a constant period of rotation $`T`$.

To look for an angle of repose we need to measure the angle that the interface makes with the horizontal (that is, the normal to gravity). The most direct way to do this is to fit the interface of the pile with a line, and measure the angle that this line makes with the horizontal. This is indeed one of our probes: we define as our interface those points which are in the pile ($`\rho >0.5`$) which have a nearest neighbor outside of the pile, and, to reduce boundary effects, we throw out those points that are not within one half radius from the center.

This method should be ideal, but because of the finite size of our system, this measure ends up depending on only a few points which makes it sensitive to various local perturbations. For this reason, we introduce another measure which we call the “bulk angle” (where the first might be called the “surface angle”). This latter measure is simply the angle that a line passing through the center of the container and the center of mass makes with the vertical. This measure depends on the entire system and so is less sensitive to noise. It will be equal to the surface angle in the event that the pile has reflection symmetry across the vertical axis; it will differ from the surface angle by some fixed amount if the shape of the pile during rotation is asymmetric but constant. Figure 8 shows a plot of the two measures for a typical run.

### VII.2 Static Angle of Repose

With the surface angle and bulk angle measures in place, let us look systematically at our data. First, we consider the behavior of the pile in the first moments of rotation, with three different periods $`T`$ (Figure 9).

In all three cases the behavior is similar: as the container turns, the interface’s slope increases, until its angle reaches some maximum and begins to decline due to one or more avalanches. Clearly, there is a static angle of repose here, which depends on the rotation speed. Notice, however, that if there were no activity in the sandpile before the first maximum in the angle, then the curve would follow the straight line denoted in the figure as “solid” until turning downward. Instead, it seems that the surface is losing slope, lagging behind the container, from a very early time: there is some preliminary flow in the pile, even before the first major avalanche.

Figure 10 contains six snapshots of the flow during this initial period, with $`T=200`$. All of the pictures here are in the lab frame (that is, gravity’s frame), so the first picture ($`t=4`$) shows the pile moving in unison with the container, as if it were a solid and fastened to the walls. At $`t=8`$, however, things are beginning to change: the center of rotation seems to have moved to the right and downward, and the sand coming up on the right side is beginning to curve over to the left. The result of this action is seen in the next picture, where there is a definite flow along the surface of the pile. Notice that, even though there is surface flow, it is not enough to prevent the angle of the interface from climbing, according to Figure 9. The next picture, at $`t=16`$, shows the surface flow continuing, and also that sand is beginning to climb up the side of the right wall, due to our no-slip boundary conditions. At $`t=20`$, a major avalanche begins, so that the angle of the interface begins to fall (see Figure 9; this corresponds to 36 on that figure’s horizontal axis). The last picture ($`t=30`$) shows a much larger flow of material than seen earlier.

In short, there is clearly an initial period where the sand moves with the container and where there is no surface flow, evidence for a nonzero static angle of repose. Recall from section II that this effect is due to the need of the pile to dilate before it can flow. In our model, the close-packed sites are the ones which are restricted in their ability to dilate (because of the barrier in the free energy). Thus, if the pile is dilating during this flow, we expect that the number of close-packed sites in the pile must be decreasing. Figure 11 shows that this is indeed so.

### VII.3 Dynamic Angle of Repose

We next consider how the pile behaves under further rotation with three or more turns of the system. Figure 8 above shows the behavior of the bulk and surface angles over three complete revolutions, for $`T=200`$. Figure 12 below compares the bulk angles between $`T=200`$ and $`T=1600`$ over three revolutions.

Most notable in Figure 12 is that the average slope of the interface is larger for a faster rate of rotation. Also, the slower system deviates much less from a constant angle than does the faster one, with more jaggedness suggesting frequent small avalanches. These points suggest an inertial effect: in the faster case the sand does not have as much time to avalanche and level its interface, while the slower case has a number of short avalanches that help to keep its interface closer to the horizontal. However, the slow line is occasionally punctuated by large avalanches that dip into the negative: curiously, the interface seems to tilt in the direction opposite that of rotation every once in a while. Figure 13 looks at one of these dips more closely. Apparently, even with the constant small avalanches that we see in the bulk angle in the slowly rotating case, there is still a build-up of potential energy that must be released by these larger events.

All of our tests so far show a nonzero dynamic angle of repose, but to rule out the possibility that these angles are due to viscosity or inertia, it is best to find how the average angle (bulk or surface) depends on the period of rotation $`T`$, and from this relationship extrapolate to $`T\mathrm{}`$. With this in mind we measure the average bulk and surface angles over time for several rotation speeds, and, in Figure 14, plot them versus the rotation speed $`1/T`$.

The first thing to note in these plots is that neither angle is going to zero in the limit of small angular velocities: there is a definite non-zero angle of repose in our system. This angle, about 2.5, is quite small compared to experiment where typically one finds angles of repose on the order of 30JLN ; Rajchenbach .

It is interesting to fit the angles in Figure 14 to a power law

$$\theta =\theta _0+cT^\alpha .$$
(26)
Such fits are shown in the figure, with the parameters given in Table 1. The intercepts $`\theta _0`$ are positive, confirming a non-zero angle of repose. The fit to the surface angle is not quantitative, but for the average bulk angle the power law fits well, with an exponent of $`1.4`$.

These data points come from single runs, and averaging over an ensemble of initial conditions may make the fits more quantitative.

Experimentally, Rajchenbach Rajchenbach  finds that the angular velocity $`\mathrm{\Omega }`$ goes as

$$\mathrm{\Omega }(\theta \theta _d)^m$$
(27)
where $`m=0.5\pm 0.1`$. This corresponds to

$$\theta T^{2\pm 0.4}$$
(28)
Of our two measures, the bulk angle comes closer to matching experiment, but there is still some discrepancy that needs to be addressed.

Finally, to compare the surface and bulk angles with each other. For high-speed rotations, the bulk angle is higher on average than the surface angle, reflecting the plume of material that creeps up the side of the container (which is accounted for in a bulk calculation, but specifically excluded from the surface angle). The surface angle actually seems to level off at high speeds. At lower speeds, the bulk angle actually dips below the surface angle, suggesting that the plume has disappeared (as is seen in Figure 13). Notice that the fluctuations in the surface angle are larger than those in the bulk angle (the former depending on fewer lattice sites and thus more volatile), and that in both angles these fluctuations get smaller as one reduces the speed of rotation (but do not seem to go to zero).

## VIII Discussion

We present evidence that one can create a nonlinear hydrodynamical model for granular materials that depends only on density and momentum fields. Qualitatively, our model demonstrates many key features of sand, including a sharp interface, a relatively uniform density, a nonzero angle of repose, and a metastable structure. The method allows us to follow the pile from creation forward, in static and dynamic situations, and the model generally behaves like a sand pile.

There are several ways in which the model could be improved. The angles of repose seen here are too small when compared to experiment, and the way that the angles scale with rotation speed is at odds with Rajchenbach’s findings. This may not be a problem since we have not yet looked at the variation of the angle of repose with the parameters characterizing the model. An explanation for the difference may also lie in the fact that our rotation probe differs from the typical experimental method of rotating a sandpile. In our simulations, where we rotate gravity, every particle feels the external force directly and immediately. In experiment, where the container is rotated, the external force must be transmitted inward from the boundaries. This difference may be enough to account for the discrepancy between simulational and experimental outcomes. One may be able to mimic the rotation of the container in our model by introducing centrifugal forces into the system. Such forces would have a magnitude of

$$F_c=\rho r\left(\frac{2\pi }{T}\right)^2,$$

where $`r`$ is the distance from the center of the container. Our round container has a diameter of 100 units, so $`r50`$, while $`\rho 1`$ and the fastest rotation speed we use is $`T=200`$; thus these centrifugal forces would have a magnitude of $`0.05`$, which is small (though not negligible) compared to the main acting force of gravity, which is of magnitude roughly equal to $`1`$. Thus centrifugal forces may provide some quantitative effect, but in our initial, qualitative presentation here we deemed it unnecessary to include them.

Another way of improving the model is to remove the constraint that the loose-packed regions have a single fixed density. There are many non-optimal, metastable ways to pack particles together, having a range of different densities. One solution may be to allow the position of the loose-packed well to fluctuate slightly at random through the pile. There are several plausible ways of implementing this idea, which we intend to pursue.

We have only begun to extract information using our model; there are other probes we can use to perturb our system. Shaking can cause the pile to slowly settle into a denser state; we hope to find the logarithmic time dependence seen in compactification experiments Knight . Applying pressure to the system may allow us to investigate force propagation in the pile, and determine the nature of stress chains. It should be possible to modify the model to depict a pile made up of two or more types of particles, to investigate the phenomena of unmixing and the Brazil nut effect in a hydrodynamical setting. The flexibility of the fluctuating nonlinear hydrodynamical approach gives us a wide range of avenues to investigate.

## Appendix: Specification of the Free Energy Density

The first potential in Figure 1 is made up of four terms:

$$f_a(\rho )=f_{\text{clump}}(\rho )+f_{\text{large}}(\rho )+f_{\text{negative}}(\rho )+f_{\text{entropy}}(\rho )$$
(29)
where

$`f_{\text{clump}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}u_0\rho ^2`$ (30)
$`f_{\text{large}}`$ $`=`$ $`Be^{\lambda (\rho ^21)}`$ (31)
$`f_{\text{negative}}`$ $`=`$ $`U_1e^{\lambda _1\rho }.`$ (32)
The first term models the clumping behavior of sand mentioned in property 1 in section II. The second term models the ultimate incompressibility of the sand grains: $`B`$ is chosen so that the minimum of the well is at $`\rho =1`$. The third term prevents the densities from becoming negative by putting a barrier at zero density.

We have not given the definition for $`f_{\text{entropy}}`$ yet. Our original intent was to have this term model the behavior of the very dilute gas that exists above our sand pile. So that the low-density regions show a Boltzmann distribution, we used the standard gas entropy term

$$f_{\text{entropy}}^{}=A(\rho \mathrm{ln}(\rho /\rho _0)\rho )$$
(33)
with $`\rho _0=0.05`$. However, this was a source of numerical instability, and since our focus was the pile and the high-density regions, we replaced this with a simpler term,

$$f_{\text{entropy}}=A\rho .$$
(34)
The reason we did not eliminate the term entirely was that the original term created a shallow minimum around $`\rho =0.05`$, and for consistency we decided to keep that minimum there which the linear term allows us to do.

Clearly, this model has a lot of parameters, and one part of our future work will be to find optimal values for these parameters. Our current choices were selected because they gave realistic results: we set $`u_0=2`$, $`\lambda =40`$, $`U_1=0.6`$, $`\lambda _1=400`$, and $`A=0.2`$. To satisfy the requirement that the minimum of the potential be at $`\rho =1`$, we set $`B=\frac{1}{2\lambda }(u_0A)`$.

The barrier introduced in Figure 1b is described by

$$f_{\text{barrier}}=u_be^{\frac{1}{2}k(\rho \rho _b)^2},$$
(35)
where $`u_b=1.5`$, $`\rho _b=0.99`$, and $`k=10^5`$. The two wells in the final potential are described by

$$f_{\text{well}\text{n}}=u_ne^{\frac{1}{2}k(\rho \rho _n)^2},$$
(36)
where $`\rho _1=0.98`$ and $`\rho _2=1.01`$, $`u_1=0.25`$ and $`u_2=0.4`$, and $`k`$ is the same as in the expression for the barrier.

###### Acknowledgements.

We would like to thank Professors Heinrich Jaeger and Sidney Nagel for helpful conversations. This work was supported by the Materials Research Science and Engineering Center through Grant No. NSF DMR 9808595.