text:
# Characteristics of the fragments produced in central collisions of 129Xe+natSn from 32 to 50 AMeV

## I Introduction

Multiple intermediate mass fragment (IMF) production in central heavy-ion collisions is related to the properties of nuclear matter under extreme conditions. Many different models have been proposed in order to explain the observed fragment production, and both theoretically and experimentally the situation is not clear. Models with widely differing basic hypotheses can be equally good at describing the same data such as charge distributions, mean energies and angular distributions. In order to gain further understanding it is therefore necessary to have more detailed information on the multifragmentation process.

One aspect of the reactions for which different models give very different predictions are the excitation energies of what we will call from now on the ‘primary fragments’: in other words the nuclei present around $``$100 fm/c after the collision, which are not necessarily the same as those arriving in the detectors a few 10’s of ns later. In Quantum Molecular Dynamics Aichelin (1991); Nebauer and Aichelin (1999); Nebauer et al. (1999); Tirel (1998) simulations or Microcanonical Metropolis Monte-Carlo model Gross (1990) calculations the primary fragments are rather cold i.e. they are almost unaffected by subsequent secondary decays and arrive unchanged in the detectors. In the former case, the (lack of) excitation energy in the nascent fragments is determined by the collision dynamics, whereas in the latter case it is an assumption of the model when calculating the statistical weights of the partitions. On the other hand, Antisymmetrised Molecular Dynamics Ono et al. (1992); Ono and Horiuchi (1996); Ono (1999) and Stochastic Mean Field Chomaz et al. (1994); Guarnera et al. (1997) simulations both predict moderately “hot” primary fragments in reactions around the Fermi energy, with $`E_{pr}^{}23`$ AMeV Hudan (2001); Hudan et al. (in preparation); Frankland et al. (2001a). Finally, the Statistical Multifragmentation Model (SMM) Bondorf et al. (1995) and the microcanonical multifragmentation model of Raduta and Raduta (1997a, b) allow primary fragments to be excited, the actual value in any given calculation being determined by energy conservation and the statistical weight given by the associated level density parametrisation. This latter may or may not take into account the level density limitation in isolated nuclei at high excitation Dean and Mosel (1985), equivalent to excluding from the primary partitions levels with very short lifetimes or introducing an effective limiting (maximum) temperature for hot nuclei Levit and Bonche (1985); Koonin and Randrup (1987).

Our previous experimental work Marie et al. (1998) has shown that the reconstruction of the average size and excitation energy of the primary fragments is possible by means of fragment-light charged particles (IMF-LCP) relative velocity correlation functions. A constant value of the excitation energy of the primary fragments has been deduced at about 3 AMeV for the Xe + Sn system at 50 AMeV. It was also possible to deduce the multiplicities of the secondary particles evaporated by the primary fragments. More recently analogous results and conclusions have been obtained for central collisions of Kr + Nb at 45 AMeV Staszel et al. (2001). An important question arises from these studies : what is the evolution of the fragment excitation energies and secondary LCP multiplicities as a function of incident energy ? The experimental answer to this question may permit to distinguish between different scenarios and assumptions made by different models. It should give a strong test of the validity of some of their basic hypotheses.

In this paper we extend the previous study Marie et al. (1998) to a wider incident energy range, from 32 to 50 AMeV for central collisions of the Xe + Sn system measured with the 4$`\pi `$ INDRA detector Pouthas et al. (1995a); Steckmeyer et al. (1995); Pouthas et al. (1995b). Excitation functions for the fragment excitation energy and the fraction of secondary emitted LCP correlated to the fragments will be shown. We will give in section 2 a brief description of the detector, the way we select the events and an overview of the fragment production. We will describe in section 3 the method employed to extract the LCP’s correlated to each fragment. The method used for the decorrelation in this work is different from the previous one Marie et al. (1998) but gives almost the same results. The experimental results are then given in section 4. In section 5 a comparison of the deduced primary excitation energy and secondary LCP multiplicities to SMM calculations is given. We then discuss the results in section 6.

## II Experiment

### II.1 Experimental set-up

The experiment was performed at GANIL with the multidetector INDRA Pouthas et al. (1995a); Steckmeyer et al. (1995); Pouthas et al. (1995b). This charged product detector covers about 90% of the 4$`\pi `$ solid angle. The total number of detection cells is 336 arranged according to 17 rings centred on the beam axis. The first ring ($`2^o`$-$`3^o`$) is made of fast NE102/NE115 phoswich detectors. Rings 2 to 9 cover the angular range from $`3^o`$ to $`45^o`$ and are made of three detector layers : a low pressure gas-ionisation chamber, a 300 $`\mu m`$ thick silicon detector and a 14 to 10 cm thick CsI(Tl) scintillator. The remaining 8 rings cover the angular range from $`45^o`$ to $`176^o`$ and have two detection layers : ionisation chamber and 7.6 to 5 cm thick CsI(Tl) scintillators. For the studied system Xe + Sn, fragments with Z up to 54 are identified in the forward region. Beyond $`45^o`$, the charge resolution is one unit up to Z=16 and few charges above. Over the whole angular range, a very good isotope identification is obtained for Z=1 to Z=3, except for particles with low energies where ambiguities are unresolved.

The energy calibration of the CsI(Tl) scintillators was obtained for light charged particles (LCP) by means of the elastic and inelastic scattering of secondary LCP beams ($`p,d,t`$,$`{}_{}{}^{3}He`$,$`{}_{}{}^{4}He`$) produced by the fragmentation of a 95 AMeV 16O beam on a thick C target. These particles were then momentum selected by the “alpha magnetic spectrometer” of GANIL and scattered in a C or Ta target installed in the INDRA reaction chamber. For $`Z3`$ fragments, the energy calibration was made by using the $`\mathrm{\Delta }E`$/E technique. A typical energy resolution was about $`4\%`$. The energy threshold was a few 100 keV for light particles, 0.7 AMeV for Z=3 and 1.4 AMeV for Z=35. A complete technical description of INDRA, its calibration and its electronics can be found in Pouthas et al. (1995a); Steckmeyer et al. (1995); Pouthas et al. (1995b); Tăbăcaru et al. (1999); Pârlog et al. (2002a, b).

### II.2 Selection of central collisions

Two selections have been made to isolate central collisions. The first one is the requirement of quasi-complete events by accepting in the off-line analysis only events having total detected charge $`(Z_{tot})`$ $``$ $`80\%`$ of the initial total charge of the system. The second is the use of the flow angle ($`\theta _{flow}`$) selection Lecolley et al. (1994). This angle is a global observable defined as the angle between the beam axis and the main direction of emission of matter in each event as determined by the energy tensor calculated from fragment ($`Z3`$) c.m. momenta Cugnon and L’Hote (1983). It has been shown for heavy ion reactions in the Fermi energy range Lecolley et al. (1994); Marie et al. (1997); Frankland et al. (2001b) that events with small $`\theta _{flow}`$ are dominated by binary dissipative collisions. On the other hand events with little or no memory of the entrance channel should be isotropic, thus favouring large $`\theta _{flow}`$ ($`P(\theta _{flow})\mathrm{sin}\theta _{flow}`$). Quasi-complete events having $`\theta _{flow}45^{}`$ for 50 AMeV bombarding energy and $`\theta _{flow}60^{}`$ for the three other systems correspond to an isotropic emission of the IMF in the centre of mass of the whole system. These events are compatible with decay of a compact object which could take place after fast emission of a direct light particle component. Indeed the velocity of the fragments are evenly distributed around the centre of mass velocity Salou (1997). By taking into account the detection efficiency and other biases due to the selection we have estimated the cross sections for ‘isotropic central collisions’ to decrease from 115$`\pm `$20 mb at 32 AMeV to 85$`\pm `$10 mb at 50 AMeV. More details about this event selection for Xe + Sn collisions at 32–50 AMeV incident energy and the extraction of the cross sections can be found in Salou (1997).

### II.3 Overview of fragment production in central collisions

Before determining the characteristics of the fragments, let us first show an overview of their production in central collisions of Xe + Sn from 32 to 50 AMeV. Figure 1 shows their charge distributions normalised to the number of events so that the four bombarding energies can

be compared. The production of small fragments ($`Z10`$) increases with incident energy. For the charge range 10 to 15 the four distributions exhibit a kind of “plateau”. In this range the fragment production rates are roughly equivalent whatever the incident energy is. Finally the charge distributions evolve from a broad shape at lower incident energy, where residues up to the size of the projectile are observed, toward an almost exponential form at 50 AMeV, favouring the production of lighter fragments. Moreover figure 2, where the distributions

of the heaviest fragment in the event are shown, confirms this behaviour. Here again the distribution at 32 AMeV is very broad, its average value is $`<Z_{max}>=25`$, it decreases to smaller $`<Z_{max}>=15`$ at 50 AMeV. It is important to notice that, even with this strong evolution in the charge distribution, the mean fragment multiplicity does not change too much with the incident energy. It evolves from 5 to 7 fragments with $`Z3`$ only.

Concerning the kinematic characteristics of the fragments, figure 3 shows an example of the fragment angle-integrated centre of mass kinetic

energy spectra for Li, O and P nuclei produced in central collisions of Xe + Sn at 32 and 50 AMeV. The distributions are broad; they are broader for the heavier elements. Comparing the spectra obtained at 32 and 50 AMeV, we observe easily that their shape, in particular the slopes of their exponential tails, are different. The distributions are broader and harder at 50 than at 32 AMeV.

We finally present in figure 4, for the four incident energies, the mean centre of mass kinetic energy of the fragments as a function of their atomic number. It increases with the charge $`Z`$ and then saturates beyond $`Z=15`$. It also increases with the bombarding energy but very little. We wondered whether this observation is true for central collisions in general, or is rather dependent on our selection. In fact it is the mean kinetic energy of the heaviest fragment which “saturates” while that of the other fragments increases monotonously with $`Z`$. The $`\theta _{flow}`$ selection we use is derived from fragment kinetic properties and therefore its effect on observables such as e.g. fragment energies and angular distributions must be taken into account in events selected in this way. Nevertheless this selection has little influence on the study of individual fragment characteristics such as excitation energy and secondary decay, whatever the mechanism of their formation.

## III Extraction of secondary evaporated light charged particles

The main aim of this work is to extract the intrinsic properties of the fragments independently of the mechanism responsible for their formation. Are they excited ? If so what are the associated LCP evaporated from the parents ? Reconstructing the primary fragments assumes that we are able experimentally to isolate the secondary contribution. This is possible if the fragments formed are not too excited, so that the time scale associated with their decay is much greater than the time scale of their production. The origin of the fragments is still an open question but is not the subject of this paper.

### III.1 Correlation functions

In the previous section it was shown that on average about 6 fragments are produced in central collisions of Xe + Sn at different energies. However the production of LCP is much more important, on average their number reaches 28 particles for the 50 AMeV beam energy. There are at least three different stages to produce these particles : i) in the early stage of the collision, in this case we call them primary particles; ii) at the same time as the formation of the fragments; iii) they can be emitted from the excited primary fragments, we call those the secondary particles. Correlation functions are a powerful tool for extracting small signals. This is the method we used to extract, on the average, the LCP emitted from each fragment. With the help of simulations we have developed a correlation technique to extract possible signals Marie et al. (1998); Hudan et al. (2000); Hudan (2001).

Figure 5 shows the relative velocity distributions : i) for P-$`\alpha `$ pairs taken from the same events, ii) for the uncorrelated events obtained by taking the fragment from a given event and the light particle from another event, iii) the correlation function defined as the ratio of the correlated and uncorrelated relative velocity distributions, iv) the difference correlation function defined as the difference between the correlated and uncorrelated distributions. In this work to decorrelate the relative velocity between the fragment and the LCP pairs we used the event mixing procedure Drijard et al. (1984). In this example, for each phosphorus found in an event having a number of alphas $`N_a`$ we take randomly $`N_a`$ alphas emitted in $`N_a`$ other events.

This technique is different from the one reported in ref. Marie et al. (1998) where Li nuclei were used to decorrelate the events. The problem with such a technique is that the Li can be the product of the known resonance of $`{}_{}{}^{7}Be`$ which decays to $`{}_{}{}^{6}Li`$ \+ $`p`$ and increases the background, thus decreasing the yield of true correlated protons. However the final result is almost the same (within the error bars) as the old method of decorrelation of events based on Li.

As we can see, the example presented in figure 5, exhibits a bump around 2.5 cm/ns relative velocity in the correlation function and difference function which may be related to the evaporation of an $`\alpha `$ particle from a parent of phosphorus. The behaviour of this correlation encourages us to make such an analysis. However it is necessary to simulate the background in order to extract the signal.

### III.2 Simulation of the background shape

The objective of this simulation is not to reproduce the data, it is more to have an idea about the shape of the background. We used a modified version of the SIMON event generator Durand et al. (1992) to simulate a scenario deduced from BNV Bonasera et al. (1994) calculations. Two steps are assumed in these simulations. The first step is the cooling of the initial fused system through a sequential light particles (LP) emission process (primary LP), the second is the fragmentation of the smaller remaining source where the remaining excitation energy is shared between a fixed number of primary fragments (typically 6 to 7 fragments). Then the primary fragments decay sequentially while moving apart under Coulomb forces plus an initial radial velocity. This simulation reproduces reasonably well the global experimental features. In particular the kinematic observables are well reproduced (see for example ref. Marie et al. (1997).)

The calculated relative velocities are shown in Figure 6.a (thick lines) for Mg-p pairs and for input parameters which reproduce data for the 32 AMeV Xe + Sn central collisions. Since in this version of SIMON we know which particle is emitted from which fragment, we plotted in the same figure the different contributions : the primary contribution (dotted histogram) that we call contribution 1, the evaporated protons from all other fragments except the parents of magnesium (dashed histogram) that we call contribution 2 and finally the protons emitted from the parents of detected magnesium fragments (hatched-dashed histogram) called contribution 3. As expected, the latter contribution is very small, it represents the protons truly correlated to the magnesium nucleus that we must extract from the data. Figure 6.b shows the uncorrelated relative velocity for Mg-p pairs reconstructed by mixing the calculated events. Figure 6.c and 6.d show the Mg-p correlation function (the ratio of the correlated and uncorrelated relative velocity distributions of fig.6.a and fig.6.b) and the difference function (the difference of the latter distributions), respectively. In the same figures are plotted the associated true backgrounds (dashed histograms) calculated by dividing (subtracting) the sum of contributions 1 and 2 by the uncorrelated distribution (of fig.6.b). The hatched areas represent the contribution of secondary emission from the parents of magnesium (contribution 3). The shape of the background shown in figure 6.c is well fitted by the function :

$$R(V_{rel})=A\frac{1}{BV_{rel}+C}$$
(1)
where A, B and C are parameters which differ for each fragment-LCP pair. In fact only 3 coordinates are needed to solve this equation, we then used particular points from figure 6.d to do so. The first one corresponds to the first point at which the difference function is equal to zero (at small relative velocity). The second point used is the local minimum seen at small relative velocity (around 2.5 cm/ns) in the difference function (fig.6.d) which corresponds to the minimum relative velocity allowed by the Coulomb barrier. The third one corresponds to the first point where the difference function is equal to zero just after the second minimum, in this region the secondary evaporation vanishes.

In order to validate the method employed to estimate the background several tests have been made. We summarize the two most important tests that we already reported in ref. Marie et al. (1997):

a) We compared the number of protons deduced by subtracting from the difference function (6.d) the real background and the background evaluated by the parametrisation of (Eq.1). We recover 91% of the evaporated protons from Mg and 84% of evaporated protons from all prefragments.

b) The second check is related to the possible upper limit of the method. We performed SIMON simulations assuming higher excitation energies in the primary fragments. For 7.5 AMeV excitation energies we recovered 81% of evaporated protons. This result indicates that the fraction of all evaporated protons recovered by this method is rather insensitive to the excitation energies of the primary fragments.

Because the experimental shape of the correlation function as well as the difference function (Fig.5) have the same behaviour as those in our simulation, we applied the same method to the experimental data to remove the background. From this simulation and method developed above we are able to isolate the LCP evaporated by the primary fragments.

### III.3 Application to the data

Figure 7 shows the experimental correlation function, the difference function and the velocity distribution of $`\alpha `$ correlated to phosphorus fragments for the central collisions of Xe + Sn at 32 AMeV. In the same figure are plotted the corresponding background calculated with Eq.1 by using three points taken from the experimental distributions as described in the above section. Therefore the $`\alpha `$ velocity spectrum is deduced by subtracting the background (the curve in fig.7 upper right panel ) from the difference function. This contribution represents the spectrum of $`\alpha `$ particles emitted by the parent of $`P`$ fragment. From the mean value of the distribution we can deduce the average kinetic energy of alpha. Its integral normalised to the total number of phosphorus nuclei provides the average multiplicity of $`\alpha `$ particles evaporated from parents of P fragments.

The uncertainties of the extracted quantities are mainly related to the uncertainty of taking the three points which define the background. In practice the first minimum in the difference function is easy to locate: the corresponding error is small (see fig.7.b). The two other points are more difficult to extract, with the possibility of significant uncertainties. We then decided to take intervals around each point which are divided into a number of bins. Considering all possible combinations of one bin in the first interval and another in the second leads to a distribution of multiplicities. This distribution has a narrow gaussian shape. We then consider the mean value of this distribution as the average multiplicity and its half-width as the error due to the method.

## IV Experimental results

### IV.1 Average multiplicities and kinetic energies of the LCP correlated to the fragments

We applied the method described above for all fragment-LCP pairs made by combining LCP isotopes ($`p,d,t,^3He`$ and $`\alpha `$) and a range of fragments emitted in central collisions between Xe and Sn at four incident energies, 32, 39, 45 and 50 AMeV. However due to a small cross section for heavy fragment production which implies a low statistics (see figs. 1 and 2) we performed these analyses for a limited range of fragment charges depending on the beam energy. Thus the maximum fragment charge we studied at 32 AMeV was 30, 27 at 39 AMeV, 22 at 45 AMeV and 20 for 50 AMeV.

The extracted average LCP multiplicities and their average kinetic energy are given in figures 8 and 9 as a function of the charge, $`Z_{IMF}`$,

of the detected fragments and for the four bombarding energies. The average multiplicities increase with the fragment size. The multiplicities are low and do not exceed a value of 1.5 which implies that the excitation energy of the corresponding primary fragments is moderate. For a given light charged particle, the multiplicity seems not to change with the beam energy. From the spectra of the LCP evaporated from the parents of the detected fragments we can extract the mean kinetic energy. This is shown in figure 9. It increases slightly with the charge of the fragment for the four incident energies and in particular for proton and alpha particles. Notice that the kinetic energies of $`{}_{}{}^{3}He`$ are high compared to the values of the other particles. The observed effect may be due to the higher identification threshold energy for $`{}_{}{}^{3}He`$.

### IV.2 Reconstruction of the size and excitation energy of the primary fragments

To reconstruct the charge of the primary fragments we used the LCP multiplicities correlated to each fragment as described in the last paragraph. Therefore the average charge of the primary fragment, $`<Z_{pr}>`$, is given by the sum of the detected fragment and all evaporated LCP’s charge weighted by their corresponding multiplicities. $`<Z_{pr}>`$ is then given by the relationship :

$$<Z_{pr}>=Z_{IMF}+z_i<M_i>$$
(2)
where $`Z_{IMF}`$ is the detected fragment charge, $`z_i`$ and $`<M_i>`$ are the charge and the average multiplicity of the evaporated particle $`i=p,d,t,^3He`$ and $`\alpha `$.

In order to reconstruct the mass of the primary fragments, a quantity needed to deduce the excitation energy, we made two extreme assumptions : the first one is that the primary fragments are produced in the valley of stability, the second assumes that they are produced with the same N/Z ratio as the composite initial system (N/Z conservation assumption). However, as mentioned above the INDRA detector does not resolve the fragment isotopes, we therefore made an additional assumption which supposes that the $`Z`$-identified detected fragments have a mass corresponding to their valley of stability isotope. In the framework of these assumptions we deduce from the primary fragment masses the number of neutrons evaporated from the primary fragments.

Figures 10 and 11 show the result of this reconstruction for the four incident energies. The values of the primary charge (fig.10, upper panel) obtained vary from 1 to 5 charge units larger than the detected fragment. The mass of the primary fragment depends on the assumption (fig.10, down panel). The average neutron multiplicities are deduced from the mass conservation, knowing the mass of the primary fragment, the detected fragment and the mass of the secondary light charged particle contribution. Figure 11 shows for the two assumptions the evolution of the number of neutrons for the four systems as a function of the deduced primary fragment atomic number. Whatever the beam energy, the multiplicity of neutrons reaches quite high values, up to 7 neutrons for the N/Z ratio conservation assumption. This is due to our assumption that detected fragments have their valley of stability mass. Clearly, when we also assume that the primary fragments are produced in the valley of stability, the deduced neutron multiplicity cannot be very high. Conversely, imposing an N/Z of 1.39 for nuclei with Z = 3-30 means that primary fragments have large neutron excess compared to the (valley of stability) detected fragments.

At this stage, the calorimetric procedure can be applied to reconstruct the average excitation energy of the primary fragments ($`<E_{pr}^{}>`$). It is given by the relationship :

$$<E_{pr}^{}>=<M_{LCP}><E_{LCP}>+<M_n><E_n>Q$$
(3)
where $`<E_{LCP}>`$ and $`<E_n>`$ are the average kinetic energies in the frame of the source (fragment) of the measured evaporated LCP’s and the deduced neutrons with the average multiplicity $`<M_n>`$. The neutron kinetic energy $`<E_n>`$ is taken as the proton kinetic energy minus the proton Coulomb barrier. $`Q`$ is the mass balance of the reaction.

Figure 12 shows the result of this procedure for the two scenarios and at the four bombarding energies. As expected from the deduced multiplicities (see paragraph IV.1), the excitation energy increases with the size of the primary fragment for all bombarding energies and for the two assumptions. However, for the 32 AMeV system, $`<E_{pr}^{}>`$ seems to saturate for the heavier fragments. We could wonder if this is due to limitations of the method. However, as we already mentioned in paragraph III.2, simulations have been performed at much higher excitation energy into the primary fragments showing that we recover more than 80% of the evaporated protons.

To decide between the scenarios for primary fragment mass, valley of stability or with the $`N/Z`$ conservation assumption, extensive statistical calculations have been performed using the GEMINI Charity et al. (1988) code, for the 50 AMeV system. In these calculations the input to the code was the experimental deduced primary charge, the fragment masses with the two assumptions and the associated excitation energies. The comparison to the experimental LCP multiplicities and kinetic energies suggests that the $`N/Z`$ conservation assumption is the most reasonable scenario. Details of these calculations are given in ref. Marie et al. (1998).

The linear trend of the $`<E_{pr}^{}>`$ with the primary charge indicates that the average excitation energy per nucleon, $`<e_{pr}^{}>`$ in MeV/nucleon, is constant whatever the size of the primary fragment. In figure 13 we verified the latter characteristic by plotting this variable. The horizontal lines in this figure represent the average over the whole set of primary fragments. Besides a few small charges the data points, within the error bars, lie on this straight horizontal line. Figure 14 shows the evolution of this value as a function of the bombarding energy. The vertical bars are the standard deviations from the mean values. They are small and do not exceed 1 AMeV, which supports the constancy of the value of $`<e_{pr}^{}>`$. For the N/Z conservation assumption the excitation energy per nucleon increases from 2.2 AMeV at 32 AMeV and saturates at 3 AMeV beyond 39 AMeV. For the valley of stability assumption, $`<e_{pr}^{}>`$ saturates also but at a lower value.

The constancy of the fragment excitation energy per nucleon, $`<e_{pr}^{}>`$ for different fragment masses, seen in figure 13 has been interpreted in Marie et al. (1998) as meaning that, on the average, thermodynamical equilibrium was achieved at the disassembly stage of the system. Only one bombarding energy (50 AMeV) was available in the previous work. On the other hand the saturation of $`<e_{pr}^{}>`$ beyond 39 AMeV beam energy (fig.14) may indicate that the fragments reach their excitation energy limit (or limiting temperature) Levit and Bonche (1985); Koonin and Randrup (1987).

### IV.3 Proportion of the evaporated light charged particles

In paragraph IV.1 we have extracted the average multiplicities of the secondary evaporated light charged particles for a given fragment. It is interesting to use this information in order to study the characteristics of the multifragmentation events. Indeed the LCP multiplicity per event can be another pertinent observable. Table 1 shows the secondary LCP multiplicities per event, the total LCP measured per event and the ratio of both quantities, for the four beam energies. The secondary LCP multiplicities per event are defined as the sum of the secondary evaporated LCP’s per fragment, extracted by the method described above, weighted by the measured fragment multiplicity per event, $`M_{IMF}`$. These values are plotted in figure 15 as a function of the beam energy.

The fraction of helium isotopes evaporated in the decay of the primary fragments is higher than for those of hydrogen. This difference is more pronounced

at lower beam energy. We observe also that the maximum proportion of evaporated particles does not exceed on the average 35% of the total number of produced light charged particles. The proportion of secondary particles increases between 32 and 39 AMeV, which reflects the increasing of the excitation energy of the primary fragments as it is seen in figure 14. Then this fraction decreases for higher incident energies, it reaches 23% at 50 AMeV, while $`<e_{pr}^{}>`$ saturates.

It should be noticed that the proportion of the secondary evaporated particles given is a lower limit, because we did not consider the contribution that can originate from the decay of unstable nuclei Nayak et al. (1992) such as $`{}_{}{}^{8}Be`$, $`{}_{}{}^{5}Li`$ etc. and the decay of short-lived excited sates. We finally have to stress that the results we obtained with the method described above are given on the average.

## V Comparison with a statistical model

An application of the experimental estimation of this secondary statistical component is to constrain the statistical multifragmentation models Bondorf et al. (1995); Botvina et al. (1987); Gross (1990); Raduta and Raduta (1997a, b); Parvan et al. (1999); Pratt and Das Gupta (2000). The comparison of the extracted quantities with these models provides a crucial test of some of their basic assumptions. Since in the MMMC Gross (1990) model the primary fragments undergo instantaneous decay with neutron emission only, it can not be used for comparison with the data.

We have chosen to compare our data, with more details, to the SMM model using input source parameters very close to the ones already optimised in previous works Bougault et al. (1997); Salou (1997); Leneindre (1999). As shown there SMM provides a very good description of experimental fragment partitions. In the present paper we aim to analyse the general behaviour of excitation energy of primary fragments, therefore, for simplicity, the size of the initial source has been fixed to be Z=83 and A=198 for the four incident energies. This corresponds to N/Z = 1.39 which is the same ratio as the initial system. This choice is justified by some dynamical calculations of source parameters in this energy range Botvina et al. (1995); Tan et al. (2001). Although the $`N/Z`$ ratio of the SMM primary fragments increases slightly with increasing $`Z`$ of the fragments, it remains very close to the $`N/Z`$ ratio of the initial source Botvina and Mishustin (2001). Therefore we will compare the results of these calculations to the extracted experimental results using the N/Z conservation hypothesis. The freeze-out volume has been fixed to three times the normal volume. Finally for each incident energy we used the excitation energy of the initial source as a free parameter. The thermal excitation energy values which reproduce best the charge distributions of the detected fragments are given in table 2.

### V.1 Characteristics of the primary fragments

We have used a version of SMM where we have access to the freeze-out configuration, i.e. to primary fragments’ characteristics before secondary decay and Coulomb propagation. This standard version is described in Bondorf et al. (1995).

The results of the SMM calculation, extracted directly from the freeze-out volume, are compared to the data in figure 16. The excitation energy of the primary fragments are globally well reproduced for the four incident energies.

Small deviations are, however, observed for large primary fragment charges in particular for the 32 AMeV case. The experimental saturation of the excitation energy is not reproduced. Quantitative comparisons with the experimental excitation energy per nucleon averaged over the charge range of the measured fragments, are presented in table 2.

The values of the calculated $`<e_{pr}^{}>`$ show smooth increase with the beam energy while the data seems to saturate at 3 AMeV above 39 A.MeV.

### V.2 Evaporated light charged particles

The contribution of the secondary evaporated LCP reflects the excitation energy of the primary fragments discussed in the previous paragraph. How do the small differences between the data and the calculation for the excitation energy affect the predicted LCP multiplicities ? We compare in figure 17 the charge contribution of total evaporated LCP resulting from SMM to the data.

The values extracted from the calculation are of the same order of magnitude as the experimental ones. However the trend is not reproduced, SMM overestimates the evaporative contribution.

At 32 AMeV, the discrepancy can be due to the limited charge range considered experimentally. Indeed, due to low statistics as we already mentioned, we do not take into account evaporation from heavier fragments, which are however more excited than lighter ones. To check this point we extracted the total evaporated particles by using the correlation functions of reduced velocities instead of the relative velocities ( $`V_{red}=V_{rel}/\sqrt{Z1+Z2}`$). This variable has the advantage to eliminate the charge dependence of the fragment-LCP relative velocity correlation functions Kim et al. (1992). By doing this procedure, taking all fragments into account, the total charge of evaporated particles increases significantly to be in agreement with the calculated value for the 32 AMeV case.

For the 50 AMeV case, where the limitation in the charge range is more important than for the other beam energies, the evaporated contribution changes very little and fails to increase the value of total evaporated charge from $`Z_{LCP}`$ = 9 to $`Z_{LCP}`$ = 14 predicted by SMM (fig. 17).

The discrepancy is real, though, partly, it is caused by the thermal source size, which should decrease with the beam energy. The SMM calculations do predict the decrease of the number of evaporated LCP, because of decreasing size and number of IMF at very large excitation energies, in the “falling” part of the “rise and fall” of multifragmentation. However, in the experiment this effect is observed when the maximum of multifragmentation is not yet reached. In the calculations this behaviour takes place because the number of evaporated alpha particles increases, contrary to the experimental result. This could be a consequence of the secondary deexcitation prescription employed in SMM Botvina et al. (1987). An other possible reason would be an overexcitation of light primary IMF’s predicted by SMM. The decay of these IMF’s contributes considerably into LCP production and their share increases with the thermal energy.

The decrease of the experimental evaporated component $`Z_{LCP}`$ at high energy could be alternatively understood if we consider the increasing effect of the collision dynamics. The direct emission of LCP increases with the incident energy while the proportion of the thermal contribution decreases. This could be mocked up in the SMM calculations by decreasing the thermal source size, but can in no case be predicted by SMM.

It is worth noting the contribution of light charged particles produced at freeze-out as predicted by SMM. Figure 17 shows that this contribution increases with the beam energy more rapidly than that of the evaporated particles.

## VI Discussion of results

In this work we have directly measured the saturation of the thermal excitation energy deposited in fragments produced in central heavy-ion collisions between 32 and 50 AMeV, by associating with each detected cold fragment the light charged particles evaporated by the primary excited parent nucleus. This saturation at excitation energies of around 3 AMeV observed in section III (see figure 14) is accompanied by a saturation of the number of evaporated LCP, that leads to a decrease in the proportion of evaporated to all detected LCP, with increasing incident energy (see figure 15 or figure 17).

A similar saturation has been observed in an earlier work by Jiang et al.Jiang et al. (1989) using a completely different experimental method, based on the measurement of neutron multiplicities. The authors claimed the saturation of the thermal energy deposited in hot nuclei formed in collisions of Ar + Au and Ar + Th in the energy range 27-77 AMeV. Their claim was based on the observation of a saturation of the multiplicity of evaporated neutrons, as well as that of the light charged particles detected in coincidence at backward angles, in central collisions at increasing beam energies. The neutron multiplicity saturates for the system Ar + Th around 30 AMeV at $`<M_n>=35`$. Let us note in passing that we estimate the neutron multiplicity per event evaporated by the system Xe + Sn to be $`<M_n>=23`$ at 39 AMeV.

In Jiang et al. (1989) the authors concluded that the observed saturation was due to the increasing inefficiency of the reaction mechanism to deposit thermal energy in to the hot nuclei it produced, rather than it being related to reaching the limits of excitation energy or temperature that a nucleus may support. In discussing the results of the present work we must ask ourselves the same question, but the situation is complicated by the fact that here we are dealing with several heated nuclei per event which may themselves result from the break-up of some other heavy, hot system. Here we will present some elements which may help to find an answer.

Although the excitation energies of primary fragments remain constant for incident energies $``$39 AMeV, the detected fragment partitions continue to evolve, becoming steeper with increasing bombarding energy (figure 1) while the average charge of the largest fragment varies from 25 at 32 AMeV to 15 at 50 AMeV (figure 2). This suggests that with increasing energy, above 39AMeV, the average number and temperature of primary fragments produced in the reactions does not change, whereas their average size decreases. Moreover, in order to conserve the total mass of the system the number of light particles produced prior to secondary evaporation from fragments must also increase with increasing energy. If indeed the mechanism for thermal energy deposition saturates, then the energy not used in forming and heating fragments has to be evacuated by some other means, for example direct particle production. This would lead to such an increase in non-evaporated particle multiplicity. Some energy may also be locked up as kinetic energy of fragments due to some kind of collective motion, either isotropic (compression-expansion effects) or anisotropic (incomplete stopping).

It should be recalled that the decrease of fragment excited state lifetimes with the excitation energy can limit the mean excitation values obtained in this paper. However, simulations we have performed indicate that the effect of shorter lifetimes on the efficiency of the method is quite small for primary fragment excitation energies up to 7.5 AMeV.

It is interesting to compare our results with a recent compilation of limiting temperatures extracted from different experimental measurements Natowitz et al. (2002). It suggests that $`T_{lim}`$ decreases with increasing nuclear mass, in good agreement with calculations Zhang et al. (1996, 1999). The primary fragments considered in the present work (figure 10)have, at the very most, masses $`A=80`$, while most of them have masses in the region $`A=10`$–50. The corresponding limiting temperature from Natowitz et al. (2002) is $`T_{lim}=9`$ MeV or $`E^{}/A=7.5`$ AMeV. As these values are much higher than the 3 AMeV maximum excitation energy we find in our primary fragments, this would imply that the observed saturation is due to reaction mechanism and not related to $`T_{lim}`$.

However, in the same compilation limiting excitation energies $`3`$ AMeV are found for the heaviest nuclei with masses in the $`A=150`$–200 range. If we suppose that fragments are produced by the break-up of some heavy composite system formed in the reaction (as in SMM calculations) then it is possible that the observed saturation of primary fragment excitation energies is due to the saturation of the excitation energy of the initial system, which attains its (mass-dependent) $`T_{lim}`$.

## VII Conclusion

We have presented in this paper the experimental results of the intrinsic properties of the fragments produced in the central collisions of Xe + Sn from 32 to 50 AMeV bombarding energy. Quantitative experimental determination of the size and excitation energy of the primary fragments produced at such collisions before their decay are given for the four beam energies. The comparison of these extracted quantities with models provides a crucial test of some of their basic assumptions.

The experimental methods used in this work are based on the relative velocity correlation functions between the detected fragment and light charged particles. Thus we have extracted the average multiplicity of the evaporated particles and their average kinetic energies in the centre of mass of the fragments. These two variables have been used in order to reconstruct the average charge, mass and excitation energy of the primary fragments.

Our results show that for a given beam energy, the excitation energy per nucleon is almost constant over the whole studied range of fragment charge. The statistical multifragmentation model, SMM, reproduces very well the internal excitation energy of the primary fragments. The average value of this quantity increase from 2.3 AMeV for a beam energy 32 AMeV to saturate around 3 AMeV for 39 AMeV and above.

We also deduced the proportion of evaporated light charged particles per event, amounting to 30% of the total measured LCP for the 32 AMeV reaction, increasing to 35% at 39 AMeV and decreasing down to 23% for 50 AMeV. Therefore the majority of light charged particles are not evaporated by excited primary fragments in these reactions. Neither the absolute values of this proportion nor its evolution are reproduced by SMM calculations assuming a constant size for the multifragmenting system.

The two last results may indicate either i) that the system which disassembles into fragments is not able to sustain more than 3 AMeV thermal excitation energy and the excitation energy of the fragments reflects the temperature limit of that system, or/and ii) that the mechanism of dissipation of beam energy in to thermal energy of hot nuclei saturates above $`32`$AMeV, and the kinetic energy in excess is evacuated via direct particle production and collective motion of the fragments.

###### Acknowledgements.

We thank A. Botvina for fruitful discussions and for performing the SMM calculations. We thank also the staff of the GANIL Accelerator facility for their support during the experiment. This work was supported by Le Commissariat à l’Energie Atomique, Le Centre National de la Recherche Scientifique, Le Ministère de l’Education Nationale, and le Conseil Regional de Basse Normandie.


format: text

text:
# Shear viscosity for a heated granular binary mixture at low-density

## I Introduction

The macroscopic behavior of rapid granular flows can be described through hydrodynamic equations accounting for dissipation among the interacting particles. A basis for the derivation of the hydrodynamic equations and explicit expressions for the transport coefficients appearing in them is provided by the corresponding Boltzmann kinetic theory in the low-density regime. In the simplest model the grains are taken to be smooth hard spheres with inelastic collisions. Assuming the existence of a normal solution for sufficiently long space and time scales, the Chapman-Enskog method can be applied to get the velocity distribution function in terms of the hydrodynamic fields and their spatial gradients. This method must be conveniently adapted to inelastic collisions due to the new time dependence of temperature resulting from collisional energy loss. In the case of a monocomponent gas, the Navier-Stokes transport coefficients have been obtained in terms of the restitution coefficient without limitation on the degree of inelasticity. This analysis for a monocomponent system has been also extended to dense gases in the context of the Enskog equation.

Similar studies for multicomponent granular gases are more scarce and most of them limited to the asymptotically weak dissipation limit. Although these studies permit in principle different temperatures for each species, they usually assume energy equipartition and so the partial temperatures $`T_i`$ are made equal to the mixture temperature $`T`$. However, some recent results obtained in molecular dynamics simulations  as well as in real experiments of vibrated mixtures in three  and two  dimensions clearly show the breakdown of energy equipartition. A more recent kinetic theory calculation which takes into account temperature differences has been carried out by Garzó and Dufty . They solved the set of Boltzmann coupled equations of the binary mixture by means of the Chapman-Enskog expansion for states near the local homogeneous cooling state. The mass, heat, and momentum fluxes were determined to first order in the gradients of the hydrodynamic fields and the associated transport coefficients were explicitly identified. As in the case of elastic collisions, these transport coefficients verify a set of coupled linear integral equations which are solved approximately by using the leading terms in a Sonine polynomial expansion. The results derived by Garzó and Dufty  provide a description of hydrodynamics in binary granular mixtures valid a priori over the broadest parameter range and not limited to the quasielastic regime. In particular, the consequences of the temperature differences on the transport coefficients were shown to be quite significant.

In the case of molecular fluid mixtures, it is known that the leading order truncation is quite accurate, except for extreme mass ratios (e.g., electron–proton systems). Much less is known in the case of inelastic collisions, although some comparisons with computer simulations for homogeneous states indicate that the accuracy is similar to that for elastic collisions. The objective here is to compare the kinetic theory predictions for the shear viscosity with those obtained from a numerical solution of the Boltzmann equation by means of the Direct Simulation Monte Carlo (DSMC) method. Specifically, the simulations are performed for a granular mixture undergoing uniform shear flow (USF), namely, a macroscopic state characterized by constant partial densities $`n_i`$ ($`i=1,2`$), uniform temperature $`T`$ and a linear flow velocity profile $`u_{i,x}=u_x=ay`$, $`a`$ being the constant shear rate. In a molecular fluid under USF, the temperature increases in time due to viscous heating. As a consequence, the average collision frequency $`\nu `$ (which is proportional to $`T^{1/2}`$ for hard spheres) increases with time and the reduced shear rate $`a^{}=a/\nu `$ (which is the relevant nonequilibrium parameter of the problem) tends to zero in the long time limit. This implies that for sufficiently long times the system reaches a regime described by linear hydrodynamics and the Navier-Stokes shear viscosity can be identified . For granular fluids, the inelasticity of collisions introduces an energy sink in the balance equation for the temperature. Thus, the relationship between the temperature and the shear viscosity is not as simple as for molecular fluids since there is a competition between viscous heating and collisional cooling. However, if the granular fluid is externally excited by an external energy source that exactly compensates for the collisional energy loss, the viscous heating effect is still able to heat the system (as in the elastic case) and one can identify the shear viscosity in the limit of small shear rate relative to the collision frequency $`\nu `$ (i.e., $`a^{}0`$). Although there are several choices for the external driving force, here we consider an external thermostat proportional to the peculiar velocity. This thermostat has been frequently used in nonequilibrium molecular dynamics simulations of molecular fluids and has the advantage that, in the absence of shear, it does not affect the dynamics of the system at all since it is formally equivalent to a rescaling of the velocities .

The motivation of our study is twofold. First, in light of some doubts about the validity of a hydrodynamic description for granular flow, the comparison with simulation allows us to test the Chapman-Enskog solution obtained by assuming the existence of a normal or hydrodynamic regime. Since the parameter space here is quite large the tests of the theory and concepts are quite stringent. Second, as said above, we can also assess the degree of reliability of the approximate solution (first Sonine polynomial approximation) to the resulting integral equation over a wide range of the parameter space. With respect to the driving external force used in our analysis, we do not claim that it is the most suited one to model any real experiment. However, it has the advantage that it can be incorporated into kinetic theory and computer simulations very easily and it allows to check the assumptions of the Chapman-Enskog method.

The plan of the paper is as follows. In Section II, we review the Boltzmann equation and associated macroscopic conservation laws in the presence of the Gaussian thermostat. The Chapman-Enskog method is applied in Section III to get all the transport coefficients of the mixture, with special emphasis in the shear viscosity coefficient. The details of the derivation are displayed in Appendix A. Section IV deals with the application of the DSMC method of the Boltzmann equation to USF with thermostat. The Chapman-Enskog and simulation results are compared in Section V at the level of the shear viscosity showing a good agreement. We close the paper in Section VI with a discussion of the results presented.

## II The Boltzmann equation and transport coefficients

We consider a binary mixture of smooth hard spheres of masses $`m_1`$ and $`m_2`$, and diameters $`\sigma _1`$ and $`\sigma _2`$. The inelasticity of collisions among all pairs is characterized by three independent constant coefficients of normal restitution $`\alpha _{11}`$, $`\alpha _{22}`$, and $`\alpha _{12}=\alpha _{21}`$, where $`\alpha _{ij}`$ is the restitution coefficient for collisions between particles of species $`i`$ and $`j`$. Due to the intrinsic dissipative character of collisions, an energy supply is requested to fluidize a granular gas. For simplicity, here the fluidization is driven by the action of a non-conservative external force, frequently referred to as the Gaussian thermostat. In this case, the mixture is heated by an “anti-drag” force, linear in the peculiar velocity $`𝐕`$ and chosen to exactly compensate for collisional cooling. As said in the Introduction, this deterministic thermostat has been widely used in computer simulations of molecular fluids . Under these conditions and in the low-density regime, the distribution functions $`f_i(𝐫,𝐯;t)`$ $`(i=1,2)`$ for the two species are determined from the set of nonlinear Boltzmann equations

$$\left(_t+𝐯_1\right)f_i+\frac{1}{2}\xi \frac{}{𝐯_1}\left(𝐕_1f_i\right)=\underset{j}{}J_{ij}\left[𝐯_1|f_i(t),f_j(t)\right],$$
(1)
where the constant $`\xi `$ is taken to be the same for each species . Here, $`𝐕_1𝐯_1𝐮`$, $`𝐮`$ being the flow velocity. The Boltzmann collision operator $`J_{ij}\left[𝐯_1|f_i,f_j\right]`$ describing the scattering of pairs of particles is

$`J_{ij}\left[𝐯_1|f_i,f_j\right]`$ $`=`$ $`\sigma _{ij}^2{\displaystyle 𝑑𝐯_2𝑑\widehat{𝝈}\mathrm{\Theta }(\widehat{𝝈}𝐠_{12})(\widehat{𝝈}𝐠_{12})}`$ (3)
$`\times \left[\alpha _{ij}^2f_i(𝐫,𝐯_1^{},t)f_j(𝐫,𝐯_2^{},t)f_i(𝐫,𝐯_1,t)f_j(𝐫,𝐯_2,t)\right],`$
where $`\sigma _{ij}=\left(\sigma _i+\sigma _j\right)/2`$, $`\widehat{𝝈}`$ is a unit vector along their line of centers, $`\mathrm{\Theta }`$ is the Heaviside step function, and $`𝐠_{12}=𝐯_1𝐯_2`$. The primes on the velocities denote the initial values $`\{𝐯_1^{},𝐯_2^{}\}`$ that lead to $`\{𝐯_1,𝐯_2\}`$ following a binary collision:

$$𝐯_1^{}=𝐯_1\mu _{ji}\left(1+\alpha _{ij}^1\right)(\widehat{𝝈}𝐠_{12})\widehat{𝝈},𝐯_2^{}=𝐯_2+\mu _{ij}\left(1+\alpha _{ij}^1\right)(\widehat{𝝈}𝐠_{12})\widehat{𝝈},$$
(4)
where $`\mu _{ij}=m_i/\left(m_i+m_j\right)`$.

The relevant hydrodynamic fields are the number densities $`n_i`$, the flow velocity $`𝐮`$, and the “granular” temperature $`T`$. They are defined in terms of moments of the distributions $`f_i`$ as

$$n_i=𝑑𝐯f_i(𝐯),\rho 𝐮=\underset{i}{}𝑑𝐯m_i𝐯f_i(𝐯),$$
(5)
$$nT=p=\underset{i}{}𝑑𝐯\frac{m_i}{3}V^2f_i(𝐯),$$
(6)
where $`n=n_1+n_2`$ is the total number density, $`\rho =m_1n_1+m_2n_2`$ is the total mass density, and $`p`$ is the hydrostatic pressure. At a kinetic level, it is convenient to introduce the kinetic temperatures $`T_i`$ for each species defined as

$$\frac{3}{2}n_iT_i=𝑑𝐯\frac{m_i}{2}V^2f_i.$$
(7)
The collision operators conserve the particle number of each species and the total momentum but the total energy is not conserved:

$$𝑑𝐯J_{ij}[𝐯|f_i,f_j]=0,\underset{i,j}{}𝑑𝐯m_i𝐯J_{ij}[𝐯|f_i,f_j]=0,$$
(8)
$$\underset{i,j}{}𝑑𝐯\frac{1}{2}m_i𝐯^2J_{ij}[𝐯|f_i,f_j]=\frac{3}{2}nT\zeta ,$$
(9)
where $`\zeta `$ is identified as the cooling rate due to inelastic collisions among all species. The macroscopic balance equations follow from the Boltzmann equation (1) and Eqs. (8) and (9). They are given by

$$D_tn_i+n_i𝐮+\frac{𝐣_i}{m_i}=0,$$
(10)
$$D_t𝐮+\rho ^1𝖯=0,$$
(11)
$$D_tT\frac{T}{n}\underset{i}{}\frac{𝐣_i}{m_i}+\frac{2}{3n}(𝐪+𝖯:𝐮)=(\zeta \xi )T.$$
(12)
In the above equations, $`D_t=_t+𝐮`$ is the material derivative,

$$𝐣_i=m_i𝑑𝐯𝐕f_i(𝐯)$$
(13)
is the mass flux for species $`i`$ relative to the local flow,

$$𝖯=\underset{i}{}𝑑𝐯m_i\mathrm{𝐕𝐕}f_i(𝐯)$$
(14)
is the total pressure tensor, and

$$𝐪=\underset{i}{}𝑑𝐯\frac{1}{2}m_iV^2𝐕f_i(𝐯)$$
(15)
is the total heat flux.

The energy balance equation (12) shows that the existence of a driving with the choice $`\xi =\zeta `$ compensates for the cooling effect due to the inelasticity of collisions. In that case, the macroscopic balance equations look like those of a conventional mixture with elastic collisions. Nevertheless, the transport coefficients entering in the constitutive equations are in general different from those of a gas of elastic particles. Furthermore, for systems with elastic collisions, the specific set of gradients contributing to each flux is restricted by fluid symmetry, time reversal invariance (Onsager relations), and the form of the entropy porduction. In the case of inelastic collisions only fluid symmetry applies and so there is more flexibility in representing the fluxes and identifying the corresponding transport coefficients. It follows from fluid symmetry that the pressure tensor has the same form to first order in the gradients as for the monocomponent gas. In the case of heat and mass fluxes, several different (but equivalent) choices of hydrodynamic fields can be used and some care is required in comparing transport coefficients in the different representations. Here, as done in the unforced case , we take the gradients of the mole fraction $`x_1=n_1/n`$, the pressure $`p`$, the temperature $`T`$ and the flow velocity $`𝐮`$ as the relevant ones. Thus, in this representation, the phenomenological constitutive relations for the fluxes in the low-density regime have the forms

$$𝐣_1=\frac{m_1m_2n}{\rho }Dx_1\frac{\rho }{p}D_pp\frac{\rho }{T}D^{}T,𝐣_2=𝐣_1,$$
(16)
$$𝐪=T^2D^{\prime \prime }x_1Lp\lambda T,$$
(17)
$$P_{\alpha \beta }=p\delta _{\alpha \beta }\eta \left(_\beta u_\alpha +_\alpha u_\beta \frac{2}{3}\delta _{\alpha \beta }𝐮\right).$$
(18)
The transport coefficients are the diffusion coefficient $`D`$, the thermal diffusion coefficient $`D^{}`$, the pressure diffusion coefficient $`D_p`$, the Dufour coefficient $`D^{\prime \prime }`$, the thermal conductivity $`\lambda `$, the pressure energy coefficient $`L`$, and the shear viscosity $`\eta `$. The Chapman-Enskog method  generalized to inelastic collisions allows one to get explicit expressions for these transport coefficients as functions of the restitution coefficients and the parameters of the mixture.

## III Shear viscosity of a heated granular mixture

The Chapman-Enskog method assumes the existence of a normal solution in which all space and time dependence of the distribution function occurs through a functional dependence on the hydrodynamic fields

$$f_i(𝐫,𝐯_1,t)=f_i[𝐯_1|x_1(𝐫,t),p(𝐫,t),T(𝐫,t),𝐮(𝐫,t)].$$
(19)
This functional dependence can be made local in space and time by means of an expansion in gradients of the fields. Thus, we write $`f_i`$ as a series expansion in a formal parameter $`ϵ`$ measuring the nonuniformity of the system,

$$f_i=f_i^{(0)}+ϵf_i^{(1)}+ϵ^2f_i^{(2)}+\mathrm{},$$
(20)
where each factor of $`ϵ`$ means an implicit gradient of a hydrodynamic field. The local reference state $`f_i^{(0)}`$ is chosen to give the same first moments as the exact distribution $`f_i`$. The time derivatives of the fields are also expanded as $`_t=_t^{(0)}+ϵ_t^{(1)}+\mathrm{}`$. The coefficients of the time derivative expansion are identified from the balance equations (10)–(12) after expanding the fluxes, and the cooling rate $`\zeta `$ in a similar series as (20). More details on the Chapman-Enskog method adapted to inelastic collisions can be found in Ref. . Now, the main difference with respect to the free cooling case  is that the sink term in the energy equation is zero (when one takes $`\xi =\zeta `$), so that the terms coming from the time derivative $`_t^{(0)}`$ vanish.

In the zeroth order, $`f_i^{(0)}`$ obeys the kinetic equation

$$\frac{1}{2}\zeta ^{(0)}\frac{}{𝐕}\left(𝐕f_i^{(0)}\right)=\underset{j}{}J_{ij}[f_i^{(0)},f_j^{(0)}],$$
(21)
where use has been made of the choice $`\xi ^{(0)}=\zeta ^{(0)}`$. Here, $`\zeta ^{(0)}`$ is determined by Eq. (9) to zeroth order. With this choice, Eq. (21) is identical to the the one obtained in the unforced case , and there is an exact correspondence between the homogeneous cooling state and this type of driven steady state. This is one of the advantages of the Gaussian thermostat. Dimensional analysis requires that $`f_i^{(0)}(𝐕)`$ must be of the form

$$f_i^{(0)}(V)=n_iv_0^3\mathrm{\Phi }_i(V/v_0)$$
(22)
where

$$v_0=\sqrt{2T\frac{m_1+m_2}{m_1m_2}}$$
(23)
is a thermal velocity defined in terms of the temperature $`T`$ of the mixture. So far, the exact form of $`\mathrm{\Phi }_i`$ has not been found, although a good approximation for thermal velocities can be obtained from an expansion in Sonine polynomials. In the leading order, $`\mathrm{\Phi }_i`$ is given by

$$\mathrm{\Phi }_i(V^{})\left(\frac{\theta _i}{\pi }\right)^{3/2}e^{\theta _iV^2}\left[1+\frac{c_i}{4}\left(\theta _iV^45\theta _iV^2+\frac{15}{4}\right)\right],$$
(24)
where $`V^{}=V/v_0`$, $`\theta _i=(\mu _{ji}\gamma _i)^1`$, and $`\gamma _i=T_i/T`$. The coefficients $`c_i`$ (which measure the deviation of $`\mathrm{\Phi }_i`$ from the reference Maxwellian) are determined consistenly from the Boltzmann equation. The approximation (24) provides detailed predictions for the cooling rate $`\zeta ^{(0)}`$, the temperature ratio $`T_1/T_2`$ and the cumulants $`c_i`$ as functions of the mass ratio, size ratio, composition and restitution coefficients. Recently, the accuracy of this approximate solution has been confirmed by Monte Carlo simulation of the Boltzmann equation over a wide range of the parameter space  .

The analysis to first order in $`ϵ`$ is similar to the one worked out in Ref.  for the free cooling case. Some details of the derivation of the transport coefficients as well as their final expressions are given in Appendix A. Here, given that the theoretical predictions for the shear viscosity coefficient $`\eta `$ will be compared with those obtained from Monte Carlo simulations, we focus our attention on the explicit final expression for $`\eta `$. According to Eq. (A25), the shear viscosity can be written as

$$\eta =\frac{nT}{\nu }\eta ^{},$$
(25)
where

$$\nu =\sqrt{\pi }n\sigma _{12}^2v_0,$$
(26)
is an effective collision frequency and the reduced shear viscosity coefficient $`\eta ^{}`$ is

$$\eta ^{}=x_1\gamma _1^2d_{1,1}^{}+x_2\gamma _2^2d_{2,1}^{}$$
(27)
with

$$d_{1,1}^{}=\frac{\gamma _2(\tau _{22}^{}\zeta ^{})\gamma _1\tau _{12}^{}}{\mathrm{\Delta }},d_{2,1}^{}=\frac{\gamma _1(\tau _{11}^{}\zeta ^{})\gamma _2\tau _{21}^{}}{\mathrm{\Delta }}.$$
(28)
Here, $`\zeta ^{}=\zeta ^{(0)}/\nu `$,

$$\mathrm{\Delta }=\gamma _1\gamma _2\left[\zeta ^2\zeta ^{}(\tau _{11}^{}+\tau _{22}^{})+\tau _{11}^{}\tau _{22}^{}\tau _{12}^{}\tau _{21}^{}\right],$$
(29)
and the dimensionless quantities $`\tau _{ij}^{}`$ are given by

$`\tau _{11}^{}`$ $`=`$ $`{\displaystyle \frac{16}{5\sqrt{2}}}x_1\left({\displaystyle \frac{\sigma _1}{\sigma _{12}}}\right)^2\theta _1^{1/2}\left[1{\displaystyle \frac{1}{4}}(1\alpha _{11})^2\right]\left(1{\displaystyle \frac{c_1}{64}}\right)`$ (33)
$`+{\displaystyle \frac{8}{15}}x_2\mu _{21}(1+\alpha _{12})\theta _1^{3/2}\theta _2^{1/2}[6\theta _1^2(\mu _{12}\theta _2\mu _{21}\theta _1)(\theta _1+\theta _2)^{1/2}`$
$`+{\displaystyle \frac{3}{2}}\mu _{21}\theta _1^2(\theta _1+\theta _2)^{1/2}(3\alpha _{12})+5\theta _1^1(\theta _1+\theta _2)^{1/2}`$
$`+{\displaystyle \frac{c_2}{16}}{\displaystyle \frac{2\theta _2(12\mu _{21}+9\mu _{12}10)\theta _1(56\mu _{21})\frac{3}{2}\mu _{21}(3\alpha _{12})(\theta _1+\theta _2)}{(\theta _1+\theta _2)^{5/2}}}],`$
$`\tau _{12}^{}`$ $`=`$ $`{\displaystyle \frac{8}{15}}x_2{\displaystyle \frac{\mu _{21}^2}{\mu _{12}}}(1+\alpha _{12})\theta _1^{3/2}\theta _2^{1/2}[6\theta _2^2(\mu _{12}\theta _2\mu _{21}\theta _1)(\theta _1+\theta _2)^{1/2}`$ (36)
$`+{\displaystyle \frac{3}{2}}\mu _{21}\theta _2^2(\theta _1+\theta _2)^{1/2}(3\alpha _{12})5\theta _2^1(\theta _1+\theta _2)^{1/2}`$
$`+{\displaystyle \frac{c_1}{16}}{\displaystyle \frac{2\theta _1(1012\mu _{12}9\mu _{21})+\theta _2(56\mu _{12})\frac{3}{2}\mu _{21}(3\alpha _{12})(\theta _1+\theta _2)}{(\theta _1+\theta _2)^{5/2}}}].`$
The corresponding expressions for $`\tau _{22}^{}`$ and $`\tau _{21}^{}`$ can be inferred from Eqs. (33) and (36) by interchanging $`12`$.

Equations (25)–(36) provide the explicit expression for the shear viscosity $`\eta `$ of a heated granular binary mixture in the first Sonine approximation. It is apparent that the reduced viscosity $`\eta ^{}`$ presents a complex nonlinear dependence on the restitution coefficients $`\alpha _{11}`$, $`\alpha _{22}`$, and $`\alpha _{12}`$ and the parameters of the mixture $`m_1/m_2`$, $`\sigma _1/\sigma _2`$, and $`n_1/n_2`$. The quality of the expression for $`\eta ^{}`$ will be later assessed by comparison with Monte Carlo simulations in the uniform shear flow (USF) problem with a thermostat. Before studying the general dependence of $`\eta ^{}`$ on the parameter space, it is instructive to consider some special limit cases. In the elastic limit, $`\alpha _{11}=\alpha _{22}=\alpha _{12}=1`$, $`\zeta ^{}=0`$, $`\gamma _i=1`$, and $`c_1=c_2=0`$, so that the expression (27) becomes

$$\eta ^{}=\frac{x_1^2R_1+x_2^2R_2+x_1x_2R_{12}}{x_1^2S_1+x_2^2S_2+x_1x_2S_{12}},$$
(37)
where

$$R_1=\frac{2}{3}+\frac{2\mu }{5},R_2=\frac{2}{3}+\frac{2}{5\mu },$$
(38)
$$R_{12}=\frac{8}{15}+\frac{\sqrt{2}}{5}\left[\left(\frac{\sigma _1}{\sigma _{12}}\right)^2\mu _{21}^{1/2}\mu _{12}^1+\left(\frac{\sigma _2}{\sigma _{12}}\right)^2\mu _{12}^{1/2}\mu _{21}^1\right],$$
(39)
$$S_1=\frac{16}{5\sqrt{2}}R_1\left(\frac{\sigma _1}{\sigma _{12}}\right)^2\mu _{21}^{1/2},S_2=\frac{16}{5\sqrt{2}}R_2\left(\frac{\sigma _2}{\sigma _{12}}\right)^2\mu _{12}^{1/2},$$
(40)
$$S_{12}=\frac{32}{15}+\frac{16}{25}\left(\frac{\sigma _1}{\sigma _{12}}\right)^2\left(\frac{\sigma _2}{\sigma _{12}}\right)^2\left(\mu _{12}\mu _{21}\right)^{1/2}.$$
(41)
Here, $`\mu =\mu _{12}/\mu _{21}=m_1/m_2`$ is the mass ratio. Equation (37) agrees with the results obtained in the first Sonine approximation to the coefficient of viscosity of a molecular gas-mixture of hard spheres . In the case of mechanically equivalent particles ($`m_1=m_2`$, $`\alpha _{11}=\alpha _{22}=\alpha _{12}\alpha `$, $`\sigma _1=\sigma _2\sigma `$), $`\gamma _i=1`$, $`\zeta ^{}=(2/3)(1\alpha ^2)(1+3c/32)`$, and

$$c_1=c_2=c=\frac{32(1\alpha )(12\alpha ^2)}{8117\alpha +30\alpha ^2(1\alpha )}$$
(42)
In this case, one gets

$$\eta ^{}=\frac{15}{4}\left((1+\alpha )\left[(2+\alpha )+\frac{c}{128}(33\alpha 39)\right]\right)^1.$$
(43)
This expression coincides with the one recently obtained for a heated granular monocomponent gas . All this shows the self-consistency of the present description.

## IV Monte Carlo simulation for uniform shear flow with thermostat

The USF is a nonequilibrium state characterized by constant partial densities $`n_i`$, a linear velocity profile $`𝐮=𝐮_i=𝖺𝐫`$, where the elements of the tensor $`𝖺`$ are $`a_k\mathrm{}=a\delta _{kx}\delta _\mathrm{}y`$, $`a`$ being the constant shear rate. In addition, the granular temperature $`T`$ and the pressure tensor $`𝖯`$ are uniform, while the mass and heat fluxes vanish by symmetry reasons. This special state is generated by Lees-Edwards boundary conditions which are simple boundary periodic boundary conditions in the local Lagrangian frame $`𝐑𝐫𝖺𝐫t`$ and $`𝐕𝐯𝖺𝐫`$. In terms of these variables the velocity distribution functions are uniform

$$f_i(𝐫,𝐯;t)=f_i(𝐕;t).$$
(44)
In the case of elastic collisions ($`\zeta =0`$) and in the absence of a thermostatting force, the energy balance equation (12) yields the heating equation

$$_tT=\frac{2}{3n}aP_{xy}.$$
(45)
Since the granular temperature $`T`$ increases in time, so does the collision frequency $`\nu (t)`$ according to Eq. (26). As a consequence, the reduced shear rate $`a^{}(t)=a/\nu (t)`$ (which is the relevant uniformity parameter) monotonically decreases with time and the system asymptotically tends towards that of (local) equilibrium. This implies that for sufficiently long times (which means here $`a^{}1`$) the system reaches a regime described by linear hydrodynamics and the Navier-Stokes shear viscosity $`\eta `$ can be identified as

$$\frac{\nu }{nT}\eta =\underset{t\mathrm{}}{lim}\frac{P_{xy}^{}}{a^{}},$$
(46)
where $`P_{xy}^{}=P_{xy}/nT`$. This route has been shown to be quite efficient to measure the Navier-Stokes shear viscosity coefficient for dilute and dense gases.

For a granular mixture, unless a thermostat is introduced, the energy balance equation (12) leads to a steady state when the viscous heating effect is exactly balanced by the collisional cooling. However, if the granular mixture is excited by the Gaussian thermostat

$$𝐅_i^{\text{th}}=\frac{1}{2}m_i\zeta 𝐕,$$
(47)
that exactly compensates for the collisional energy loss, the viscous heating still heats the system and Eq. (45) remains valid. Consequently, the linear relationship (46) allows one to determine the shear viscosity coefficient in the long time limit. Recently, this idea has been used to measure the shear viscosity of a heated granular monocomponent gas. The comparison with kinetic theory shows an excellent agreement over a wide range of values of the restitution coefficient. It must be noted that here $`\eta `$ represents the shear viscosity of an excited granular mixture and thus it does not necessarily coincide with the Navier-Stokes shear viscosity obtained in the unforced case. As a matter of fact, the results obtained in Sec. III indicate that the transport properties are affected by the Gaussian thermostat and the expression (27) for the (reduced) shear viscosity differs from the one derived in the freely cooling case. The use of thermostats to control collisional cooling in undriven systems is quite common. Usually, the motivation is to produce a steady state while here is to remove the steady state in favor of one whose dynamics determines the viscosity.

The Boltzmann equation for a mixture of inelastic hard spheres under USF and subject to the external Gaussian force (47) reads:

$$_tf_iaV_y\frac{}{V_x}f_i+\frac{1}{2}\zeta \frac{}{𝐕}\left(𝐕f_i\right)=\underset{j}{}J_{ij}\left[𝐕|f_i(t),f_j(t)\right].$$
(48)
The second term on the left-hand side represents an inertial force of the form $`𝐅_i^{\text{in}}=m_i𝖺𝐕`$, while the third term represents the thermostat force $`𝐅_i^{\text{th}}`$ given by Eq. (47). Thus, in this frame, the system is in a homogenous state subjected to the action of the (total) force $`𝐅_i^{\text{in}}+𝐅_i^{\text{th}}`$. We have numerically solved Eq. (48) by means of the Direct Simulation Monte Carlo (DSMC) method. This method was devised to mimic the processes involved in the Boltzmann collision term and its extension to deal with inelastic collisions is straightforward. In addition, since the USF is spatially homogeneous in the Lagrangian frame, the simulation method is easy to carry out and only the (peculiar) velocities of the particles need to be stored. The restriction to this homogeneous state prevents us from studying the possible formation of particle clusters (microstructure).

Technical details of the DSMC method and its application to the USF state can be found in Refs.  and . In our simulations we have typically taken a total number of particles $`N=10^5`$, a number of replicas $`𝒩=5`$, and a time step $`\mathrm{\Delta }t=3\times 10^3\mathrm{}_{11}/v_{01}`$, where $`\mathrm{}_{11}=\left(\sqrt{2}\pi n_1\sigma _1^2\right)^1`$ is the mean free path for collisions $`11`$ and $`v_{01}=\sqrt{2T/m_1}`$.

At given values of the shear rate $`a`$, the restitution coefficients $`\alpha _{ij}`$, and the parameters of the mixture, the system is initially prepared in a local equilibrium state with a temperature $`T(0)=T_0`$ such that the initial value of the reduced shear rate is $`a_0^{}=a/\nu (T_0)`$. As the system evolves, we monitor the time evolution of the reduced shear rate $`a^{}(t)=a/\nu (T(t))`$ and the reduced $`xy`$ element of the pressure tensor $`P_{xy}^{}(t)=P_{xy}(t)/nT(t)`$. We observe that in all the cases, after a transient period, the ratio $`\eta ^{}P_{xy}^{}/a^{}`$ reaches a constant value that is independent of the shear rate and time. This allows us to measure the corresponding shear viscosity coefficient $`\eta `$ as

$$\eta (t)=\frac{nT(t)}{\nu (t)}\eta ^{},$$
(49)
where the dimensionless shear viscosity $`\eta ^{}`$ is independent of time but depends on dissipation and the parameters of the mixture (masses, sizes and concentrations).

The theoretical prediction for $`\eta ^{}`$ can be obtained from the Chapman-Enskog solution to Eq. (48) up to first order in the shear rate $`a`$. In Appendix B it is easily proved that the first order solution to (48) leads to the same expression for the shear viscosity as the one obtained in Sec. III from the general Chapman-Enskog method specialized to USF. Thus, in the first Sonine approximation, the theoretical prediction of $`\eta ^{}`$ is given by Eqs. (27)–(36).

Before analyzing the dependence of the dimensionless shear viscosity coefficient $`\eta ^{}`$ on the parameters of the problem, it is instructive to test the consistency of the simulation method in the limit $`a^{}0`$ (which corresponds here to $`t\nu 1`$) . For long times and for given values of $`\alpha `$, $`m_1/m_2`$, $`\sigma _1/\sigma _2`$ , and $`n_1/n_2`$, the reduced viscosity $`\eta ^{}`$ must reach a value independent of the inital preparation of the system. In Fig. 1, we plot the shear-rate dependent viscosities $`\eta ^{}(a^{})`$ measured in the simulation, relative to its Navier-Stokes value $`\eta _\text{B}^{}`$ given by the Boltzmann theory \[Eqs. (27)–(36)\] for three different choices of the initial shear rate $`a_0^{}`$: 0.2, 0.3, and 0.4. Here, the restitution coefficient is $`\alpha =0.9`$, the mass ratio is $`m_1/m_2=4`$, the concentration ratio is $`n_1/n_2=1`$ and the size ratio is $`\sigma _1/\sigma _2=3`$. After a transient regime of a few mean free times, we observe that the curves corresponding to the three different initial conditions practically coincide. This means that a hydrodynamic regime independent of the inital conditions has been achieved. In addition, for very small values of $`a^2`$, the ratio $`\eta ^{}(a^{})/\eta _\text{B}^{}`$ fluctuate around 1 showing that in this regime the viscosity coefficient measured in the simulation is consistent with the value obtained from the Boltzmann kinetic theory. The same behavior has been found for other values of the restitution coefficient as well as of the parameters of the mixture. Notice that the limit $`a^{}0`$ is strictly unattainable in the USF because it requires an infinite amount of time. Also, the signal-to-noise ratio decreases in that limit so that the fluctuations increase.

## V Comparison between theory and simulation

Once the consistency of the simulation method has been tested, we focus our attention on the study of transport properties in the Navier-Stokes regime. In this Section we compare the predictions of the Sonine approximation with the results obtained from the DSMC method. A complete presentation of the results is complex due to the high dimensionality of the parameter space: $`\{\alpha _{11},\alpha _{22},\alpha _{12},m_1/m_2,\sigma _1/\sigma _2,n_1/n_2\}`$. For the sake of concreteness, henceforth we will assume that the spheres are made of the same material, i.e., $`\alpha _{11}=\alpha _{22}=\alpha _{12}\alpha `$. This reduces the parameter space to four quantities.

Apart from the shear viscosity coefficient, another interesting quantity at this level of description is the temperature ratio $`T_1/T_2`$. This ratio measures the breakdown of the energy equipartition. The analysis of the temperature differences has been a subject of growing interest in the past few years among both theorists  and experimentalists. As was previously found from the Boltzmann kinetic theory, except for mechanically equivalent particles, the partial temperatures $`T_i`$ are different. For the sake of illustration, Fig. 2 shows the dependence of the temperature ratio on the size ratio $`\sigma _1/\sigma _2`$ for an equimolar mixture ($`n_1/n_2=1`$) and three different values of the restitution coefficient $`\alpha =0.9,0.8`$, and $`0.7`$. We consider a binary mixture of constant density and so, $`m_1/m_2=(\sigma _1/\sigma _2)^3`$. We observe that for large size ratios the temperature differences are quite important, even for moderate dissipation. It is also apparent that an excellent agreement between the theory (given by the first Sonine correction) and Monte Carlo simulations (symbols) is found over the entire range of values of size and mass ratios considered.

Next, we explore the influence of dissipation on the reduced shear viscosity $`\eta ^{}(\alpha )`$ for different values of the mass ratio, the size ratio, and the concentration ratio. Three different values of the (common) restitution coefficient are considered: $`\alpha =0.9,0.8`$, and $`0.7`$. In Fig. 3, we plot the ratio $`\eta ^{}(\alpha )/\eta ^{}(1)`$ versus the mass ratio $`m_1/m_2`$ for $`\sigma _1/\sigma _2=n_1/n_2=1`$. Here, $`\eta ^{}(1)`$ refers to the elastic value for the shear viscosity coefficient. Again, the symbols represent the simulation data while the lines refer to the theoretical results obtained from the Boltzmann equation in the first Sonine approximation. We see that in general the deviation of $`\eta ^{}(\alpha )`$ from its functional form for elastic collisions is quite important. This tendency becomes more significant as the mass disparity increases. The agreement between the first Sonine approximation and simulation is seen to be in general excellent. This agreement is similar to the one previously found in the monocomponent case. At a quantitative level, the discrepancies between theory and simulation tend to increase as the restitution coefficient decreases, although these differences are quite small (say, for instance, around 2% at $`\alpha =0.7`$ in the disparate mass case $`m_1/m_2=10`$).

The influence of the size ratio on the shear viscosity is shown in Fig. 4 for $`m_1/m_2=4`$ and $`n_1/n_2=1`$. We observe again a strong dependence of the shear viscosity on dissipation. However, for a given value of $`\alpha `$, the influence of $`\sigma _1/\sigma _2`$ on $`\eta ^{}`$ is weaker than the one found before in Fig. 3 for the mass ratio. The agreement for both $`\alpha =0.9`$ and $`\alpha =0.8`$ is quite good, except for the largest size ratio at $`\alpha =0.8`$. These discrepancies become more significant as the dissipation increases (say $`\alpha =0.7`$), especially for mixtures of particles of very different sizes. Finally, Fig. 5 shows the dependence of $`\eta ^{}(\alpha )/\eta ^{}(1)`$ on the concentration ratio for $`m_1/m_2=4`$ and $`\sigma _1/\sigma _2=1`$. We observe that both the theory and simulation predicts a very weak influence of composition on the shear viscosity. With respect to the influence of dissipation, the trends are similar to those of Figs. 3 and 4: the main effect of inelasticity in collisions is to enhance the momentum transport with respect to the case of elastic collisions. The agreement now between theory and simulation is very good, even for disparate values of the concentration ratio and/or strong dissipation. Therefore, according to the comparison carried out in Figs. 3, 4, and 5, we can conclude that the agreement extends over a wide range values of the restitution coefficient, indicating the reliability of the first Sonine approximation for describing granular flows beyond the quasielastic limit.

## VI Discussion

Although the utility of a hydrodynamic description for granular media under rapid flow conditions has been recognized for many years, its domain of validity as well as the forms of the transport coefficients remain a topic of interest and controversy. In this context, there are some doubts about the possibility of going from a kinetic theory to a hydrodynamic level of description by using a Chapman-Enskog expansion around the homogenous cooling state. Given that the search for exact solutions of the Boltzmann equation is far beyond the present perspectives, an alternative to get some insight into the above question is to numerically solve the kinetic equation and compare these results with the corresponding solution obtained by assuming the validity of a hydrodynamic description. In this paper, we have performed such a comparison at the level of the shear viscosity coefficient of a heated granular mixture. The system is heated by the action of a thermostatting external force which exactly compensates for cooling effects associated with the inelasticity of collisions. Although some previous works  have compared kinetic theory predictions for transport coefficients with computer simulations in the case of a monocomponent gas, studies for multicomponent granular gases are more scarce. Very recently, a seemingly similar analysis for the shear viscosity $`\eta `$ of a dense mixture has been given in Ref. . Nevertheless, the above kinetic theory only holds for nearly elastic particles and the expression of $`\eta `$ in the first Sonine approximation coincides with the one obtained in the elastic case.

As a first step in our issue, in Sec. III we have derived the general hydrodynamic equations of a heated binary mixture of smooth inelastic spheres from the Boltzmann kinetic equation by using the Chapman-Enskog method. The corresponding transport coefficients have been expressed in terms of the solution to integral equations, which are then solved approximately (first Sonine polynomial approximation) just as in the case of elastic collisions. The explicit expressions for the set of relevant transport coefficients $`\{D,D_p,D^{},D^{\prime \prime },L,\lambda ,\eta \}`$ are displayed in Appendix A. In contrast to previous works, our results are not limited a priori to weak inelasticity and they take into account the effect of the temperature differences on the transport coefficients. On the other hand, the results obtained here for the transport coefficients slightly differ from those obtained in the freely cooling case , showing that in general the introduction of a thermostat affects the transport properties of the gas . The Chapman-Enskog results obtained for the mixture have been then specialized to the hydrodynamic state of transverse shear. In this state, the shear viscosity coefficient $`\eta `$ is the relevant transport coefficient of the problem. The explicit form of $`\eta `$ is given by Eqs. (27)–(36) in terms of the restitution coefficients $`\alpha _{ij}`$ and the parameters of the mixture (masses, diameters, and concentrations).

To test the assumptions of the Chapman-Enskog method and the approximate Sonine solution to the resulting integral equation, the DSMC method has been used to solve the Boltzmann equation in the uniform shear flow state. In the absence of a thermostat, in a granular fluid there is a competition between two opposite effects: viscous heating and collisional cooling. In that case, when both effects exactly cancel each other, a steady state is reached after a transient period. In this steady state, due to the coupling between dissipation and the shear rate, the system is far away from the Navier-Stokes regime, except when $`\alpha 1`$ . However, if the external thermostat is adjusted to compensate for the energy lost in collisions, the shearing work still heats the system. As a consequence, as the system evolves, the reduced shear rate $`a^{}(t)`$ goes to zero and the system achieves a regime described by linear hydrodynamics. In this regime, the Navier-Stokes shear viscosity coefficient can be measured from simulations. In this paper, the thermostat is used to remove the steady state in favor of a time-dependent state whose dynamics allows one to get the Navier-Stokes shear viscosity just as for the case of elastic collisions.

The dependence of the viscosity $`\eta `$ on the full parameter space has been explored. Specifically, the parameter space over which our solution has been verified is the mass ratio $`m_1/m_2`$, the concentration ratio $`n_1/n_2`$, the ratio of diameters $`\sigma _1/\sigma _2`$, and the (common) restitution coefficient $`\alpha \alpha _{11}=\alpha _{11}=\alpha _{22}=\alpha _{12}`$. The theory and simulation clearly show how in general, the influence of dissipation on momentum transport is quite important since there is a relevant dependence of the viscosity $`\eta (\alpha )`$ on the restitution coefficient $`\alpha `$. At a given value of the restitution coefficient, the dependence of $`\eta (\alpha )/\eta (1)`$ on the mass ratio is more significant than the one found on the composition and diameters. This feature has been also found for the temperature ratio in the experiments recently carried out in vibrated mixtures , although experimental confirmation of the trends observed here for the viscosity is still lacking. With respect to the accuracy of the theory predictions we see that, in general, the Chapman-Enskog results in the first Sonine approximation exhibit an excellent agreement with the simulation data. This supports the idea that the Sonine polynomial approximation for granular fluids has an accuracy comparable to that for elastic collisions. Exceptions to this agreement are extreme mass ratios and strong dissipation. These discrepancies are due basically to the approximations introduced in applying the Chapman-Enskog method, and more specifically in using the first Sonine approximation.

One of the main limitations of the results obtained here from the Boltzmann equation is its restriction to the low-density regime. In this regime, the collisional transfer contributions to the fluxes are negligible and only their kinetic contributions are taken into account. Possible extension in both aspects, theory and simulation, of the present simple hydrodynamic state to higher densities can be carried out in the context of the revised Enskog theory. In this case, many of the phenomena appearing in dense granular fluids (such as spontaneous formation of dense clusters surrounded by regions of low-density ) could be studied. On the other hand, although the comparison performed here has been made undergoing uniform shear flow without paying attention to the possible formation of particle clusters, our Chapman-Enskog results apply for general inhomogeneous situations. The only restriction is that they provide the irreversible parts of the mass, heat, and momentum fluxes to leading order in the spatial gradients of the hydrodynamic fields. In this context, the results derived in this paper can be used to analyze the behavior of granular mixtures in a lot of physical situations. Thus, for instance, the knowledge of the complete hydrodynamic equations for a binary mixture allows one to say whether the mixture hydrodynamics is more or less unstable to long-wavelength perturbations than that of the one-component case, and what are the mechanisms involved in phenomena very often observed in nature and experiments such as phase separation or segregation. We hope that the present results give some insight into the understanding of these interesting and complex problems.

###### Acknowledgements.

V. G. acknowledges partial support from the Ministerio de Ciencia y Tecnología (Spain) through Grant No. BFM2001-0718.

## A Chapman-Enskog expansion

In this Appendix, the expressions of the transport coefficients for a heated granular mixture are obtained. The derivation follows similar steps as those made in Ref.  in the free cooling case. Here, we will use the same notation as in Ref. . In the first order, the distribution function $`f_i^{(1)}`$ verifies the kinetic equation

$$\left(_i+\frac{1}{2}\zeta ^{(0)}\frac{}{𝐕}𝐕\right)f_i^{(1)}+_if_j^{(1)}=\left(D_t^{(1)}+𝐕\right)f_i^{(1)},$$
(A1)
where $`D_t^{(1)}=_t^{(1)}+𝐮`$, and

$$_if_i^{(1)}=\left(J_{ii}[f_i^{(0)},f_i^{(1)}]+J_{ii}[f_i^{(1)},f_i^{(0)}]+J_{ij}[f_i^{(1)},f_j^{(0)}]\right),$$
(A2)
$$_if_j^{(1)}=J_{ij}[f_i^{(0)},f_j^{(1)}].$$
(A3)
In these equations, it is understood that $`ij`$ and use has been made of the fact that $`_t^{(0)}T=0`$ and the results $`𝐣_i^{(0)}=𝐪^{(0)}=\zeta ^{(1)}=0`$. The last equality follows from the fact that the cooling rate is a scalar, and so $`\zeta ^{(1)}`$ should be proportional to $`𝐮`$. However, as shown later, there is no contribution to $`f_i^{(1)}`$ proportional to the divergence of the flow field so that $`\zeta ^{(1)}=0`$ by symmetry. This property is special to the low density Boltzmann kinetic theory and such terms occur at higher densities . The macroscopic balance equations to first order are

$$D_t^{(1)}x_1=0,\frac{3}{5}D_t^{(1)}\mathrm{ln}p=\frac{3}{2}D_t^{(1)}\mathrm{ln}T=𝐮,D_t^{(1)}𝐮=\rho ^1p.$$
(A4)
Use of these in (A1) yields

$$\left(_i+\frac{1}{2}\zeta ^{(0)}\frac{}{𝐕}𝐕\right)f_i^{(1)}+_if_j^{(1)}=𝐀_ix_1+𝐁_ip+𝐂_iT+D_{i,\alpha \beta }_\alpha u_\beta ,$$
(A5)
where

$$𝐀_i(𝐕)=\left(\frac{}{x_1}f_i^{(0)}\right)_{p,T}𝐕,$$
(A6)
$$𝐁_i(𝐕)=\frac{1}{p}\left[f_i^{(0)}𝐕+\frac{nT}{\rho }\left(\frac{}{𝐕}f_i^{(0)}\right)\right],$$
(A7)
$$𝐂_i(𝐕)=\frac{1}{T}\left[f_i^{(0)}+\frac{1}{2}\frac{}{𝐕}\left(𝐕f_i^{(0)}\right)\right]𝐕,$$
(A8)
$$D_{i,\alpha \beta }(𝐕)=V_\alpha \frac{}{V_\beta }f_i^{(0)}\frac{1}{3}\delta _{\alpha \beta }𝐕\frac{}{𝐕}f_i^{(0)}.$$
(A9)
The solutions to Eqs. (A5) are of the form

$$f_i^{(1)}=𝒜_ix_1+_ip+𝒞_iT+𝒟_{i,\alpha \beta }_\alpha u_\beta .$$
(A10)
The coefficients $`𝒜_i,_i,𝒞_i,`$ and $`𝒟_{i,\alpha \beta }`$ are functions of the peculiar velocity $`𝐕`$ and the hydrodynamic fields. The cooling rate depends on space through its dependence on $`x_1`$, $`p`$, and $`T`$. The integral equations for the unknowns are easily identified as coefficients of the independent gradients in (A10). The result is

$$\left(_i+\frac{1}{2}\zeta ^{(0)}\frac{}{𝐕}𝐕\right)\left(\begin{array}{c}𝒜_i\\ _i\\ 𝒞_i\\ 𝒟_{i,\alpha \beta }\end{array}\right)+_i\left(\begin{array}{c}𝒜_j\\ _j\\ 𝒞_j\\ 𝒟_{j,\alpha \beta }\end{array}\right)=\left(\begin{array}{c}𝐀_i\\ 𝐁_i\\ 𝐂_i\\ 𝐃_{i,\alpha \beta }\end{array}\right).$$
(A11)
Note that, in contrast to what happens in the free cooling case, here each one of the quantities $`𝒜_i`$, $`_i`$, $`𝒞_i`$, and $`𝒟_{i,\alpha \beta }`$ obey closed integral equations. The solution to Eq. (A11) provides the expression for the transport coefficients. In the case of the mass flux $`𝐣_1`$, these coefficients are identified as

$$D=\frac{\rho }{3m_2n}𝑑𝐯𝐕𝒜_1,$$
(A12)
$$D_p=\frac{m_1p}{3\rho }𝑑𝐯𝐕_1,$$
(A13)
$$D^{}=\frac{m_1T}{3\rho }𝑑𝐯𝐕𝒞_1.$$
(A14)
The transport coefficients for the heat flux are

$$D^{^{\prime \prime }}=\frac{1}{3T^2}\underset{i=1}{\overset{2}{}}𝑑𝐯\frac{1}{2}m_iV^2𝐕𝒜_i,$$
(A15)
$$L=\frac{1}{3}\underset{i=1}{\overset{2}{}}𝑑𝐯\frac{1}{2}m_iV^2𝐕_i,$$
(A16)
$$\lambda =\frac{1}{3}\underset{i=1}{\overset{2}{}}𝑑𝐯\frac{1}{2}m_iV^2𝐕𝒞_i.$$
(A17)
Finally, the shear viscosity is given by

$$\eta =\frac{1}{10}\underset{i=1}{\overset{2}{}}𝑑𝐯m_iV_\alpha V_\beta 𝒟_{i,\alpha \beta }.$$
(A18)
Accurate aproximations to the solutions to the integral equations for $`(𝒜_i,_i,𝒞_i,𝒟_{i,\alpha \beta })`$ may be obtained using low order truncation of expansions in a series of Sonine polynomials. In the case of the mass flux, we consider the leading Sonine approximation (lowest degree polynomial)

$$\{𝒜_i,_i,𝒞_i\}f_{i,M}𝐕\{a_{i,1},b_{i,1},c_{i,1}\},f_{i,M}(𝐕)=n_i(m_i/2\pi T_i)^{3/2}\mathrm{exp}(m_iV^2/2T_i),$$
(A19)
where $`a_{i,1}=(m_1m_2n/\rho n_iT_i)D`$, $`b_{i,1}=(\rho /pn_iT_i)D_p`$, and $`c_{i,1}=(\rho /Tn_iT_i)D^{}`$. The coefficients $`a_{i,1}`$, $`b_{i,1}`$, and $`c_{i,1}`$ are determined by multiplying the three first equations of Eq. (A11) by $`m_i𝐕`$ and integrating over the velocity. The result is

$$a_{1,1}=\left(\nu _D\frac{1}{2}\zeta ^{(0)}\right)^1\left(\frac{}{x_1}\mathrm{ln}n_1T_1\right)_{p,T},$$
(A20)
$$b_{1,1}=\left(\nu _D\frac{1}{2}\zeta ^{(0)}\right)^1\frac{n_1T_1}{p}\left(1\frac{m_1nT}{\rho T_1}\right),$$
(A21)
$$c_{1,1}=0.$$
(A22)
Here, the collision frequency $`\nu _D`$ is given by Eq. (73) of Ref. .

In the case of the pressure tensor, the leading Sonine approximation for the function $`𝒟_{i,\alpha \beta }`$ is

$$𝒟_{i,\alpha \beta }f_{i,M}d_{i,1}R_{i,\alpha \beta },R_{i,\alpha \beta }=m_i(V_\alpha V_\beta \frac{1}{3}V^2\delta _{\alpha \beta }).$$
(A23)
The shear viscosity coefficient is given by

$$\eta =nT^2\left(x_1\gamma _1^2d_{1,1}+x_2\gamma _2^2d_{2,1}\right).$$
(A24)
The coefficients $`d_{i,1}`$ can be determined by multiplying the fourth equation of Eq. (A11) by $`R_{i,\alpha \beta }`$ and integrating over the velocity to get the coupled set of equations

$$\left(\begin{array}{cc}\tau _{11}\zeta ^{(0)}& \tau _{12}\\ \tau _{21}& \tau _{22}\zeta ^{(0)}\end{array}\right)\left(\begin{array}{c}d_{1,1}\\ d_{2,1}\end{array}\right)=\left(\begin{array}{c}T_1^1\\ T_2^1\end{array}\right).$$
(A25)
The collision frequencies $`\tau _{ij}=\tau _{ij}^{}\nu `$, where $`\tau _{ij}^{}`$ are given by Eqs. (33) and (36). From Eq. (A25) one easily gets the expression (25) for the shear viscosity given in the main text.

The calculations for the heat flux are similar to those previously made for the other fluxes. As in the unforced case, this requires going to the second Sonine approximation. In this case, the transport coefficients defining the heat flux (17) are given by

$$D^{\prime \prime }=\frac{5}{2}T\left(\frac{n_1\gamma _1^3}{m_1}a_{1,2}+\frac{n_2\gamma _2^3}{m_2}a_{2,2}\right)+\frac{5}{2}\frac{nm_1m_2}{\rho T}\left(\frac{\gamma _1}{m_1}\frac{\gamma _2}{m_2}\right)D,$$
(A26)
$$L=\frac{5}{2}T^3\left(\frac{n_1\gamma _1^3}{m_1}b_{1,2}+\frac{n_2\gamma _2^3}{m_2}b_{2,2}\right)+\frac{5}{2}\frac{\rho }{n}\left(\frac{\gamma _1}{m_1}\frac{\gamma _2}{m_2}\right)D_p,$$
(A27)
$$\lambda =\frac{5}{2}T^3\left(\frac{n_1\gamma _1^3}{m_1}c_{1,2}+\frac{n_2\gamma _2^3}{m_2}c_{2,2}\right)+\frac{5}{2}\rho \left(\frac{\gamma _1}{m_1}\frac{\gamma _2}{m_2}\right)D^{}.$$
(A28)
The coefficients $`a_{i,2}`$, $`b_{i,2}`$ and $`c_{i,2}`$ obey the equations

$$\left(\begin{array}{cc}\nu _{11}\frac{3}{2}\zeta ^{(0)}& \nu _{12}\\ \nu _{21}& \nu _{22}\frac{3}{2}\zeta ^{(0)}\end{array}\right)\left(\begin{array}{c}a_{1,2}\\ a_{2,2}\end{array}\right)=\left(\begin{array}{c}X_1\\ X_2\end{array}\right),$$
(A29)
$$\left(\begin{array}{cc}\nu _{11}\frac{3}{2}\zeta ^{(0)}& \nu _{12}\\ \nu _{21}& \nu _{22}\frac{3}{2}\zeta ^{(0)}\end{array}\right)\left(\begin{array}{c}b_{1,2}\\ b_{2,2}\end{array}\right)=\left(\begin{array}{c}Y_1\\ Y_2\end{array}\right),$$
(A30)
$$\left(\begin{array}{cc}\nu _{11}\frac{3}{2}\zeta ^{(0)}& \nu _{12}\\ \nu _{21}& \nu _{22}\frac{3}{2}\zeta ^{(0)}\end{array}\right)\left(\begin{array}{c}c_{1,2}\\ c_{2,2}\end{array}\right)=\left(\begin{array}{c}Z_1\\ Z_2\end{array}\right),$$
(A31)
where

$`X_1`$ $`=`$ $`{\displaystyle \frac{\zeta ^{(0)}m_1m_2nD}{\rho n_1T_1^2}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{T^2}{n_1T_1^3}}{\displaystyle \frac{}{x_1}}\left(n_1\gamma _1^2c_1\right)+{\displaystyle \frac{2}{15}}{\displaystyle \frac{m_1^2m_2nD}{\rho n_1^2T_1^4}}[{\displaystyle }d𝐯_1𝐒_1_1(f_{1,M}𝐕_1)`$ (A33)
$`\delta \gamma {\displaystyle }d𝐯_1𝐒_1_1(f_{2,M}𝐕_2)],`$
$`Y_1`$ $`=`$ $`{\displaystyle \frac{\zeta ^{(0)}\rho D_p}{pn_1T_1^2}}{\displaystyle \frac{1}{2}}{\displaystyle \frac{c_1}{pT_1}}+{\displaystyle \frac{2}{15}}{\displaystyle \frac{m_1\rho D_p}{pn_1^2T_1^4}}[{\displaystyle }d𝐯_1𝐒_1_1(f_{1,M}𝐕_1)`$ (A35)
$`\delta \gamma {\displaystyle }d𝐯_1𝐒_1_1(f_{2,M}𝐕_2)],`$
$`Z_1`$ $`=`$ $`{\displaystyle \frac{\zeta ^{(0)}\rho D^{}}{Tn_1T_1^2}}{\displaystyle \frac{2+c_1}{2TT_1}}+{\displaystyle \frac{2}{15}}{\displaystyle \frac{m_1\rho D^{}}{Tn_1^2T_1^4}}[{\displaystyle }d𝐯_1𝐒_1_1(f_{1,M}𝐕_1)`$ (A37)
$`\delta \gamma {\displaystyle }d𝐯_1𝐒_1_1(f_{2,M}𝐕_2)],`$
and

$$c_i=\frac{8}{15}\left[\frac{m_i^2}{4n_iT_i^2}𝑑𝐯_1V_1^4f_i^{(0)}\frac{15}{4}\right].$$
(A38)
The corresponding expressions of the elements $`X_2`$, $`Y_2`$ and $`Z_2`$ can be deduced from Eqs. (A33), (A35) and (A37), respectively, by interchanging $`12`$ and setting $`DD`$, $`D_pD_p`$ and $`D^{}D^{}`$. The frequencies $`\nu _{ij}`$ and the collision integrals appearing in Eqs. (A31), (A33), (A35) and (A37) were explicitly evaluated in the Appendix D of Ref. . Thus, the transport coefficients $`D^{\prime \prime }`$, $`L`$, and $`\lambda `$ are completely determined.

## B First order solution to the USF

In this Appendix we get the solution to Eq. (48) in the first order in the shear rate $`a`$. The normal solution to Eq. (48) is provided by the Chapman-Enskog method, i.e., a solution given as a power series in $`a`$:

$$f_i=f_i^{(0)}+f_i^{(1)}+\mathrm{}$$
(B1)
The zeroth-order solution $`f_i^{(0)}`$ verifies Eq. (21) and it corresponds to the homogeneous cooling state distribution in the local Lagrangian frame. Its first Sonine approximation is given by Eq. (24). Inserting the expansion (B1) into Eq. (21) leads to the following integral equation for $`f_i^{(1)}`$:

$$_tf_i^{(1)}aV_y\frac{}{V_x}f_i^{(0)}+\frac{1}{2}\zeta ^{(0)}\frac{}{𝐕}\left(𝐕f_i^{(1)}\right)+\frac{1}{2}\zeta ^{(1)}\frac{}{𝐕}\left(𝐕f_i^{(0)}\right)=_if_i^{(1)}_if_j^{(1)},$$
(B2)
where the operators $`_i`$ and $`_i`$ are defined by Eqs. (A2) and (A3), respectively. Since $`f_i^{(1)}`$ depends on time only through the temperature, Eq. (45) implies that $`_tf_i^{(1)}=𝒪(a^2)`$ and so the first term on the left hand side of Eq. (B1) vanishes in the first order. Further, $`\zeta ^{(1)}=0`$ by symmetry because $`𝐮=0`$ in the USF. Taking into acount the above properties, Eq. (B2) reduces to

$$\left(_i+\frac{1}{2}\zeta ^{(0)}\frac{}{𝐕}𝐕\right)f_i^{(1)}+_if_j^{(1)}=aV_y\frac{}{V_x}f_i^{(0)}.$$
(B3)
This integral equation is identical to Eq. (A5) when one particularizes the latter one to USF. Therefore, the expression for the shear viscosity obtained from (B3) is given by Eqs. (25)–(36).