text:
# 65~Cybele in the thermal infrared: Multiple observations and thermophysical analysisPartly based on observations with ISO, an ESA project with instruments funded by ESA Member States (especially the PI countries: France, Germany, the Netherlands and the United Kingdom) and with the participation of ISAS and NASA.

## 1 Introduction

The asteroid 65~Cybele is a main-belt asteroid and the main representative of a group of asteroids with semi-major axis in the range 3.27$`<`$a$`<`$3.70 AU, eccentricities e$``$0.30 and inclinations of i$``$25 (Gradie et al. gradie89 (1989)). 65~Cybele is generally considered as one of the 10 largest asteroids, with a published diameter between 230 (Dobrovol’skij dobrovol80 (1980); Taylor taylor81 (1981)) and about 330 km (Green et al. green85 (1985)). Although discovered in 1863, it took until 1980 before its rotational properties were determined as one of the last remaining large minor planets (Schober et al. schober80 (1980)). Recent analysis of many lightcurve observations lead to a rotation period of about 6 hours, a unique spin vector solution pair with a retrograde sense of rotation and an elongated body shape (Erikson erikson00 (2000)).

Based on its low albedo and the flat to slightly reddish spectrum in the 0.3 to 1.1 $`\mu `$m range, it was classified as P-type object (Tholen tholen89 (1989)) and due to a weak ultraviolet absorption feature short ward of 0.4 $`\mu `$m as C0-type (Barucci et al. barucci87 (1987)). Recently Bus (bus99 (1999)) classified 65~Cybele as Xc type, based on reflectance spectra covering the wavelength interval from 0.44 to 0.92 $`\mu `$m. These taxonomic types indicate hydrated silicates, carbon and organics on the asteroids surface (Gaffey et al. gaffey89 (1989)).

Another interesting aspect comes from the occultation measurement: Dobrovol’skij (dobrovol80 (1980)) reported a second event during the occultation of AGK3 +19 599 which was attributed to a satellite of 11 km diameter, located 917 km from the center of 65~Cybele (in the case of a central occultation). So far, no confirmation of the satellite through other techniques was published. It was one of the targets in the HST search for binaries ‘Imaging Snapshots of Asteroids’ Storrs et al. (storrs99 (1999)) & Storrs (priv. comm.), but with unknown outcome so far. No radar measurements are available yet.

Tedesco & Veeder (tedesco92 (1992)) and Tedesco et al. (tedesco02 (2002)) derived from IRAS observations diameter and albedo values of 237.26$`\pm `$4.2 km and 0.0706$`\pm `$0.003, respectively. The only direct size information from the occultation event in 1979 resulted in a diameter of 230$`\pm `$16 km, with a largest measured chord of 245 km (Dobrovol’skij dobrovol80 (1980); Taylor taylor81 (1981); Xiu-Yi xiuyi79 (1979)). Both values, the radiometric and the occultation diameter, agree nicely. But ground-based N-band observations (Morrison & Chapman morrison76 (1976); Morrison morrison77 (1977); Bowell et al. bowell79 (1979); Green et al. green85 (1985)) lead to diameters larger than 300 km. The published radiometric diameters and albedos were all based on the concept of the Standard Thermal Model (STM; Lebofsky et al. lebofsky86 (1986); Lebofsky & Spencer lebofsky89 (1989)), but using different values for the infrared beaming parameter $`\eta `$, the phase angle correction and the emissivity. Using for example the STM $`\eta `$-value of 0.756 (Lebofsky et al. lebofsky86 (1986)) would lower some of the earlier 300 km diameter values to about 250 km, but the Green et al. (green85 (1985)) values would still be close to 300 km. The higher phase angle correction of Bowell et al. (bowell79 (1979)) would again increase the diamter of all measurements taken at large phase angles.

The diameter is the most fundamental parameter, but in case of 65~Cybele it has very large uncertainties, much larger than for other asteroids of comparable size. Where is this uncertainty coming from? Are the data sets too different or some of them unreliable? What role play the STM concept, the simplification in the shape, the missing thermal and rotation effects or the various other model input parameters?

The first thermal spectrum of 65~Cybele was published by Green et al. (green85 (1985)). They did not see any pronounced emission features in the 8 to 13 $`\mu `$m wavelength range. In the years 1996–1998 65~Cybele was extensively observed by ISO with the scientific goal of providing new insights into the chemical composition and the mineralogy for a sample of main-belt asteroids. Are there any spectral features visible? If yes, what are the meteoritic and mineralogic counterparts?

With the availability of a recently established thermophysical model (TPM) by Lagerros (lagerros96 (1996); lagerros97 (1997); lagerros98 (1998)) we are now in a position to combine the information on the asteroid, its rotational behaviour, different observing and illumination geometries with the thermal measurements taken at different wavelengths and different times.

In the following, we show the results of the re-analysis and re-calibration of the existing thermal observations from ground- and satellite-based instruments, including photometric and spectroscopic measurements (Sect. 2). In Sect. 3 we present the basic parameters of our TPM analysis and derive radiometric diameters and albedos for all photometrically reliable data sets (Sect. 3.2). In the spectroscopic analysis (Sect. 3.3) we focused on the modelling of the shape of the thermal continuum emission which is strongly connected to thermal effects on a regolith covered surface. The TPM analysis section is followed by extensive discussions of radiometric and spectroscopic results (Sect. 4). We used a large sample of infrared observations, spanning a wide range of phase angles and wavelengths, to demonstrate the possibilities and limitations of the thermophysical modelling of asteroids (Sect. 5).

## 2 Observations in the Thermal Infrared

### 2.1 Re-analysis of previous measurements

#### TRIAD Observations

The original TRIAD111Tucson revised Index of Asteroid Data fluxes were based on three 65~Cybele observations from March 18, 1975, 08:24 UT, with an average N-band magnitude of $`+`$0.85$`\pm `$0.05 mag (Morrison & Chapman morrison76 (1976)). We used the Rieke et al. (rieke85 (1985)) absolute photometric system for the standard stars to re-calibrate this flux. This led to a new N-band value of $`+0.92`$ mag, corresponding to a colour corrected flux at 10.6 $`\mu `$m of 14.55$`\pm `$0.87 Jy (measurement and calibration errors have been added quadratically).

#### IR Spectroscopy from UCL

Green et al. (green85 (1985)) presented the first 8- to 13 $`\mu `$m spectrum of 65~Cybele. It was taken on April 18, 1982 and calibrated against $`\alpha `$ Boo. We re-calibrated the 10.0 $`\mu `$m flux value to 25.5$`\pm `$2.9 Jy with the Cohen et al. (cohen95 (1995)) model for $`\alpha `$ Boo (the Green et al. value for $`\alpha `$ Boo is 15.0 % higher than the Cohen et al. value at 10.0 $`\mu `$m). The re-calibration resulted in flux values of 14.4$`\pm `$2.0 Jy at 8.5 $`\mu `$m, 25.5$`\pm `$2.9 Jy at 10.0 $`\mu `$m and 50.9$`\pm `$4.4 Jy at 12.0 $`\mu `$m for the given epoch.

#### IRAS Observations

IRAS (Beichman et al. beichman88 (1988)) scanned 6 times over 65~Cybele, each time with all 4 wavelength bands. The latest flux values and the corresponding analysis were presented by Tedesco et al. (tedesco02 (2002)). We used the colour corrected monochromatic IRAS fluxes (without applying the STM-related band-to-band corrections) for our reanalysis (see also Müller & Lagerros mueller98 (1998)).

#### IRTF Observations

We included in our analysis 3 IRTF measurements (D. Osip, priv. communications, see Table 1):

No observational or calibration errors were given, but the data were considered as a high accuracy dataset (Osip, priv. comm.). We assumed a 10 % flux uncertainty for the N-band and 15 % for the Q-band photometry.

#### UKIRT Observations

The UKIRT mid-IR N- and Q-band data of 65~Cybele have been obtained in August 1996 with the MPIA 5–25 $`\mu `$m camera MAX (Robberto & Herbst robberto98 (1998)). Observations, data reduction and calibration are described in Müller & Lagerros (mueller98 (1998)).

#### ISOPHOT photometric data

65~Cybele was measured at wavelengths between 65 and 200 $`\mu `$m with ISOPHOT (Lemke et al. lemke96 (1996)) at 3 different observing epochs. Observations, data reduction and calibration are discussed in Müller & Lagerros (mueller98 (1998); mueller02 (2002)).

### 2.2 ISOPHT PHT-S Measurement

The observational details of the PHT-S measurements are given in Table 2. The PHT-S short wavelength part (2.47–4.87 $`\mu `$m) covers the transition between reflected light and thermal emission. It has been excluded in our thermophysical analysis. The data reduction of the 5.84–11.62 $`\mu `$m part was done in a standard way, using PIA222ISOPHOT Interactive Analysis (for details of the PHT-S data reduction steps see Dotto et al. dotto02 (2002)). There was no dedicated background measurement, we therefore used COBE-DIRBE measurement results (Hauser et al. hauser98 (1998)) as replacement. The color corrected, weekly maps from the COBE-DIRBE data sets were adjusted for the PHT-S beam size. The two DIRBE values at 4.9 and 12 $`\mu `$m were then connected using a 266 K blackbody (Ábrahám et al. abraham99 (1999)) and subtracted from the measured spectrum of 65~Cybele. The background contribution was between 4 and 10 % in the PHT-S wavelength range.

Figure 1 shows the reduced and background subtracted result of the PHT-S observation together with the corresponding TPM prediction. Typical signal processing uncertainties were between 2 and 6 %. The absolute calibration error for bright point-like sources in staring mode is less than 10 % (Laureijs et al. laureijs02 (2002)). For our analysis we added quadratically the individual signal processing errors and an absolute calibration error of 5 %. The lower part of Fig. 1 shows the ratio between the observed spectrum and the TPM prediction (see Sect. 3).

### 2.3 ISOCAM CVF & Filter Observations

On Christmas day 1997, a for- and backward scan was performed in the 5 to 17 $`\mu `$m wavelength range with the ISOCAM Circular Variable Filter (CVF) (Cesarsky et al. cesarsky96 (1996)). The scan was performed skipping one CVF position per CVF step, corresponding to steps of approximately 0.2 $`\mu `$m for the 9-17 $`\mu `$m and 0.1 $`\mu `$m for the shorter wavelength range. The pixel field of view was 3′′ leading to a total field of view of 96′′. The data reduction was done within the CAM Interactive Analysis package (Ott et al. ott97 (1997)) and followed the same strategy as was done for the standard stars observations, used for the calibration of the Spectral Response Function (Blommaert et al. blommaert01 (2001)). The flux calibration was done using the OLP v10 SRF333Off Line Processing version 10, Spectral Response Function. Aperture photometry was applied avoiding the contribution from reflected light on the detector (Blommaert et al. blommaert01 (2001)). The strategy of combining a for- and backward scan allows to check observed spectral features and also allows to a certain extent the correction of the ISOCAM transient response. We applied the Fouks-Schubert method (Coulais & Abergel coulais00 (1996)) to correct the transient of the signal. The CVF scan starting at 5 $`\mu `$m after the fixed filter measurements suffers strongly from a remnant which is not corrected for by the transient correction method. For this reason we judge the data for 5 - 7 $`\mu `$m in the up-scan as less reliable.

In the time between the two CVF scans nine fixed filter measurements were performed in the 6–15 $`\mu `$m wavelength range. The zodiacal background becomes important at wavelengths longer than 9 $`\mu `$m. Therefore, dedicated background measurements were done in the long-wavelength filters. They were performed before 65~Cybele itself was observed to avoid remnants of the 65~Cybele measurement because of the transients. The on-source measurements were long enough to allow stabilisation of the signal so that no correction had to be applied (Blommaert et al. blommaert00 (2000)). The obtained flux densities have been colour corrected using the average CVF as input spectrum.

The resulting monochromatic flux densities are listed in Table 3. These values are also plotted together with error bars in Fig. 2. The lower part of the figure shows the ratio between the observed CAM filter values and the corresponding TPM prediction (see Sect. 3).

### 2.4 ESO-TIMMI2 N- and Q-band Observations

We obtained additional N- and Q-band observations of 65~Cybele using the TIMMI2 camera (Reimann et al. reimann00 (2000)), mounted on the ESO-La Silla 3.6 m telescope. Table 4 summarises the observational results444Based on observations collected at the European Southern Observatory, Chile; ESO, No. 69.C-0090.

We did a standard data reduction for the chopping-nodding imaging mode. Asteroid and calibration star measurements have been taken within 30 min, all N1 integration time were 144 sec, the Q1 integration times were 108 sec. The relatively large uncertainties in table 4 are due to the variable atmospheric conditions. The flux errors are estimates based on the monitoring of the reference stars. No colour or airmass correction was performed, both effects are only on a 1-2 % level and much smaller than the photometric errors.

## 3 Thermophysical Model Analysis

We applied the thermophysical model (TPM) developed by Lagerros (lagerros96 (1996); lagerros97 (1997); lagerros98 (1998)) and extensively used and refined by Müller & Lagerros (mueller98 (1998); mueller02 (2002)). The TPM includes shape, rotational behaviour and thermal behaviour of a given asteroid at the specific observing and illumination geometry. It allows to derive radiometric sizes and albedos from thermal infrared observations, taking shape information, geometry and thermal aspects into account.

### 3.1 Model input values

The $`H`$555Absolute magnitude in the H-G magnitude system and $`G`$666V-band slope parameter in the H-G magnitude system values are taken from Lagerkvist & Magnusson (lagerkvist90 (1990)), with $`H=6.70\pm 0.15`$ as the weighted mean value of four apparitions and $`G=0.09\pm 0.05`$. Erikson (erikson00 (2000)) presented a good quality solution pair for the spin vector based on a retrograde sense of rotation and a synodic rotation period of 6.1 hours (see Table 5). However, due to noisy and under-sampled input data, the spin vector pair is still considered only as class II, indicating that higher quality data might improve the situation.

The sense of rotation can in principle be determined from the anlysis of thermal observations before and after opposition (e.g. Müller mueller02a (2002)). In cases where two spin vector solutions with different senses of rotation are given, the thermophysical analysis allows to determine the correct one. But to distinguish between two spin vector solutions with the same sense of rotation is very difficult on basis of infrared data. The before-/after opposition effect is then not very pronounced. Our TPM analysis for 65~Cybele gave almost identical diameter and albedo values for both spin vector solutions. For the first solution (first line in Table 5) we obtained a slightly smaller standard deviation around the mean diameter and albedo values. Therefore we accepted solution 1 for all following calculations. The corresponding ellipsoidal shape model solution is closely connected to this spin vector solution.

The determination of thermophysical quantities such as thermal inertia, beaming parameters or emissivity require large samples of infrared measurements (Müller & Lagerros mueller98 (1998)), but for regolith covered, large main-belt asteroids, default thermal values seem to work fine (Müller et al. mueller99 (1999)). They found a thermal inertia value of $`\mathrm{\Gamma }=15\mathrm{J}\mathrm{m}^2\mathrm{s}^{0.5}\mathrm{K}^1`$ for 1~Ceres which they also used successfully for other main-belt asteroids (Müller & Lagerros mueller02 (2002)). We adopted this $`\mathrm{\Gamma }`$-value, as well as the derived beaming model parameters $`f=0.6`$ (fraction of surface covered by craters) and $`\rho =0.7`$ (r.m.s. of the surface slopes). These beaming values correspond roughly to the “default” STM beaming parameter (Lebofsky et al. lebofsky86 (1986)) if the heat conduction is neglected (Lagerros lagerros98 (1998)). For the emissivity we used a constant value of $`ϵ=0.9`$ at all wavelengths between the near-IR and 100 $`\mu `$m. At longer wavelengths we used a wavelength dependent emissivity with decreasing values from 0.9 at 100 $`\mu `$m to 0.8 at 400 $`\mu `$m.

### 3.2 Radiometric analysis

Based on TPM parameters from Sect. 3.1, we derived radiometric effective diameter and albedo values from all photometric and spectroscopic measurements in Sect. 2. The results are given in Table 6.

The typical observational and calibrational errors of individual measurements are between 5 and about 20 % in flux density. The translation into diameter leads to 3 to 10 % uncertainties for individual radiometric diameters. A 10 % change in radiometric diameter for 65~Cybele requires an albedo modification of 0.010 (a larger diameter requires a lower albedo and vice versa). The $`D_{eff}`$ values per instrument are within $`\pm `$8 % of the weighted mean, only the UKIRT value is significantly off for unknown reasons. The ISOPHOT data were taken under the same aspect angles as the UKIRT data, but do not show such high $`D_{eff}`$ values. The highest quality photometric data are the CAM filter measurements. The derived diameter and albedo values are very close to the weighted mean values of all 10 sets of independent observations.

The effective diameters are connected to the shape model via the equal volume sphere formula $`D_{eff}=2(abc)^{1/3}`$ (see also Müller & Lagerros mueller98 (1998)). Shape, object orientation (represented by spin vector and aspect angle777Angle between the rotation axis of 65~Cybele and the radius vector to Earth), illumination, rotation and thermal effects are taken into account in the calculation of $`D_{eff}`$ and $`p_V`$. The calculated weighted mean effective diameter of 273.0 km corresponds to an absolute ellipsoidal size of $`2a\times 2b\times 2c=`$ 301.9$`\times `$290.3$`\times `$232.2 km for the given spin vector solution (first line in Table 5).

### 3.3 Spectroscopic Analysis

The original idea of the spectroscopic observations was to get new insides into the chemical composition and mineralogy of this main-belt asteroid. Mainly the nature of silicate features can be studied in the wavelength range between 8 and 12 $`\mu `$m (Dotto et al. dotto00 (2000)). Green et al. (green85 (1985)) found no prominent emission features in the N-band spectrum of 65~Cybele.

Dotto et al. (dotto02 (2002)) and Barucci et al. (barucci02 (2002)) analysed the spectroscopic features by dividing the observed spectrum by the corresponding TPM prediction under the assumption of a perfect emissivity at all wavelength ($`ϵ=1.0`$). Here, we follow this procedure. Figure 3 shows the PHT-S and the CAM-CVF measured spectrum divided by the corresponding TPM with $`ϵ=1.0`$. Both spectra agree in spectral shape within the given observational and instrumental uncertainties.

The PHT-S spectrum (Fig. 1) was derived from 64 individual pixels which have a completely independent calibration for each pixel. The PHT-S data are therefore not ideal for an analysis of very low level spectroscopic features. The CAM-CVF scans are much better suited for such an analysis. The source is always seen by the same pixels and with the variable filters the source spectrum is scanned. Despite the disadvantages of the PHT-S spectrum for the analysis of shallow features, a 10 element smoothed PHT-S spectrum resembles the CVF spectrum very well.

## 4 Discussion

### 4.1 Physical parameters

Our reanalysis of 10 independent data sets (see Sect. 2) give now very consistent radiometric result through the TPM. We determined the effective diameter to Deff=273.0$`\pm `$12.0 and the albedo to pV=0.050$`\pm `$0.005. The discrepancy to the Tedesco et al. (tedesco02 (2002)) results is mainly related to the spherical shape simplification of the STM. Thermal effects, the beaming model, phase angle corrections and band-to-band corrections also contribute to the deviations. The occultation diameter was derived from a circle fitted to 3 measured chords. The uncertainties of the measurements and the limited information from a 3-chord occultation event also allow for a fit with a larger slightly elongated ellipse, corresponding to the projected cross section of our derived radiometric size values. The fact that different data sets, taken at different wavelengths, different observing and illumination geometries give now one consistent solution is certainly in favor of our larger diameter. And, all IRAS observations fit nicely into that picture, including the 100 $`\mu `$m data which were not used for the Tedesco et al. (tedesco02 (2002)) calculations. Here, we would like to emphasize that we did not apply any band-to-band corrections or any other wavelengths-, phase angle or beaming corrections, like it was done for the STM IRAS applications.

Despite the usage of best available ellipsoidal shape model and the nice agreement on average between the TPM and the observations, we believe that the unknowns of the true shape limit the accuracy of TPM predictions for individual epochs. Deviations between the true object cross section and the projected model cross section are easily possible, but so far no improved shape model (e.g. from radar or direct imaging methods) is available.

The shape model is only useful together with the spin vector solution. Erikson (erikson00 (2000)) gives 2 solutions as result of lightcurve inversion methods (see Table 5), but both solutions have the same sense of rotation. The determination of the correct solution on basis of radiometric analysis is therefore very difficult (see Sect. 3.1 for details). Here, we accepted solution 1 for all calculations.

### 4.2 Thermal parameters

We used a flat emissivity model with $`ϵ=0.9`$ over all wavelengths out to 100 $`\mu `$m, at longer wavelengths the model emissivity is slowly decreasing towards values of 0.8 at sub-mm wavelength. Previously published emissivity values of 1.0 at around 20 $`\mu `$m were based on ground-based Q-band observations and IRAS 25 $`\mu `$m fluxes together with very uncertain beaming values (Müller & Lagerros mueller98 (1998)). With our modification of taking a constant 0.9 value at mid-IR wavelengths we investigate the emissivity behaviour in this critical region where beaming effects play an important role (Dotto et al. dotto00 (2000)). In Fig. 4 there is no trend with wavelength visible out to about 150 $`\mu `$m. The constant ‘observation over TPM’ ratios in Figs. 1 and 2 confirm our assumption of a constant $`ϵ=0.9`$ at mid-IR wavelengths. The 5 ISOPHOT measurements at wavelength beyond 150 $`\mu `$m (calibrated against a few stars and the planets Uranus and Neptune) show fluxes which are about 10 % lower than the TPM predictions (Fig. 4). Müller & Lagerros (mueller02 (2002)) see nice agreement with this current emissivity model for 1~Ceres, 2~Pallas, 4~Vesta and 10~Hygiea out to 200 $`\mu `$m, although for 2~Pallas and 4~Vesta there were indications for a similar trend at the very longest wavelengths. The observed far-IR fluxes are to a good approximation directly proportional to the hemispherical model emissivity (Müller & Lagerros 1998). Instead of a slow decrease from 0.9 at 100 $`\mu `$m to 0.8 at 400 $`\mu `$m, the 65~Cybele observations now indicate a more abrupt decrease in emissivity to values around 0.7-0.8 in the wavelength range 150 to 200 $`\mu `$m.

Redman et al. (redman92 (1992); redman98 (1998)) found effective emissivities between 0.5 and 0.9 in the sub-mm range for 7 large asteroids. They attribute differences in the emissivity behaviour to different optical depth of the warm surface material and its density. Scattering processes by grains within the regolith reduce the emissivity in a wavelength dependent fashion. For the Moon the far-IR emission is also suppressed by scattering processes (Simpson et al. simpson81 (1981)). Scattering processes are most effective if wavelength and grain sizes have similar dimensions. Following this interpretation, our far-IR low emissivities indicate that the dominant particle size in the 65~Cybele’s regolith (and maybe also for 2~Pallas and 4~Vesta) would be around 200 $`\mu `$m. More observations are needed to identify the precise emissivity behaviour. Well-calibrated spectroscopic observations in the far-IR and sub-millimetre (e.g. with HERSCHEL) will then allow a detailed description of particle size distributions within the regolith.

The measurements before and after opposition come from different instruments and are too inhomogeneous to identify clear trends with phase angle, especially as the measurements before opposition (negative phase angle) have too large error bars. We varied the thermal inertia between 0 and 50 $`\mathrm{J}\mathrm{m}^2\mathrm{s}^{0.5}\mathrm{K}^1`$ to study the systematic changes in the data representation of Fig. 5. The zero thermal inertia produces a significantly larger scatter in the observation-to-model ratios. The $`\mathrm{\Gamma }=50`$ calculations showed a very pronounced asymmetry between before and after opposition ratios. Our default value of 15 $`\mathrm{J}\mathrm{m}^2\mathrm{s}^{0.5}\mathrm{K}^1`$ seems to work fine, but somewhat smaller or larger values cannot be excluded. The much larger scatter in the before opposition ratios is pure coincidence. E.g. the before opposition IRAS data cover a much larger range in flux and the most deviating IRAS ratios on the left side of Fig. 5 belong to the lowest measured fluxes.

The beaming model with $`\rho =0.7`$ and $`f=0.6`$ influences strongly the slope of the observation to model ratios in the mid-IR (Dotto et al. dotto00 (2000)). The flat curves in Figs. 1 and 2 confirm our “default beaming model”. About 10-20 % modifications of the r.m.s. slope value ($`\rho `$) and the fraction of surface covered by craters (f) would still be in agreement with the PHT-S and the CAM filter measurements. Müller (mueller02a (2002)) demonstrated in different model simulations the influences of the thermal parameters with wavelength and phase angles. But thermal inertia, beaming effects and the low level mineralogic emission features play together and influence the slope of the thermal emission in a complicated way and the different contributions are not easy to disentangle.

### 4.3 Emissivity analysis

Barucci et al. (barucci87 (1987)) classified 65~Cybele as C0-type asteroid on basis of visual multi-color and IRAS data. They associated this taxonomic type with carbonaceous chondrites. The asteroids 10~Hygiea and 511~Davida, also C0-types, were found to show similarities in near-IR and 3 $`\mu `$m spectroscopy with carbonaceous chondrite meteorites of CI/CM subclass (Hiroi et al. hiroi96 (1996)). Therefore, we concentrated our search for meteoritic analogues to carbonaceous chondrites, which represent, chemically and mineralogically, the most primitive material in the solar system.

In Fig. 6 we show the average CVF scan (divided by a TPM prediction with $`ϵ=1.0`$) with the error bars derived from the two CVF scans ($`\sigma _{mean}/\sqrt{2}`$). For comparison, we plotted the emissivity of different powdered stony meteorite samples, representing 17 carbonaceous chondrite meteorites at small particle sizes of 0–75 $`\mu `$m (Salibury et al. 1991a ; 1991b , ASTER spectral library888http://speclib.jpl.nasa.gov). This library contains 60 meteorite samples all at grain sizes of 0-75 $`\mu `$m. The sample includes five meteorites of subclass CI, five of subclass CM, five of subclass CO and two of subclass CO3. The emissivity spectra (1 – measured reflectivity) are offset for clarity.

The 65~Cybele emissivity shows no prominent features, but an increase in emissivity between 8.0 and 9.5 $`\mu `$m is seen in both down and upward CVF scans and is thus a reliable feature, contrary to the part below 7.0 $`\mu `$m where the source becomes weak and transient problems hamper the measurement (see Sect. 2.3). Between 12.5 and 14 $`\mu `$m we recognise a very shallow feature which is also present in both scans.

Salibury et al. (1991a ) noted that a rather steep increase in emittance, the so-called Christiansen feature (CF), occurs between the wavelength region dominated by volume scattering (roughly below 7.5 $`\mu `$m) and the wavelength region dominated by surface scattering and reststrahlen bands (roughly beyond 8.5 $`\mu `$m). They placed the CF wavelengths of different types of carbonaceous chondrites within the region between about 8.7 and 9.4 $`\mu `$m. The CO and CO3 types of the carbonaceous (Fig. 6) come closest in CF similarity with the 65~Cybele emittance peak between 9.0 and 9.4 $`\mu `$m, but the data quality is not sufficient for a conclusive identification. The two CO3 stony meteorites in Fig. 6 are Isna (upper curve) and ALHA 77003 (lower curve), the CO meteorites are Lance (2$`\times `$), ALH 83108, Ornans and Warrenton (Salibury et al. 1991a ). The CI and CM subclasses show the CF at slightly shorter wavelengths.

Our results from the emissivity analysis (Sect. 4.2) pointed at predominantly larger particles as cause for the far-IR emissivity drop. But the Salibury et al. (1991a ) library provides only spectra of small grain sizes. Le Bras & Erard lebras03 (2003) found that the grain size effect can shift the CF in wavelengths and change also the contrast. Meteoritic samples with larger grain sizes are therefore needed to proove direct connections between the measured spectroscopic features on asteroids and meteoritic samples.

For the feature between 12.5 and 14 $`\mu `$m we found no counterpart yet, but phyllosilicates, like Kaolinite, Lepidolite or Serpentine, show emissivity increases at around 12.5 $`\mu `$m, similar to what we see here. The fine-grained phyllosilicates are one main component of CI/CM carbonaceous chondrites, while the CO types are anhydrous materials. This would be an argument in favour of similarities with CI or CM types, but again, laboratory studies at larger grain sizes are needed to clarify the situation.

Figure 6 clearly reveals the difficulties of identifying mineralogic features based on mid-IR spectroscopy. Laboratory spectra of different meteorites at various grain sizes would be necessary. Concerning the observations of asteroids, a high S/N ratio over a large wavelengths range is necessary to detect the low level emissivity changes. These features would be very difficult to detect from ground-based N-band spectroscopy (8 to 13 $`\mu `$m) where atmospheric effects limit the wavelength range and degrade the data quality. Additional complications arise from the strong slope changes around the thermal continuum energy peak at around 14 $`\mu `$m. Poor spectral information together with wrongly modelled continuum slopes can then easily lead to misinterpretation at the band borders. The modelling of the surface temperature distribution, the thermal inertia as well as the beaming parameters $`\rho `$ and $`f`$ influence the slope of the thermal emission. Our thermal values (Sec. 3.1) worked perfectly for the PHT-S and the CAM-CVF data sets, but the spectroscopic interpretation still remains difficult.

### 4.4 Comparison with the STM

The shape and spin vector information is very important for a consistent analysis of the combined observations. Although we only used a simple ellipsoidal shape model, one can see the shape influence by comparing Fig. 7 with Fig. 4 in Sect. 4.1. For Fig. 7 we calculated the model values via the STM (Lebofsky et al. lebofsky86 (1986)), using a spherical shape model. The scatter between data sets taken under different aspect angles is much larger. Some measurements which agree in Fig. 4 are now offset by large amounts (e.g. the IRAS 12 and 25 $`\mu `$m observations and the CAM filter measurements). Assuming a spherical shape (with the given spin vector and thermal behaviour) leads to a weighted mean diameter of 280.0 km and an albedo of 0.047. But now the standard deviations would increase by about 60 % for the diameter and about 50 % for the albedo! This result clearly demonstrates the importance of using spin vector and shape information for combining observations from different epochs. Modifying the STM beaming value might improve the “Observation/STM” ratios for some observations, but different observations would require then different STM beaming values. And, for some data sets the STM would require a wavelength dependent $`\eta `$ to account for the trends in Fig. 7.

We used the most reliable data set –the CAM filter measurements– for a comparison between TPM and STM radiometric results. The TPM analysis leads to a constant diameter and albedo value for all CAM filter measurements, while the STM produces wavelength dependent values (Fig. 8). The albedo values show the opposite behaviour. The derived diameters show no trend with wavelengths (see also Fig. 2, lower part for a comparison of the observation to TPM ratios over wavelengths) and agree within 2-3 % (about 6 % for the albedo). This also confirms our simplified constant emissivity values in this wavelength range close to the thermal emission peak. Simulations of the 65~Cybele constellation during the CAM filter epoch show that this wavelength dependent trend continues out to about 30 $`\mu `$m. At longer wavelengths also the STM would produce an almost constant diameter. At 30 $`\mu `$m the STM diameter would be about 10 % larger. This wavelength dependent effect of the STM is mainly produced by the STM correction factor for infrared beaming. Detailed studies of the beaming influence with wavelengths and phase angles can be found in Müller (mueller02a (2002)).

Figures 4, 5 and 7 reveal additional problems of the STM. The scatter in the observed fluxes divided by model calculations is much larger for the STM (Fig. 7). The main reason for the larger scatter is the assumption of a spherical shape. But also the simplified beaming, emissivity and phase angle corrections can be recognized in Figs. 4 and 7. The residual beaming effects cause wavelengths dependent ratios in Fig. 7, easily visible for the IRAS data ($``$-symbol) or the CAM filter data ($`\mathrm{}`$-symbol). The offset between these 2 sets is caused by the wrong cross section due to the spherical shape model of STM. The STM calculation of the ratios of Fig. 5 shows a strong asymmetry between the ratios before and after opposition. Additionally, the scatter between the different data sets would be much larger. Unfortunately, the overlaying shape effects of the inhomogeneous data did not allow to compare the 0.01 mag/ phase angle correction factor of the STM with the more sophisticated calculations of the true illumination and observing geometry by the TPM.

## 5 Conclusions

The recently developed TPM allows to combine thermal infrared observations with information from visual photometry, lightcurves, close-up observations or direct measurements. Using 65~Cybele as a “text book example”, we applied this powerful TPM in several different ways. First, we used re-calibrated observations from literature and recently obtained thermal infrared measurements to derive the size and albedo of 65~Cybele, taking also shape and spin vector information into account. Our derived effective diameter (273.0$`\pm `$11.9 km) and albedo values (0.050$`\pm `$0.005) deviate significantly from the IRAS based results of 237.3$`\pm `$4.2 km and 0.071$`\pm `$0.003 (Tedesco et al. tedesco02 (2002)). We attribute these differences to the limitations of the STM (Sects. 4.1, 4.2 and 4.4).

The spectroscopic analysis confirmed the proposed “default thermal behaviour” (Müller et al. mueller99 (1999)) for large main-belt asteroids which is dominated by a surface regolith with very low thermal inertia. The related beaming model and thermal inertia values are the key elements to match even the high accuracy ISO measurements at the critical mid-IR range. In return, reliable diameter and albedo values can be determined from individual filter measurements at any thermal infrared wavelength. No band-to-band or wavelength dependent correction factors are necessary anymore. Additionally, the thermal emission modelling provides the basis for the analysis of mineralogic features and allows to investigate relations to meteoritic samples.

###### Acknowledgements.

We would like to thank U. Klaas and P. Abraham for their support in the PHT-S data analysis and D. Osip for the provision of the IRTF observations. R. Siebenmorgen and M. Sterzik supported the ESO observing run and the data analysis of the TIMMI2 measurements. The ISOPHOT data presented in this paper were reduced using PIA, which is a joint development by the ESA Astrophysics Division and the ISOPHOT Consortium with the collaboration of the Infrared Processing and Analysis Center (IPAC). Contributing ISOPHOT Consortium institutes are DIAS, RAL, AIP, MPIK, and MPIA. The ISOCAM data presented in this paper were analysed using ‘CIA’, a joint development by the ESA Astrophysics Division and the ISOCAM Consortium. The ISOCAM Consortium is led by the ISOCAM PI, C. Cesarsky. We also thank the referee Dr. Delbo for his careful review.


format: text

text:
# Hard Disks in Narrow Channels

## I Introduction

The thermodynamic behavior of hard-disk systems has been a subject of considerable interest in recent years . In particular, the freezing transition taking place in these systems has been analysed by computer simulations for cases where the aspect ratio of the simulation box is close to one. Furthermore, the dynamical stability of such systems has been studied in detail and Lyapunov spectra have been computed. The maximum Lyapunov exponent was found to exhibit a maximum in the density regime characteristic for the fluid-to-solid transition. In a restricted geometry, such as a narrow channel, one expects other interesting features to emerge when the width of the channel becomes comparable to the disk’s size. Studying narrow channels also enables one to examine long-wavelength perturbed states associated with the small positive Lyapunpov exponents found in the tangent-space dynamics of hard-disk and hard-sphere systems . This provided the motivation for us to re-investigate the thermodynamic and dynamical properties of a hard-disk gas in a narrow channel.

It is found that for narrow channels the pressure exhibits a singularity at a channel width equal to twice the diameter of the disks. Narrower channels do not allow particles to pass each other. At the same width a singularity is also observed for the collision frequency. In addition, the maximum Lyapunov exponent and the Kolmogorov-Sinai entropy display a non-trivial dependence on the channel width in this regime.

The model we study consists of $`N`$ hard disks in a two-dimensional box with side lengths $`L_x,L_y`$ and aspect ratio $`A=L_y/L_x`$. We vary the width of the box, keeping the volume $`V=L_xL_y`$ and the particle density $`\rho =N/V`$ constant. Very small aspect ratios $`A1`$ are considered, for which the box resembles a narrow channel. For most of our work we use periodic boundary conditions in both directions. Exact thermodynamic properties are known for extremely narrow channels in this case $`L_y<\sqrt{3}\sigma `$ . Recently, the transport properties of two particles in a square periodic box, $`A=1`$, have been studied by computer simulation . The case of two particles with reflecting boundary conditions in both directions has also received attention, both by molecular dynamics  and analytical approaches . It was found that for certain aspect ratios the compressibility becomes negative. We find analogous results also for partly-reflecting boundaries, which are reflecting for the walls parallel to the long direction of the box, $`x`$, and periodic for the walls parallel to $`y`$.

Throughout, reduced units are used for which the disk diameter $`\sigma `$, the particle mass $`m`$ and the kinetic energy per particle, $`K/N`$, are unity. There is no potential energy in this case, and the total energy $`E`$ is identical to the kinetic energy $`K`$. As a consequence, the temperature $`T`$ is just an irrelevant parameter and is fixed to unity. Also Boltzmann’s constant $`k`$ is taken unity.

The paper is organized as follows. In Section II we consider the low-density limit. In Section II A we present an analytical derivation of the pressure and show that the singularity taking place at $`L_y=2`$ already exists in this limit. In addition, we obtain an analytical expression for the collision frequency which shows a behavior very similar to that of the pressure. Results of numerical simulations of these quantities and of the Lyapunov spectra are presented in Section II B. A numerical study of the high-density case is presented in Section III. Full Lyapunov spectra for periodic systems with periodic boundaries are discussed in Section IV. We close with a brief discussion of the results in Section V.

## II Low-density limit for systems with periodic boundaries

### A Virial approach to the equation of state and collision frequency

The pressure in the low-density limit is obtained by the leading terms in the virial expansion as applied to the narrow-channel geometry. Let $`q(L_y)`$ be the excluded volume of a single disk in the channel. The phase-space volume associated with placing $`N`$ disks in the channel is given by

$$\mathrm{\Omega }(N)=V(Vq)\mathrm{}(V(N1)q),$$
(1)
which, to leading order in the density, is reduced to

$$\mathrm{\Omega }(N)V^N\left(1\frac{N(N1)}{2}\frac{q}{V}\right)V^N\left(1\frac{q}{2v}\right)^{N1}.$$
(2)
Here, $`v=V/N`$ is the specific volume. For low density the specific entropy of the gas is given by

$$s\frac{S}{N}=k\mathrm{ln}\left(v\frac{q(L_y)}{2}\right),$$
(3)
where $`k`$ is Boltzmann’s constant.

For $`L_y>2`$ the excluded area is simply given by the area of a disk of radius $`1`$, namely $`q(L_y>2)=\pi `$. To leading order in the density, the pressure $`P`$ is given by

$$\frac{Pv}{kT}=v\frac{s}{v}1+\frac{\pi }{2v},$$
(4)
and is thus independent of the channel width $`L_y`$. This is the expected leading-order expression of the virial expansion of a gas of hard disks.

For $`1<L_y<2`$, however, the excluded area is reduced due to the interference with the disk image resulting from the periodic boundary condition in the $`y`$ direction (see Fig. 1). Simple geometrical considerations yield

$$q(L_y)=2\theta +\mathrm{sin}(2\theta ),$$
(5)
where

$$L_y=2\mathrm{sin}(\theta ).$$
(6)
The excluded volume $`q(L_y)`$ depends on $`L_y`$. Taking the derivative of the entropy with respect to the volume and keeping the aspect ratio $`AL_y/L_x`$ constant, one finds to leading order in the density

$$\frac{Pv}{kT}=1+\frac{q(L_y)}{2v}\frac{L_y}{4v}\sqrt{4L_y^2}.$$
(7)
Note that the pressure curve exhibits a square root singularity as one approaches $`L_y=2`$ from below.

The singularity at $`L_y=2`$ does not signify a phase transition in the usual sense as it is not a collective phenomenon. In fact, the singularity exists even for finite $`N`$ and it explicitly shows up in the second term of the virial expansion. This is different from the usual liquid-solid transition in the bulk, which is obtained only after summing over all orders of that expansion. The singularity in the narrow channel is a consequence of the fact that the available volume for the disks is a singular function of the control parameter $`L_y`$ at $`L_y=2`$. Higher order terms in the virial expansion are expected to exhibit a singularity at other integer channel widths, $`L_y=3,4\mathrm{}`$ . Thus singularities in the pressure are expected to show up at these widths at higher densities.

We now turn to the evaluation of the collision frequency $`\nu _2`$. The virial expression for the pressure, appropriately modified for two-dimensional hard disks, is given by

$$\frac{PV}{NkT}=1+\frac{1}{2NkT\tau }\underset{c}{}𝐫_c𝐯_c,$$
(8)
where the sum is over all collisions taking place during the time interval $`\tau `$. Here, the relative position vector of two colliding particles $`i`$ and $`j`$ at the time of the collision, $`𝐫_c𝐫_i𝐫_j`$, satisfies $`|𝐫_c|=\sigma `$, and $`𝐯_c𝐯_i𝐯_j`$ is the corresponding relative pre-collision velocity. This may be written as

$$\frac{PV}{NkT}1=\frac{1}{2kT}\nu _2\mathrm{\Delta }p$$
(9)
where $`\mathrm{\Delta }p`$ is the average momentum transfer in a collision. For a large aspect ratio $`A`$ the collision frequency is given by

$$\nu _2=2\pi ^{1/2}(kT)^{1/2}(\sigma /v)g(\sigma ),$$
(10)
where $`g(r)`$ is the pair distribution function. At contact, $`r=\sigma `$, and in the low-density limit one has $`g(\sigma )=1`$. Combining this result with the Eqs. (4) and (9), one finds for the average momentum transfer

$$\mathrm{\Delta }p=(\sqrt{\pi }/2)\sqrt{kT}.$$
(11)
This expression may also be derived directly from the Maxwell-Boltzmann velocity distribution. It is expected to hold also for narrow channels. The reason is that in the narrow channel the velocity distribution is still isotropic and is given by the Maxwell-Boltzmann distribution (see Section II B). Therefore, one can use Eq. (11), together with Eq. (9), to obtain the collision frequency for narrow channels:

$$\frac{\nu _2}{\sqrt{kT}}=\frac{4}{\sqrt{\pi }}\left(\frac{Pv}{kT}1\right)$$
(12)
This expression is valid both below and above $`L_y=2`$.

It is expected that in such a highly anisotropic system the pressure tensor is not isotropic. In the low-density limit one can evaluate the diagonal components of the pressure tensor, $`P_{xx}`$ and $`P_{yy}`$, using again the expression for the entropy in Eq. (3). One finds

$`P_{xx}`$ $`=`$ $`L_x\left({\displaystyle \frac{s}{L_x}}\right)_{L_y}={\displaystyle \frac{v}{v(1/2)q(L_y)}},`$ (13)
$`P_{yy}`$ $`=`$ $`L_y\left({\displaystyle \frac{s}{L_y}}\right)_{L_x}={\displaystyle \frac{v(1/2)L_y\sqrt{4L_y^2}}{v(1/2)q}}.`$ (14)
Furthermore, one has

$$v\left(\frac{s}{v}\right)_A=\frac{1}{2}L_x\left(\frac{s}{L_x}\right)_{L_y}+\frac{1}{2}L_y\left(\frac{s}{L_y}\right)_{L_x},$$

and, hence, $`P=(P_{xx}+P_{yy})/2`$.

It is evident that in the low-density limit the pressure tensor is isotropic for $`L_y>2`$. It becomes anisotropic for $`L_y<2`$. At higher densities the two pressure-tensor components differ even above $`L_y=2`$ (see Section III). The square-root singularity of the total pressure $`P`$ at $`L_y=2`$ originates from $`P_{yy}`$. The other component, $`P_{xx}`$, exhibits a weaker singularity.

### B Numerical simulations

Here we carry out numerical simulations for the computaton of the pressure, the collision frequency and the Lyapunov exponents. The pressure is evaluated from the impulsive version of the virial theorem Eq. (8) . The collision frequency of a single particle, $`\nu _2`$, is obtained from the total number of collisions per unit time, divided by $`N/2`$. For the Lyapunov spectrum we use the algorithm outlined in Refs.  and . For the periodic-boundary case the linear momentum is conserved both in $`x`$ and $`y`$ directions, and, as expected, altogether 6 Lyapunov exponents vanish due to the conservation of energy, momentum, center of mass (only in tangent space), and the regularity of the dynamics in the phase-flow direction. For the partly-reflecting boundary case only four of the exponents vanish due to a lack of momentum and center-of-mass conservation in the $`y`$ direction. The Kolmogorov-Sinai (KS) entropy, $`h_{KS}`$, which measures the rate of information gain on the initial conditions due to the time-reversible dynamics, is given by the sum of the positive exponents. Since its dependence on the channel width is qualitatively similar to that of the maximum Lyapunov exponent $`\lambda _1`$ we concentrate on the latter in the following.

In Figure 2 we compare the theoretical curves for the pressure and collision frequency, Eqs. (4), (7) and (12), with the numerical results. In this low-density regime an excellent agreement is obtained. Also shown in this figures is the channel-width dependence of the maximum Lyapunov exponent. It shows a shoulder at the critical width $`L_y=2`$.

The tangent-space perturbation associated with $`\lambda _1`$ has been shown to be localized in space . This exponent is closely connected to the collision frequency and, as a consequence, has a similar qualitative behavior. This is particularly pronounced in the high density regime as dealt with in the next section.

In Fig. 3 we compare the theoretical curves for the pressure-tensor components $`P_{xx}`$ and $`P_{yy}`$ with computer simulations for a density $`\rho =0.01`$ and find excellent agreement. As noted above, it is evident that the square-root singularity of $`P`$ originates from $`P_{yy}`$.

It is interesting to note that in the narrow-channel regime, $`L_y<2`$, the velocity distribution is still isotropic and is given by the Maxwell-Boltzmann distribution as is shown in Fig. 4. This supports the use of Eq. (11) for the average momentum transfer at a collision even for $`L_y<2`$.

The singular behavior of the pressure is due to a sharp increase of the collision frequency as one approaches $`L_y=2`$ from below. As can be seen in Fig. 2 the maximum Lyapunov exponent is reminiscent of this singularity as it is closely related to the collision frequency. At this low density this singularity is not clearly pronounced. However as will be shown below, at high densities the existence of a cusp at $`L_y=2`$ is evident. The increase of the collision frequency is due to arresting configurations in which a pair of disks become trapped due to the boundaries.

## III High density regime

For hard-disk systems at high density one has to rely only on numerical simulations. In Fig. 5 we present the pressure, the collision frquency and the maximum Lyapunov exponent as a function of $`L_y`$ at a density $`\rho =0.4`$.

It is observed that the singularity taking place at $`L_y=2`$ persists at high density in all three curves. According to Eq. (12) the curves for the collision freuqency and for $`Pv/T1`$ are proportional with a propertionailty constant independent of the density. This is evidently confirmed in Fig. 5. For the special case of only $`N=2`$ particles in a square periodic box, the singularity at $`L_y=2`$ appears at the critical density $`\rho _c=0.5`$ and has recently been observed .

When the density is increased beyond 0.5, a broad peak emerges at some $`L_y<2`$ in the curves for the pressure and the maximum Lyapunov exponent. This is demonstrated in Fig. 6. This feature is attributed to the

pronounced short-range crystalline-like order which is induced by the narrow-channel boundary conditions. Note that this order is not the natural triangular lattice of the system in the bulk. A typical microscopic configuration corresponding to this structure for a density 0.9 and $`L_y=1.5`$ is shown in Fig. 7. The particles are arrested and cannot travel across the system, which has been referred to as the localized regime .

We have computed also the two pressure-tensor components $`P_{xx}`$ and $`P_{yy}`$ for a system with a density $`\rho =0.8`$. As can be seen in Fig. 8, the two components are significantly different both above and below $`Ł_y=2`$. The singularity of the pressure at $`L_y=2`$ is evidently caused by $`P_{yy}`$, as is also the case for low densities.

All the examples so far are for the case of periodic boundaries in both the $`x`$ and $`y`$ directions. A qualitatively similar result is obtained, if the long sides of the channel are elastically reflective, whereas the short sides remain periodic. This is demonstrated in Fig. 9 for $`N=100`$ particles at a density of $`N/V=0.4`$, where $`P`$, $`\lambda _1`$ and $`\nu `$ are shown as a function of the box width $`L_y`$. The maxima are not as sharp as in Fig. 5 for the periodic case, but are as pronounced.

In order to get some insight into the dynamics of such systems it is usually useful to study the diffusion coefficient. It is well known that this coefficient does not exist in two-dimensional systems in the thermodynamic limit . However, it is possible to compute the mean-squared displacement for a finite system and to extract effective diffusion constants $`D_x`$ and $`D_y`$ in $`x`$ and $`y`$ directions, respectively. They are obtained from fits to the

linear growth of $`\mathrm{\Delta }x^2`$ and $`\mathrm{\Delta }y^2`$ for times $`t>40`$. The results for a density $`\rho =0.8`$ are shown in Figure 10. As expected $`D_x`$ vanishes for $`L_y<2`$, while $`D_y`$ is non-zero. In fact, $`D_y`$ increases with $`L_y`$ approaching unity as a result of the periodic-boundary restrictions disfavoring momentum exchange in the $`y`$ direction. Another intersting feature of Fig. 10 is the vanishing of $`D_x`$ and $`D_y`$ for $`L_y>2.7`$. This is due solid-like triangular configurations characteristic of the bulk at this high density.

## IV Some remarks on Lyapunov spectra

In Fig. 11 we plot the Lyapunov spectra for a 100-disk system in a narrow periodic box. All spectra are for a density $`\rho =0.4`$, which corresponds to a fluid. Since the volume, $`L_xL_y`$, is fixed, an increase of the channel width $`L_y`$, varied in the figure in steps of 0.4, decreases $`L_x`$ accordingly. $`L_y`$ is indicated by the labels. Each spectrum consists of $`4N`$ Lyapunov exponents, of which only the positive branch, $`\{\lambda _l0;l=1,2,\mathrm{},2N\}`$, is shown. The spectra are only defined for integer values of the index, $`l`$, which labels the exponents according to size. The lines are only drawn for clarity. According to the equilibrium formulation of the conjugate-pairing rule, the negative branch of the spectra, $`\{\lambda _l0;l=2N+1,\mathrm{},4N\}`$, is the mirror image of the positive branch .

The inset of Fig. 11 provides a magnified view of the small exponents. They turn out to be degenerate, which gives rise to a step-like appearance of the spectra . The tangent-space perturbations associated with these exponents are collective, wave-like fields defined over the simulation box, and are referred to as Lyapunov modes . The multiplicity is determined by the intrinsic symmetries of the Hamiltonian and the boundary conditions, which give rise to the conserved quantities, energy and momentum. A complete classification in terms of transverse (T), longitudinal (L), and momentum (P) modes is given in Ref. . Theoretical attempts have been made to interpret the modes in terms of fluctuating hydrodynamics .

We observe that a particular exponent, say the smallest positive, increases with $`L_y`$ and, hence, with the wave number, $`k=2\pi /L_x=2\pi L_y/V`$, belonging to this mode. This is referred to as a “dispersion relation” . Thus, $`\lambda _{2N3}`$ increases monotonously with $`L_y`$. However, as Fig. 11 shows, this proportionality does not hold for the large exponents. There is a crossover of the spectra reversing the sequence of $`\lambda `$ near $`L_y=2`$. This is mainly a consequence of the collision frequency. The perturbations for the large exponents are found to be localized in space with only a small fraction of the particles contributing to the large exponents at any instant of time. The active zone moves around in space, restoring homogeneity on average.

The smooth lines in Fig. 12 show Lyapunov spectra for analogous 100-disk systems for a density $`\rho =0.9`$, where the channel width is indicated by the labels. The case $`L_y=1.5`$ corresponds to the arrested structure of Fig. 7 and is of particular interest. In Fig. 12 we compare it with the dash-dotted spectrum of $`100`$ disks in a periodic channel with an aspect ratio $`L_y/L_x=\sqrt{3}/(N/2)`$. For $`L_y=1.8258`$ and $`N=100`$, the density is $`\rho =1.039`$, and the particles form a triangular lattice characteristic of a solid in the bulk. Although the nearest neighbor separations and, hence, the collision frequencies, are nearly the same for both systems, the Lyapunov spectra are distinctly different. The maximum exponent, $`\lambda _1`$, for the arrested structure is significantly larger than for the triangular system, whereas the behavior of the smaller exponents is just the reverse. This example demonstrates that the Lyapunov spectra are sensitive to structural details and are potentially a useful tool. But it is fair to say that no theory exists at present to interpret such spectra in full detail.

## V Discussion

In this work the thermodynamic and dynamical behaviour of hard disks in a narrow channel is analysed using analytical and numerical approaches. The main results are for periodic boundary conditions in both $`x`$ and $`y`$ directions. Related studies of reflecting boundaries parallel to the channel axis reveal that similar features exist there as well. It is found that the pressure, collision frequency, maximum Lyapunov exponent, and Kolmogorov-Sinai entropy curves exhibit a singularity for a channel width equal to twice the disk diameter. For low densities this singularity for $`P`$ and $`\nu _2`$ is well understood by a virial expansion as it already shows up in the second virial coefficient. The singularity is not the result of a collective behavior and, in fact, is present for systems with a finite number of particles. Thus, it is not a genuine phase transition. It is rather a result of a rapid change in the available phase space taking place close to $`L_y=2`$. Similar but less pronounced singularities are expected to show up at $`L_y=3`$ and larger integers at higher virial coefficients.

The singularity for the maximum exponent and for the KS-entropy is most pronounced for higher densities. This is a consequence of the fact that these quantities are closely related to the collision frequency . It has been established that in bulk systems these quantities exhibit a maximum at a phase transition . The maximum Lyapunov exponent is a measure of the fastest dynamical events taking place in a system which are localized processes in space . Thus, its behavior near the transition point is not related to a diverging length scale at the transition. Whereas in an ordinary phase transition in the bulk the enhanced collision frequency - and, thus, of $`\lambda _1`$ \- is due to the emergence of a new structure, here it is a consequence of the constraints imposed by the boundaries.

We expect the singularity found for $`L_y=2`$ to be present also for soft disks. In this case, however, the singularity is likely to be weaker than the square-root singularity found in the pressure curve of hard disks. Similar features found in the present study should be present in three-dimensional channels as well. It would be of interest to study these problems in more detail.

###### Acknowledgements.

We thank the Austrian Science Foundation (FWF), Grant P15348-PHY, and the Israel Science Foundation for support. We also acknowledge the Einstein Center and the Österreichische Gesellschaft der Freunde des Weizmann Institute of Science for support during mutual visits.