format: text

text:
# Modes of counterion density-fluctuations and counterion-mediated attractions between like-charged fluid membranes

## Abstract

Counterion-mediated attractions between like-charged fluid membranes are long-ranged and non-pairwise additive at high temperatures. At zero temperature, however, they are pairwise additive and decay exponentially with the membrane separation. Here we show that the nature of these attractions is determined by the dominant modes of fluctuations in the density of counterions. While the non-pairwise additive interactions arise from long-wavelength fluctuations and vanish at zero temperature, the short-ranged pairwise additive interactions arise from short-wavelength fluctuations and are stronger at low temperatures.

87.22.Bt,82.65Dp

Counterion-mediated attractions play a significant role in many physical and biological phenomena . These attractions are responsible for fitting DNA inside a small biological container such as a viral capsid or a nuclear envelope  and can also be crucial in promoting adhesion and fusion of biological membranes . Accordingly, significant effort has been expended in developing a practical way of investigating the nature of these counterion-mediated attractions. Besides an integral equation method , there have emerged two distinct approaches. The first approach , based on a charge fluctuation picture, suggests that these attractions are mediated by correlated fluctuations of ion clouds of counterions. This approach is consistent with our conception of counterions as fluctuating objects and thus merits significant consideration . In the second approach, based on a zero temperature picture , the appearance of attractions between like-charged molecules is attributed to the strong charge correlations that drive the systems, together with counterions, into an ionic crystal. At first glance these two approaches appear to be contradictory to one another, but there is an evidence that they can, in fact, be complimentary . Despite this, there still remain fundamental discrepancies between the two that have yet to be resolved. For the case of two planar surfaces a distance $`h`$ apart, the zero-temperature picture leads to an attractive force that decays expenentially with $`h`$ . In the charge-fluctuation approach, however, the attractive force decays algebraically as $`h^3`$ , as long as $`h`$ is sufficiently large. Most recently, the apparent discrepancy between these two approaches has been examined  by Lau et al. Lau et al. have shown that the zero-temperature quantum fluctuations give rise to a long-ranged attraction which varies as $`h^{7/2}`$. Their theory, however, does not explain the crossover between the long-ranged interactions driven by thermal fluctuations and the exponentially-decaying interactions expected at low temperature. At finite temperatures, the long-ranged interactions supported by the charge-fluctuation approach should exist and constitute a dominant contribution to the plate-plate interactions. Furthermore, when applied to many rod systems , the charge-fluctuation approach suggests that interactions between rods are not pairwise additive while the exponentially decaying interactions between plates as implied by the zero-temperature approaches are pairwise additive.

In this Letter, we present a simple theoretical approach to counterion-mediated interactions between fluid membranes in order to bridge the gap between the two existing theories. In particular, we show that the nature of these attractions is controlled by the dominant modes of fluctuations in the density of counterions. At high temperature, the membrane interactions are dominated by long-wavelength fluctuations in the counterion density, ı.e., fluctuations at large length scales. The resulting long-wavelength interactions are shown to be long-ranged and non-pairwise additive. As the temperature decreases, the high-temperature behavior of the membrane interactions crosses over to the behavior expected at low temperatures. At low temperatures, the membrane interactions are mainly determined by the short-wavelength charge fluctuations. We find that the resulting interactions decay exponentially with the separation between the membranes and are approximately pairwise additive. Finally, we obtain a phase diagram to depict the two distinctive regimes characterized by the corresponding dominant modes, and the crossover boundaries between the two.

The system we consider here consists of negatively charged parallel membranes surrounded by neutralizing counterions of opposite charge. In the following derivation, we assume that counterions are localized in the plane of the membrane. Thus our approach presented here is relevant to the strongly charged case. More precisely, the separation between plates is much larger than the Gouy-Chapman length, a typical length scale within which counterions are confined. The main purpose of the present work is to study the crossover from the high-temperature results for the membrane attractions to the behavior expected at low temperatures. Since this crossover occurs at low temperaures or at high densities of counterions, the assumption of localized counterions is reasonable.

The charge distribution on a layer $`j`$ is described by the local surface charge density, $`\widehat{\sigma }_j(𝐫_{})=e\sigma _0+emZ`$, where $`e`$ is the electronic charge, $`\sigma _0`$ is the average counterion number density, $`m=0,1,2,3,`$ etc., is the number of counterions per unit area at $`𝐫_{}(x,y)`$, and $`Z`$ is the counterion valency. In order to calculate the free energy and the corresponding charge correlations (cf. Eq. 6), we use two-dimensional Debye-Hückel (DH) theory. This approximation, however, fails to capture the strong charge-correlations at low temperatures. This defect in the DH theory has been corrected in an approximate way by taking into account short-ranged charge correlations over the size of ions . We thus implement DH theory with the counterion size via the two-dimensional form factor $`g(𝐫_{},𝐫_{}^{})=\mathrm{\Theta }(|𝐫_{}𝐫_{}^{}|D)/\pi D^2`$ where $`D`$ is the diameter of the counterions. We find that the charge-fluctuation contribution to the free energy per area is

$$\frac{\mathrm{\Delta }F_N}{k_BT}=\frac{1}{2}\frac{d𝐤_{}}{(2\pi )^2}\left\{\mathrm{log}\left[detQ(𝐤_{})\right]N\frac{g(𝐤_{})}{\lambda k_{}}\right\},$$
(1)
where $`\lambda ^1=2\pi Z\mathrm{}_B\sigma _0`$ is the inverse Gouy-Chapman length and $`\mathrm{}_B=e^2/ϵk_BT`$ is the Bjerrum length, ı.e, the length scale at which the electrostatic energy between two charges is comparable to the thermal energy. The matrix $`Q(𝐤_{})`$ is defined by the matrix elements

$$Q_{ij}(𝐤_{})=\delta _{ij}+\frac{g_{ij}(𝐤_{})}{\lambda k_{}}\mathrm{e}^{k_{}h_{ij}}$$
(2)
where $`g_{ij}(𝐤_{})`$ is $`g(𝐤_{})\frac{2J_1(k_{}D)}{k_{}D}`$ if $`i=j`$ and 1 otherwise, and $`J_1(x)`$ is the first-order Bessel function of the first kind. Finally, $`h_{ij}`$ is the separation between plates $`i`$ and $`j`$

First note that the free energy of an $`N`$-plate system is not simply a pairwise sum of the corresponding two-plate results over all pairs of plates. Thus the pairwise additivity is not always satisfied as will be detailed later. For the $`N=2`$ case, our result in Eq. 1 reduces to the previous two-plate result \[see Eq. 3 of Ref. \], if $`D`$ is set to zero. When $`D=0`$, the free energy has a single minimum at a nonzero value $`k_{}=k_{}^<1\AA ^1`$ for all values of $`\lambda `$. The dominance of the long-wavelength charge fluctuations is responsible for the breakdown of the pairwise additivity of electrostatic interactions between macroions. It has been shown that the pairwise additivity for the case of charged rods breaks down if the expansion of the corresponding interaction free energy in powers of $`\mathrm{}_B`$ diverges . For rod systems, the free energy is dominated by the zero-$`k`$ mode and thus this expansion converges only when the charge fluctuation along the rods is vanishingly small. For the two-dimensional case, the convergence of the $`\mathrm{}_B`$-expansion can be tested by estimating $`\delta (k_{}^<)^1\lambda ^1`$; the $`\mathrm{}_B`$-expansion is convergent if $`\delta <1`$. We, however, find that $`\delta `$ is smaller for smaller values of $`\lambda ^1`$ and is comparable to unity if $`\lambda ^1<10^3`$. As in the rodlike systems, the pairwise additivity is violated unless the inverse of the Gouy-Chapman length is very small. Thus the pairwise additivity can easily be violated even in two-dimensional systems.

The above analysis, however, cannot be pursued when the temperature is low or when the density of ions is high. In this case, it is crucial to incorporate $`D0`$ ; the ion size constitutes an important length scale at low temperatures or at high density of ions. Having understood that the dominance of long-wavelength charge fluctuations for the case $`D=0`$ leads to the breakdown of pairwise additivity, we now examine the low-$`T`$ behavior of the free energy (thus with $`D`$ set to a finite value) in Fourier space. Let us first consider the case of two plates seperated by a distance $`h`$. For convenience, we consider the following quantity: $`𝒬(𝐤_{})k_{}\mathrm{log}[detQ(𝐤_{})]`$, ı.e., the first term in $`\{\mathrm{}\}`$ of Eq. 1 mutiplied by $`k_{}`$. In Fig. 1, we plot this quantity $`𝒬(𝐤_{})`$ as a function of $`k_{}`$ for several different values of the Gouy-Chapman length $`\lambda `$. We have chosen $`D=5\AA `$ and $`h=5\AA `$. When $`\lambda ^1=1\AA ^1`$, $`𝒬(𝐤_{})`$ has a single minimum at $`k_{}=k_{}^<1`$. This implies that the free energy is dominated by long-wavelength charge fluctuations as in the previous case of point charges. We find that the function $`𝒬(𝐤_{})`$ has two minima at $`k_{}=k_{}^<1\AA ^1`$ and at $`k_{}=k_{}^>=𝒪(1\AA ^1)`$, respectively, for $`\lambda ^1>4.2\AA ^1`$ (not shown in the figure). As $`\lambda ^1`$ changes, the minimum at $`k_{}=k_{}^>`$ varies monotonically and is deeper for larger $`\lambda ^1`$. The minimum at $`k_{}=k_{}^<`$, however, is roughly independent of $`\lambda `$. When $`\lambda ^15.9\AA ^1`$, the two minima are comparable in magnitude. When $`\lambda ^1=7\AA ^1`$, the function $`𝒬(𝐤_{})`$ is overwhelmingly dominated by the second minimum at $`k_{}=k_{}^>`$, as shown in the figure. At $`\lambda ^1=\lambda _X^1=7.2\AA ^1`$, the second minimum diverges. This is suggestive of the onset of crystallization of counterions as will be detailed later. Even though the region $`\lambda ^1\lambda _X^1`$ is certainly beyond the validity of our theory, it is neverthless interesting to see what our theory impiles for that region. Notably, the free energy curve corresponding to $`\lambda ^1=10\AA ^1`$ has two local minima at large $`k_{}=𝒪(1\AA ^1)`$. The existence of multiple minima at large $`k_{}`$ assures that the system is in a solidlike phase.

Our results in Fig. 1 imply that there are two distinct contributions to the free energy: long-wavelength fluctuations and short-wavelength fluctuations in the density of counterions. They also imply that the short-wavelength fluctuation contribution to the free energy has a much narrower peak if $`\lambda ^15.9\AA ^1`$. This enables us to separate the short-wavelength contribution from the long-wavelength contribution. By noting that $`\mathrm{e}^{k_{}h}`$ does not change appreciably over the region inside the peak at $`k_{}=k_{}^>`$, we find, up to $`h`$-independent terms

$`\mathrm{\Delta }F_2{\displaystyle \frac{k_BT}{16\pi }}{\displaystyle \frac{\zeta (3)}{h^2}}`$ (3)
$`{\displaystyle \frac{k_BT}{8\pi ^2}}{\displaystyle \frac{\mathrm{e}^{2k_{}^>h}}{\lambda ^2}}{\displaystyle _{k_{}k_{}^>}}{\displaystyle \frac{k_{}^1dk_{}}{[1+k_{}^1\lambda ^1g(k_{})]^2}},`$ (4)
where $`\zeta (x)`$ is the zeta function (thus $`\zeta (3)/16\pi 0.024`$). The first term denoted by $`F_{LW}`$ is the free energy calculated with $`D`$ set to zero  and is the long-wavelength free energy. Our previous analysis on $`𝒬(𝐤_{})`$ implies that the short-wavelength free energy denoted by $`F_{SW}`$, ı.e., the second term in Eq. 3, is dominant over the long-wavelength free energy $`F_{LW}`$ at low temperatures ($`\lambda ^15.9\AA ^1`$) and decays exponentially in space. The exponentially decaying interaction between two plates is consistent with the zero-$`T`$ analysis in Refs. . At high temperatures corresponding to $`\lambda ^15.9\AA ^1`$, however, the free energy is mainly determined by $`F_{LW}`$. In this case, the fluctuation contribution to the pressure between the two plates scales as $`h^3`$ . Also note that this long-wavelength contribution vanishes at $`T=0`$. This follows from the fact that $`𝒬(k_{}^<)`$ is roughly independent of $`\lambda `$. To estimate the temperature dependence of $`F_{SW}`$, note that the prefactor of this term varies as $`T^1`$. The $`k_{}`$-integral of this term depends on the depth and width of the second minimum of $`𝒬(k_{})`$ at $`k_{}=k_{}^>`$. While the width is roughly independent of $`\lambda `$ for given $`h`$, the minimum becomes deeper with the increasing $`\lambda ^1`$ (or decreasing $`T`$). This proves that $`F_{SW}`$ is more negative at low $`T`$ than at high $`T`$, as opposed to $`F_{LW}`$.

For a given value of $`h`$, there exists a special value $`\lambda _{cr}`$ at which the crossover between the two distinctive behaviors of the plate interaction ($`F_{LW}`$ and $`F_{SW}`$) takes place. By requiring $`\frac{}{h}\left(F_{LW}F_{SW}\right)=0`$, we have the following transcendental equation for $`\lambda _{cr}`$:

$$\left(\frac{\pi }{2}\right)\frac{\zeta (3)}{k_{}^>}\frac{\mathrm{e}^{2k_{}^>h}}{h^3}=\frac{1}{\lambda _{cr}^2}_{k_{}k_{}^>}\frac{k_{}^1dk_{}}{[1+k_{}^1\lambda _{cr}^1g(k_{})]^2}.$$
(5)
To solve the trascendental equation, we have chosen $`D=5\AA `$. Fig. 2 describes distinct regimes characterized by the corresponding dominant modes of fluctuations, and the crossover boundaries between them; the regimes where the long-wavelength (LW) and short-wavelength (SW) fluctuations dominate are denoted by LW and SW, respectively. When $`\lambda ^1`$ is smaller than $`2.7\AA ^1`$, the plate interaction is solely determined by the LW fluctuations for the whole range of $`h`$. At $`\lambda ^12.7\AA ^1`$, marked by the vertical dotted line on the left, the SW fluctuations start to contribute to the plate interaction. When $`2.7\AA ^1\lambda ^17.6\AA ^1`$, however, the plate interaction is determined by the competition between the two; the crossover from the LW to SW regime takes places for larger value of $`h`$ at low temperatures (corresponding to larger $`\lambda ^1`$). At $`\lambda ^1=\lambda _X^17.6\AA ^1`$ , marked by the vertical dotted line on the right, the crossover occurs only when $`h\mathrm{}`$. Beyond $`\lambda _X^1`$, the SW fluctuations solely determine the plate interaction.

The appearance of two distinctive competing interactions, ı.e., $`F_{LW}`$ and $`F_{SW}`$, can also be understood in terms of in-plane charge correlations for a single plate: $`G_1(𝐫_{},𝐫_{}^{})=\widehat{\sigma }(𝐫_{})\widehat{\sigma }(𝐫_{}^{})\widehat{\sigma }(𝐫_{})\widehat{\sigma }(𝐫_{}^{})`$. The long-wavelength contribution to the charge correlation was shown to scale as $`G_{LW}(𝐫_{},𝐫_{}^{})k_BT/|𝐫_{}𝐫_{}^{}|^3`$ . Note that this correlation vanishes as $`T0`$. Our result in Eq. 1 implies that the short-wavelength charge correlation function is given by

$`G_{SW}(𝐫_{},𝐫_{}^{})g(𝐫_{},𝐫_{}^{})`$ (6)
$`{\displaystyle \frac{\sigma _0e^2}{2\pi }}{\displaystyle _{k_{}0}}{\displaystyle \frac{k_{}dk_{}g(𝐤_{})}{1+\frac{\lambda k_{}}{g(𝐤_{})}}}J_0(k_{}|𝐫_{}𝐫_{}^{}|),`$ (7)
where $`J_0(x)`$ is the zeroth-order Bessel function of the first kind. Unlike the long-wavelength correlation, the behavior of the short-wavelength correlation is determined by the nature of the poles of $`\left[1+\lambda k_{}/g(𝐤_{})\right]^1`$. We find that $`G_{SW}`$ shows an oscillatory decay, and that the amplitude of $`G_{SW}`$ varies as $`T^1\left[1+(k_{}^>\lambda )^1g(k_{}^>)\right]^1`$ and is larger at low temperatures. For the two-plate case, this oscillatory charge correlation essentially gives rise to the SW interaction.

At high temperatures, we can consider each plate to consist of large domains, ı.e., counterion-rich and counterion-poor domains. The size of the domains is on the order of $`(k_{}^<)^11\AA `$, and thus these domains can form huge dipoles, resulting in a long-ranged attraction. The long-wavelength fluctuations couple over many plates, leading to breakdown of pairwise additivity  in the $`N`$-plate case. We also find that the power-law behavior of the in-plane charge correlation crosses over to an exponential decaying form as $`N\mathrm{}`$ . This is clearly suggestive of the breakdown of pairwise additivity. As the temperature decreases, however, the plate interaction driven by long-wavelength fluctuations crosses over to a distinct behaviour controlled by non-monotonically decaying charge correlations. At low temperatures, each domain becomes overall charge neutral and thus the distiction between domains is meaningless. In this case, the local correlation between a counterion and a backbone charge in its neighborhood dominates the free energy. There is thus strong cancellation of repulsions (between like charges) with attractions (between opposite charges). This results in an exponentially decaying, short-ranged attraction between the plates.

In conclusion, we have presented a systematic approach to study the nature of counterion-mediated attractions between fluid membranes. We have shown that the nature of these attractions is determined by the dominant modes of fluctuations in the density of counterions. At high temperatures, fluctuations at large length scales determine the membrane interactions; the resulting interactions are long-ranged and not pairwise additive. Charge densities of biomembranes range from $`0.03`$ to $`0.24/\mathrm{nm}^2`$, corresponding to LW regimes at room temperature. It is thus clear that many-body, non-pairwise additive interactions operate in biomembrane systems at room temperature. At low temperatures, however, the membrane interactions are dominated by SW charge fluctuations and are exponentially decaying with the membrane spacing. In this case, only the nearest pairs of plates couple strongly with each other. Surprisingly, this implies that the pairwise additivity is restored at very low temperatures (the non-pairwise additive interaction becomes smaller, and eventually vanishes, as $`T0`$). The approach presented here allows one to systematically study the crossover, as the temperature decreases, from the high-temperature, long-ranged attractions to the behaviors expected at low temperatures.

We have benifited from illuminating discussions with A.W.C. Lau and A.J. Liu. We are also grateful to W.M. Gelbart, H. Schiessel, and R. Menes for scientific stimulation. We thank M. Howard for reading the manuscript carefully, and D. Boal, M. Wortis, and M. Plischke for providing an excellent research enviornment. This work was in part supported by the NSERC of Canada.


text:
# Commensurate and Incommensurate Structure of the Neutron Cross Section in LaSrCu and YBaCuO

## Abstract

We study the evolution of the $`d`$-wave neutron cross-section with variable frequency $`\omega `$ and fixed $`T`$ (below and above $`T_c`$) in two different cuprate families. The evolution from incommensurate to commensurate to incommensurate peaks is rather generic within an RPA-like scheme. This behavior seems to be in reasonable accord with experiments, and may help distinguish between this and the “stripe” scenario.

The goal of the present paper is to address the frequency evolution of the neutron cross section, over the entire range of $`\omega ,𝐪`$, using a scheme which we have previously applied to the normal and $`d`$-wave superconducting states. This work is viewed as significant because it leads to a fairly generic frequency evolution, which seems to be observed experimentally in two cuprate families. These calculations, which have no adjustable parameters (besides those which were used to fit the normal state), can help establish whether the details of the fermiology plus $`d`$-wave pairing can account for the observed incommensurabilities and their evolution with frequency, or whether (by default) some new phenomenon such as stripe  or other exotic phases may be required.

Our starting point is a three band, large $`U`$ calculation which yields a dynamical susceptibility $`\chi (q,\omega )=\chi ^o(q,\omega )/[1+J(q)\chi ^o(q,\omega )]`$ where the Lindhard function $`\chi ^o`$ is appropriate to the ($`d`$-wave) superconducting state and the underlying Fermi surface. Here the residual exchange $`J(q)=J_o[\mathrm{cos}q_x+\mathrm{cos}q_y]`$ arises from Cu-Cu interactions via the mediating oxygen band. While the YBaCuO system is a two layer material, our past experience has shown that most of the peak structures associated with the neutron cross section are captured by an effective one layer band calculation, which we will investigate here. For definiteness, we fix the temperature at $`4`$ K and vary frequency in increments of a few meV. We take the electronic excitation gap to be described by an ideal $`d`$-wave, $`\mathrm{\Delta }(q)=\mathrm{\Delta }(\mathrm{cos}q_x+\mathrm{cos}q_y)`$, where at $`T=4`$ K, $`\mathrm{\Delta }`$ is taken to be $`17`$ meV in YBaCuO6.6 and $`8`$ meV in optimally doped LaSrCuO.

In Figure 1a we show the frequency evolution of an YBaCuO6.6 sample. A node-to-node peak appears at low frequencies but the magnitude is small compared to that of all of the other features shown. Incommensurate peaks at $`(\pi ,\pi \pm \delta ),(\pi \pm \delta ,\pi )`$ are first seen at around $`\omega \mathrm{\Delta }`$, albeit their magnitude is somewhat smaller than found experimentally. As frequency increases there is a clear trend: the incommensurability decreases continuously; this decrease is most apparent, in the immediate vicinity of the resonance frequency. At $`27`$ meV the resonant, commensurate peak is now well established. As can be seen, there is some fine structure near $`(\pi ,\pi )`$ which is a remanent of the incommensurate peaks. And there is a pronounced evolution in the peak shape and height above resonance in the underdoped YBaCuO sample. The $`(\pi ,\pi )`$ regime is a flat-topped, possibly weakly incommensurate peak, just after resonance. It then broadens and remains structureless (as shown by our $`34`$ meV plot) between $`3448`$ meV. Finally, above $`48`$ meV, clear incommensurate structure appears.

Figure 1b shows the frequency evolution for the LaSrCuO family, here shown at optimal doping. Node-to-node structures are seen at low frequencies with very small amplitude. A four peaked structure is seen at frequencies just around the gap frequency, $`8`$ meV. These peaks are sharper, but in roughly the same position as their counterparts in the normal state. These peaks persist (with slightly growing amplitude) until around $`12`$ meV, at which point the incommensurability seems to decrease while their overall magnitude increases. In between $`14`$ meV and $`2\mathrm{\Delta }=16`$ meV is an interesting (evidently “resonant”) structure, but there have been, thus far, no reports of this resonance. However, it should be stressed, that if there is a commensurate feature in LaSrCuO, it should be seen only over a very narrow ( $`2`$ meV) frequency window. Finally, just beyond $`2\mathrm{\Delta }`$, as shown by the last panel at $`\omega =19`$ meV, the cross section becomes very similar in shape to its normal state counterpart, although the magnitude is larger.

This evolution from incommensurate to commensurate peaks and then back to incommensurate behavior with increasing $`\omega `$ can, thus, be seen to be the case for both cuprates. There are claims for these effects in recent experiments on YBaCuO and possibly in BiSSCO. What appears to be different between our observations and these particular experiments is that we do not find distinct two energy scales $`E_c`$, where the incommensurate peaks merge and $`E_r`$, where the resonance occurs. Our calculations for the reduced YBaCuO case, suggest that the incommensurate peaks probably never merge, but rather that at resonance the $`(\pi ,\pi )`$ feature fills in the gap between the two incommensurate features. In this sense $`E_c`$ and $`E_r`$ are the same frequency, although it should be stressed that here we have incorporated no resolution limiting effects.

This work was supported by the NSF under awards No. DMR-91-20000 (through STCS) and No. DMR-9808595 (through MRSEC).


text:
# Damped Ly𝛼 Systems in Semi-Analytic Models: Sensitivity to dynamics, disk properties, and cosmology.

## 1. Introduction

Previously (Maller et al. 1999, hereafter MSPP) we have shown that it is possible to account for the kinematic properties of damped Lyman alpha systems (DLAS) as measured by Prochaska & Wolfe (1997, 1998, hereafter PW97 and PW98) in the context of semi-analytic models (SAMs) (Somerville & Primack 1999, hereafter SP99). In these models, hierarchical structure formation is approximated by constructing a merger tree for each dark matter halo. A natural consequence is that every virialized halo may contain not only a central galaxy, but also a number of satellite galaxies as determined by its merging history. Thus the kinematics of the DLAS arise from the combined effects of the internal rotation of gas disks and the motions between gas disks within a common halo. Here we investigate the sensitivity of this model to some of the assumptions made in MSPP, including the modeling of satellite dynamics, the scale height of the gas, and the cosmology.

## 2. Satellite Dynamics

In the SAMs, a merger tree represents a halo’s growth through the mergers of smaller halos (see Somerville 1997). When halos merge, the central galaxy of the largest progenitor halo becomes the central galaxy of the new halo, and all other galaxies become satellites. These satellites then fall in towards the central galaxy due to dynamical friction, and eventually merge with it (see SP99). Because the treatment of the dynamics of the satellites is necessarily simplified in the usual semi-analytic spirit, and since the kinematics of the DLAS arise from both the rotation of disks and the motions of satellite galaxies in the common halo, it is important to test whether our results are sensitive to the details of our modeling.

In the models presented in MSPP, all satellites were assumed to start falling in from the virial radius of the newly formed halo, and all satellites were assumed to be on circular orbits. These assumptions will tend to maximize the dynamical friction timescale of the satellites, and correspond to the assumptions made in some earlier versions of the SAMs (e.g. Somerville 1997). In the models presented here, we have modified these assumptions in two ways. We initially place the satellites at one half of the virial radius, and assign the orbits from a random distribution as observed in N-body simulations by Navarro, Frenk, & White (1995). They found that the “circularity” parameter $`ϵ`$ has a flat distribution, where $`ϵJ/J_c`$ is defined as the ratio of the angular momentum of the satellite to that of a circular orbit with the same energy. The dynamical friction time scale is then scaled by a factor $`ϵ^{0.78}`$ (Lacey & Cole 1993). This causes some of the satellites to fall in faster than in our previous modeling. In addition, we have changed from the SCDM cosmology used in MSPP to a more fashionable $`\mathrm{\Omega }_0=0.4`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.6`$ cosmology. These ingredients are compatible with the models that were shown to produce good agreement with both local galaxy observations (SP99) and the high redshift Lyman-break galaxies (Somerville, Primack, & Faber 1999).

We find that the combined effect of these changes has a negligible effect on our results. Although we do see fewer satellites within a halo of a given mass, the number of satellites in the inner part of the halo is similar. These are the satellites most likely to give rise to multiple hits, which as we argued in MSPP is the crucial factor in matching the observed kinematics of the DLAS. As in MSPP, we still find that the gaseous disks must have very large radial extents in order to match the observations.

## 3. Disk Thickness

Another important feature of the models is the assumed distribution of the gas within the disks. In MSPP, we investigated several radial profiles and obtained the best results with an assumed gas distribution of the form

$$n(R,z)=\frac{N_t}{2h_z}\frac{R_t}{R}\mathrm{exp}\frac{|z|}{h_z}(R<R_t)$$
(1)
where the truncation density $`N_t`$ is an adjustable parameter. The vertical scale height of the gas disks is another uncertainty. In MSPP we assumed it to be one tenth of the stellar disk scale radius, i.e. $`h_z=0.1R_{}`$ where $`R_{}`$ is the stellar disk scale length as given by the SAMs (see SP99). Because the gas disks in the successful models tend to have a large radial extent, much larger then the stellar disk, this leads to very thin gaseous disks. When the scale height is increased to half the stellar scale radius ($`h_z=0.5R_{}`$), we are able to use more physically plausible values for the truncation density $`N_t`$ (see Table 1), and the gas disks are now typically contained within the truncation radius of the dark matter sub-halos surrounding the satellites. This model therefore seems to be more physical, yet still produces good agreement with all four of the diagnostic statistics of PW97 (see Table 1). Note that PW98 have shown that the effect of including a warp in the gas disk is the same as increasing its thickness, so our exploration of thicker disks can also be thought of as including warps in the disk.

Figure 2 shows the distribution of thicknesses and inclinations for the disks that contribute to DLAS in the two models. We find that the disks that produce damped systems in model B are more likely to have a face-on geometry. With thinner disks, more of the cross section comes from an inclined geometry as $`nh_x^1`$. The total cross section in the two models is roughly the same; in model A, more of it comes from extended highly inclined disks, while in model B, denser, face-on disks are more important.

Thus the need for gas disks with large radial extent can be reduced by using thicker disks, though not by a large amount. Increasing the thickness by a factor of five only reduces the radial truncation value by $`30\%`$. The gaseous disks in model B are still quite large compared to the stellar component.

## 4. Conclusions

In MSPP we argued that the observed kinematics of DLAS can be reproduced in hierarchical models if a significant fraction of the lines of sight pass through multiple disks orbiting within a common dark matter halo. However, we found that in order to obtain a high enough cross section for multiple hits, we had to assume very large radial extents for the gaseous disks in our models. We have tested the robustness of these conclusions by modifying several of the uncertain ingredients of our models. We find that modifying the dynamical friction timescale of the satellites by assuming a different intial radius or orbit has a small effect on our results. Similarly, our results do not seem to be sensitive to the assumed cosmology. Increasing the vertical scale height of the disks has a larger effect and leads to models that are more physically plausible and still produce good agreement with the diagnostic statistics of PW. Our results suggest that in order to reconcile the observed kinematics of DLAS with hierarchical theories of structure formation, gaseous disks at high redshift must be large in radial extent and thickened or warped. The physical cause of these properties remains obscure, however we speculate that tidal encounters or outflows could be responsible.


text:
# Distortion of leaves in product foliations

## 1. Introduction

It is a basic problem, given a foliation of a Riemannian manifold, to compare the extrinsic and the intrinsic geometry of leaves of the foliation. For leaves which are compact, this amounts to studying the distortion of $`\pi _1`$ of the leaf in $`\pi _1`$ of the ambient manifold; much is known about this problem. For a partial survey, see .

For taut codimension one foliations of $`3`$–manifolds, one knows by a result of Sullivan in  that one can choose a metric on the ambient manifold such that the leaves are minimal surfaces. One could think of this as saying that taut foliations “measure area well”. However, leaves of taut foliations are far from being quasi–isometrically embedded, in general. For instance, a theorem of Fenley in  states that for an $``$–covered foliation $``$ of a hyperbolic $`3`$–manifold $`M`$, leaves of the pulled–back foliation $`\stackrel{~}{}`$ of $`\stackrel{~}{M}`$ limit to the entire sphere at infinity of $`\stackrel{~}{M}=^3`$. A foliation is $``$–covered if the pulled–back foliation $`\stackrel{~}{}`$ is topologically the product foliation of $`^3`$ by horizontal $`^2`$’s. Fenley’s proof is attractive but somewhat lengthy, and therefore we think it is useful to give a simple proof of this fact here:

###### Theorem 1.1.

For $``$ an $``$–covered foliation of $`M`$, a finite volume hyperbolic $`3`$–manifold, every leaf $`\lambda `$ of $`\stackrel{~}{}`$ limits to the entire sphere at infinity $`S_{\mathrm{}}^2`$ of $`\stackrel{~}{M}=^3`$.

Proof: For a leaf $`\lambda `$ of $`\stackrel{~}{}`$, let $`\lambda _{\mathrm{}}=\overline{\lambda }\lambda S_{\mathrm{}}^2`$ be the set of points in $`S_{\mathrm{}}^2`$ in the closure of $`\lambda `$. We claim that $`\lambda _{\mathrm{}}=S_{\mathrm{}}^2`$ for each $`\lambda `$.

Suppose some $`\lambda _{\mathrm{}}`$ omits some point $`pS_{\mathrm{}}^2`$, and therefore an open disk $`U`$ containing $`p`$. Let $`U^{}`$ and $`U^{\prime \prime }`$ be two open disks in $`S_{\mathrm{}}^2`$ whose closures are disjoint and contained in $`U`$. Since $`M`$ has finite volume, there are elements $`\alpha ,\beta \pi _1(M)`$ such that $`\alpha (S_{\mathrm{}}^2U)U^{}`$ and $`\beta (S_{\mathrm{}}^2U)U^{\prime \prime }`$. Let $`\lambda ^{}=\alpha (\lambda )`$ and $`\lambda ^{\prime \prime }=\beta (\lambda )`$. Since the leaf space of $`\stackrel{~}{}`$ is $``$, After relabeling the leaves if necessary, we can assume $`\lambda `$ separates $`\lambda ^{}`$ from $`\lambda ^{\prime \prime }`$. However, there is clearly an arc $`\gamma `$ in $`S_{\mathrm{}}^2`$ joining $`\lambda _{\mathrm{}}^{}`$ to $`\lambda _{\mathrm{}}^{\prime \prime }`$ which avoids $`\lambda _{\mathrm{}}`$. $`\gamma `$ is a Hausdorff limit (in $`^3S_{\mathrm{}}^2`$) of arcs $`\gamma _i^3`$ running between $`\lambda ^{}`$ and $`\lambda ^{\prime \prime }`$. Each $`\gamma _i`$ must intersect $`\lambda `$ in some point $`q_i`$, and by extracting a subsequence, we find $`q_iq\gamma `$. But by construction, $`q\lambda _{\mathrm{}}`$, giving us a contradiction.

If we want to study the distortion of leaves in taut foliations, then $``$–covered foliations are the most delicate case. For, if $``$ is a taut foliation which is not $``$–covered, then the leaf space of $`\stackrel{~}{}`$ is non–Hausdorff. That is, there are a sequence of pairs of points $`p_i,q_i\stackrel{~}{M}`$ with $`p_ip`$ and $`q_iq`$ such that $`p_i,q_i\lambda _i`$ but $`p`$ and $`q`$ lie on different leaves of $`\stackrel{~}{}`$. It follows that the distance between $`p_i`$ and $`q_i`$ as measured in $``$ is going to infinity; that is, $`d_{\lambda _i}(p_i,q_i)\mathrm{}`$. On the other hand, the distance between $`p_i`$ and $`q_i`$ as measured in $`\stackrel{~}{M}`$ is bounded: $`d_{\stackrel{~}{M}}(p_i,q_i)d_{\stackrel{~}{M}}(p,q)`$.

By contrast, leaves of $``$–covered foliations are never infinitely distorted. In fact in  we show the following:

###### Theorem 1.2.

Let $``$ be an $``$ covered foliation of an atoroidal $`3`$–manifold $`M^3`$. Then leaves $`\lambda `$ of $`\stackrel{~}{}`$ are quasi–isometrically embedded in their $`t`$–neighborhoods, for any $`t`$. That is, there is a uniform $`K_t,ϵ_t`$ such that for any leaf $`\lambda `$ of $`\stackrel{~}{}`$, the embedding $`\lambda N_t(\lambda )`$ is a $`(K_t,ϵ_t)`$ quasi–isometry.

A heuristic scheme to measure the distortion of a leaf $`\lambda `$ runs as follows.

We assume that $``$ is co–oriented, so that in the universal cover there is a well–defined notion of the space “above” and the space “below” a given leaf. Since leaves are minimal surfaces, the mean curvature is zero, so there will be well–defined approximate directions at every point in which the leaf “bends upwards” and approximate directions in which it “bends downwards”. For an arc $`\alpha `$ between points $`p,q`$ of a leaf $`\lambda `$ which is roughly parallel to a direction of positive extrinsic curvature, the geodesic $`\alpha ^+`$ in $`\stackrel{~}{M}`$ running between $`p`$ and $`q`$ will lie mostly above $`\lambda `$. We can cap off the circle $`\alpha \alpha ^+`$ with a disk $`D_\alpha `$ of minimal area, and consider its projection to $`M`$. The intersection of $`D_\alpha `$ with $`\stackrel{~}{}`$ gives a foliation of $`D_\alpha `$ by arcs parallel to the principal directions of positive curvature in the leaves of $`\stackrel{~}{}`$. Thus the projection of $`D_\alpha `$ to $`M`$ should be close to an embedding, since its self–intersections cannot have a large dihedral angle. If we consider longer and longer leafwise geodesics $`\alpha `$, the geometric limit of the disks $`D_\alpha `$ in $`M`$ should be a geodesic lamination $`\mathrm{\Lambda }^+`$ in $`M`$ which describes the “eigendirections” of positive curvature of $``$. If we consider the principal directions of negative curvature, we should get a complementary lamination $`\mathrm{\Lambda }^{}`$, transverse to $``$ and to $`\mathrm{\Lambda }^+`$, which describes the “eigendirections” of negative curvature of $``$.

To study the distortion of leaves of $``$ in such a setup, one need only study, for a typical leaf $`l`$ of $`\mathrm{\Lambda }^+`$ say, how distorted the induced one–dimensional foliation $`l`$ is in $`l`$. Such a foliation has bounded geometry — its extrinsic curvature is bounded above and below by some constant — because of the compactness of $`M`$. Moreover it is topologically a product. The point of this paper is to show that these two conditions in no way allow one to establish any bound on the distortion function of $``$.

One remark worth making is that a pair of laminations $`\mathrm{\Lambda }^\pm `$ transverse to an $``$–covered foliation of an atoroidal $`3`$–manifold $`M`$ are constructed in  and independently by Sérgio Fenley (). The interpretation of these foliations as eigendirections of extrinsic curvature is still conjectural, however.

## 2. Non–recursive distortion

###### Definition 2.1.

The distortion function $`𝒟`$ of a foliation $``$ in a Riemannian manifold $`M`$ is defined as follows: if $`d_{}`$ denotes the intrinsic distance function in a leaf of $``$ thought of as a geodesic metric space, and $`d_{\text{amb}}`$ denotes the Riemannian distance function in the ambient space, then

$$𝒟(t)=\underset{x,y\lambda |d_{\text{amb}}(x,y)=t}{sup}d_{}(x,y)$$

where the supremum is taken over all pairs of points in all leaves of $``$. We usually expect that $`M`$ is simply connected.

###### Theorem 2.1.

Let $`\lambda `$ be a smooth oriented bi–infinite ray properly immersed in $`^2`$.

1. If $`1\kappa 1`$ everywhere, then the distortion function $`𝒟`$ is at most exponential.
2. If $`\kappa 1`$ everywhere and $`>1`$ somewhere, then $`\lambda `$ has a self–intersection.
3. If $`\kappa =1`$ everywhere then $`\lambda `$ is a horocircle and therefore has exponential distortion.

Proof: The proof of 1. is a standard comparison argument, and amounts to showing that $`\lambda `$ makes at least as much progress as a horocircle.

To see that 2. holds, we consider the progress of the osculating horocircle to $`\lambda `$. At a point $`x`$ on $`\lambda `$, the osculating horocircle $`H_x`$ is the unique horocircle on the positive side of $`\lambda `$ which is tangent to $`\lambda `$ at $`x`$. Let $`B_xS_{\mathrm{}}^1`$ be the basepoint of $`H_x`$. Then the curvature condition on $`\lambda `$ implies that as one moves along $`\lambda `$ in the positive direction, $`B_x`$ always moves anticlockwise. Since $`\kappa >1`$ somewhere, the derivative of $`B_x`$ is nonzero somewhere. If we truncate $`\lambda `$ on some sufficiently big compact piece and fill in the remaining segments with horocircles, we get a properly immersed arc $`\lambda ^{}`$ in $`B^2`$ with positive winding number, which consequently has a self–intersection somewhere. If this intersection is in the piece agreeing with $`\lambda `$, we are done. Otherwise, the curvature condition implies that $`\lambda `$ must have an intersection in the region enclosed by $`\lambda ^{}`$.

In light of this theorem, it is perhaps surprising that we can make the following construction:

###### Theorem 2.2.

For any $`ϵ>0`$ there is a foliation $``$ of $`^2`$ which is topologically the standard foliation of $`^2`$ by horizontal lines with the following properties:

* Leaves are smooth.
* The extrinsic curvature of any leaf at any point is bounded between $`1ϵ`$ and $`1+ϵ`$.
* The distortion function $`𝒟`$ grows faster than any recursive function.

Proof: We consider the upper half-space model of $`^2`$ and let $`\lambda `$ be the leaf passing through the point $`i`$. $`\lambda `$ will be the graph of a function $`r=\varphi (\theta )`$ in polar co–ordinates for $`\theta (0,\pi )`$. Choose some very small $`\delta `$. Then for $`\pi \delta >\theta >\delta `$ we let $`\lambda `$ agree with the horocircle with “center” at $`\mathrm{}`$ passing through $`i`$. Let $`r_n`$ be a sequence of positive real numbers which grows faster than any recursive function. Then we define

$$\varphi \left(\frac{\delta }{2^n}\right)=r_n$$

Since $`r_n`$ grows so very fast, one can easily choose $`\varphi `$ to interpolate between $`\frac{\delta }{2^n}`$ and $`\frac{\delta }{2^{n+1}}`$ so that $`\lambda `$ is very close to a radial (Euclidean) line of very small slope in this range. The extrinsic curvature of such a line in $`^2`$ is very close to $`1`$. Define $`\varphi `$ in the range $`\theta >\pi \delta `$ by $`\varphi (\theta )=\varphi (\pi \theta )`$. Then the extrinsic distance between $`(r,\theta )`$ and $`(r,\pi \theta )`$ depends only on $`\theta `$. It follows that for $`\theta `$ between $`\frac{\delta }{2^n}`$ and $`\frac{\delta }{2^{n+1}}`$ the extrinsic distance between $`(r,\theta )`$ and $`(r,\pi \theta )`$ is bounded above by some recursive function in $`n`$. However, the length of the path in $`\lambda `$ between $`(r,\theta )`$ and $`(r,\pi \theta )`$ is growing approximately like $`\mathrm{ln}(r_n)`$. The distortion function $`𝒟`$ is therefore non–recursive.

Let $`\lambda _1=\lambda `$ and for $`t>0`$, define $`\lambda _t`$ to be the graph of the function $`r=t\varphi (\theta )`$ in polar co-ordinates. Then there is a hyperbolic isometry taking $`\lambda _t`$ to $`\lambda _s`$ for any $`s,t`$ and therefore each of the $`\lambda _t`$ satisfies the same curvature bounds as $`\lambda `$. Moreover, the union of the $`\lambda _t`$ gives a product foliation of $`^2`$, as required.

###### Remark 2.1.

For those unfamiliar with the concept, it should be noted that it is easy to produce a function (even an integer valued function) such as $`r_n`$ which grows faster than any recursive function. For, enumerate all the recursive functions somehow as $`\varphi _n`$. Then define $`r(n)=\mathrm{max}_{mn}\varphi _m(n)`$. Such an $`r`$ grows (eventually) at least as fast as any recursive function, and therefore faster than any recursive function.

Essentially the same construction works for Euclidean space, a fact which we now establish.

###### Theorem 2.3.

For any $`ϵ>0`$ there is a foliation $``$ of $`𝔼^2`$ which is topologically the standard foliation of $`^2`$ by horizontal lines with the following properties:

* Leaves are smooth.
* The extrinsic curvature of any leaf at any point is bounded by $`ϵ`$
* The distortion function $`𝒟`$ grows faster than any recursive function.

Proof: The leaf $`\lambda `$ of $``$ passing through the origin will be the graph of an even function $`y=\varphi (x)`$ and every other leaf will be a translate of $`\lambda `$. After choosing a sufficiently large $`K`$ and small $`\delta `$, we make $`\varphi (x)=\delta x^2`$ for $`|x|<K`$, so that $`\lambda `$ has curvature bounded below $`ϵ`$, and is very nearly vertical at $`x=K`$. Then let $`r_i`$ be a sequence growing faster than any recursive function, as before, and define $`\varphi (K+n)=r_n`$. Then we can choose $`\varphi `$ to interpolate between $`K+n`$ and $`K+n+1`$ for each $`n`$ to make it very smooth and almost straight, since $`\lambda `$ will be almost vertical in this region.

## 3. Minimality

###### Definition 3.1.

A decorated metric space is a metric space in the usual sense with some auxiliary structure. For instance, this structure could consist of a basepoint, some collection of submanifolds, a foliation or lamination, etc. together with a notion of convergence of such structures in the geometric topology, in the sense of Gromov. One says a decorated metric space $`M`$ has bounded geometry if the metric space itself has bounded geometry, and if for every $`t`$, the decorated metric spaces obtained as restrictions of $`M`$ to the balls of radius $`t`$ form a precompact family.

For instance, we might have a family of foliated Riemannian manifolds $`M_i,_i`$. The geometric topology requires a choice of basepoints $`p_i`$ in each $`M_i`$. Then we say the family $`(M_i,_i,p_i)`$ is a Cauchy sequence if there are a sequence of radii $`r_i`$ and $`ϵ_i0`$ so that the Gromov–Hausdorff distance between the ball of radius $`r_i`$ in $`M_i`$ and $`M_j`$ about $`p_i,p_j`$ for $`j>i`$ is at most $`ϵ_i`$, and that such “near–isometries” can be chosen in a way that the leaves of $`_i`$ can be taken $`ϵ_i`$–close to the leaves of $`_j`$. From such a Cauchy sequence we can extract a limit $`(M,,p)`$.

###### Definition 3.2.

A decorated metric space $`M`$ with bounded geometry is minimal if for every limit $`(M^+,p)`$ of the pointed decorated metric spaces $`(M,p_i)`$ and every point $`qM`$, there is a sequence $`q_iM^+`$ such that the pointed decorated metric spaces $`(M^+,q_i)(M,q)`$.

The foliations constructed in the last section are certainly not minimal. We ask the question: does the assumption of minimality allow one to make some estimate on the distortion function for a product foliation?


text:
# Rotations and Abundances of Blue Horizontal-Branch Stars in Globular Cluster M15Based in large part on observations obtained at the W.M. Keck Observatory, which is operated jointly by the California Institute of Technology and the University of California

## 1 Introduction

M15 (NGC 7078) is one of the most metal-poor globulars known, with a metallicity $`[\mathrm{Fe}/\mathrm{H}]2.4`$ dex measured from red giant abundances (Cohen 1979, Sneden et al. 1997). Like many other such clusters, M15’s horizontal branch lies predominantly bluewards of the instability strip, and color-magnitude diagrams (Buonanno et al. 1983; Durrell & Harris 1993) show an extended “blue tail” reaching $`T_{\mathrm{eff}}`$ as high as $`20000\mathrm{K}`$, which is separated from the horizontal part of the HB by a “gap” in the distribution of stars along the HB (Moehler et al. 1995). Similar gaps appear in the CMDs of M13, M80, NGC 6752, NGC 288, and other clusters, but are difficult to explain via standard models of RGB mass loss or HB evolution.

Recently, attention has focused on atmospheric effects as a possible explanation for these photometric features (Caloi 1999; Grundahl et al. 1998). At $`T_{\mathrm{eff}}10000\mathrm{K}`$, they suggest, the stellar atmospheres become susceptible to diffusion effects, and thus develop chemical peculiarities similar to those that appear in main-sequence CP stars. The resulting changes in atmospheric opacity alter the emitted flux distributions and thus the measured photometry of the hotter stars, giving rise to the gaps. These claims have been bolstered by measurements of large photospheric abundance anomalies among hotter BHB stars in M13 (Behr et al. 1999) and NGC 6752 (Moehler et al. 1999), which show 30 to 50 times the iron abundance expected for these metal-poor clusters. In M13, the transition from normal-metallicity cooler stars to metal-enhanced hotter stars is remarkably abrupt, and coincides with the location of the gap, further supporting the surface-effect hypothesis.

It appears, however, that stellar rotation also plays some role in this process. Theoretical treatments of the diffusion mechanisms (Michaud et al. 1983) suggest that circulation currents induced by higher rotation velocities can easily prevent abundance variations from appearing. This prediction is borne out by measurements of $`v\mathrm{sin}i`$ for the M13 stars (Behr et al. 2000), which show that although the cooler stars exhibit a range of $`v\mathrm{sin}i`$, some as high as $`40\mathrm{km}\mathrm{s}^1`$, the hot metal-enhanced stars all show very low rotations, $`v\mathrm{sin}i<10\mathrm{km}\mathrm{s}^1`$. This correlation suggests that slow rotation may be required in order for the metal enhancement and helium depletion appear in the photosphere.

Although these results from M13 and NGC 6752 imply that we are on the right track towards explaining the photometric peculiarities, BHB stars in many other clusters will have to be analyzed in a similar fashion before we can make any firm claims. In particular, since the radiative levitation that causes the observed metal enhancements is thought to depend strongly on the intrinsic metallicity of the atmosphere, we should study clusters spanning a range of \[Fe/H\], to see whether the onset and magnitude of the enhancements vary. To this end, we have observed BHB stars in five clusters in addition to M13. In this Letter, we describe the results for both rotation and photospheric abundances in M15. The data for the other clusters — M92, M68, NGC 288, and M3 — are being analyzed currently, and will be presented in a future publication.

## 2 Observations and Reduction

The eighteen stars in our sample were selected from Buonanno et al. (1983), and are listed in Table 1. We acquired supplementary Strömgren photometry of the target stars using the Palomar 60-inch telescope, in order to better constrain the effective temeratures. The stars generally lie in the cluster outskirts, where crowding and confusion are less of a problem than towards the core, and our CCD imaging confirmed an absence of faint companions.

The spectra were collected using the HIRES spectrograph (Vogt et al. 1994) on the Keck I telescope, during four observing runs on 1997 August 01–03, 1997 August 26–27, 1998 June 27, and 1999 August 14–17. A 0.86-arcsec slit width yielded $`R=45000`$ ($`v=6.7\mathrm{km}\mathrm{s}^1`$) per 3-pixel resolution element. We limited frame exposure times to 1200 seconds, to minimize susceptibility to cosmic ray accumulation, and then coadded four frames per star. $`S/N`$ ratios were on the order of $`3060`$ per resolution element.

We used a suite of routines developed by J.K. McCarthy (1988) for the FIGARO data analysis package (Shortridge 1988) to reduce the HIRES echellograms to 1-dimensional spectra. Frames were bias-subtracted, flat-fielded against exposures of HIRES’ internal quartz incandescent lamps (thereby removing much of the blaze profile from each order), cleaned of cosmic ray hits, and coadded. A thorium-argon arc lamp provided wavelength calibration. Sky background was negligible, and 1-D spectra were extracted via simple pixel summation. A 10th-order polynomial fit to line-free continuum regions completed the normalization of the spectrum to unity.

## 3 Analysis

The resulting spectra show a few to over 140 metal absorption lines each, with the hottest stars showing the largest number of lines. Line broadening from stellar rotation is evident in a few of the stars, but even in the most extreme cases, the line profiles were close to Gaussian, so line equivalent widths ($`W_\lambda `$) were measured by least-square fitting of Gaussian profiles to the data. Equivalent widths as small as 10 mÅ were measured reliably, and errors in $`W_\lambda `$ (estimated from the fit $`\chi ^2`$) were typically 5 mÅ or less. Lines were then matched to the atomic line lists of Kurucz & Bell (1995), and those that were attributed to a single species (i.e. unblended) were used to determine radial velocity $`v_r`$ for each of the stars. On the basis of $`v_r`$, all of the targets appear to be cluster members.

To derive photospheric parameters $`T_{\mathrm{eff}}`$ and $`\mathrm{log}g`$, we compared the published photometry and our own Strömgren indices to synthesized colors from ATLAS9 (Kurucz 1997), adopting a cluster reddening of $`E(BV)=0.09`$. Temperatures are well-determined ($`\pm `$ a few hundred K) for the cooler stars, but are somewhat less firm (as much at $`\pm 1000\mathrm{K}`$) for the hotter stars. This situation will improve with a more sophisticated treatment of the various Strömgren colors, and by using transitions with different excitation potentials $`\chi `$ to constrain $`T_{\mathrm{eff}}`$. We estimated surface gravities using the $`AB_\nu `$ flux method (Oke & Gunn 1984), which relates $`\mathrm{log}g`$, stellar mass $`M_{}`$, distance $`d`$, and photospheric Eddington flux $`H_\nu `$ at 5480 Å. Since $`AB_\nu (5480\AA )=V`$ magnitude, we can derive

$$\mathrm{log}g=3.68+\mathrm{log}(M_{}/M_{})+\mathrm{log}(H_\nu )(Mm)_V+0.4V_0,$$

where unextincted magnitude $`V_0=V+A_V=V+0.30`$ for M15. We assumed $`M_{}=0.6M_{}`$ as a representative BHB star mass, used a distance modulus of 15.26 (Silbermann & Smith 1995), and drew the $`H_\nu `$ values from the ATLAS9 model atmospheres (Kurucz 1997), iterating until $`\mathrm{log}g`$ converged. The resulting $`\mathrm{log}g`$ values agree well with model ZAHB tracks (Dorman 1993), except for the hotter stars, which are “overluminous” as described by Moehler et al. (1995). Table 1 lists the final photospheric parameters used for each of the target stars, as well as the heliocentric radial velocities.

For the chemical abundance analyses, we use the LINFOR/LINFIT line formation analysis package, developed at Kiel, based on earlier codes by Baschek, Traving, and Holweger (1966), with subsequent modifications by M. Lemke. Our spectra of these very metal-poor stars are sufficiently uncrowded that we can simply compute abundances from equivalent widths, instead of performing a full spectral synthesis fit. Only lines attributed to a single chemical species were considered; potentially blended lines are ignored in this analysis. Upper bounds on \[Ti/H\] were determined for two of the hotter stars, which did not show any titanium lines, by assigning an equivalent width of 20 mÅ (double the strength of the weakest metal lines actually measured in those spectra) to the six strongest Ti II transitions and calculating the implied abundance. Microturbulent velocity $`\xi `$ was chosen such that the abundance derived for a single species (Fe II for most cases) was invariant with $`W_\lambda `$. For those stars with an insufficient number of lines to utilize this technique, we adopted a typical value of $`\xi =2\mathrm{km}\mathrm{s}^1`$. We assumed a cluster metallicity of \[Fe/H\] $`=2.4`$ dex in computing the initial model atmospheres, and then adjusted as necessary for those stars that turned out to be considerably more metal-rich, although these adjustments to the atmospheric input were found to have only modest effects ($`<0.2`$ dex) on the abundances of individual elements.

Measurement of $`v\mathrm{sin}i`$ broadening traditionally entails cross-correlation of the target spectrum with a rotational reference star of similar spectral type, but this approach assumes that the template star is truly at $`v\mathrm{sin}i=0`$, which is rare. Furthermore, given the abundance peculiarities that many of these BHB stars exhibit, it is difficult to find an appropriate spectral analog. Since we are able to resolve the line profiles of our stars, we instead chose to measure $`v\mathrm{sin}i`$ by fitting the profiles directly, taking into account other non-rotational broadening mechanisms. We used bright unsaturated arc lines to construct an instrumental profile, and then combined it with an estimated thermal Doppler broadening of $`3\mathrm{km}\mathrm{s}^1`$ FWHM and the previously-determined $`\xi `$. This profile was convolved with hemicircular rotation profiles for various $`v\mathrm{sin}i`$ to create the final theoretical line profiles. Each observed line in a spectrum was fit to the theoretical profile using an iterative least-squares algorithm, solving for $`v\mathrm{sin}i`$ and line depth. The values for $`v\mathrm{sin}i`$ from different spectral lines generally agree quite well, once we removed helium lines and blended lines such as the Mg II 4481 triplet.

## 4 Results

In Figure 1, abundance determinations for key chemical species are plotted as a function of stellar $`T_{\mathrm{eff}}`$. The values \[X/H\] represent logarithmic offsets from the solar values of Anders & Grevesse (1993). The error bars incorporate the scatter among multiple lines of the same species, plus the uncertainties in $`T_{\mathrm{eff}}`$, $`\mathrm{log}g`$, $`\xi `$, $`W_\lambda `$ for each line, and \[Fe/H\] of the input atmosphere.

In the cooler stars ($`T_{\mathrm{eff}}<10000\mathrm{K}`$), the compositions are largely as expected. The iron abundance \[Fe/H\] averages $`2.5`$, slightly below the value of $`2.4`$ expected for this metal-poor cluster. Magnesium and titanium appear at \[Mg/H\] $`2.2`$ and \[Ti/H\] $`1.8`$, respectively, which are also reasonable for this environment.

As we move to the hotter stars, however, the iron abundances change radically. Six of the eight stars at $`T_{\mathrm{eff}}10000\mathrm{K}`$ show solar iron abundances, \[Fe/H\] $`0`$, an enhancement of a factor of 300 from the cooler stars. (The other two hot stars show the same \[Fe/H\] as the cool stars, for reasons which will be discussed shortly.) Titanium, in a similar fashion, rises to \[Ti/H\] $`>0`$, although there are hints of a monatonic increase with $`T_{\mathrm{eff}}`$ rather than an abrupt jump. Magnesium, on the other hand, appears to be unaltered, maintaining a roughly constant metal-poor level across the entire temperature range. The hotter stars also start to show helium lines, providing evidence of helium depletion, since \[He/H\] $`0`$ at $`11000\mathrm{K}`$ but it then drops by factors of 10 to 30 for $`T_{\mathrm{eff}}=12000`$$`13000\mathrm{K}`$.

Figure 1 also charts the values of $`v\mathrm{sin}i`$ derived for the target stars. Among the cooler stars, we find a range of rates, with most of the stars rotating at $`15\mathrm{km}\mathrm{s}^1`$ or less, except for two stars at $`29\mathrm{km}\mathrm{s}^1`$ and $`36\mathrm{km}\mathrm{s}^1`$, respectively. This sort of distribution of $`v\mathrm{sin}i`$ is not what one would expect given a single intrinsic rotation speed $`v`$ and random orientation of the rotation axes, since large $`\mathrm{sin}i`$ are more likely than small $`\mathrm{sin}i`$. Instead, the cool end of the HB appears to contain two rotational populations, one with $`v35\mathrm{km}\mathrm{s}^1`$, and another with $`v15\mathrm{km}\mathrm{s}^1`$, much like in M13 (Peterson et al. 1995).

For the hotter stars, there also appears to be a bimodal distribution in $`v`$, although the difference is less pronounced. Six of these eight stars show $`v\mathrm{sin}i<7\mathrm{km}\mathrm{s}^1`$, while the other two have $`v\mathrm{sin}i12\mathrm{km}\mathrm{s}^1`$. Interestingly, these two faster-rotating hot stars are the same stars that show “normal” (metal-poor) iron abundances.

## 5 Discussion

These BHB stars in M15 exhibit abundance and rotation characteristics very similar to those reported for BHB stars in M13 (Behr et al. 1999 & 2000). The enhancements of metals except for magnesium, the depletion of helium, and the difference in maximum $`v\mathrm{sin}i`$ between hotter and cooler stars are shared by both clusters. The abundance anomalies in M15 are therefore likely to be due to the same diffusion processes that were invoked for the prior study — radiative levitation of metals, and gravitational settling of helium, in the stable non-convective atmospheres of the hotter, higher-gravity stars, as hypothesized by Michaud et al. Unfortunately, the stars selected in M15 do not sample the immediate vicinity of the photometric gap as well as those in M13, so the association between the onset of diffusion-driven abundance variations and the gap is not as clear-cut, but the general trend still supports this association.

There are, however, two notable differences between the results for M15 and those for M13. First, the magnitudes of the metal enhancements are somewhat different. Iron and titanium are each enhanced by $`2.5\mathrm{dex}`$ in M15, while in M13, iron increases by only $`2\mathrm{dex}`$, and titanium by $`1.5\mathrm{dex}`$. Despite this difference, the enhancement mechanism yields the same final metallicities for the hot stars in both clusters: \[Fe/H\] $`0`$ and \[Ti/H\] $`0`$. This correspondence suggests that the radiative levitation mechanism reaches equilibrium with gravity at or near solar metallicity, independent of the initial metal content of the atmosphere, as the metal lines become saturated and are thus unable to support further enhancements.

Second, there is the issue of the two hotter M15 stars which do not show metal enhancement, denoted by circles in Figure 1. Their derived temperatures associate them with the hotter population, as does their photometry, which places them bluewards of the photometric gap. Their iron abudances, however, are $`<2.5\mathrm{dex}`$, like the cooler unenhanced stars in the cluster, and stringent upper bounds on their titanium abundances again suggest that they are metal-poor. These two stars are also distinguished by having $`v\mathrm{sin}i12\mathrm{km}\mathrm{s}^1`$, nearly twice as large as any of the other hot stars. Although their temperatures and gravities are high enough to support radiative levitation, it appears that their faster rotations induce meridonal circulation, which keeps the atmosphere well-mixed and prevents the metal enhancements from appearing. Such sensitive dependence upon rotation speed was mentioned by Michaud et al., but these observations provide direct evidence that rotationally-driven mixing can and does directly influence the atmospheric composition.

With these results from M15, we add to the growing body of evidence that element diffusion, regulated by stellar rotation, are at least partially responsible for the observed photometric morphology of globular cluster HBs. The findings from this cluster corroborate the prior work on M13 and NGC 6752, while adding some new twists which further illuminate the diffusion mechanisms. Analysis of many other clusters, spanning a range of metallicity and HB morphology, will be necessary before all of the implications of these effects can be known.

These observations would not have been feasible without the HIRES spectrograph and the Keck I telescope. We are indebted to Jerry Nelson, Gerry Smith, Steve Vogt, and many others for making such marvelous machines, to the W. M. Keck Foundation for making it happen, and to a bevy of Keck observing assistants for making them work. J.G.C. and B.B.B. acknowledge support from NSF Grant AST-9819614.


text:
# Nonlinear analysis of the shearing instability in granular gases

## I Introduction

The understanding of the dynamics of granular fluids is crucial for various industrial processes. This has led to many investigations where theory, experiments and simulations are used in order to construct a predictive theory . There is hope that for moderate densities and slightly dissipative grains, grain dynamics may be described on large scales by using fluid hydrodynamics with slight modifications.

One way to guess which changes to make in standard hydrodynamics is to use the tools developed by statistical mechanics in order to derive hydrodynamic equations. This is justified by the fact that a simple grain model, the inelastic hard sphere model, has been shown to reproduce most of the phenomena occuring in granular systems: in some sense, it has proven to contain the essential ingredients necessary to predict the peculiar physics observed .

The scheme used in the study of nonequilibrium fluids can then be extended to granular fluids: “microscopic” simulations permit to compute the equation of state and values for transport coefficients which are then fed into the guessed macroscopic equations. Comparison between direct nonequilibrium simulations of microscopic and macroscopic models allows then to test the validity of the proposed macroscopic equations. This approach was used recently by two of us to investigate heat transport (heat being identified with the kinetic energy associated with the grains’ motion) and it has been shown that Fourier’s law has to be generalized with a density gradient term appearing in the expression for the heat flux .

In the Inelastic Hard Sphere model (IHS), grains are spherical hard particles with only translational degrees of freedom. The energy dissipation is included through a restitution coefficient $`r`$ lower than one. As for hard spheres, the collision is an instantaneous event and the grain velocities after a collision are given by

$`𝐯_1^{}`$ $`=`$ $`𝐯_1+{\displaystyle \frac{1}{2}}(1+r)\left[\widehat{𝐧}(𝐯_2𝐯_1)\right]\widehat{𝐧}`$ (1)
$`𝐯_2^{}`$ $`=`$ $`𝐯_2{\displaystyle \frac{1}{2}}(1+r)\left[\widehat{𝐧}(𝐯_2𝐯_1)\right]\widehat{𝐧}`$ (2)
where $`\widehat{𝐧}`$ is the unit vector pointing from the center of particle 1 towards the center of particle 2. It is convenient to define the dissipation coefficient $`q=(1r)/2`$ which vanishes when collisions are elastic. In what follows units are chosen such that the disk diameter $`\sigma `$ and the particles masses are set to one.

The IHS model, like the elastic hard sphere model, does not have an intrinsic energy scale. This means that two systems with the same initial configuration and with the particle velocities of one system being equal to those of the other multiplied by one scaling factor will follow the same trajectory but at different speeds. This lack of an intrinsic energy scale also implies a simple temperature dependence for all hydrodynamic quantities which may be found by simple dimensional analysis.

When a fluidized granular medium is let to evolve freely in a box with periodic boundary conditions, the energy decreases continuously in time and the system remains homogeneous. This non-steady state is called the Homogeneous Cooling State (HCS) and it is the reference state from which perturbations are studied: it is analogous to the thermodynamic equilibrium state for elastic systems . In the IHS model this state is particularly simple and the energy decreases obeying the Haff’s law

$$E(t)=\frac{E(0)}{\left(1+t/t_0\right)^2}$$
(3)
In order to avoid the cooling of the system, the simulations are made at constant energy: at every collision between the grains, the dissipated energy is redistributed by scaling all the velocities. This is equivalent to a time rescaling and can be treated as such in the equations for the continuous model. Using appropriate hydrodynamic equations, this homogeneous state is predicted to be unstable under certain conditions of density, system size and dissipativity . Considering the dissipativity coefficient q as a bifurcation parameter, while increasing q at constant density and number of grains, the system first develops an instability characterized by two counterflows, the shearing instability, and then either a clustering regime in which the density becomes inhomogeneous or a vortex state where many small vortices develop throughout the system.

For given values of total number of grains $`N`$ and number density $`\rho =N/V`$ the shearing instability appears when the dissipation coefficient is larger than a critical value. In the low density and large system limit,the critical value is given (in two dimensions) by

$$\widehat{q}=\frac{\pi }{N\rho }$$
(4)
Note that in the thermodynamic limit the system is always unstable for any finite dissipation.

In this article we will study the shearing instability using a nonlinear hydrodynamic approach. In Sec. II we will develop a formalism that allows us to treat the HCS as a non equilibrium steady state, and we will present a stability analysis around the steady state. Next, in Sec. III we will present the nonlinear analysis of the instability, obtaining expressions for the hydrodynamic fields beyond the instability. It will also be shown that the presence of the instability modifies the collision rate and energy dissipation rate. Finally, in Sec. IV we compare the predictions of the continuous model with molecular dynamics simulations of inelastic hard disks. Conclusions are presented in Sec. V.

In what follows we will treat the two dimensional case, but the extension to three dimensions is direct.

## II Rescaled time formalism

Consider a two dimensional system composed of $`N`$ grains interacting with the IHS model in a square box of size $`L`$ that has periodic boundary conditions in both directions. Units are chosen such that Boltmann’s constant and the particle masses are set to one. Granular temperature is defined analogously to the kinetic definition for classical fluids

$$T=\frac{1}{N}\underset{i}{}\frac{1}{2}(𝐯_i𝐯)^2$$
(5)
where $`𝐯`$ is the hydrodynamic velocity.

The shearing instability has been predicted by linear analysis of the hydrodynamic equations for the IHS. Also an analysis of the first stages of the nonlinear regime has been done in Ref. . As we shall see, the relative simplicity of the model allows for a complete nonlinear analysis as well.

The hydrodynamic equations for the low-density IHS are similar to the usual hydrodynamic equations for fluids except that there is an energy sink term and the heat flux has contributions from the density gradient . The equations read

$`{\displaystyle \frac{\rho }{t}}+(\rho 𝐯)=0`$ (6)
$`\rho \left({\displaystyle \frac{𝐯}{t}}+(𝐯)𝐯\right)=IP`$ (7)
$`\rho \left({\displaystyle \frac{T}{t}}+(𝐯)T\right)=𝐉IP:𝐯\omega `$ (8)
with the following constitutive equations

$`IP_{ij}`$ $`=`$ $`\rho T\delta _{ij}\eta _0\sqrt{T}\left({\displaystyle \frac{v_i}{x_j}}+{\displaystyle \frac{v_j}{x_i}}(𝐯)\delta _{ij}\right)`$ (9)
$`𝐉`$ $`=`$ $`k_0\sqrt{T}T\mu _0{\displaystyle \frac{T^{3/2}}{\rho }}\rho `$ (10)
$`\omega `$ $`=`$ $`\omega _0\rho ^2T^{3/2}`$ (11)
where $`\eta _0`$, $`k_0`$, $`\mu _0`$, and $`\omega _0`$ do not depend on density or temperature but on the dissipation coefficient $`q`$. In particular, $`\omega _0`$ and $`\mu _0`$ vanish with $`q`$.

The total energy dissipation rate can be computed from the hydrodynamic equations, obtaining

$`{\displaystyle \frac{dE}{dt}}`$ $`=`$ $`{\displaystyle \frac{d}{dt}}{\displaystyle (\rho T+\rho v^2/2)𝑑V}`$ (12)
$`=`$ $`{\displaystyle \omega 𝑑V}`$ (13)
In the HCS, where density and temperature are homogeneous and there is no velocity field, the energy dissipation rate is

$$\frac{dE}{dt}=V\omega _0\rho _0^2T^{3/2}$$
(14)
This continuous cooling down has the disadvantage that any perturbation analysis must be done with respect to this non-steady state. To overcome this difficulty we propose a differential time rescaling $`ds=\gamma dt`$ such that in this rescaled time, energy is conserved. The transformation corresponds to a continuous rescaling of all the particle velocities such that the kinetic energy remains constant. This rescaling does not, however, introduce new phenomena since as we have mentioned the IHS does not have an intrinsic time scale, and thus a time rescaling gives rise to the same phenomena viewed at different speed.

In order to keep the rescaled energy constant the appropriate value of $`\gamma `$ is given by

$$\gamma =\sqrt{\frac{E(t)}{E(0)}}$$
(15)
The rescaled time hydrodynamic fields transform as

$`𝐯`$ $`=`$ $`\gamma \stackrel{~}{𝐯}`$ (16)
$`T`$ $`=`$ $`\gamma ^2\stackrel{~}{T}`$ (17)
$`IP_{ij}`$ $`=`$ $`\gamma ^2\stackrel{~}{IP_{ij}}`$ (18)
$`𝐉`$ $`=`$ $`\gamma ^3\stackrel{~}{𝐉}`$ (19)
$`\omega `$ $`=`$ $`\gamma ^3\stackrel{~}{\omega }`$ (20)
where the tilde denotes the rescaled variable.

Due to the non-constant character of the time rescaling, extra source terms appear in the hydrodynamic equations, which can be simplified using the relation

$$\frac{1}{\gamma }\frac{d\gamma }{ds}=\frac{1}{2E(0)}\stackrel{~}{\omega }𝑑V$$
(21)
Suppressing the tildes everywhere, the equations now read:

$`{\displaystyle \frac{\rho }{s}}+(\rho 𝐯)=0`$ (22)
$`\rho \left({\displaystyle \frac{𝐯}{s}}+(𝐯)𝐯\right)=IP+{\displaystyle \frac{\rho 𝐯}{2E(0)}}{\displaystyle \omega 𝑑V}`$ (23)
$`\rho \left({\displaystyle \frac{T}{s}}+(𝐯)T\right)=𝐉IP:𝐯\omega +{\displaystyle \frac{\rho T}{E(0)}}{\displaystyle \omega 𝑑V}`$ (24)
Note that the constitutive relations (9) remain unchanged under time rescaling.

In the rescaled time, the HCS reduces to a non-equilibrium steady state with a continuous energy supply that compensates the energy dissipation. As there is no energy scale, we fix the reference temperature to be one, so that $`E(0)=N`$. The HCS is then characterized by $`\{\rho =\rho _H,T=1,𝐯=0\}`$. To study the stability of this state, we first introduced the change of variables:

$`\rho `$ $`=`$ $`\rho _H+\delta \rho `$ (25)
$`𝐯`$ $`=`$ $`\delta 𝐯`$ (26)
$`T`$ $`=`$ $`1+\delta T`$ (27)
Taking the (discrete) Fourier transform of the linearized hydrodynamic equations around HCS, it is easy to check that the transverse velocity perturbation decouples from the rest and satisfies the equation

$$\rho _H\frac{\delta 𝐯_k_{}}{s}=\lambda _k\delta 𝐯_k_{}$$
(28)
with

$$\lambda _k\frac{4k^2\pi ^2}{L^2}\eta _0+\frac{\rho _H^2\omega _0}{2};\{k_x,k_y\}=\mathrm{\hspace{0.17em}0},\pm 1,\pm 2,\mathrm{}$$
(29)
For small values of $`\omega _0`$ (proportional to the dissipation) $`\lambda _k<0`$, so that the perturbations $`\delta 𝐯_k_{}`$ decay exponentially (the case $`k=0`$ should not be taken into account since the center of mass velocity is strictly zero). But there exists a critical value of $`\omega _0`$ for which $`\lambda _k`$ vanishes, thus indicating the stability limit of the corresponding mode. The first modes to become unstable correspond to $`|𝐤|=1`$ (i.e, $`k_x=\pm 1,k_y=0`$ and $`k_x=0,k_y=\pm 1`$). The instability threshold for $`\omega _0`$ is then given by

$$\widehat{\omega }_0=\frac{8\pi ^2\eta _0}{\rho _H^2L^2}$$
(30)
The stability for the other modes has been studied previously . In this last reference, it was shown that for, low dissipation, the first instability that arises is indeed the transverse velocity instability.

The origin of the instability can be understood also in terms of the real time hydrodynamics. In fact, the relaxation of the transverse hydrodynamic velocity is basically due to the viscous diffusion, which depends on the system size, whereas the cooling process is governed by local dissipative collisions. There exists therefore a system length beyond which the dissipation of thermal energy is faster that the relaxation of the transverse hydrodynamic velocity. The latter will then increase, when observed on the scale of thermal motion, thus producing the shearing pattern.

## III Nonlinear analysis of the instability

In this section we propose to work out the explicit form of the velocity field beyond the instability. The calculations are tedious and quite lengthy, so that, here, we only describe explicitely the basic steps. We start by taking the Fourier transform of the full nonlinear hydrodynamic equations, obtaining a set of coupled nonlinear equations for the modes { $`\delta \rho _k,\delta T_k,\delta 𝐯_k_{},\delta 𝐯_k_{}`$ }, where $`\delta 𝐯_k_{}`$ represents the longitudinal component of the velocity field. Close to the instability ($`\omega _0\widehat{\omega }_0`$, $`|𝐤|=1`$), the modes $`\delta 𝐯_1_{}`$ exhibit a critical slowing down since $`\lambda _10`$, whereas the other hydrodynamic modes decay exponentially (cfr. ref. ). On this slow time scale, i.e. $`sO(\lambda _1^1)`$, the “fast” modes {$`\delta \rho _k;\delta T_k;\delta 𝐯_k_{};\delta 𝐯_k_{},k1`$} can be considered as stationary, their time dependence arising mainly through $`\delta 𝐯_1_{}`$. Setting the time derivatives of these fast modes to zero one can express them in terms of the slow modes $`\delta 𝐯_1_{}`$. If now one inserts the so-obtained expressions into the evolution equation for the slow modes, one gets a set of closed nonlinear equations for $`\delta 𝐯_1_{}`$ (adiabatic elimination ), usually referred to as normal form or amplitude equations \[G. Nicolis, Introduction to Nonlinear Science, Cambridge University Press (1995) \]. Note that such a calculation is only possible close to the instability threshold, where the amplitude of $`\delta 𝐯_1_{}`$ approaches zero as $`\omega \widehat{\omega }_0`$. On the other hand, in this limit the only Fourier modes that will have large amplitudes are the transverse velocity modes with wavevector equal to $`\pm \mathrm{\hspace{0.17em}2}\pi /L`$ (i.e. $`|𝐤|=1`$). There are four such modes:

$$\delta 𝐯_k_{}=\{\begin{array}{cc}A_1\hfill & 𝐤=(1,0)\hfill \\ A_1^{}\hfill & 𝐤=(1,0)\hfill \\ A_2\hfill & 𝐤=(0,1)\hfill \\ A_2^{}\hfill & 𝐤=(0,1)\hfill \end{array}$$
(31)
After some tedious algebra, one finds

$`\rho _H{\displaystyle \frac{dA_1}{ds}}`$ $`=`$ $`\lambda _1A_14\pi ^2L^2\eta _0\left(C_1\left|A_1\right|^2+C_2\left|A_2\right|^2\right)A_1`$ (32)
$`\rho _H{\displaystyle \frac{dA_2}{ds}}`$ $`=`$ $`\lambda _1A_24\pi ^2L^2\eta _0\left(C_1\left|A_2\right|^2+C_2\left|A_1\right|^2\right)A_2`$ (33)
where

$`C_1`$ $`=`$ $`{\displaystyle \frac{8(k_0\mu _0)3\eta _0}{8(k_0\mu _0)2\eta _0}}`$ (34)
$`C_2`$ $`=`$ $`{\displaystyle \frac{2(k_0\mu _0)+\eta _0}{2(k_0\mu _0)\eta _0}}`$ (35)
The amplitude equations (32) and (33) admit three different stationary solutions:

$`(a)`$ $`A_1=0,A_2=0`$ (36)
$`(b)`$ $`\left|A_1\right|=\widehat{A},A_2=0`$ (38)
$`A_1=0,\left|A_2\right|=\widehat{A}`$
$`(c)`$ $`\left|A_1\right|=\left|A_2\right|=\widehat{A}\sqrt{{\displaystyle \frac{C_1}{C_1+C_2}}}`$ (39)
where we have set

$$\widehat{A}=\frac{L}{2\pi }\sqrt{\frac{\lambda _1}{\eta _0C_1}}$$
(40)
Note that the phases of the above stationary solutions are arbitrary (recall that $`A_1`$ and $`A_1`$ are complex conjugate).

The trivial solution (a) corresponds to a motionless fluid, whereas the solutions (b) are shearing states with the corresponding fluxes oriented either in the $`y`$ or $`x`$ direction. The mixed mode solution (c) represents a vortex state with two counter-rotating vortices in the box.

Below the critical point ($`\lambda _1<0`$), the trivial solution (a) is the only stable one. As we cross the critical point this solution becomes unstable. A linear stability analysis of the Eqs. 33 shows that the shearing states are stable provided

$$\frac{C_2}{C_1}>1$$
(41)
while the mixed mode solution is unstable. Satisfying Eq.41 depends on the values of the transport coefficients. For small dissipativity it is always fulfilled, provided the number density remains relatively low. In fact, a mixed mode state has been observed recently in a highly dense system .

The occurence of either shearing states depends on the initial state. In a statistical sense, they are equally probable. For example, in the case that the system chooses the solution {$`A_1=\widehat{A}`$, $`A_2=0`$}, the velocity field reads

$`𝐯=2\widehat{A}\mathrm{cos}(2\pi x/L)\widehat{𝐲}`$ (42)
where $`\widehat{𝐲}`$ represents the unit vector in $`Y`$ direction. The temperature and density perturbations are then given by

$`\delta T`$ $`=`$ $`(\widehat{A})^2\left[1+(1C_1)\mathrm{cos}(4\pi x/L)\right]`$ (43)
$`\delta \rho `$ $`=`$ $`\rho _H(\widehat{A})^2(1C_1)\mathrm{cos}(4\pi x/L)`$ (44)
The instability of the transverse velocity field thus gives rise to modifications of the temperature and density fields. The temperature decreases globally, since part of the energy is taken by the convective motion. Moreover, because of the viscous heating, the temperature profile exhibits a spatial modulation, i.e. it is higher where the viscous heating is higher. The density profile also shows a spatial modulation that keeps the pressure homogeneous (recall that for a two dimensional low density gas, $`\delta p\rho _h\delta T+\delta \rho `$, in system units). This modulation of the density, that was first observed by Goldhirsh and Zanetti  (Fig. 2 of their article) in the study of the clustering instability, is stable and does not lead to clustering.

Using the above expressions for the hydrodynamic fields, the energy density profile reads

$$e=\rho _H\left(1+\widehat{A}^2\mathrm{cos}(4\pi x/L)\right)$$
(45)
Assuming the molecular chaos hypothesis, the mean collision rate, $`\overline{\nu }`$, and dissipation rate, $`\overline{\omega }`$, can be computed as well:

$`\overline{\nu }`$ $`=`$ $`{\displaystyle \frac{2\sqrt{\pi }}{V\rho _H}}{\displaystyle \rho ^2\sqrt{T}𝑑V}=\nu _0\rho _H\left(1{\displaystyle \frac{\widehat{A}^2}{2}}\right)`$ (46)
$`\overline{\omega }`$ $`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle \omega 𝑑V}=\omega _0\rho _H^2\left(1{\displaystyle \frac{3\widehat{A}^2}{2}}\right)`$ (47)
These relations show that the global decrease of the temperature leads to corresponding decreases of the collision frequency and dissipation rate.

It is important to note that the origin of the nonlinear coupling of slow modes lies in the viscous heating term and the state-dependence of the transport coefficients, and not in the usual convective derivatives ($`(𝐯)𝐯`$). In fact, as we have shown above, the shearing state produces a variation in the density and temperature fields that modify locally the value of the transport coefficients. This effect, which is negligible in classical fluids, can become very important in granular fluids, mainly because of the the lack of scale separation of the kinetic and hydrodynamic regimes. Contrary to normal fluids, here, the convective energy is comparable to the thermal energy. In fact, as will be shown in the MD simulations, in a well developed shearing state up to half of the total kinetic energy corresponds to the convective motion. The microscopic source of this phenomenon lies in the fact that only the relative energy is dissipated in binary collisions, i.e. the center of mass energy is conserved. In other words, the thermal energy is dissipated but the convective one is conserved. This asymmetry in the energy dissipation mechanism is at the very origin of the shearing instability.

## IV Molecular dynamics simulations

For the molecular dynamic simulations we have considered a system made of $`N=10000`$ hard disks with a global number density $`\rho _H=0.005`$. Inelastic collision rules are adopted, with a dissipativity varying from $`q=0.0`$ to $`q=0.12`$. The boundary conditions are periodic in all directions. A spatially homogeneous initial condition is adopted, with velocities sorted from an equilibrium (zero mean velocity) Maxwellian distribution. We note that the density is low enough so that the system remains within the low-density regime.

The simulations have been performed in the rescaled time. Computationally, this is achieved by doing a normal IHS simulation (event driven molecular dynamics ), but at each collision the value of the kinetic energy is updated according to the energy dissipated. The instantaneous value of $`\gamma `$ is computed from the kinetic energy, allowing to evaluate all the rescaled quantities. Note that the $`s`$-time can be integrated in the simulation because $`\gamma `$ is a piecewise constant function, thus allowing to make periodic measurements in the system. Finally, to avoid roundoff errors, a real velocity rescaling is performed whenever the kinetic energy decreases by a given amount (typically $`10^7`$ of the initial value).

In each simulation, the collision frequency and temperature dissipation rate are computed with respect to the $`s`$-time, after the system has reached a stationary regime. In Fig. 1, the collision frequency and the temperature dissipation rate are presented as a function of the dissipation coefficient $`q`$. As expected, the shearing instability is associated with an abrupt decrease of these functions. It must be noted that the decreasing of the collision frequency is more that $`30\%`$, which corresponds to a global decrease of the temperature of more than $`50\%`$. This means that when the shearing is fully developed, about half the total kinetic energy is taken by the macroscopic motion. This phenomenon is typical of granular media and has no counterpart in classical fluids.

To measure the critical point, $`q_0`$, we fit the collision frequency according to the following piecewise function

$$\overline{\nu }=\{\begin{array}{cc}a_0\hfill & qq_0\hfill \\ a_0+a_1(qq_0)+a_2(qq_0)^2\hfill & q>q_0\hfill \end{array}$$
(48)
obtaining

$`q_0`$ $`=`$ $`0.0686`$ (49)
$`a_0`$ $`=`$ $`0.0178`$ (50)
$`a_1`$ $`=`$ $`0.167`$ (51)
$`a_2`$ $`=`$ $`1.03`$ (52)
Similarly, the dissipation rate is fitted according to

$$\overline{\omega }=\{\begin{array}{cc}b_0q\hfill & qq_0\hfill \\ b_0q+b_1(qq_0)q+b_2(qq_0)^2q\hfill & q>q_0\hfill \end{array}$$
(53)
Using $`q_0=0.0686`$, one finds

$`b_0`$ $`=`$ $`0.000166`$ (54)
$`b_1`$ $`=`$ $`0.00466`$ (55)
$`b_2`$ $`=`$ $`0.0546`$ (56)
To compare these results with our theoretical predictions we need the explicit form of the transport coefficients, up to critical dissipativity. Unfortunately, there are no known expressions for them in the case the $`2d`$ IHS model in the low-density regime. However, as the critical dissipativity is small we can use the the quasielastic approximation for the transport coefficients (that is, taking the first non-trivial order in $`q`$)

$`\omega _0`$ $`=`$ $`4\sqrt{\pi }q`$ (57)
$`\eta _0`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{\pi }}}`$ (58)
$`k_0`$ $`=`$ $`2/\sqrt{\pi }`$ (59)
$`\mu _0`$ $`=`$ $`0`$ (60)
The critical dissipativity and the amplitude of the shearing state are given by

$`\widehat{q}`$ $`=`$ $`{\displaystyle \frac{\pi }{\rho _H^2L^2}}`$ (61)
$`\widehat{A}`$ $`=`$ $`\rho _HL\sqrt{{\displaystyle \frac{30}{29\pi }}\delta q}`$ (62)
$`=`$ $`\sqrt{{\displaystyle \frac{30}{29}}{\displaystyle \frac{\delta q}{\widehat{q}}}}`$ (63)
For the presented simulation, the predicted critical dissipativity is

$$\widehat{q}=0.0628$$
(64)
which shows a discrepancy of $`8\%`$ with the observed value. This difference is consistent with the adopted approximations.

The predicted values for $`a_0`$, $`a_1`$, $`b_0`$, $`b_1`$ (cfr. Eqs. 46, 47, 61, and 63) are

$`a_0`$ $`=`$ $`0.0177`$ (65)
$`a_1`$ $`=`$ $`0.146`$ (66)
$`b_0`$ $`=`$ $`0.000177`$ (67)
$`b_1`$ $`=`$ $`0.00438`$ (68)
which are also consistent with the adopted approximations.

Since the system is periodic, the developed convective pattern can diffuse in the direction perpendicular to the flow (the phases of the complex amplitudes $`A_i`$ are arbitrary due to Galilean invariance). As a result, the average hydrodynamic fields remain vanishingly small, mainly because of “destructive” interference. To overcome this difficulty we have performed another series of simulations, keeping periodic boundary conditions in the vertical direction, while introducing a pair of stress-free and perfectly insulating parallel walls in the horizontal direction (in a collision with a wall the tangential velocity is conserved whereas the normal one is inverted). As a consequence, the total vertical momentum is conserved, which will be simply set to zero initially.

The nonlinear analysis for this case is similar to the periodic one, except that, here, the direction of the flow pattern remains always parrallel to the walls. Furthermore, the unstable wavevector is now $`k=\pi /L`$, because of to the fixed boundary conditions. As a result all the previous predictions remain valid, except that everywhere $`L`$ must be replaced by $`2L`$.

We have used the very same number of particles and density for this series of simulations, but, of course, the different boundary conditions produce a new critical dissipativity. Performing the same analysis as before, the measured critical point and fit parameters turn out to be

$`q_0`$ $`=`$ $`0.0163`$ (69)
$`a_0`$ $`=`$ $`0.0178`$ (70)
$`a_1`$ $`=`$ $`0.562`$ (71)
$`a_2`$ $`=`$ $`8.93`$ (72)
$`b_0`$ $`=`$ $`0.000175`$ (73)
$`b_1`$ $`=`$ $`0.0179`$ (74)
$`b_2`$ $`=`$ $`1.06`$ (75)
while the predicted ones are

$`\widehat{q}`$ $`=`$ $`0.0157`$ (76)
$`a_0`$ $`=`$ $`0.0177`$ (77)
$`a_1`$ $`=`$ $`0.583`$ (78)
$`b_0`$ $`=`$ $`0.000177`$ (79)
$`b_1`$ $`=`$ $`0.0175`$ (80)
Eqs. 42, 44, and 45, once $`L`$ is replaced by $`2L`$, indicate that the perturbation of the transverse momentum density ($`𝐣=\rho 𝐯`$) has a wavevector equal to $`k_x=\pi /L`$, while the density and energy density have wavevectors $`k_x=2\pi /L`$. In the simulations we computed the amplitudes of these Fourier modes using the microscopic definitions for the particle, momentum, and energy densities. Figs. 2, 3, and 4 show the predicted and computed Fourier modes amplitudes. The predictions are in good agreement with the simulations in the neighborhood of the critical point, showing that not only the average quantities like the collision frequency are well predicted, but also the whole hydrodynamic picture is correct.

## V Conclusions

Taking advantage of the lack of energy scale in the IHS model, a rescaled time formalism was introduced that allows us to study the homogeneous cooling state as a nonequilibrium steady state. Using a hydrodynamic description for granular media written with a rescaled time variable, the shearing instability has been studied in the nonlinear regime. It has been shown that the shearing state is the stable solution and its amplitude has been computed. The appearance of the velocity field produces that part of the kinetic energy goes from the kinetic to the hydrodynamic scale. In usual fluids this redistribution of the energy is negligible, but in granular fluids it can represent an important fraction of the total energy. This phenomenon is a manifestation of a global property of granular fluids: there is not in general a clear separation between the kinetic and the hydrodynamic regimes. This could lead to put into question the validity of a hydrodynamic description. Nevertheless, at small values of the dissipativity coefficient, predictions based on the nonlinear hydrodynamic equations are in excellent agreement with molecular dynamics simulations. Both the value of the critical dissipativity and the behavior after the instability has developed are well predicted. This is a remarkable result which shows again how robust the hydrodynamic fluid equations are when they are tested at time and length scales where their validity could be questioned.

## ACKNOWLEDGMENTS

This work is supported by a European Commission DG 12 Grant PSS\*1045 and by a grant from FNRS Belgium. One of us (R.S.) acknowledges the grant from MIDEPLAN.