text:
# Hydromagnetic instabilities in proto–neutron stars

## 1 Introduction

Neutron stars are born in the aftermath of core collapse Supernovae explosions, the fate of massive stars within a range of 8 to 30$`M_{}`$. In the final stages of stellar evolution an iron core grows, until electron captures or photodesintegration of heavy nuclei trigger its gravitational collapse. The collapse proceeds until nuclear density is reached and the stiffening of the equation of state (EOS) provokes the bounce of the infalling material. A shock wave forms and propagates outwards, leaving behind a newly born, hot, lepton–rich PNS. These objects are prodigious emitters of neutrinos of all types, which dominate the energetics of core collapse Supernovae and play a crucial role in the explosion mechanism by carrying away the binding energy from the PNS and depositing a part of it in the outer regions.

Whether or not the energy relocated by neutrinos is able to relaunch the stalled shock is still a controversial issue (Burrows, et al., 2000; Rampp & Janka, 2000; Janka, 2000; Liebendorfer et al., 2001) that has led to speculation about alternative ways the explosion is powered such as convection, rotation, and magneto-hydrodynamics effects. In particular, it has been recognized by a number of authors that convective energy transport in the newly born PNS might play an important role by increasing the neutrino luminosity required to ignite the Supernova explosion (Burrows, & Fryxell, 1993; Janka & Müller, 1996; Bruenn & Dineva, 1996) Convection in PNSs can be driven not only by the lepton gradient, as originally suggested by Epstein (1979), but also by the development of negative entropy gradients, which is common in many simulations of supernovae models (Bruenn & Mezzacappa, 1994; Bruenn, Mezzacappa, & Dineva, 1995; Rampp & Janka, 2000; Liebendorfer et al., 2001) and evolutionary models of PNSs (Burrows & Lattimer, 1986; Keil & Janka, 1995; Keil, Janka & Müller, 1996; Pons et al., 1999, 2001a, 2001b) despite differences in the equation of state and neutrino transport.

Although stellar core collapse may also have dramatic consequences on the magnetic field strength, the effect of magnetic fields was not considered in the above mentioned works, mainly due to the tough technical difficulties inherent to multidimensional magnetohydrodynamical simulations. It is known that the magnetic field in a collapsing star can be amplified by many orders of magnitude due to the conservation of the magnetic flux. This idea is often used to explain the origin of strong magnetic fields in neutron stars (Ginzburg, 1964; Ginzburg & Ozernoy, 1964; Woltjer, 1964; Zeldovich & Novikov, 1971; Shapiro & Teukolsky, 1983). The conductivity of plasma inside the star is so large that the magnetic decay timescale is much shorter than the collapse time for a lengthscale of the magnetic field comparable to the stellar radius. Therefore, the field is frozen-in, and the magnetic flux is conserved during the collapse stage. In the case of a frozen-in magnetic field, the average field strength increases approximately as $`B\rho ^{2/3}`$ where $`\rho `$ is the density. By assuming that the magnetic field strength of the core of a massive star is comparable to that of strongly magnetized white dwarfs ($`10^9`$ G), at the end of core collapse, the neutron star may have a magnetic field $`10^{13}10^{14}`$ G, comparable to that observed in young radio-pulsars. Following the same arguments, if we consider that the magnetic flux is conserved during the evolution of a massive star from the main sequence, starting from an initial magnetic field of $`10^310^4`$ G, we end up with a magnetic field of $`10^{13}10^{14}`$ G at the time the core of the star has become a neutron star. Similar arguments are often used to estimate the magnetic field of neutron stars formed from magnetized white dwarfs by the accretion induced collapse mechanism. Note that these estimates give the order of magnitude of a poloidal field. The toroidal field, however, might be much stronger (Ardelyan, Bisnovatyi-Kogan, & Popov, 1980).

Under certain conditions, convective motions in PNSs can amplify the magnetic field via dynamo action. The most optimistic estimates lead some authors (Thompson & Duncan, 1993; Thompson & Murray, 2001) to the conclusion that turbulent dynamo action could be responsible for the formation of ultra-magnetized neutron stars (magnetars), with fields as strong as $`10^{16}`$ G at a very early evolutionary stage. Certainly, such strong magnetic fields influence magneto-hydrodynamic processes in PNSs, particularly their stability properties. At this point, the fact that the convective overturn is intrinsically a 3D phenomenon makes fully consistent numerical simulations a very laborious and computationally expensive way to address the problem, and simplified semi-analytical approaches are extremely useful to evaluate the necessity of more detailed studies. In this respect, some attempts have been made to establish definite criteria as to the different types of instability and their growth times (Bruenn & Dineva, 1996; Miralles, Pons, & Urpin, 2000, MPU in the following). In MPU we derived the general criteria for instability, taking dissipative processes into account such as neutrino transport and viscosity. In the present paper, we consider the stability properties of magnetic PNSs in detail. It is well known (see, e.g., Chandrasekhar 1961) that the magnetic field may stabilize a fluid against convection causing an additional effective viscosity due to Lorentz force. On the other hand, the field can lead to new branches of instability associated, for instance, with the Alfvén waves. In this paper, we present the main results of a general stability analysis in magnetized PNSs. The stability properties will be considered for a wide range of magnetic fields.

The paper is organized as follows. In §2, we obtain the dispersion equation which determines the condition of linear instability from the equations governing magneto-hydrodynamics of PNSs and neutrino transport. In §3, we discuss the stability properties of both non-magnetic (§3.1) and magnetic (§3.2) PNSs, appling the criteria of hydromagnetic stability to the results of numerical simulations to identify the unstable zones. We also calculate the critical magnetic field that stabilizes a PNS. Finally, in §4 we summarize and discuss our main findings.

## 2 The dispersion equation

Consider the stability properties on adopting a plane-parallel geometry. We assume that within the layer between $`z=0`$ and $`z=d`$ the gravity $`𝐠`$ is directed in the negative $`z`$-direction and neglect a non-uniformity of $`𝐠`$ as well as general relativistic corrections to the equations of hydrodynamics. The characteristic cooling time scale of a PNS is assumed to be much longer than the growth time of instability thus the latter can be treated in a quasi-stationary approximation. Since the convective velocities are typically much smaller than the speed of sound, one can describe instability by making use of the standard Boussinesq approximation (see, e.g., Landau & Lifshitz 1987). We consider the linear instability regime when the equations governing small perturbations can be obtained by the linearization of hydrodynamic equations. In the following, small perturbations of hydrodynamic quantities will be marked by a subscript “1”. For simplicity, the unperturbed magnetic field is assumed to be vertical. The linearized momentum, continuity and induction equations read

$$\rho \dot{𝐯}_1=p_1+𝐠\rho _1+\frac{1}{4\pi }(\times 𝐁_1)\times 𝐁+\rho \nu \mathrm{\Delta }𝐯_1,$$
(1)
$$𝐯_1=0,$$
(2)
$$\dot{𝐁}_1=\times (𝐯_1\times 𝐁)+\nu _m\mathrm{\Delta }𝐁_1,$$
(3)
$$𝐁_1=0,$$
(4)
where $`p`$ and $`\rho `$ are the pressure and density, respectively, $`\nu `$ is the viscosity and $`\nu _m`$ is the magnetic diffusivity. We assume that the matter inside a PNS is in chemical equilibrium, thus the density can generally be considered as a function of the pressure $`p`$, temperature $`T`$ and lepton fraction $`Y`$, $`Y=(n_e+n_\nu )/n`$, with $`n_e`$ and $`n_\nu `$ being the net (particles minus antiparticles) number densities of electrons and neutrinos, respectively, and $`n=n_p+n_n`$ is the number density of baryons. Then, the perturbations of density, entropy per baryon ($`s`$) and neutrino chemical potential ($`\mu `$) can be expressed in terms of $`p_1`$, $`T_1`$ and $`Y_1`$. In the Boussinesq approximation, the perturbations of pressure are negligible because the fluid motions are assumed to be slow and the moving fluid elements are nearly in pressure equilibrium with their surroundings. Therefore, we have

$$\rho _1\rho \left(\beta \frac{T_1}{T}+\delta Y_1\right),$$
(5)
$$s_1m_pc_p\frac{T_1}{T}+\sigma Y_1,$$
(6)
where $`m_p`$ is the proton mass (we neglect the mass difference between protons and neutrons), $`\beta `$ and $`\delta `$ are the coefficient of thermal and chemical expansion, respectively, $`\beta =(\mathrm{ln}\rho /\mathrm{ln}T)_{pY}`$, $`\delta =(\mathrm{ln}\rho /Y)_{pT}`$; $`\sigma =(s/Y)_{pT}`$; and $`c_p=(T/m_p)(s/T)_{pY}`$ is the specific heat at constant pressure.

The above equations should be complemented by the equation driving the evolution of chemical composition and heat balance. We employ the equilibrium diffusion approximation, which is sufficiently accurate and reliable during the early stage of PNS evolution, when the mean free path of neutrino is short compared to the density and temperature length scales. In this approximation, the diffusion and thermal balance equations can be linearized (see MPU for details) and written as

$$\dot{Y}_1+𝐯_1Y=\lambda _T\frac{\mathrm{\Delta }T_1}{T}+\lambda _Y\mathrm{\Delta }Y_1,$$
(7)
$$\frac{\dot{T}_1}{T}𝐯_1\frac{\mathrm{\Delta }T}{T}=\kappa _T\frac{\mathrm{\Delta }T_1}{T}+\kappa _Y\mathrm{\Delta }Y_1,$$
(8)
where

$$\mathrm{\Delta }T\frac{T}{m_pc_p}(s\sigma Y)=\left(\frac{T}{p}\right)_{s,Y}pT$$
(9)
is the superadiabatic temperature gradient, i.e., the difference between the temperature gradient of a fluid with constant entropy and composition and the actual temperature gradient, and the transport coefficients $`\lambda _T,\kappa _T,\lambda _Y,\kappa _Y`$ have been defined in MPU.

Equations (1)–(4), (7) and (8) together with the corresponding boundary conditions determine the behaviour of small perturbations. For simplicity, we consider the case of vanishing perturbations at the boundaries $`z=0`$ and $`z=d`$. Note that other boundary conditions cannot change the main conclusions of our analysis qualitatively.

Without loss of generality, the dependence of all perturbations on time and on the horizontal coordinate can be chosen as $`\mathrm{exp}(\gamma tikx)`$, where $`k`$ is the horizontal wavevector and $`\gamma `$ is a constant which can be complex and its real part gives the inverse growth (or decay) timescale. Then, equations (1)-(4), (12), and (13) can be reduced to only one equation of higher order. Solving for $`v_{1z}`$, we obtain

$`\mathrm{\Delta }[(\gamma \kappa _T\mathrm{\Delta })(\gamma \lambda _Y\mathrm{\Delta })\lambda _T\kappa _Y\mathrm{\Delta }^2]`$
$`\left[(\gamma \nu \mathrm{\Delta })(\gamma \nu _m\mathrm{\Delta })c_A^2{\displaystyle \frac{d^2}{dz^2}}\right]v_{1z}=`$
$`gk^2\{{\displaystyle \frac{dY}{dz}}[\delta (\gamma \kappa _T\mathrm{\Delta })+\beta \kappa _Y\mathrm{\Delta }]`$
$`{\displaystyle \frac{\mathrm{\Delta }T}{T}}[\beta (\gamma \lambda _Y\mathrm{\Delta })+\delta \lambda _T\mathrm{\Delta }]\}(\gamma \nu _m\mathrm{\Delta })v_{1z},`$ (10)
where $`c_A=B/\sqrt{4\pi \rho }`$ is the Alfvén velocity and $`\mathrm{\Delta }=\frac{d^2}{dz^2}k^2.`$ The coefficients of this equation are constant in our simplified approach. The solution for the fundamental mode with “zero boundary conditions” has the form $`v_{1z}\mathrm{sin}(\pi z/d)`$. Thus the dispersion equation for the fundamental mode is

$`[(\gamma +\omega _\nu )(\gamma +\omega _m)+\omega _A^2][(\gamma +\omega _T)(\gamma +\omega _Y)\omega _{TY}\omega _{YT}]`$ $`=`$
$`(\gamma +\omega _m)\left[\omega _g^2(\gamma +\omega _Y\alpha \omega _{YT})+\omega _L^2\left(\gamma +\omega _T{\displaystyle \frac{\omega _{TY}}{\alpha }}\right)\right],`$ (11)
where $`\alpha =\delta /\beta `$. In this expression, we introduced the characteristic frequencies

$`\omega _\nu =\nu Q^2,\omega _m=\nu _mQ^2,\omega _T=\kappa _TQ^2,\omega _Y=\lambda _YQ^2,`$
$`\omega _{YT}=\lambda _TQ^2,\omega _{TY}=\kappa _YQ^2,\omega _A=c_A\pi /a,`$ (12)
$`\omega _g^2={\displaystyle \frac{\beta gk^2\mathrm{\Delta }T}{TQ^2}},\omega _L^2={\displaystyle \frac{\delta gk^2}{Q^2}}{\displaystyle \frac{dY}{dz}},`$
with $`Q^2=(\pi /d)^2+k^2`$. The quantities $`\omega _\nu `$, $`\omega _m`$, $`\omega _T`$, and $`\omega _Y`$ are the inverse time scales of dissipation of the perturbations due to viscosity, magnetic diffusivity, thermal conductivity and chemical diffusivity, respectively; $`\omega _{YT}`$ characterizes the rate of diffusion caused by the temperature inhomogeneity (thermodiffusion), and $`\omega _{TY}`$ describes the influence of chemical inhomogeneities on the rate of thermal evolution; $`\omega _g`$ is the frequency (or, in the case of instability, the inverse growth time) of the buoyant wave; $`\omega _L`$ characterizes the dynamical time scale of the processes associated with the lepton gradient, and $`\omega _A`$ is the Alfvén frequency.

## 3 Stability of proto-neutron stars

### 3.1 Non-magnetic PNSs

Consider initially the case when $`𝐁=\mathrm{𝟎}`$. Then, equation (11) is reduced to a cubic equation with the roots given by

$$\gamma ^3+a_2\gamma ^2+a_1\gamma +a_0=0,$$
(13)
where

$`a_2`$ $`=`$ $`\omega _\nu +\omega _T+\omega _Y,`$
$`a_1`$ $`=`$ $`\omega _T\omega _Y\omega _{TY}\omega _{YT}+\omega _\nu (\omega _T+\omega _Y)(\omega _g^2+\omega _L^2),`$
$`a_0`$ $`=`$ $`\omega _\nu (\omega _T\omega _Y\omega _{TY}\omega _{YT})`$ (14)
$``$ $`\omega _g^2(\omega _Y\alpha \omega _{YT})\omega _L^2\left(\omega _T{\displaystyle \frac{1}{\alpha }}\omega _{TY}\right).`$
Equation (13) describes three essentially different modes, corresponding to the three roots of the equation, which can generally exist in a chemically inhomogeneous fluid. The condition that at least one of the roots has a positive real part (unstable mode) is equivalent to one of the following inequalities

$$a_2<0,a_0<0,a_1a_2<a_0$$
(15)
being fulfilled (see, e.g., Aleksandrov, Kolmogorov, & Laurentiev 1985; DiStefano III, Stubberud, & Williams 1994). Since $`\kappa _T`$ and $`\lambda _Y`$ are positive defined quantities, the first condition $`a_2<0`$ will never apply, and only the two other conditions determine the stability of non-magnetic PNSs.

Generally, the last two conditions (15) are rather complex and depend on the horizontal wavevector of perturbations, $`k`$. The critical temperature and lepton gradients which determine the transition from stability to instability can be quite different for perturbations with different $`k`$. In order to obtain the true criteria we have to find the minimal values of these gradients which can cause instability. From the properties of the kinetic coefficients we have $`(\kappa _T\lambda _Y\kappa _Y\lambda _T)>0`$, therefore the gradients in the last two conditions (15) are minimal for perturbations with the wavevector $`k`$ minimizing the quantity $`Q^6/k^2`$. This happens when $`k^2=(\pi /d)^2/2`$. For such a horizontal wavevector, the instability criteria read

$`{\displaystyle \frac{dY}{dz}}(\delta \kappa _T`$ $``$ $`\beta \kappa _Y){\displaystyle \frac{\mathrm{\Delta }T}{T}}(\beta \lambda _Y\delta \lambda _T)`$ (16)
$`<`$ $`{\displaystyle \frac{27\pi ^4}{4gd^4}}\nu (\kappa _T\lambda _Y\kappa _Y\lambda _T),`$
$`{\displaystyle \frac{dY}{dz}}[\delta (\nu +\lambda _Y)+\beta \kappa _Y]`$ $``$ $`{\displaystyle \frac{\mathrm{\Delta }T}{T}}[\beta (\nu +\kappa _T)+\delta \lambda _T]`$ (17)
$`<`$ $`{\displaystyle \frac{27\pi ^4}{4gd^4}}(\kappa _T+\lambda _Y)[(\nu +\kappa _T)(\nu +\lambda _Y)\lambda _T\kappa _Y].`$
Conditions (16)–(17) divide the $`\mathrm{\Delta }TY`$ plane into four regions, which are characterized by different stability properties. The size and configuration of these regions depend on the thermodynamic and kinetic properties of nuclear matter. The detailed analysis of stability of non-magnetic PNSs based on criteria (16) and (17) has been done in MPU. In that paper, as in the present one, we use results from numerical simulations of PNS evolution by Pons et al. (1999) to calculate the different thermodynamical derivatives, diffusion coefficients, and conductivities appearing in the stability criteria. To illustrate the large variety of situations encountered during the Kelvin-Helmholtz phase of a PNS, we have chosen three different cases. The physical conditions and the different transport parameters are summarized in Table 1. Case A corresponds to a layer near the center of the star, with an enclosed mass of 0.05 $`M_{}`$, at an evolutionary time of $`t=2`$ s, when the PNS is in the deleptonization stage and the core is being heated. Cases B and C correspond to two different layers at a time of $`t=20`$ s, during the cooling stage, mainly driven by thermal neutrinos since an important part of the lepton content of the star has already been radiated away. The layer corresponding to case B is located at the same depth as A (0.05 $`M_{}`$) star while the layer in case C has an enclosed baryonic mass of 1 $`M_{}`$. In our calculations we have used the value of $`d`$ given by the pressure height $`d=\mathrm{}d\mathrm{ln}p/dr^1`$.

In Figure 1 we plot the lines of critical stability given by equations (16) and (17) (solid and dotted line respectively) in a $`\mathrm{\Delta }TY`$ plane. The asterisks indicate the particular values of the gradients for each of the three cases considered, the parameters of which are summarized in Table 1. The shaded area covers the region where all the roots of the cubic equation (13) are real. We denote by neutron fingers a region where condition (16) is fulfilled but not condition (17); by convection we denote those regions where both conditions (16) and (17) are satisfied; we denote by semiconvection the case when condition (17) is satisfied but condition (16) is not.

By applying Routh criterium (DiStefano III, Stubberud, & Williams 1994), which states that the number of unstable modes of a cubic equation is given by the number of changes of sign in the sequence

$$\{1,a_2,\frac{a_2a_1a_0}{a_2},a_0\},$$
(18)
we deduce that the number of unstable modes in the regions labelled by convection or neutron finger is always one. Thus, the root corresponding to this mode is real and positive. Notice that, while in the convection case all roots are real, in the neutron finger case, the two stable modes can be either real or complex conjugates. In the region labelled semiconvection we have two unstable modes, which may be either real or complex conjugates. Semiconvection can be, therefore, oscillatory (complex conjugate unstable modes) or non-oscillatory (two real unstable modes). We emphasize that the convention we use to divide the plane $`\mathrm{\Delta }TY`$ according to the instability criteria (16) and (17) differs from that of Bruenn & Dineva (1996). According to their convention, the convectively unstable region is extended to contain that part of our neutron finger and semiconvection regions where all roots are real. In general, the classification in different types of instability is somehow arbitrary, because there is no clear qualitative distinction between neutron fingers and convection: in both cases the unstable root is real. Our classification simply responds to the fact that, in the limit of vanishing diffusivity, criteria (16) and (17) reduce to the classical Rayleigh-Taylor and Schwarzschild criteria, and certainly is not the only valid choice. The important idea to bear in mind is that, when all the kinetic coefficients are of similar magnitude, the distinction among the different types of instabilities is less evident and only the concept of overall stability is meaningful.

### 3.2 Magnetic PNSs

In the general case, $`𝐁0`$, equation (18) can be rewritten as

$$\gamma ^4+q_3\gamma ^3+q_2\gamma ^2+q_1\gamma +q_0=0,$$
(19)
where

$`q_3`$ $`=`$ $`a_2+\omega _m,`$
$`q_2`$ $`=`$ $`a_1+\omega _ma_2+\omega _A^2,`$
$`q_1`$ $`=`$ $`a_0+\omega _ma_1+\omega _A^2(\omega _T+\omega _Y)`$
$`q_0`$ $`=`$ $`\omega _ma_0+\omega _A^2(\omega _T\omega _Y\omega _{TY}\omega _{YT}).`$ (20)
Equation (19) describes four essentially different modes which can generally exist in chemically inhomogeneous fluid in the presence of the magnetic field. The Routh criterium applied to a quartic equation states that the number of roots with positive real part (unstable modes) is given by the number of changes of sign in the sequence

$$\{1,q_3,\frac{q_3q_2q_1}{q_3},q_1\frac{q_3^2q_0}{q_3q_2q_1},q_0\}.$$
(21)
Since $`q_3`$ is always positive, there is, at least, one stable mode. The criteria for hydro–magnetic instability is equivalent to requiring that any one of the following inequalities be satisfied:

$`q_0`$ $`<`$ $`0,`$ (22)
$`q_3q_2q_1`$ $`<`$ $`0,`$ (23)
$`q_1{\displaystyle \frac{q_3^2q_0}{(q_3q_2q_1)}}`$ $`<`$ $`0.`$ (24)
In general, the magnetic field has a stabilizing effect on all three conditions. Consequently, there will be a critical value of $`𝐁`$ above which none of the instability conditions is fulfilled and all modes become stable.

To warm up, we begin analyzing the behaviour of the roots with increasing magnetic field. In Figure 2 we draw a sketch, plotting the real part ($`\gamma `$) of the four modes as a function of $`\omega _A^2`$. Modes with positive $`\gamma `$ are unstable. In order to show more clearly the qualitative effect, we have assigned values to the dissipative frequencies which do not correspond to real values encountered in the early stages of a neutron star’s life. The qualitative behaviour using real coefficients is very similar, although the different scales would not allow one to capture the trends in a simple plot. In the limit of a vanishing magnetic field ($`\omega _A0`$) there are three non-magnetic modes (solid lines), discussed in the previous section, which remain essentially unaltered, and a stable magnetic mode with a value equal to $`\omega _m`$ (dashed line). The situation shown in Figure 2 would correspond to the convection case (1 unstable and 2 stable non-magnetic modes when $`𝐁=0`$). By switching on the magnetic field, the real part of the new magnetic mode (dashes) increases, and becomes positive in the point labelled $`a`$, which corresponds to the physical conditions such that $`q_0=0`$. At that point, the original unstable mode (upper solid line) has been partially damped. Further increasing the magnetic field, the real parts of the originally unstable mode and the new unstable magnetic mode merge, to form a pair of complex conjugate modes (point $`b`$). Convection becomes, thus, oscillatory, due to the stabilizing action of the magnetic field. At a critical magnetic field, given by condition (24), the real part of the unstable modes vanishes (point $`c`$), and for higher magnetic field all modes are stable. The originally stable modes are not seriously affected.

The general trends discussed above (namely: damping of the unstable mode, a new magnetic mode becoming unstable, merging of modes and change to oscillatory convection, and final suppression of instabilities) are also found when realistic kinetic coefficients from numerical simulations are used. To quantify the characteristic values of the magnetic field causing the effects described above, we have taken the parameters corresponding to case B from Table 1. Those modes with a typical wavelength of the order of the pressure height (see below for further discussion), suffer the transitions $`a`$, $`b`$, $`c`$ in Figure 2 when $`B=2.5\times 10^{11}`$ G, $`1.8\times 10^{16}`$ G and $`2.09\times 10^{16}`$ G, respectively. Thus, a field larger than $`2.09\times 10^{16}`$ G entirely suppresses all types of instability in the considered region.

In order to understand the effect of the magnetic field on each instability condition, it is more convenient to rewrite equations (2224) separating the terms containing the Alfvén frequency and magnetic dissipation:

$$a_0<\frac{\omega _A^2}{\omega _m}(\omega _T\omega _Y\omega _{TY}\omega _{YT}),$$
(25)
$$a_2a_1a_0<\omega _ma_2(a_2+\omega _m)\omega _A^2(\omega _\nu +\omega _m),$$
(26)
$`[a_0+\omega _ma_1+\omega _A^2(\omega _T+\omega _Y)][a_2a_1a_0+\omega _ma_2(a_2+\omega _m)+\omega _A^2(\omega _\nu +\omega _m)]`$
$`<(a_2+\omega _m)^2\times [(\omega _ma_0+\omega _A^2(\omega _T\omega _Y\omega _{TY}\omega _{YT})].`$ (27)
These inequalities can be further simplified if we take into account that the frequencies entering these conditions are of different orders of magnitude in PNSs. The dissipative frequencies, $`\omega _\nu `$, $`\omega _T`$, $`\omega _Y`$, $`\omega _{TY}`$, and $`\omega _{YT}`$ are approximately comparable to each other and range between $`1010^3`$ s<sup>-1</sup>. The magnetic dissipative frequency, $`\omega _m`$, is much smaller. The electrical resistivity of hot nuclear matter ($``$) is relatively low, $`6\times 10^{29}T_8^2`$ s (Yakovlev & Shalybkov, 1991), thus the characteristic time scale of dissipation of the magnetic field is extremely long compared to other dissipative time scales. Taking into account that $`\nu _m=c^2/4\pi `$, we can estimate $`\omega _m10^{11}`$ s<sup>-1</sup> for $`k\pi /\mathrm{}`$. Therefore, we can neglect $`\omega _m`$ when it enters conditions (23) and (24) in a row with other dissipative frequencies. In addition, due to the small magnetic diffusivity, it also turns out that $`\omega _m\omega _\nu \omega _A^2`$ and equation (22) can be written as

$`\omega _g^2(\omega _Y\alpha \omega _{YT})\omega _L^2\left(\omega _T{\displaystyle \frac{1}{\alpha }}\omega _{TY}\right)+{\displaystyle \frac{\omega _A^2}{\omega _m}}(\omega _T\omega _Y\omega _{TY}\omega _{YT})<0.`$ (28)
All three conditions depend on the wavevector $`k`$. It is well known (Chandrasekhar, 1961) that in the inviscid limit the minimal values of gradients which can cause instability of a magnetized fluid correspond to $`k\mathrm{}`$. Taking account of viscosity increases the most unstable horizontal wavelength and the critical gradients are reached for a large, though finite, value of $`k`$. In the case of low magnetic dissipation, $`k`$ is much larger than $`\pi /d`$. Substituting $`kQ`$ into equation (35), we have

$`{\displaystyle \frac{dY}{dz}}(\delta \kappa _T\beta \kappa _Y){\displaystyle \frac{\mathrm{\Delta }T}{T}}(\beta \lambda _Y\delta \lambda _T)<{\displaystyle \frac{\pi B^2}{4g\rho d^2\nu _m}}(\kappa _T\lambda _Y\kappa _Y\lambda _T).`$ (29)
Due to the particular physical conditions in PNSs ($`\omega _m1`$), even a relatively low magnetic field can make this condition not satisfied. We must point out, however, that this criterium is associated to the change of sign of new the magnetic mode (dashed line in Figure 2), which is irrelevant, because the dominant unstable mode is always that associated to one of the other two conditions.

Let us turn back to the other two conditions, (26) and (27). The dissipative frequencies are usually small compared to the dynamical frequencies, $`\omega _g`$ and $`\omega _L`$, for typical lengthscales in a PNS. For example, the contribution of terms proportional to the cube of dissipative frequencies in equations (26)–(27) does not exceed 10%. Therefore, we can neglect the terms containing only dissipative frequencies compared to those containing the dynamical ones (including the Alfvén frequency). This latter approximation is justified if the magnetic field is large enough for $`\omega _A`$ to be larger than any dissipative frequency, $`\omega _A>\omega _{\mathrm{diss}}=\mathrm{max}(\omega _\nu ,\omega _T,\omega _Y,\omega _{TY},\omega _{YT})`$. Since $`\omega _{\mathrm{diss}}0.11`$ s<sup>-1</sup>, the condition $`\omega _A>\omega _{\mathrm{diss}}`$ is fulfilled if the magnetic field exceeds some characteristic value which is of the order of $`10^{12}10^{13}`$ G. For $`B<10^{13}`$ G the influence of the magnetic field is insignificant and conditions (26)–(27) are exactly equivalent to conditions (16)–(17). For $`B>10^{13}`$ G, the instability conditions can be simplified as follows:

$`\omega _g^2(\omega _Y\alpha \omega _{YT})\omega _L^2\left(\omega _T{\displaystyle \frac{\omega _{TY}}{\alpha }}\right)+\omega _A^2(\omega _T+\omega _Y)<0,`$ (30)
$`\omega _g^2(\omega _T+\omega _\nu +\alpha \omega _{YT})\omega _L^2\left(\omega _Y+\omega _\nu +{\displaystyle \frac{\omega _{TY}}{\alpha }}\right)+\omega _\nu \omega _A^2<0.`$ (31)
As for equation (28), the minimal values of the gradients which cause the instability correspond to $`kQ`$ (short wavelengths). Substituting $`k=Q`$ into equations (30) and (31), we have

$`{\displaystyle \frac{dY}{dz}}(\delta \kappa _T\beta \kappa _Y){\displaystyle \frac{\mathrm{\Delta }T}{T}}(\beta \lambda _Y\delta \lambda _T)<{\displaystyle \frac{\pi B^2}{4g\rho d^2}}(\kappa _T+\lambda _Y),`$ (32)
$`{\displaystyle \frac{dY}{dz}}[\delta (\nu +\lambda _Y)+\beta \kappa _Y]{\displaystyle \frac{\mathrm{\Delta }T}{T}}[\beta (\nu +\kappa _T)+\delta \lambda _T]<{\displaystyle \frac{\pi \nu B^2}{4g\rho d^2}}.`$ (33)
Note that these criteria are analogous to the instability criteria for non-magnetic PNSs (16)–(17), but replacing

$`\nu (\kappa _T\lambda _Y\kappa _Y\lambda _T){\displaystyle \frac{d^2B^2}{27\pi ^3\rho }}(\kappa _T+\lambda _Y),`$ (34)
$`(\kappa _T+\lambda _Y)[(\nu +\kappa _T)(\nu +\lambda _Y)\lambda _T\kappa _Y]{\displaystyle \frac{d^2B^2}{27\pi ^3\rho }}\nu ,`$ (35)
and the magnetic field plays the role of an efficient dissipative process. In magnetic PNSs the dissipative effect caused by the Lorentz force is much greater than neutrino transport or neutrino viscosity effects.

In Figure 3 we show the unstable zones at each time of the evolution of a young PNS from one of the numerical simulations in Pons et al. (1999). Top and bottom panels correspond to criteria (32) and (33), respectively. Criterion (33) depends on the neutrino viscosity, for which available estimates are sometimes contradictory (van den Horn & van Weert, 1981; Goodwin & Pethick, 1982; Thompson & Duncan, 1993). In our simulations we used the value $`\nu =\eta ^2\kappa _T`$ (van den Horn & van Weert, 1981). The contours indicate the convectively unstable region for the following values of a passive magnetic field: $`10^{14}`$ G (solid), $`5\times 10^{15}`$ G (dotted), $`10^{16}`$ G (dashed), $`3\times 10^{16}`$ G (dash-dot) and $`5\times 10^{16}`$ (dash-dot-dot-dot). It turns out that, with this viscosity, the critical magnetic fields for criteria (32) and (33) are of the same order of magnitude. The magnetic correction on (33) is somewhat more restrictive and, generally, its associated instability is suppressed with weaker magnetic fields. The difference is caused mainly by the different combinations of kinetic coefficients that appear in the equations. Both criteria require a rather strong magnetic field, $`B(15)\times 10^{16}`$G, in order to stabilize the fluid completely, although the region where condition (33) is fulfilled is typically smaller than that given by condition (32). Note that moderate magnetic fields in the order of $`10^{14}`$ G leave the unstable regions essentially unchanged, and the effects are visible only for $`B10^{15}`$ G.

In Figure 4 we plot the critical magnetic field ($`B_{\mathrm{st}}`$) which completely stabilizes a certain layer of the PNS (corresponding to enclosed baryonic masses of 0.2, 0.5, 1.0, and 1.4 $`M_{}`$) as a function of time. Layers with $`M_\mathrm{B}<1M_{}`$ are stable at the beginning of the evolution. Later on, when the corresponding layer becomes unstable, $`B_{\mathrm{st}}`$ grows rapidly, reaching some maximum value, and then decreases monotonously as the temperature and lepton gradients are smeared out via neutrino diffusion. The maximum value of $`B_{\mathrm{st}}`$ at different points varies in a very narrow range ($`25\times 10^{16}`$ G), but this maximum is reached at different instants in time. Note that $`B_{\mathrm{st}}`$ provides stabilization in the largest unstable lengthscale. Stabilization at smaller scales is accomplished with weaker magnetic fields. The instability seems to be especially persistent in the central region with $`M_\mathrm{B}<0.5M_{}`$, where convective motions might last for a long time (Keil, Janka & Müller, 1996) and $`B_{\mathrm{st}}5\times 10^{16}`$.

## 4 Discussion

In this paper we have considered the effect of magnetic field (consistently coupled to other dissipative processes, such as neutrino transport and neutrino viscosity) on the convective instability of PNSs. Our results show that transport phenomena play an important role in determining the size of the region of the star which is subject to instability. The particular values of the transport and kinetic coefficients also determine the kind of instabilities which can be present in an unstable region (convection, neutron fingers, semiconvection).

For low magnetic fields ($`B<10^{13}`$ G) the region where instabilities are present remains unaffected when compared to the non-magnetic case. That region covers an important part of the PNS and lasts for almost the whole Kelvin-Helmholtz phase of the PNS evolution. Consequently, due to the existence of such an extended region subject to instabilities, convective motions are likely to appear and turbulent dynamo action can significantly increase, locally, the magnetic field strength.

Once the magnetic field has a significant strength ($`10^{13}`$ G $`<B<10^{15}`$ G), its effect is twofold. First, it is responsible for the appearance of a new mode which can generally be unstable, although in this case, its growth time is always longer than that of the other unstable modes. Secondly, the magnetic field acts as a stabilizing agent for all types of instability studied in the non-magnetic case, progressively reducing the size of the convective regions. Our calculations show, however, that only a relatively strong magnetic field can influence significantly the stability properties of PNSs, particularly in the central region with $`M_\mathrm{B}<0.5M_{}`$. For instance, if the dynamo action is sufficiently efficient, the duration of the convective epoch can be noticeably shorter as the magnetic field increases up to $`10^{16}`$G. We find that the critical magnetic field that suppresses all sorts of instability happens to be in a narrow range, at about $`5\times 10^{16}`$ G. This can be considered as an upper limit to the poloidal magnetic field generated by convective dynamo.

In our simplified approach, we have considered the case of a vertical magnetic field. It is known (Chandrasekhar, 1961) that this magnetic configuration requires the maximal field strength for stabilization. In real conditions, the field geometry may be more complex with a substantial component perpendicular to the temperature and composition gradients. In PNSs, this field component can correspond, for example, to a toroidal magnetic field. For this geometry, the stabilizing magnetic field is typically several times lower. Therefore, stabilization against hydromagnetic instabilities can really take place at a slightly weaker magnetic field than that obtained in our calculations. In the presence of differential rotation (a very likely situation after the core collapse), the main contribution to the total field strength would be given by the toroidal magnetic field ($`B_\mathrm{t}`$) and the poloidal component ($`B_\mathrm{p}`$) would be less important. In a quasi-stationary regime, the rate of stretching of the toroidal field lines from the poloidal configuration by differential rotation is compensated by turbulent dissipation (Parker, 1979). By equating both rates one can write the following equation:

$$\frac{B_\mathrm{p}}{\tau }\frac{v_T\mathrm{}_T}{3}^2B_\mathrm{t},$$
(36)
where $`\tau `$ is the characteristic stretching time, $`v_T`$ the turbulent velocity and $`\mathrm{}_T`$ the turbulent lengthscale. Since the stretching is due to differential rotation, $`\tau `$ might be estimated as $`\tau =\mathrm{\Delta }\mathrm{\Omega }^1`$, $`\mathrm{\Delta }\mathrm{\Omega }`$ being the difference between the maximum and minimum angular velocity. An estimation of $`^2B_\mathrm{t}`$ is $`B_\mathrm{t}/R^2`$, $`R`$ being the typical radius of the star. Using the following values, $`R20`$ km, $`v_T10^8`$ cm/s and $`\mathrm{}_T10^5`$cm, we obtain

$$\frac{B_t}{B_p}\frac{\mathrm{\Delta }\mathrm{\Omega }}{1\mathrm{s}^1}.$$
(37)
Assuming that neutron stars are born with large angular velocities (say $`\mathrm{\Omega }=10^3`$ s<sup>-1</sup>), and that $`\mathrm{\Delta }\mathrm{\Omega }/\mathrm{\Omega }`$ is comparable to that of the Sun ($`0.10.3`$), we can estimate the toroidal field to be $`100300`$ times stronger than the poloidal component. Since the total magnitude of the magnetic field generated by dynamo action should not exceed, in any case, $`10^{16}`$ G we have to conclude that the poloidal component which may represent the observed pulsar field must be in the range $`3\times 10^{13}10^{14}`$ G. This estimate is in agreement with the observed fields of young pulsars in supernova remnants. The toroidal field could be, however, much stronger, which might lead to important consequences on the structure of neutron stars (Bocquet et al., 1995; Cardall, Prakash, & Lattimer, 2001).

This work was supported by the Spanish Ministerio de Ciencia y Tecnología (grant AYA 2001-3490- C02-02), the EU Program ‘Improving the Human Research Potential and Socio-Economic Knowledge Base’ (Research Training Network contract HPRN-CT-2000-00137), the Russian Foundation of Basic Research and Deutsche Forschungsgemeinschaft (grant 00-02-04011). V.U. also thanks the Ministerio de Educacion, Cultura y Deporte of Spain for financial support under the grant SAB1999-0222. J.P. is supported by the Marie Curie Fellowship No. HPMF-CT-2001-01217. The authors acknowledge to an anonymous referee for useful suggestions.

Stability plane for the three thermodynamic states detailed in Table 1. Four different regions in the plane are separated by criteria (16)-(17) (solid and dotted lines, respectively). The shaded area indicates the region where the three roots of equation (13) are real. The particular values of the gradients in each studied case (A,B,C) are indicated by asterisks. The scale factor $`\mathrm{}`$ has been taken as the pressure height of each case.

Sketch showing the behaviour of the real part of the roots of equation (19) as a function of the square of the Alfvén frequency, for fixed values of the dynamical and dissipative frequencies. Dashed line corresponds to the mode decoupled in the non-magnetic case (referred to as magnetic mode in the text). Units are arbitrary and the points $`a`$, $`b`$ and $`c`$ indicate changes of sign and character (real–complex) of the modes. See §3.1 for further details.

Contour plots showing the variation of the different convectively unstable zones, assuming different strengths of the magnetic field. The contour lines correspond to the following values of the magnetic field: $`10^{14}`$ G (solid), $`5\times 10^{15}`$ G (dotted), $`10^{16}`$ G (dashed), $`3\times 10^{16}`$ G (dash-dot) and $`5\times 10^{16}`$ (dash-dot-dot-dot). The top panel refers to the unstable region according to criteria (32) and the bottom panel refers to (33). The x-axis indicates the evolution time and the y-axis the enclosed baryonic mass of the layer. The innermost part of the star remains unstable at intermediate times even when magnetic fields as large as 10<sup>16</sup> G are present. The asterisks correspond to the cases A, B and C described in the text. This PNS model was taken from numerical simulations performed in Pons et al. (1999).

Critical magnetic field strength that completely stabilizes the fluid as a function of the evolution time. The four lines stand for four different layers, with enclosed baryonic masses of 0.2 (solid), 0.5 (dotted), 1.0 (dashed) and 1.4 (dash-dot) $`M_{}`$, respectively. The PNS model is the same as in Figure 3. In general, the deeper is the layer, the larger magnetic field is required for unconditional stabilization.


format: text

text:
# Forcing and Velocity Correlations in a Vibrated Granular Monolayer

## Abstract

The role of forcing on the dynamics of a vertically shaken granular monolayer is investigated. Using a flat plate, surprising negative velocity correlations are measured. A mechanism for this anti-correlation is proposed with support from both experimental results and molecular dynamics simulations. Using a rough plate, velocity correlations are positive, and the velocity distribution evolves from a gaussian at very low densities to a broader distribution at high densities. These results are interpreted as a balance between stochastic forcing, interparticle collisions, and friction with the plate.

Granular gases, systems of large numbers of macroscopic grains in rapid motion and interacting through dissipative collisions, appear in a wide range of industrial applications and natural phenomena. Energy must be supplied externally to compensate for the inelastic collisions, so granular gases are necessarily systems out of equilibrium. As a result, they may display dramatic non-equilibrium effects such as non-gaussian velocity distributions OlafsenUrbach ; Losert99 ; Rouyer00 ; VanNoije98 ; ChrisBizon ; NieBenNaimChen  and long-range spatial velocity correlations VanNoije99 ; Pagonabarraga01 ; MoonShattuckSwiftPRE2001 ; BlairKudrolliRPRE2001 . Non-gaussian velocity distributions are a direct demonstration of the inapplicability of the Gibbs distribution, and significant velocity correlations indicate the absence of ‘molecular chaos’, which is a crucial approximation normally used to solve the Boltzmann equation and to calculate other fundamental quantities in kinetic theory. Recent theoretical work has focused on the non-equilibrium steady state obtained when the energy supplied by spatially homogeneous random external forcing is balanced by the dissipation due to the collisions such that the average energy of the system remains constant VanNoije98 ; VanNoije99 ; Pagonabarraga01 . Non-gaussian velocity distributions and algebraically decaying velocity correlations arise as a direct consequence of the energy injection. Non-gaussian velocity distributions and velocity correlations have been observed in a number of experiments OlafsenUrbach ; Losert99 ; Rouyer00 ; BlairKudrolliRPRE2001 , but in each case the forcing is sufficiently different from that of the theoretical models that a direct comparison is difficult. In this Letter, we provide a direct demonstration of the determining role that the forcing plays on the spatial velocity correlations in a homogeneously forced granular gas. We will also show how the framework of kinetic theory can provide a coherent description of the origin of the observed velocity distributions and correlations.

We have investigated a quasi-2D granular system consisting of a layer of a large number of spherical particles partially covering a vertically driven horizontal plate. Two plates with different surface properties were used in the experiments: a smooth circular plate (20 cm in diameter) made of black anodized aluminum, and a rough hexagonally shaped plate (30 cm between opposite corners). The roughness of the latter is provided by a close-packed lattice of blackened steel balls (1.19 mm diameter) glued to a flat plate. The granular gas is made of uniform stainless steel spheres with diameter $`\sigma =1.59`$ mm (smooth plate), or $`\sigma =3.97`$ mm (rough plate). On the smooth plate, the layer is constrained from above by an anti-static coated plexiglas lid, 1.7 $`\sigma `$ above the plate. Using an electromagnetic shaker, the plate is driven sinusoidally, and the granular layer is brought to a non-equilibrium steady state. The strength of the shaking $`\mathrm{\Gamma }=A\omega ^2/g`$ is directly measured with a fast response accelerometer ($`A`$ is the amplitude and $`\omega `$ the frequency of the plate oscillation, and $`g`$ the acceleration due to gravity).

A high resolution camera (Pulnix TM1040) placed above the plate records the bright spots on the tops of the balls produced by a stroboscopic LED array. Instantaneous horizontal velocities $`𝐯_i`$ are obtained from the displacements between strobe pulses. The strobe is synchronized with the camera so that the first pulse occurs at the end of the exposure of one frame, and the next pulse occurs at the beginning of the subsequent frame. The time interval between pulses ($``$ 1 ms) is chosen to be significantly smaller than the mean collision time. The width of both pulses is typically less than 0.05 ms. Displacements of the centers of the bright spots are measured with 0.1 pixel accuracy. To eliminate systematic errors in the velocity correlations, snapshots are taken at a fixed phase relative to the driving signal corresponding to the maximum average position of the layer.

The general behavior of the granular layer on the smooth plate has been reported previously OlafsenUrbach  and a rich phase diagram depending on both $`\mathrm{\Gamma }`$ and $`\omega `$ was found. Here we focus mainly on the fluidized regime, obtained for $`\mathrm{\Gamma }1`$. All results presented in this Letter were obtained with $`\mathrm{\Gamma }=1.5`$ and $`\omega /2\pi =60`$ Hz.

The longitudinal and transverse velocity correlations, $`C_{||}`$ and $`C_{}`$, respectively, are calculated by

$`C_{||,}(r)={\displaystyle \underset{ij}{\overset{N_r}{}}}v_i^{||,}v_j^{||,}/N_r,`$
where the sum runs over the $`N_r`$ pairs of particles separated by a distance $`r`$, $`v_i^{||}`$ is the projection of $`𝐯_i`$ along the line connecting the centers of particles $`i`$ and $`j`$, and $`v_i^{}`$ is the projection perpendicular to that line.

Figure 1 shows $`C_{||}`$ (inset: $`C_{}`$) normalized by the granular temperature $`T=𝐯_i^2`$ obtained on the smooth plate for $`\rho =0.4`$ and 0.5, and $`\mathrm{\Gamma }=1.5`$ coverage . Surprisingly, the longitudinal velocity correlations are strongly negative for all $`r`$ and all accessible densities: the velocities are anti-correlated. $`C_{||}`$ decays slower than exponentially, and in the range $`r=[2\sigma ,16\sigma ]`$ can be reasonably well described by $`C_{||}(r)r^2`$. The structure visible between $`1\sigma `$ and $`2\sigma `$ is stronger at high density, and its origin is not understood. The transverse component $`C_{}`$ is much smaller and has a shorter range than $`C_{||}`$ (note the difference in vertical scale). Nearly identical correlations were observed without the lid in place. We have performed molecular dynamics simulations that closely match the conditions of the experiment computer  and find qualitatively similar correlations.

The existence of long-range velocity correlations in a randomly forced granular gas has been explained theoretically using a hydrodynamic approach VanNoije99 . The velocity correlations result from the balance between the excitation of hydrodynamic modes by the random external forcing and their dissipation through diffusion. At shorter range, velocity correlations arise from an increased probability of re-collision due to the forcing Pagonabarraga01 . This suggests that the correlations could be quite sensitive to the particular mechanism of energy injection, and we believe that this is the origin of the dramatic difference between the negative correlations observed in our experiment and the positive correlations observed in the randomly forced model systems VanNoije98 ; VanNoije99 ; Pagonabarraga01 ; MoonShattuckSwiftPRE2001 . In those systems, energy is injected through white noise forcing of each particle independently, while on the smooth plate energy is injected only into the vertical motion of the spheres and is then transferred to the horizontal motion through inter-particle collisions. As a result of this transfer, collisions will result, on average, in an increase in the magnitude of the relative horizontal velocities $`𝐯_r`$ of pairs of particles. This effect can be measured directly in the molecular dynamics simulations, and at $`\rho =0.5`$ we find that the average value of $`𝐯_r^2`$ for pairs of particles leaving a collision is 13% higher than for pairs entering a collision. Experimental support for this scenario can be found in Fig. 2, which shows a calculation of $`𝐯(0)𝐯(r)`$, where the average includes only particles restricted to a narrow band along the direction of motion of the particle at the origin, either ahead (for positive $`r`$) or behind (for negative $`r`$). There is a clear asymmetry with respect to zero, showing that balls moving away from each other contribute significantly more to the velocity anti-correlation than do particles moving towards each other. The velocity anti-correlations would presumably dissipate through diffusion as the momentum gets transferred to surrounding particles through collisions, thus it is not surprising that the observed correlations extend over several mean free paths. This mechanism for the generation of anti-correlated velocities may also operate in inhomogeneously forced granular media when the average kinetic energy of the grains is anisotropic.

To test our hypothesis about the origin of the correlations, we have changed the forcing by using a rough surface, so that the vibrating plate also injects energy directly into the horizontal motion of the spheres. This energy input mechanism is closer to the white noise forcing of Ref. VanNoije99 , and allows for a more direct comparison between experiment and theory chaotic .

The velocity correlations measured on the rough plate are shown on Fig. 3 for $`\mathrm{\Gamma }=1.5`$ Cluster , and two extremal densities. We find that both $`C_{||}`$ and $`C_{}`$ are positive for all measured $`\rho `$. $`𝐯(0)𝐯(r)`$, calculated as described above, is positive everywhere and is roughly symmetric about zero. Molecular dynamics simulations of the layer on the rough plate computer  show similar correlations, and the effect of the collisions is to decrease the average value of $`𝐯_r^2`$ by 3.7% at $`\rho =0.5`$. Thus in the absence of large vertical to horizontal collisional energy transfer, the velocity anti-correlations are replaced by positive velocity correlations, presumably arising from the random forcing as described in Refs. VanNoije99 ; Pagonabarraga01 . The dramatic difference in the correlations between the smooth and rough plates is a clear indication that the forcing, rather than the inelastic interparticle collisions, determines the velocity correlations.

Interestingly, we find that $`C_{||}`$ and $`C_{}`$ can be roughly scaled on two different curves of the form $`C_{||,}=T\alpha (\rho )f_{||,}(r)`$ (Fig. 4). The functions $`f_{||}`$ and $`f_{}`$ have different $`r`$ dependence, but the function $`\alpha (\rho )`$ is the same for both $`C_{||}`$ and $`C_{}`$. $`f_{||}`$ and $`f_{}`$ are reasonably well described by $`f(r)e^{r/r_o}`$, with $`r_o=2\sigma `$ for $`f_{||}`$ and $`0.6\sigma `$ for $`f_{}`$. It is surprising that this decay length has little density dependence, despite the fact that the mean free path estimated from Enskog-Boltzmann kinetic theory varies from 2.3$`\sigma `$ for $`\rho =0.125`$ to 0.16$`\sigma `$ for $`\rho =0.6`$. A similar scaling relation is predicted in the randomly forced model at large $`r`$ VanNoije99 , but the scaling is the same for both $`C_{||}`$ and $`C_{}`$, and $`f(r)`$ decays algebraically.

An important difference between our experiment and the randomly forced model is revealed by the velocity distributions $`P(v)`$ of the balls on the rough plate. Figure 5 shows the normalized $`P(v_x)`$ obtained for $`\mathrm{\Gamma }=1.5`$, and $`\rho `$ ranging from $`0.014`$ to $`0.6`$. At the lowest density, $`P(v_x)`$ is well described by a gaussian. As $`\rho `$ increases, the tails of the distribution rise significantly above a gaussian. At high density, the velocity distribution resembles those reported previously on flat plates OlafsenUrbach ; Losert99 . To better understand the effect of density on $`P(v)`$, we have followed the trajectories of individual balls using a fast camera (838 frames/s, OlafsenUrbach ), and measured the time-dependent velocity autocorrelation function $`C_v(t)=𝐯(0)𝐯(t)`$ and find that $`C_v(t)e^{t/\tau (\rho )}`$ at all densities. The scattering length defined by $`l\sqrt{T}\tau (\rho )`$ varies from about 0.5$`\sigma `$ at $`\rho =0.6`$ to about 2.5$`\sigma `$ at $`\rho =0.014`$. At the lowest density, $`l`$ is much smaller than the mean free path for ball-ball scattering ($``$ 20 $`\sigma `$) obtained from kinetic theory. This suggests that it is ball-plate scattering that is responsible for the decay of the velocity auto-correlation function, and that the scattering must therefore have a viscous component. The evolution of $`P(v)`$ for a single ball on the rough plate can thus be very likely modelled with an equation including only random white noise forcing and a viscous term, such as the Fokker-Plank equation for Brownian motion, which produces a gaussian $`P(v)`$. When the random forcing is instead balanced by ball-ball scattering in a modified Boltzmann equation, there is an over-population of the tails of the distributions, $`P(v)\mathrm{exp}(|v|^{3/2})`$ VanNoije98 . This is consistent with the experimentally measured $`P(v)`$ at high densities, where the scattering is dominated by ball-ball collisions. At intermediate densities both scattering processes contribute and the evolution of $`P(v)`$ with density can likely be described with a Boltzmann equation consisting of random white noise forcing, collisional scattering, and a density independent viscous drag. This model has been recently investigated Puglisi99 , although primarily at higher inelasticities, where significant clustering is observed.

The ball-plate scattering may also explain the lack of significant density dependence of the decay length of the velocity correlations. We speculate that the density independent scattering with the plate, which unlike ball-ball collisions does not conserve momentum, is more effective at destroying the long-range velocity correlations, and therefore controls the observed decay. Support for this scenario is provided by the fact that the decay length for the exponentially decaying $`C_{||}(r)`$ at all measured densities is close to the ball-plate scattering length deduced from the individual trajectories at very low densities (2$`\sigma `$ versus $`2.5\sigma `$). However, a full quantitative understanding of the velocity correlations observed on the rough plate will require a theoretical analysis along the lines of Ref. VanNoije99 , including the effect of viscous drag. Our results provide a clear experimental demonstration that the velocity distribution and correlations in a granular gas are not determined exclusively by internal dynamics (inelastic ball-ball collisions) but by a balance between internal dynamics and external forcing.

The authors gratefully acknowledge Mickey Moore for providing us with the rough plate. This work was supported by grants from the Research Corporation, the Petroleum Research fund, DOE, and Grant No. DMR-9875529 from the NSF. J.S.U. acknowledges support from the Sloan foundation.