text:
# Small Scale Structure at High Redshift: IV. Low Ionization Gas Intersecting Three Lines of Sight to Q2237+0305 The observations were made at the W.M. Keck Observatory which is operated as a scientific partnership between the California Institute of Technology and the University of California; it was made possible by the generous support of the W.M. Keck Foundation.

## 1 Introduction

Absorption lines in the spectra of the separate images of gravitationally lensed QSOs provide a unique opportunity to investigate the detailed structure of the interstellar gas and the intergalactic medium on scales of ten parsecs to a few tens of kpc. The ultimate goal of such observations is to study the evolution of the microstructure of the gas in galaxies and in the intergalactic medium from the epoch of formation of the first luminous objects in the universe down to the present epoch. For the past few years we have undertaken a survey of lensed QSOs with the HIRES spectrograph (Vogt et al 1994) on the Keck I telescope on Mauna Kea. The first paper in the present series (Rauch, Sargent and Barlow 1999) analyzed a low ionization absorption system at $`z_{abs}=3.538`$ on a scale of only 26h$`{}_{50}{}^{}{}_{}{}^{1}`$ parsecs. The second paper (Rauch, Sargent and Barlow 2001) investigated the structure of C IV absorbing clouds, including an estimate of the rate of input of mechanical energy, while the third paper (Rauch, Sargent, Barlow and Carswell 2001) set a limit on inhomogeneities in the “Lyman alpha forest” clouds at redshift $`z3.3`$.

In the present paper we turn our attention to low ionization systems which are characterized by MgII or FeII absorption lines. This gas is likely to be denser, more metal-rich, and more closely associated with galaxies than either the average Lyman $`\alpha `$ forest or the high ionization CIV absorption systems.

At intermediate redshifts MgII absorbers have usually been interpreted as the large gaseous halos of luminous galaxies. The sizes of entire MgII complexes are relatively well constrained: Smette et al (1995), from observations of HE1104-1805A,B, estimate lower limits of $`>20h_{50}^1`$ kpc for the diameters of MgII systems (see also the high resolution study of the z=1.66 system in that QSO by Lopez et al 1999). Monier et al (1998) find an upper limit of $`20h_{50}^1`$ kpc for a strong Lyman limit absorber. HIRES spectra we have obtained of the two lines of sight to UM673 A and B (to be published) show three strong MgII absorption systems. The one at z=0.426 (transverse separation between the lines of sight: $`14.6h_{50}^1`$ kpc) shows absorption in only one of the images, whereas in the other two systems (z=0.492, sep. = $`15.5h_{50}^1`$ kpc; z=0.564, sep. = $`13.9h_{50}^1`$ kpc) there is MgII absorption of similar strength in both sight lines. Surveys for galaxies in the vicinity of MgII absorption systems and attempts to reconcile the total absorption cross section with the luminosity function of galaxies (Bergeron & Boissé 1991; Steidel 1995; Bowen, Blades & Pettini 1995; Churchill et al 1999) have shown that galaxies appear to be surrounded by halos of MgII gas with radii estimated to be on the order of 60-130 $`h_{50}^1`$ kpc.

The spatial extent of individual MgII absorption components as opposed to entire MgII systems is somewhat harder to pin down. The low ionization system discussed by Rauch, Sargent and Barlow (1999) showed that SiII and CII gas cloudlets (at least at z=3.6) can be as small as a few tens of parsecs. Petitjean et al (2000) have investigated intermediate redshift MgII systems in a (spatially unresolved) HIRES spectrum of the lensed BAL APM08279+5255 and conclude from the residual intensity in the MgII doublets that the sizes of individual clouds are less than $`1.5h_{50}^1`$ kpc. The only other independent size estimate on individual cloudlets comes from ionization calculations (e.g., Bergeron & Stasinska 1986; Steidel 1990), which, however, provide more indirect constraints. Rigby et al (2002) using such arguments, find that a subset of iron-rich MgII systems may be consistent with cloudlets as small as 10 pc in diameter.

In the present paper we address the question of the small scale structure of low and intermediate ionization gas in more detail. In particular, we investigate the absorption lines in the spectrum of the ”Einstein Cross” in which Q2237+0305 ($`z_{em}=1.69`$) is gravitationally lensed by a 15th magnitude galaxy with a redshift z = 0.0390 (Huchra et al. 1985). Subsequent work (Yee 1988) revealed that the galaxy is a barred Sb with a ring, and has a disk scale length $`r_0=9.4\mathrm{"}`$ and an effective radius for the bulge $`r_e=1.9\mathrm{"}`$. Q2237+0305 turned out to have four images roughly arranged in the form of a cross centered on the nucleus of the galaxy (Yee 1988; see his fig.4). The typical separation between images is less than 2”. Crane et al. (1991) obtained B = 17.96 magn. for component A, B = 17.82 magn. for B, 18.66 magn. for C and 18.98 magn. for image D. Later observations showed that the relative brightnesses of the images is subject to changes, in part due to microlensing (Irwin et al 1989). Note that the flux ratios obtained during the current observation (in the caption of Fig. 1) are very different indeed.

The original discovery spectra only showed self absorption in the Lyman alpha and C IV emission lines. Hintzen et al. (1990) searched for but failed to find Ca II H and K absorption due to the lensing galaxy in the spectra of the quasar images. However, they detected Mg II and Fe II absorption at $`z_{abs}=0.5664`$ and a C IV complex at $`z_{abs}=1.694`$ and 1.697.

In section 2 we describe the high resolution observations of the three brighter components of Q2237+0305 with HIRES. Section 3 presents the detection of CaII absorption from the Milky Way and Section 4 describes a CaII system arising in the lensing galaxy. Section 5 discusses the properties of a peculiar, strong MgII absorption and tries to give a tentative interpretation. Sections 6 and 7 briefly describe the other two low ionization systems towards Q2237+0305, and section 8 gives the conclusions.

## 2 Observations and data analysis

Q2237+0305 has four images roughly arranged in the form of a cross centered on the nucleus of a 15th magnitude galaxy. A picture of the lensed QSO taken off the Keck HIRES guider is shown in Fig. 1 where the four images are identified. Spectra of the three brighter images A, B and C were obtained with HIRES in October 1998. The blue cross disperser was used, giving a wavelength coverage of 3645-5204 Å which puts the C III\] 1909 emission line at the red end of the range covered. We used a 0.86” slit, giving a resolution of 6.6 km s-1. The exposures, broken into 3000 second sections, totaled 18,000 seconds on image A, 18,000 s on B and 15,000 s on C. A heliocentric correction was applied. All wavelengths quoted here are vacuum values. The individual spectra of each image were combined and scaled to the same continuum level in order to compare the details of the absorption lines. The absorption systems were fitted with the Voigt profile fitting routine VPFIT (Carswell et al 1991, 1992). Fitting components were put in until a statistically satisfying $`\chi ^2`$ value was achieved. Usually this involved obtaining a reduced $`\chi ^2`$ with a probability $`>1\%`$ to have ocurred by chance. The continuum level was treated as an additional free parameter (see paper III) to reduce systematic errors that could have crept in with the interactive continuum fitting.

The interpretation of this dataset is somewhat more challenging than usual because the QSO images closely surround (to within an arcsecond) the peaked emission from the bulge of the foreground lensing galaxy, so there is a high background level which rapidly changes in the radial direction (Fig. 1). The data show several instances of absorption lines which to the eye appear flat-bottomed and saturated while not dropping down to zero intensity in the line center. This phenomenon is most obvious in the strong absorption of the MgII complex at z=0.5656, especially in images $`B`$ and $`C`$ (see Figs. 5 and 6 and below, and the discussion in section 4.1), where there appear to be partly flat-bottomed MgII 2796 lines with a residual intensity of about 11% (in image $`B`$). At values of 1.35 and 1.16 the MgII doublet ratios for the strongest components in the $`A`$ and $`C`$ images, respectively, also suggest partial saturation. There are several possible explanations for this behaviour: (a) the absorbing gas cloud does not cover the continuum source of the background QSO and the troughs are partly filled in with light by-passing the absorber; (b) the absorption lines are very narrow and unresolved by the instrument; and (c) there is a foreground light source which partly fills in the absorption lines of systems arising between the foreground source and the background QSO. Possibility (a) is quite unlikely in that all of the low ionization systems are very far away from the QSO redshift. Case (b), unresolved narrow lines, may indeed by present, but the flat bottoms of the absorption lines and the fact that different lines often reach similar levels argues against it being the main effect: it would require a conspiracy between individual optical depths and velocity separation between the components to produce such a pattern. Possibility (c), pollution by the light of the bright lensing galaxy appears most likely: inspection of various strong lines in all three lines of sight shows that the residual intensities at the line bottoms correlate with the local flux of the QSO (they are weaker where the QSO is brighter, e.g., in the CIV emission line), as one would expect if the contribution of the QSO continuum light to the total spectrum increases relative to the light of the foreground galaxy. The relative contribution of galaxy and QSO continuum to the flux at a given wavelength, $`F_{gal}/F_{QSO}^{cont}`$, can be written in terms of the residual intensity at the bottom of an apparently saturated absorption line, $`F^{res}(\lambda )`$, as

$`{\displaystyle \frac{F_{gal}(\lambda )}{F_{QSO}^{cont}(\lambda )}}={\displaystyle \frac{F^{res}(\lambda )}{1F^{res}(\lambda )}}`$ (1)
For image $`A`$ the MgII 2796 line at z=0.5656 gives approximately $`F_{gal}/F_{QSO}^{cont}0.08`$, while the strongest CIV 1548 line at z=1.693 gives 0.04. The CIV z=1.693 line is an associated system riding almost on top of the CIV emission line and hence may be a candidate for partial continuum coverage. However, we ignore this effect here since it can only make the line less black. For image $`B`$ the MgII 2796 at z=0.5656 gives approximately $`F_{gal}/F_{QSO}^{cont}0.11`$, while the CIV 1548 line at z=1.693 gives again about 0.04. In other words the contribution of the light from the lensing galaxy to the integrated spectrum is up to about 11% near the MgII 2796 z=0.5656 system.

A naive measurement of the column densities without correcting for the light pollution by the lens would underestimate the column density for weak lines on the linear part of the curve of growth by on the order of 0 - 12% (see section 4 below), but for saturated lines they may be highly underestimated. Therefore, the column densities and measurement uncertainties in tables 1,2, 4 and 5 below (which reflect only the profile fits without correction for the residual flux from the lens) should be taken as lower limits. The exception is the $`z=0.5656`$ low ionization system (section 4) where many lines are obviously saturated, so a correction was attempted by simultaneously adjusting the zero level of the spectroscopic flux while fitting the absorption line profiles. The correction was possible because there were enough saturated regions to constrain the zero flux level, but it could not be applied meaningfully to the other absorption systems which appear to be mostly optically thin and do not have enough transitions available. Thus table 3 contains actual estimates of the absorption line properties.

## 3 CaII absorption from the Milky Way

The CaII absorption from our Galaxy can be modelled satisfactorily with a single velocity component in all three images (see the fit parameters in table 1). A pairwise comparison between any two of the three sightlines is shown in Fig. 2. The column densities and Doppler parameters are identical to within the errors, and the projected velocities along the line of sight differ by less than 2.9 $`\sigma `$ between the images. In other words, there is no significant difference in the appearance of the Galactic CaII systems in the three lines of sight. The mean weighted heliocentric velocity of the three sightlines is $`cz=13.9\pm 0.5`$ kms-1. The direction of the Q2237+0305 sight lines in galactic coordinates ($`l^{II}`$ = 71.8, $`b^{II}`$ = -46.1 deg) subtends an angle of $`63`$ deg with the heliocentric velocity vector of the local bulk flow of interstellar cloudlets within $`30`$ pc ($`v_{\mathrm{h}elio}`$ = $`28.1\pm 4.5`$ kms-1; Frisch, Grodnicki, & Welty 2002) so the projected velocity of that flow relative to the 2237+0305 line of sight would be -12.9 kms-1. This is well consistent with the idea that the CaII observed here arises from the cluster of local interstellar gas clouds within 30 pc. The maximum angle between the three lines of sight to Q2237+0305 is about 1.8 arcsec. Thus, if the Galactic CaII absorber towards 22327+0305 resides within 30pc, our observation means that there are no significant detectable differences in that particular CaII cloudlet over transverse separations smaller than 54 astronomical units.

## 4 High velocity CaII clouds related to the lensing galaxy (z=0.039)

### 4.1 Observed properties

We detect an absorption system at z=0.039, the redshift of the lensing galaxy 2237+0305. The system is visible in all three lines of sight and consists of moderately narrow ($`7.6<b<24.3`$ kms-1), presumably nonstellar Ca II H and K absorption lines. A comparison of the Ca II H and K profiles in each pair of the three images is shown in Fig. 3. Each line of sight is shown separately in Fig.4.111Note that the wavelengths of Ca II H and K quoted in Fig. 2,3 and 4 are the vacuum values. Any underlying, broad stellar Ca II absorption feature due to the lensing galaxy has been removed during the background subtraction and by drawing the continuum such as to eliminate large scale features. The velocity scale in Figs. 3 and 4 is relative to the lensing galaxy redshift $`z_{abs}=11696(\pm 40)`$ kms$`{}_{}{}^{1}/c=0.039014`$ given by Foltz et al. (1992).

Fig. 3 shows that image $`A`$ has a prominent absorption feature at +103 km s-1 and a weaker feature at -91 km s-1. 222To avoid confusion we do not show the fits on top of the spectra but give the fitting parameters redshift, Doppler parameter $`b`$ and logarithmic column density $`\mathrm{log}N`$ and their $`1\sigma `$ errors in table 2. Image $`B`$ has three (and possibly more) weak components at -288, -202, and +105 kms-1. Image $`C`$ has two strong components at -153 and -95 kms-1 and two additional weaker ones at +38 and +90 kms-1. Only the two strongest lines, the +103 kms-1 component in image $`A`$ and the -95 kms-1 component in $`C`$ are strong enough to be clearly present in both transitions of the CaII doublet, and it is only for those that we can be sure of them being CaII. Any of the other lines could in principle individually be an unidentified interloper from other transitions (there is one such case, marked by an $`x`$ in the CaII 3970 regions in Fig. 3 which is the 1550Å line of a stray CIV system at z=1.661), but the line density is much too low for this to explain the majority of components.

Using intermediate resolution spectroscopy of the integrated light from the four QSO images, Hintzen et al (1990) have searched for CaII from the lensing galaxy, deriving a $`3\sigma `$ upper limit of 72 $`m\AA `$ on the rest frame equivalent width of the CaII K line. They compared this limit with the equivalent width of CaII from the thick disk in our Galaxy and concluded that the absence of strong CaII rules out the presence of a thick disk near the center of 2237+0305. Calculating the CaII 3935 Å equivalent widths from the column densities given in table 2 we get $`W`$ = 97$`\pm 4`$, 122$`\pm 12`$, and 424$`\pm 9`$ $`m\AA `$, for the total equivalent widths of all the fitted components in images $`A`$, $`B`$, and $`C`$, respectively. The unweighted mean of the three images is 214 $`m\AA `$. These values are not corrected for the zero level problem mentioned above but if anything the real values can only be larger (probably by not more than 12%). The equivalent widths are formally inconsistent with the null detection of Hintzen et al. However, their result is an average of the CaII absorption along all four lines of sight which they observed through a 2.5” diameter circular aperture, whereas our observation was done only for the brightest three images individually with a slit 0.86” wide. The fourth line of sight (included in their average but not in ours) may not have any absorption, which would reduce the discrepancy. Another possible explanation has to do with the uncertainties due to pollution of the QSO spectra by the bright foreground galaxy. Hintzen et al used an aperture larger by more than a factor 8 than ours (assuming that we had 0.7” seeing), and the depth of the stellar H and K lines in their spectrum shows that the galactic contribution to the QSO spectrum was severe. Hintzen et al dealt with this problem by subtracting sky ”from a position 36 arcsec offset from the QSO”, and ”subtracting a suitably scaled spectrum of the M31 globular cluster K225-B280” to remove the broad stellar H and K lines. In our case, the sky background (including the stellar contribution) was subtracted in windows immediately adjacent to the QSO image. We can only speculate that given our own problems in subtracting the spatially fluctuating sky background (discussed in section 2) it is conceivable that Hintzen et al may have had similar ones. They may have subtracted less of the stellar background than we did which would have led to them underestimating the depth of the absorption and thus the CaII equivalent width.

In any case, CaII is detected and it appears that the measured total equivalent width is not dramatically different from what is seen in the Galaxy: Bowen (1991) found the typical high-latitude CaII equivalent width distribution to peak at $`W\mathrm{sin}b=110m\AA `$ for lines of sight going outward from the Galactic disk. In the present case the lensing galaxy at z=0.039 is inclined by 60 degrees (Irwin et al 1989) with respect to the plane of the sky so the equivalent width expected based on the Bowen study would be $`W=220m\AA `$. Our measurement gives the two-sided equivalent width (intersecting the whole galaxy) and translates into 214/2 = 107 $`m\AA `$ for a line of sight going out from the central plane of the galaxy. This is about half of what is seen in the Milky Way for the thick disk. However, most of the absorption components in the HIRES data formally qualify as intermediate or high velocity clouds (e.g. Wakker 2001), so a comparison of their equivalent width with a subset of higher velocity clouds in our Galaxy would be more appropriate. Galactic intermediate velocity clouds studied by Morton & Blades (1986) lead to an average of $`100\pm 50`$ $`m\AA `$ (Hintzen et al. 1990), similar to the present result, 107 $`m\AA `$.

We do not have information on the HI column densities so we cannot obtain the gas phase abundances for these clouds. If the abundances measured in the interstellar medium of the Galaxy were similar to those in the 2237+0305 galaxy, then one can use the CaII-HI abundance relation of Wakker & Mathis (2000) to determine the HI column densities. According to their formula, the CaII column densities in table 2 would roughly correspond to logarithmic HI column densities between 18.6 and 20.1, with the exception of the strongest component in the $`C`$ image (at -95 kms-1). This one formally gives a very high HI column (23.1) but it is beyond the range of validity of the Wakker-Mathis relation, and may be more indicative of a higher than usual CaII gas phase abundance.

### 4.2 Origin of the CaII velocity structure

The velocity structure of the CaII system is different from that expected from absorption arising in a thick disk:

The center of CaII absorption for all CaII components in all three sightlines (defined as the column density weighted mean velocity relative to the redshift of the lensing galaxy) occurs at -71 kms-1 with respect to the Foltz et al (1992) velocity of recession, $`11696(\pm 40)`$ kms-1, or at +4 kms-1 relative to the HI 21cm velocity, $`11621(\pm 40)`$ kms-1 given by Barnes et al (1999). Thus, within a 2-sigma deviation the column density weighted velocity centroid of the CaII absorption agrees with the centroid of the stellar and HI mass distribution. However, CaII absorption from the disk of our Galaxy seen against high latitude stars tends to have a strong component within about 10-20 kms-1 of the local standard of rest (e.g., Greenstein 1968; Blades & Morton 1983; Songaila et al 1986; D’Odorico et al 1989; Robertson et al 1991; Meyer & Roth 1991; Ho & Filippenko 1995; this paper). Similar CaII absorption near systemic galactic velocities is also seen towards other galaxies, e.g., against SN1989M, SN1993J, and SN1994I in the disks of M81 (Vladilo et al. 1994), NGC4579 (Steidel, Rich & McCarthy 1990), and NGC5194 (Ho & Filippenko 1995), respectively, and towards the nucleus of M31 (Morton & Andereck 1976). This is in contrast to the present case where there does not seem to be a zero velocity component in any of the three sightlines, no matter which of the above velocity centroids for the galaxy is adopted.

Furthermore, the spread in velocity of the CaII components is larger than expected if caused by the disk. The total unweighted velocity widths (highest minus lowest velocity component) are $`\mathrm{\Delta }v_A=193`$, $`\mathrm{\Delta }v_B=393`$, and $`\mathrm{\Delta }v_C=243`$ kms-1 for the three images $`A`$, $`B`$, and $`C`$, respectively. Even if only the relatively certain strong components at -153 kms-1 (C image) and +103 kms-1 (A image) are used the velocity difference is at least 256 kms-1. The typical line of sight from the source QSO has an impact parameter of $`1\mathrm{"}`$ or 1.1h$`{}_{50}{}^{}{}_{}{}^{1}`$ kpc. Therefore, the three lines of sight pass through the bulge of the lensing galaxy, rather than the disk, and it is clear that the rotation of the disk cannot be responsible for the large velocity widths even if it extended to the very center of the galaxy. Thus we agree with Hintzen et al, albeit for different reasons, that this case does not look like thick disk absorption.

The origin of the large velocity width of the CaII absorption is not obvious. The width of the absorption system is not unlike what is seen in sightlines towards the Magellanic clouds or in the Magellanic stream (Songaila 1981; West et al. 1985; Songaila et al 1986). It is possible that a similar configuration involving tidal gas flows (see also Bowen et al. 1994) exists here as well. In their search for HI emission in the vicinity of the lensing galaxy, Barnes et al (1999) have detected a group of galaxies which contains the lens. The two objects nearest to the lens (in projection), the dwarf irregulars NW1 and NW2 (in their nomenclature), are separated from it by about 150 kpc and 125 kpc in the plane of the sky, and in velocity of recession by -380 and -264 kms-1, respectively. It is conceivable that this group may produce a tidal feature like the Magellanic stream, but the negative velocity differences of the closest objects with respect to the lensing galaxy 2237+0305 may be inconsistent with the positive velocities of some of the CaII components. More suitable but fainter satellites of the 2237+0305 galaxy may have eluded detection.

Even with the lack of evidence for a tidal mechanism there still remains a wealth of possible causes as have been discussed in the general context of Galactic high velocity clouds (which we do not wish to enumerate here; for a review see Wakker, van Woerden & Gibson 1999, and references therein). One possibility not yet mentioned is that the CaII clouds may arise directly in the bulge of the z=0.039 galaxy. Barnes et al (1999) have argued that the stellar velocity dispersion measured by Foltz et al (1992), corrected for observational effects, agrees with various mass models in giving a line of sight velocity dispersion $`145`$ kms-1 for the bulge. The velocity differences between CaII components quoted above exceed this value by at least 110 and possibly as much as 250 kms-1 (if all CaII components are real). To account for the additional velocity dispersion one could envoke outflows (winds, SN bubbles) from the bulge, or there could be other hydrodynamic effects in the bulge that stir up the gas of an otherwise quiescent disk. Morton & Andereck (1976) have observed additional velocity components at -260 and -450 kms-1 against the nucleus of M31 and argue that they are likely to be caused by outflows. However, there is no evidence for nuclear activity in the galaxy 2237+0305.

Finally there are the size constraints from the three lines of sight. While the column density of the components is rather different between the sightlines there is a hint that at least two components of the C image (-95 and +90 kms-1 seem to correspond to weak features in the A image (at -91 and +103 kms-1). Accordingly, this particular system does have at least some coherence over 1.4 kpc, whereas the average cloud size (if defined as the distance over which the column density varies by 50% ) may be much smaller. Coherence on kpc scales may argue for clouds outside the bulge or disk, but the evidence is weak in the present case. Summarizing, the appearance of the z=0.039 system is consistent with that of high velocity clouds in the Galaxy. The velocity structure of the CaII components and the apparent hole in any CaII absorption by a thick disk argues for a kinematic connection of the CaII cloudlets with the bulge and / or locations above the plane of the lensing galaxy.

## 5 The MgII-FeII system at $`z=0.5656`$: expanding shell, collapsing sheet ?

### 5.1 Observational properties

This system is the strongest absorption complex in the spectral region covered by our observations (Figs. 5 and 6). The appearance of the MgII 2796, MgI 2853, and FeII 2600 absorption lines makes it likely that the structure is a strong Lyman limit or damped Ly$`\alpha `$ system with the hydrogen predominantly neutral: first, MgI is rather strong (e.g., $`\mathrm{log}N`$(MgI) - $`\mathrm{log}N`$(MgII) = $`1.42\pm 0.08`$ for the reddest MgI/MgII component in the A image (z=0.5669)). Photoionization models with CLOUDY (assuming an $`\alpha =1`$ powerlaw ionizing background with intensity $`J(912)=10^{21}`$ at $`z=2.5`$, blueshifted to z=0.56) show that for a strong Lyman limit system (N(HI)=$`10^{18}`$ cm-2) the total gas density would have to be rather high ($`350`$ cm-3) to get a MgI/MgII ratio as high as observed. For a damped Lyman $`\alpha `$ system with N(HI)=$`10^{20}`$ cm-3, the density would need to be only as high as $`0.5`$ cm-3. Second, strong FeII and MgII lines often indicate damped systems (Boissé et al 1998; Rao & Turnshek 2000).The total equivalent widths of MgII 2796 and FeII 2600 for the z=0.5656 system in line of sight $`A`$ are 1.20 Å and 0.73 Å respectively (after correction for the zero level). According to Rao & Turnshek, 50% of those systems with EW(MgII 2796) $`0.5\AA `$ and with EW(FeII 2600) $`0.5\AA `$ are damped (N(HI)$`2\times 10^{20}`$ cm-2).

The $`z=0.5656`$ absorption system shows an interesting degree of coherence between the three lines of sight (Figs. 5,6): in MgII 2796, two massive groups of components (near 130 and 175 kms-1 in the $`A`$ spectrum (bottom)), separated by a gap in absorption near 150 kms-1, can be traced in all three images (with the $`A`$ spectrum having an additional strong component near 255 kms-1). While the gap between the components does not shift by more than a few kms-1 between the spectra the two components seem to get further apart in velocity space as one goes from image $`A`$ to $`C`$ to $`B`$, i.e., in the order of increasing distance from $`A`$. The velocity differences $`v_A`$, $`v_C`$, and $`v_B`$, between the components are $`v_A=69.8`$ kms-1, $`v_C=81.4`$ kms-1, and $`v_B=102.5`$ kms-1). The two massive components are clearly composed of multiple components indicative of differential motion in the gas; in particular the B spectrum shows a pronounced ragged appearance. The projected relative positions of the images $`A`$, $`B`$, and $`C`$, can be expressed in terms of distances $`\overline{AB}`$ (=0.66 kpc) and $`\overline{BC}`$ (=0.50 kpc) between images $`A`$ and $`B`$,and $`B`$ and $`C`$, respectively, and the angle $`\mathrm{}(ABC)(=47.5\mathrm{deg})`$.

The fits to the absorption line profiles are listed in table 3. Experience shows that MgII, FeII and sometimes MgI have a very similar component structure in velocity space even if the relative column density levels differ. Thus for each fitting component the redshifts of MgII and FeII (and in the case of image $`A`$ also MgI) were tied together, i.e. they were forced to have the same redshift. The redshift errors are quoted in the entry for the MgII ion in table 3. Where the corresponding other ions at the same redshift have nominal errors of zero it should be understood that the MgII redshift error is actually the combined error of the tied ions.

To correct for residual flux from the lensing galaxy the zero level within the fitting regions was adjusted as an independent variable simultaneously with the profile fits. The required adjustments of the zero level relative to a unit continuum in the three regions containing the MgII doublet and MgI (region 1), the FeII 2383 line (region 2), and the FeII 2600 line (region 3), turned out to be $`7.9\pm 1.2`$%, $`0.02\pm 4.8`$%, and $`0.07\pm 7.9`$%, respectively, in the $`A`$ spectrum, $`10.6\pm 1.4`$%, $`12.2\pm 4.8`$%, and $`5.5\pm 6.6`$% in the $`B`$ image, and $`11.6\pm 0.7`$%, $`9.5\pm 1.4`$% and $`10.7\pm 1.4`$% in the $`C`$ image. These values can be used as $`F^{res}(\lambda )`$ in eqn. 1 to obtain the ratio between galaxy and QSO fluxes which made it through the spectrograph slit.

Limits on the gas temperature can be estimated from the line widths of the MgII doublet and the FeII 2600 Å line. An upper limit on the thermal width is given by the total line widths of the MgII lines in Fig. 6. Limiting ourselves to the B image (which appears to have the weakest lines, least likely to be saturated) a simultaneous fit to the MgII doublet gives a weighted mean Doppler parameter $`5.3\pm 0.4`$ kms-1 for the components in the $`z=0.5656`$ complex, corresponding to an upper limit on the temperature $`T<4.1\times 10^4`$ K. A more stringent limit may be obtained from comparison between two ions with different atomic masses, MgII and FeII. For example, in Figs. 5 and 6 the MgII 2796 spectrum of the $`B`$ image shows a group of 3 – 4 ragged components between 140 and 210 kms-1. The three strongest MgII components are well matched by corresponding FeII 2600 Å lines (see also Fig. 5). Another good case for an unblended line occurs near 250 kms-1 in the $`A`$ image. Using those 3 MgII/FeII components from image $`B`$ where $`b`$(MgII)$``$ $`b`$(FeII) together with the abovementioned line from the $`A`$ spectrum, we employ the method discussed by Rauch et al (1996) to disentangle thermal and turbulent contributions to the total Doppler width. We have omitted two other cases where FeII was broader than MgII as not providing information. The mean weighted thermal MgII Doppler parameter turns out to be $`4.0\pm 0.8`$ kms-1, corresponding to a temperature $`2.3\times 10^4`$ K. Thus there appears to be a substantial thermal contribution to the broadening (on the order of a a few times $`10^4`$ K as opposed to a few hundred Kelvin or less which would not show up observationally). Unfortunately, the large errors and the line widths being close to the resolution limit preclude any further constraints.

### 5.2 Modelling the absorbing structure as a moving shell of gas

It is impossible to say with certainty how this particular absorption system arises. However, the double component structure mentioned above, the increase of velocity separation between the two components with transverse separation between the lines of sight, and the large coherence length (Fig. 6) are suggestive of the three lines of sight intersecting two walls of a coherent, expanding (or collapsing) ”shell” (Fig. 12). Earlier we described a small low ionization absorption structure in a high redshift galaxy towards Q1422+231 which had several of the absorption characteristics expected from an expanding shell (paper I). The size of that structure, however, was at least one to two orders of magnitude smaller than the current object, which must be larger than several hundred parsecs in order to cover all three lines of sight. It is conceivable that this absorption system could be associated with galactic HI shells and supershells which are known to be at least that large (e.g., Heiles 1976,1979; Deul & den Hartog 1990; Kennicutt et al 1994; Oey & Clarke 1997). Such structures may be caused by SN explosions or wind blown bubbles (Heiles 1979) or by the infall of gas from the intergalactic medium e.g., in the form of high velocity clouds (Oort 1967). Recently, Bond et al (2001a,b), in a detailed analysis, have given plausible arguments for at least some MgII absorbers being caused by galactic ejecta, i.e., super-bubbles and -winds. The line of sight component structure of the $`z=0.5656`$ system is quite similar to that of several MgII systems in Bond et al. 2001b.

Here we attempt to model the absorption system as a spherical bubble in isotropic expansion (or contraction). Having the velocity information from three lines of sight we can infer the radius and expansion velocity of the bubble from observations of the same system intersected by multiple lines of sight. The measured input parameters for such a calculation are the three velocity differences between the main components and the three transverse distances between the lines of sight. Unfortunately, these constraints are not enough to determine both radius $`R`$ and expansion velocity $`v`$ independently, and only the relation $`R(v)`$ can be obtained. The line of sight to the fourth ($`D`$) image also would have to be observed, which was not attempted by us because of the small separation and the faintness of image $`D`$. The Appendix shows how to calculate the relation $`R(v)`$, subject to the observational constraints. The result is given in the form of a diagram in Fig. 11, and the geometry is shown in Fig. 12. There is a broad minimum in the expansion velocity leading to values below 100 kms-1 for bubble radii between a few hundred parsecs and a few kiloparsecs, with the minimum expansion velocity, $`v=55`$ kms-1, occuring at a radius $`R=0.83`$ kpc. At very small radii the diameter of the bubble approaches the transverse separation between the lines of sight. Such small bubbles require not only an increasingly unlikely alignment between the bubble walls and the lines of sight to still cover all three of them. They also lead to a near-grazing incidence of the lines of sight with respect to the bubble and they require very large expansion velocities to reproduce the velocity splitting seen between the components. At the opposite end, for very large bubbles with $`R>>1`$ kpc, the expansion velocity also would have to be very high: the spatial curvature of the bubble wall between the lines of sight would be smaller with increasing radius, and the differential velocity splitting between the major components seen in Fig. 6 can only be reproduced with an ever increasing expansion velocity.

Galactic and local extragalactic superbubbles tend to show maximum sizes of about 1-2 kpcs, and expansion velocities are typically below 150 kms-1 (Heiles 1979; Tenorio-Tagle & Bodenheimer 1988; Martin 1998). Thus the parameter space of observed superbubbles shows considerable overlap with the allowed combinations of $`R`$ and $`v`$ for our model (hatched area in Fig. 11). The observed total MgII column densities in the two main components are $`2.7\times 10^{13}`$ and $`3.3\times 10^{13}`$ cm-2 for the blue and red component in the $`A`$ image, $`1.8\times 10^{14}`$ and $`2.7\times 10^{13}`$ cm-2 for $`B`$, and $`7.8\times 10^{14}`$ and $`4.5\times 10^{13}`$ cm-2 for the $`C`$ image, respectively. We can obtain a crude upper limit on the size of the bubble by assuming that the absorption is arising entirely in a bubble shell composed of gas swept up by a wind outside the disk’s interstellar medium (assuming subsequent cooling to produce MgII as the dominant Mg ionization state). Then to produce a MgII column density as strong as observed the gas would have to come from a volume with radius

$`R=3{\displaystyle \frac{N}{n}}=2\mathrm{k}pc\left({\displaystyle \frac{N_{MgII}}{10^{14}\mathrm{c}m^2}}\right)\left({\displaystyle \frac{n_{tot}}{10^2\mathrm{c}m^3}}{\displaystyle \frac{A_{Mg}}{3.3\times 10^5}}{\displaystyle \frac{Z/Z_{}}{0.3}}{\displaystyle \frac{f_{MgII}}{0.5}}\right)^1.`$ (2)
Here $`N`$ is the total column density of the shell and $`n`$ is the background number density of the gas to be swept up, where we assume the value for the halo gas density from model A1 by Suchkov et al. (1994). $`N_{MgII}`$ is the column density, $`A_{Mg}`$ the relative elemental abundance of Mg, $`Z/Z_{}`$ the metallicity in terms of solar abundances, and $`f_{MgII}`$ the ionization fraction of MgII. The latter is quite uncertain. $`Z/Z_{}`$ is probably between 0.1 and 1, judging from the abundances of damped Ly$`\alpha `$ systems at low z (e.g. Boissé et al 1998). The radius $`R`$ is again fully consistent with observed superbubbles (if much of the interstellar medium is swept up as well the radius could be smaller, of course).

It is conceivable though that the $`z=0.5656`$ system is not produced in the bubble walls of a smallish ($`1`$ kpc) shell but rather by superwinds. Then the absorption would come from the dense shreds and filaments formed when the interstellar medium gets entrained by the hot outflowing wind (cf Bond et al 2001; Suchkov et al 1996; Heckman et al 2000, 2001; Mori, Ferrara & Madau 2002). The ragged component structure seen here may be caused by differential motion between gas clouds accelerated by the surrounding wind.

The possibility of the absorber being caused by (symmetric ?) infall of gas has also been discussed, e.g., in the form of high velocity clouds which may be sweeping up galactic material. Apparently sizes and infall velocities of known high velocity clouds are sufficient to produce large scale features in the galactic gas upon impact and may even be able to more easily satisfy the energetic requirements of the largest shells than collective supernova explosions (e.g., Oort 1967; Tenorio-Tagle 1981; Tenorio-Tagle & Bodenheimer 1988 and references therein; Rand & Stone 1996). At high redshift low ionization metal absorption systems from infalling gas are expected to arise as a natural by-product of galaxy formation (Rauch, Haehnelt & Steinmetz 1997). Once the gas makes it past the accretion shock surrounding the forming galaxies it cools and can be seen in MgII or damped Ly$`\alpha `$ absorption. However, it is not clear whether cosmological infall of pre-enriched material can reproduce all the small scale density and velocity structure seen in real metal absorption systems (paper I, II).

On the basis of the observations presented here and the rather idealized theoretical models we can currently not elucidate the origin of the z=0.5656 system any further. It would be interesting to study more MgII systems in the spectra of lensed QSOs to see whether there is a statistical difference in spatial and velocity extent between the two-component systems like the one seen here and randomly selected systems.

## 6 The z=0.82743 system

This system also has a prominent double-component structure that appears in all lines of sight, albeit shifted and with rather different column densities (Figs. 7 and 8). The velocity separation $`\mathrm{\Delta }v`$ between the two most prominent MgII lines in each of the A, C and B images is 35.0, 47.9 and 18.6 kms-1, respectively. The velocity structure is more complicated than in the z=0.5656 system in that here are also shifts of the column-density weighted velocity along the line of sight, $`\overline{v}`$, between the lines of sight, which amount to $`\overline{v_C}\overline{v_A}=32.2`$ kms-1, $`\overline{v_B}\overline{v_A}=16.5`$ kms-1, and $`\overline{v_B}\overline{v_C}=64.2`$ kms-1.

## 7 The z=0.97163 system

For this system (Figs. 9 & 10) only FeII lines are covered. Again there is some resemblance between the absorption patterns for the three lines of sight. The differences between the column density weighted velocities along the line of sight are $`\overline{v_C}\overline{v_A}=4.1`$ kms-1, $`\overline{v_B}\overline{v_A}=13.4`$ kms-1, and $`\overline{v_B}\overline{v_C}=17.4`$ kms-1. The differential motion in this and the z=0.8274 system are small enough to be attributed to rotation, .e.g., of an entire gaseous halo (e.g., Lanzetta & Bowen 1992; Steidel et al 2002), but even for the simplest rotation model there are too many free parameters to be constrained by the anecdotal evidence presented here, and a larger sample is clearly required.

## 8 Conclusions

(1) CaII absorption from the Galaxy is found in all three lines of sight to Q2237+0305 with an average heliocentric velocity $`13.9\pm 0.5`$ kms-1 and without significant differences between the lines of sight. This velocity is consistent with the absorption arising in the local bulk flow in the ISM within $`30`$ pc from the sun (Frisch et al 2002). In that case the largest angular separation between the sightlines corresponds to $`<1.8\mathrm{"}\times 30`$ pc = 54 astronomical units. Apparently, the Galactic CaII cloud seen toward Q2237+0305 is smooth on scales smaller than that.

(2) We have detected a high velocity cloud system related to the lensing galaxy 2237+0305 (z=0.039). Unlike CaII systems in our Galaxy this one does not show a strong component at the systemic velocity of the galaxy as measured either from the stellar or HI centroid. The system arises within impact parameters on the order of 1 kpc to the center of the galaxy (i.e., behind and / or in front of its bulge), so the large velocity differences among the components along individual lines of sight and between different lines of sight cannot be a result of the rotation of the galactic disk. The absorption line strengths are similar to those measured for intermediate and high velocity clouds in our Galaxy. The system may have an origin in the bulge or halo of the galaxy 2237+0305.

(3) Three MgII or FeII selected low ionization systems not associated with the lensing galaxy show coherent absorption over the transverse separation between the lines of sight, ranging from about 200 - 600 $`h_{50}^1`$ pc. This is a strict lower limit to the size of the the gaseous structures giving rise to the entire MgII system: they must be at least that large to cover all three lines of sight. The closeness of the lines of sight in the present case does not permit an upper limit to be obtained.

(4) individual MgII components are much smaller than entire MgII systems: interpreting the separation where the optical depth between two lines of sight differ by 50% as a minimum ‘size’ of the clouds, the three systems examined here show fluctuations at that level for separations at least as small as $`200300h_{50}^1`$ pc. In fact, several components do not show up at all in more than one image (Figs. 5-10). This is consistent with the Petitjean et al upper limit but it indicates that many of the cloudlets are actually much smaller than a kpc. This result appears in contrast to size estimates for CIV systems. Individual CIV clouds appear to be larger than about $`300h_{50}^1`$ pc (paper II). Ionization calculations place the MgII gas at a higher density than the more common, higher ionization, CIV gas clouds. Thus it is plausibly that the MgII cloudlets are smaller than those clouds as well.

(5) We have tentatively interpreted the low ionization system at $`z=0.5656`$ as a moving shell, possibly caused by an expanding superbubble or galactic wind. The large scale coherence of the system and the differential velocities between the absorption in the three lines of sight are consistent with supershell parameters similar to those observed in the Galaxy and in local dwarf galaxies. In this picture the two MgII components observed may correspond to the wall of the cooling shell at the interface with the intergalactic medium and / or the entrained material within the hot bubble. Observations of transitions associated with higher ionization stages, especially CIV and OVI, in the UV should be able to shed further light on the nature of this system.

We thank Crystal Martin for helpful discussions, and Bob Carswell and Gary Ferland for making VPFIT and CLOUDY, respectively, available. At Keck, Gary Puniwai, Ron Quick, Gabrelle Saurage, David Sprayberry, Wayne Wack, and Greg Wirth helped with the observations. The anonymous referee is thanked for several helpful suggestions. MR is grateful to the NSF for grant AST-0098492 and to NASA for grant AR-90213.01-A. The work of WLWS was supported by NSF grant AST-9900733.

## 9 Appendix

This is an outline of the calculation of the relation between the expansion velocity $`v`$ of a spherical shell, and its radius $`R`$, as a function of the observed velocity spreads between absorption line pairs like the one in the $`z=0.5656`$ MgII system ($`v_A=69.8`$ kms-1, $`v_B=102.5`$ kms-1, $`v_C=81.4`$ kms-1) and the geometric constraints (the projected relative positions of the images $`A`$, $`B`$, and $`C`$, as given by the distances $`\overline{AB}`$ (=0.66 kpc) and $`\overline{BC}`$ (=0.50 kpc), and the angle $`\mathrm{}(ABC)(=47.5\mathrm{deg})`$).

There are four unknowns, which can be expressed as (for example) $`v`$, $`R`$, $`b`$ (the impact parameter of line of sight $`B`$ with respect to the center of the sphere), and $`\varphi `$ (i.e., the angle between $`b`$ and $`\overline{AB}`$, the line connecting lines of sight $`A`$ and $`B`$), and three equations connecting them (see below). Thus we can only solve for a relation between any two of these. The relation $`v(R)`$ is probably the most useful one as it determines which values of size and expansion velocity are consisting with the observations. The various geometric relations between these quantities can be written as:

$`R^2\left(1\left({\displaystyle \frac{v_B}{2v}}\right)^2\right)=b^2`$ (3)
$`R^2\left(1\left({\displaystyle \frac{v_A}{2v}}\right)^2\right)=b^2+\overline{AB}^2+2b\overline{AB}\mathrm{cos}\varphi `$ (4)
$`R^2\left(1\left({\displaystyle \frac{v_C}{2v}}\right)^2\right)=b^2+\overline{BC}^2+2b\overline{BC}cos(\varphi +\mathrm{}(ABC))`$ (5)
Elimination of the unwanted parameters leads to the single equation

$`R^4(8\overline{AB}^3\overline{BC}v^2(4\overline{BC}^2v^2+R^2(v_C^2v_B^2))y+`$
$`2\overline{AB}\overline{BC}R^2(v_B^2v_A^2)(4\overline{BC}^2v^2+R^2(v_C^2v_B^2))y+`$
$`16\overline{AB}^4\overline{BC}^2v^4+\overline{BC}^2R^4(v_B^2v_A^2)^2+\overline{AB}^2(16\overline{BC}^4v^4+R^4(v_B^2v_C^2)^2`$
$`8\overline{BC}^2R^2v^2(v_C^2v_A^2y^2+8v^2z^2v_A^2z^2+v_B^2(1+y^2z^2))))`$ $`=`$ $`0,`$ (6)
where $`y=\mathrm{cos}(\mathrm{}(ABC))`$, and $`z=\mathrm{sin}(\mathrm{}(ABC))`$. It can be solved for the velocity $`v`$. Of the four solutions, two are imaginary and one describes a contracting ($`v<0`$) rather than an expanding sphere. The remaining solution is well-behaved for $`0.33`$ kpc $`<R<\mathrm{}`$:

$`v`$ $`=`$ $`1/2[(\overline{AB}\overline{BC}^3R^2(v_A^2v_B^2)y+\overline{AB}^3\overline{BC}R^2(v_C^2v_B^2)y+`$ (7)
$`\overline{AB}^2\overline{BC}^2R^2(v_B^2(1+y^2z^2)v_C^2v_A^2)+`$
$`[(\overline{AB}^2\overline{BC}^2R^4((\overline{BC}^22\overline{AB}\overline{BC}y4R^2z^2+\overline{AB}^2)(\overline{AB}^2(v_B^2v_C^2)^2`$
$`2\overline{AB}\overline{BC}(v_B^2v_A^2)(v_B^2v_C^2)y+\overline{BC}^2(v_B^2v_A^2)^2)+`$
$`(\overline{BC}^2(v_B^2v_A^2)y+\overline{AB}^2(v_B^2v_C^2)y+\overline{AB}\overline{BC}(v_C^2+v_B^2(z^21y^2)+v_A^2))^2)]^{1/2})/`$
$`(\overline{AB}^2\overline{BC}^2(\overline{BC}^22\overline{AB}\overline{BC}y4R^2z^2+\overline{AB}^2))]^{1/2}`$
Inserting the numerical values gives

$`v=0.5\sqrt{4.878\times 10^{15}\sqrt{{\displaystyle \frac{R^4(0.3959+R^2)}{(9.345684.5541R^2)^2}}}+{\displaystyle \frac{3.4712\times 10^{15}R^2}{(84.5541R^29.3456)}}},`$ (8)
where $`R`$ is now to be measured in units of kpc, and $`v`$ in cms-1. The relation $`v(R)`$ is shown in Fig. 11.


format: text

text:
# Inhomogeneous Quasi-stationary States in a Mean-field Model with Repulsive Cosine Interactions

## 1 Introduction

Mean-field models, that is models in which all the particles interact with equal intensity regardless of distance, have been the object of considerable interest (see  for a review and the references therein). One of the simplest models of this kind consists of $`N`$ particles moving on a circle coupled globally by a cosine interaction, with Hamiltonian:

$$H=\frac{1}{2}\underset{i=1}{\overset{N}{}}p_i^2+\frac{J}{2N}\underset{i,j=1}{\overset{N}{}}\mathrm{cos}(\theta _i\theta _j).$$
(1)
It is sometimes called the Hamiltonian Mean-Field model (HMF). Here the variables $`p_i`$ are the momenta conjugate to $`\theta _i`$, which is the angle describing the state of the $`i`$-th particle. As stated above, the strength of the interparticle interaction does not depend on the distance and all particles interact with all others. Alternatively, one can think of this model as representing a mean-field approximation to the classical XY model, though this is certainly not a realistic model for a spin system. Its equilibrium statistical mechanics can be treated exactly via the Hubbard-Stratonovich transformation, as shown in . The attractive case, corresponding to $`J<0`$, is shown to have a phase transition at a certain value of the inverse temperature $`\beta _c=2/J`$, above which the stable state has a clustered structure: that is, the particle density on the circle shows a non-trivial profile. Analogous results were also obtained using entropy maximization methods  and, numerically, for a time discrete version of model (1. On the other hand, in the repulsive case ($`J>0)`$, no phase transition is found . This can be heuristically justified through the following physical argument. First, note that Hamiltonian (1) can be rewritten as follows:

$$H=\frac{1}{2}\underset{i=1}{\overset{N}{}}p_i^2+\frac{J}{2N}\left[\left(\underset{i=1}{\overset{N}{}}\mathrm{cos}\theta _i\right)^2+\left(\underset{i=1}{\overset{N}{}}\mathrm{sin}\theta _i\right)^2\right].$$
(2)
From this follows immediately that, at least if $`N`$ is an even number, any state in which all spins appear in pairs $`(\theta _k,\pi +\theta _k)`$ for all $`k`$ corresponds to a ground state configuration of the system, all of such states having the same energy. Since the positions of half the angles can be chosen arbitrarily without in any way affecting the energy, uniform ground states will be entropically favoured, thus precluding the appearance of a thermodynamically stable density modulation at low temperatures, as is indeed rigorously established in Refs. .

It was therefore quite intriguing to find spatially inhomogeneous solutions in the antiferromagnetic case . Numerically, they have been found in molecular dynamics simulations starting from certain initial conditions at low temperatures. These so-called bicluster solutions show two paired clusters in the density, that is, particles appear preferentially in correlated positions $`(\theta ,\pi +\theta )`$. These biclusters were found to be extremely stable, though it is of course never easy to decide on numerical ground whether they are simply very long-lived or truly stationary.

This paper addresses the following issues: in Section 2, we give a cursory description of the numerical findings concerning the bicluster. In Section 3, we show how the $`N`$ particle system can be reduced to a system in which the particles are only linked to each other through the common interaction with a collective mode. We shall show that this phenomenon is reminiscent of the wave-particle interaction mechanism that governs a large amount of plasma physics phenomena . To make this picture even clearer and more quantitative, we show in Section 4 how this system can be further reduced to a one-particle time dependent Hamiltonian. Finally, in Section 5, we perform a detailed analysis of the one-particle model, showing that several features of the bicluster can be explained in these terms. Finally, in Section 6, we present our conclusions.

## 2 Description of the bicluster

In the following, we consider a system of $`N`$ particles interacting according to Hamiltonian (1) with a positive value of $`J`$, corresponding to repulsion among the particles. We shall consider almost exclusively the following initial conditions: the initial velocities are all identically zero, and the angles are randomly distributed according to a uniform distribution on the interval $`[0,2\pi ]`$. Since this is not in general an equilibrium position, the particles start moving and acquire typical velocities of the order of $`N^{1/2}`$. We are therefore always dealing with the low-energy dynamics of the system in the large $`N`$ limit. Under these circumstances, it is found numerically that a bicluster forms, that is, the particle density becomes inhomogeneous. The early stages of this formation were studied analytically in Ref. , using a convenient zero-temperature limit of the Vlasov equation. These involve the rapid formation of rather complex spatio-temporal inhomogeneities, so-called chevrons, which, however, disappear on a fairly rapid time scale. Since we are here primarily interested in the long-time behaviour, these shall not concern us in the following.

Once the initial stages of growth described in Ref.  are over, an apparently stationary density $`\rho (\theta )`$ arises. The time evolution of this density can be characterized by the following moments:

$$M_k(t)=\frac{1}{N}\underset{m=1}{\overset{N}{}}\mathrm{exp}\left(ik\theta _m(t)\right).$$
(3)
The second moment $`|M_2(t)|`$ is shown in Figure 1 for a large system ($`N=10^5`$ particles). After some oscillations corresponding to the chevron structures described in Ref. , the second moment eventually reaches a constant value of order one (though on the times shown in the figure, the increase is not yet quite over). Apart from quite rapid oscillations, which we always average over, the density $`\rho (\theta )`$ settles down onto a provisionally stable profile, the bicluster, which is shown in Figure 2.

It is found on numerical evidence  that the time-averaged moments $`M_k^{(0)}`$ in this quasi-stationary state are well approximated by the expression

$`\left|M_0^{(0)}\right|`$ $`=`$ $`1`$
$`\left|M_{2k}^{(0)}\right|`$ $`=`$ $`{\displaystyle \frac{1}{|k|}}(k0),`$ (4)
whereas they are zero for odd values of $`k`$. This leads to the following analytical expression for the density (again after averaging over the rapid oscillations):

$$\rho _0(\theta )=\frac{1}{2\pi }(1\mathrm{ln}|2\mathrm{sin}\theta |.)$$
(5)
Here the moments $`M_k^{(0)}`$ and the density $`\rho _0(\theta )`$ are connected by the following relation:

$$M_k^{(0)}=\frac{1}{2\pi }_0^{2\pi }\rho _0(\theta )e^{ik\theta }𝑑\theta .$$
(6)
The bicluster state does not last forever, at least for finite values of $`N`$. Indeed, we have performed additional numerical simulations of the system and could observe the decay of the second moment to its equilibrium value over a long time scale. Specifically, for $`N=50`$, the relaxation of $`|M_2(t)|`$ is displayed in Figure 3 for two different energy densities. We see that a significant relaxation to equilibrium has taken place by a time of approximately $`410^8`$ for both energies. These times rapidly increase as $`N`$ increases, so that, in some sense, one might argue that the bicluster is some kind of equilibrium state. This is related to the issue of the non-commutation of the limits $`t\mathrm{}`$ and $`N\mathrm{}`$ in mean-field or Vlasov-like systems (see e.g. references  for a numerical approach). On the other hand, the appearance of a bicluster is closely linked to the peculiar choice of the initial conditions: in particular, if instead of starting with zero velocities and random positions, one starts with random velocities of the same order as those that will eventually develop in the former initial condition, no bicluster formation is observed. Together with the fact that we are dealing with a very low-energy phenomenon in a system with long-range forces, this leads us to suspect that the bicluster is simply a phenomenon involving metastable configurations. Since it is well known that metastable states become infinitely long-lived in the mean-field limit, this point of view is in full agreement with the one involving non-commuting long-time and thermodynamic limits.

## 3 A simple equivalent model

Our purpose in this section is to map Hamiltonian (1) onto a system of the following form:

$$\stackrel{~}{H}=\frac{\omega (\mathit{\varphi })^2P^2}{2NJ}+\frac{NJX^2}{2}+\frac{1}{2}\underset{k=1}{\overset{N}{}}p_k^2.$$
(7)
Here $`P`$ and $`X`$ stand for two canonically conjugate variables. The vector $`\mathit{\varphi }`$ has coordinates $`\varphi _k`$, conjugate to the $`p_k`$ and $`\omega (\mathit{\varphi })`$ is a smooth function which will be computed explicitly in the following. Here the $`\varphi _k`$ correspond to the angles of the original particles, whereas $`P`$ and $`X`$ correspond to an explicit separation of a collective mode. One therefore has the following situation: the system consists of non-interacting particles coupled to an oscillatory excitation, in a similar way as that familiar in the wave particle models in plasma physics . Furthermore, as we shall see, the function $`\omega (\mathit{\varphi })`$ is an average over all $`\varphi _k`$, and therefore a slowly varying function. Under these circumstances, one expects a very accurate (though in general not exact) conservation of the adiabatic invariant related to the energy stored in the collective mode. As we shall see, if this invariant has an initial value much larger than its thermal average value, a bicluster forms and lasts as long as the invariant maintains its value.

To show this correspondence, we start by considering system (1) as a system of particles in two dimensions, where all the particles are additionally constrained to move on the unit circle. In order to deal with the constraint, it is convenient to go to the Lagrangian formulation. One finds for the Lagrangian corresponding to (1):

$`L`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{k=1}{\overset{N}{}}}(\dot{\overline{x}}_k^2+\dot{\overline{y}}_k^2){\displaystyle \frac{J}{2N}}[({\displaystyle \underset{k=1}{\overset{N}{}}}\overline{x}{}_{k}{}^{})^2+({\displaystyle \underset{k=1}{\overset{N}{}}}\overline{y}{}_{k}{}^{})^2]+`$ (8)
$`+{\displaystyle \underset{k=1}{\overset{N}{}}}\lambda _k(\overline{x}{}_{k}{}^{2}+\overline{y}{}_{k}{}^{2}1).`$
This has the form of a non-linearly constrained $`N`$-dimensional harmonic oscillator. We now separate the motion of the center of mass, since it plays a peculiar role. We therefore introduce:

$`X_1`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{k=1}{\overset{N}{}}}\overline{x}_kX_2={\displaystyle \frac{1}{N}}{\displaystyle \underset{k=1}{\overset{N}{}}}\overline{y}_k`$
$`x_k`$ $`=`$ $`\overline{x}_kX_1y_k=\overline{y}_kX_2.`$ (9)
In these new coordinates the Lagrangian reads

$`L`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{k=1}{\overset{N}{}}}(\dot{x}_k^2+\dot{y}_k^2)+{\displaystyle \frac{N}{2}}\left[\dot{X}_1^2+\dot{X}_2^2J(X_1^2+X_2^2)\right]+`$ (10)
$`+{\displaystyle \underset{k=1}{\overset{N}{}}}\lambda _k(x_k^2+y_k^2+2X_1x_k+2X_2y_k+X_1^2+X_2^21).`$
Going back to polar coordinates and making the following approximations, which are valid in the limit in which the formation of the bicluster is observed,

$$X_1,X_21\dot{X}_1,\dot{X}_21,$$
(11)
one finally obtains for the Lagrangian, after some manipulations shown in Appendix A,

$$L_0=\frac{\dot{X}^2}{2}\left(N+\underset{k=1}{\overset{N}{}}\mathrm{cos}^2\varphi _k\right)+\frac{N}{2}X^2\dot{\mathrm{\Phi }}^2+\underset{k=1}{\overset{N}{}}\dot{\frac{\varphi _k^2}{2}}\frac{NJ}{2}X^2.$$
(12)
Here $`X`$ and $`\mathrm{\Phi }`$ are related to $`X_1`$ and $`X_2`$ via

$$X_1=X\mathrm{cos}\mathrm{\Phi }X_2=X\mathrm{sin}\mathrm{\Phi }.$$
(13)
Let us now switch back to the Hamiltonian picture. The Hamiltonian corresponding to (12) is exactly given by

$$H=\frac{\omega (\mathit{\varphi })^2P^2}{2NJ}+\frac{NJX^2}{2}+\frac{\mathrm{\Lambda }^2}{2NX^2}+\frac{1}{2}\underset{k=1}{\overset{N}{}}p_k^2.$$
(14)
Here $`P`$ is conjugate to $`X`$ and $`\mathrm{\Lambda }`$ to $`\mathrm{\Phi }`$. As a final remark, we note that $`\mathrm{\Lambda }`$ is an exact constant of the motion, namely the total angular momentum. It can therefore always be set equal to zero without loss of generality. This therefore leads us to (7). The frequency function $`\omega (\mathit{\varphi })`$ is given by

$$\omega (\mathit{\varphi })=\sqrt{\frac{NJ}{N+_{k=1}^N\mathrm{cos}^2\varphi _k}}.$$
(15)
Hamiltonian (7) is very similar in form to those developed in the theory of wave particle interactions in plasma physics. However, it should be emphasized at this stage that the system described by (7) is still an $`N`$-dimensional system. Contrarily to the reduction proposed in  where bulk and tail (i.e. suprathermal) particles can be discriminated, our system is ‘cold’ so that all particles participate at the same level to the collective mode.

A very appealing interpretation of (7) is the following. Consider a ring of mass one on which $`N`$ beads of mass $`1/N`$ can move freely. If $`X`$ is taken to be the amplitude of the torsional vibrations of such a ring around the $`y`$ axis, and if no other motion of the ring is allowed, then (7) describes the motion both of the ring and the beads. Consequently, if we start with particles at rest and with an initial vibration entirely in the degree of freedom associated to the ring, we will have the ring oscillating first around its axis, and then driving the beads towards those regions where the ring motion is least, that is, towards the two points of the ring lying on the $`y`$ axis. Such a structure corresponds exactly to what is observed for the bicluster.

We now proceed with the calculation. Our aim is twofold. On the one hand, we show that equilibrium states do not contain a density modulation. On the other, we display a mechanism capable of leading to very long-lived inhomogeneous states.

Let us first consider the partition function of the effective Hamiltonian (7). The $`p_k`$ integrals factorize and one is left with

$$Z=_0^{2\pi }𝑑\varphi _1\mathrm{}𝑑\varphi _N\omega (\mathit{\varphi })^1𝑑P𝑑X\mathrm{exp}\left(\frac{\beta }{2}(P^2+X^2)\right)$$
(16)
This, however, is easily shown not to lead to any density modulation in the angles, since $`\omega (\mathit{\varphi })`$ is of order one. We must therefore look beyond statistical equilibrium, something already indicated by the plausibility arguments mentioned in the Introduction as well as by exact results.

To this end, we transform $`X`$ and $`P`$ locally to action angle variables, that is, we transform to

$`I`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{\omega (\mathit{\varphi })P^2}{NJ}}+{\displaystyle \frac{NJX^2}{\omega (\mathit{\varphi })}}\right)`$
$`\psi `$ $`=`$ $`\mathrm{arctan}{\displaystyle \frac{\omega (\mathit{\varphi })P}{NJX}}.`$ (17)
This transformation must be complemented by an approppriate transformation of the $`\varphi _k`$ to make it canonical. This is computed and discussed further in Appendix B. The final hamiltonian is

$$\stackrel{~}{H}=\omega (\mathit{\varphi })I(\mathit{\varphi })+\frac{1}{2}\underset{k=1}{\overset{N}{}}\left(P_k+\frac{I(\mathit{\varphi })\psi }{\omega (\mathit{\varphi })}\frac{\omega (\mathit{\varphi })}{\varphi _k}\right)^2,$$
(18)
where the $`P_k`$ are the new particle momenta after the canonical transformation. They are given by (39). Note that the interaction terms, which are linear in the $`p_k`$, all average to zero in the fast variables $`\psi `$, and the others are of order $`1/N^2`$. They can therefore be accurately taken into account via a perturbative averaging approach. The adiabatic theorem then states that the action $`I(\mathit{\varphi })`$ is conserved to a high degree of precision and over long time scales. We may therefore construct an approximately stationary state as follows: set an arbitrary value for the action $`I`$ and, using this now as an external parameter, compute the resulting equilibrium state for the $`\varphi _k`$.

It should be noted here that recently, using a completely different approach involving multiple time expansion, close results have been derived by Barré et alin Ref. . We believe that our approach complements their work by seeking to provide some intuition for the physical mechanisms involved in these processes.

If we now take $`I`$ to have given values in (18) and proceed to compute the Gibbs partition function, we can still eliminate the shifts in the $`P_k`$’s and then factor the resulting integral out. We are then left with

$$Z=_0^{2\pi }𝑑\varphi _1\mathrm{}𝑑\varphi _N\mathrm{exp}\left[\beta \omega (\mathit{\varphi })I\right].$$
(19)
Since $`\omega (\mathit{\varphi })`$ is of order one, it is necessary for $`I`$ to be of order $`N/\beta `$ (that is, $`N`$ times larger than predicted by equipartition) in order to induce a density modulation of order one in the angles. On the other hand, if such a modulation is induced, it is straightforward to verify that it will indeed be a $`\pi `$-periodic modulation directed along the direction in which the collective mode oscillates, and therefore has most of the obvious properties observed in the bicluster. Strictly speaking, however, the microcanonical ensemble should be used in this calculation, which may possibly open a broader range of possibilities, as regions of negative temperature become accessible. For a discussion of these issues, see .

As stated in Section 2, various initial conditions may or may not lead to the formation of a bicluster. Among those that are successful, there are

1. Initial conditions with equispaced angles and a sinusoidal amplitude velocity perturbation: these initial conditions yield a bicluster always, as they put in an extensive (though very small) value for the action $`I`$ of the collective oscillation.
2. Initial conditions with random angles and zero initial velocities: For these it is readily verified that the actions of the collective mode are of order one, whereas the typical velocity of a particle decays as $`1/\sqrt{N}`$. From this follows that the effective inverse temperature $`\beta `$ increases as $`N`$, thus yielding the appropriate result in (19).

## 4 Reduction to a One Particle Model

We here show how to further reduce the model to a one-particle time dependent Hamiltonian. The assumptions involved are less well controlled and the results we obtain only coincide with the original model over fairly short times. However, in many features, it will turn out that the description given by this simplified model is extremely satisfactory. To this end, we make the following remark: in the model given, say, by (7), all particles are coupled to one single collective mode. The dynamics of this mode is determined both by an external harmonic force and by the reaction forces to all the particles involved. The system is therefore not a one-particle system. In order to reduce it to such a system, we require an approximate expression for the behaviour of the macroscopic mode. Once this is given, a one-particle approximation for the whole system immediately follows. From (7) one obtains

$$\ddot{X}=\omega ^2(\mathit{\varphi })X.$$
(20)
Let us now consider the picture of the bicluster developed in Section 3. We had found there that the bicluster was a quasi-stationary state in which the action of the collective mode had a given (macroscopic) value and the $`\mathit{\varphi }`$ a corresponding stationary distribution. From this picture and (20) follows that, in the bicluster, the motion of the collective mode can be well approximated by a simple harmonic motion.

These observations are confirmed numerically in Figure 7 of  which shows the phase-points of $`X_1`$ in the dynamical regime where the bicluster forms. One can see that $`X_1`$ undergoes rapid oscillations of small amplitude around the origin, with its phase only slowly changing in time. Note that, in reality, a slow drift in the relative phase between $`X_1`$ and $`X_2`$ is actually observed. This is not, however, due to a non-zero value of the angular momentum, but rather to the fact that the description of the motion as a two-dimensional harmonic oscillator is only approximate and not valid over large times.

Armed with this knowledge concerning the dynamics of the collective mode, we can proceed to analyze the behaviour of the particles. We go back to the original $`\theta _k`$ variables. The equation of motion following from (2) can now be rewritten as follows:

$$\ddot{\theta }_k=J\left[X_1(t)\mathrm{sin}\theta _kX_2(t)\mathrm{cos}\theta _k\right].$$
(21)
From this, using the approximate description of the collective mode described above, we obtain the following one-particle approximation:

$$\ddot{\theta }_k=\epsilon \mathrm{cos}\omega t\mathrm{sin}\theta _k.$$
(22)
Here the small parameter $`\epsilon `$ is equal to $`J|X_1|`$. Note that, by hypothesis, the amplitude of the collective mode $`X_1(t)`$ is always small, so that this last approximation is in fact justified (see (11)). In the next section we analyze this model in greater detail numerically. Before we do this, however, let us quickly show how this can be treated perturbatively. The dynamics (22) is described by the following Hamiltonian

$$H=\frac{p^2}{2}\epsilon \mathrm{cos}\omega t\mathrm{cos}\theta .$$
(23)
Using standard averaging techniques, one finds in the limit of small $`\epsilon `$ that the motion is well described by the following effective Hamiltonian:

$$\overline{H}=\frac{P^2}{2}\frac{\epsilon ^2}{8\omega ^2}\mathrm{cos}2\mathrm{\Theta }.$$
(24)
Here $`\mathrm{\Theta }`$ is $`\theta `$, from which a small rapid oscillation has been substracted and $`P`$ is the canonicallly conjugate momentum. Explicitly

$$\mathrm{\Theta }=\theta +\frac{\epsilon }{2}\left(\frac{\mathrm{sin}(\omega t\theta )}{(\omega P)^2}\frac{\mathrm{sin}(\omega t+\theta )}{(\omega +P)^2}\right).$$
(25)
It follows therefore that the particle moves, at least in an average sense, in an effective potential of the form $`\mathrm{cos}2\mathrm{\Theta }`$, thereby leading to an attraction to the two points zero and $`\pi `$.

## 5 Analysis of the one-particle model

In this Section, we investigate numerically whether the dynamical features of the original model (1) with cold initial conditions can be accounted for by the minimal effective model (22) with one-and-a-half degrees of freedom derived in the preceding Section. We first observe that we can set $`\omega =1`$ without loss of generality, because this parameter can be absorbed in a trivial time rescaling. We then introduce the following area-preserving map with time-step $`\tau `$, that constitutes the lowest order symplectic algorithm that numerically solves (22),

$`\theta _{n+1}`$ $`=`$ $`\theta _n+p_{n+1}\tau `$ (26)
$`p_{n+1}`$ $`=`$ $`p_n\epsilon \tau \mathrm{cos}\left(\left(n+1\right)\tau \right)\mathrm{sin}\theta _n.`$
In order to reproduce the conditions of the simulations performed in , we consider a set of $`N`$ initial conditions featuring a spatially homogeneous distribution of particles of zero momentum. That is, we take for all $`j`$, $`1jN`$

$$\theta _0^j=j\frac{2\pi }{N}\mathrm{and}p_0^j=0.$$
(27)
We then compute the resulting $`N`$ trajectories, obtaining their superposed phase-space plots of figures 4 and 6. Since $`\omega =1`$ the period $`T`$ of the rapid oscillations of the collective variable $`X_1`$ takes the value $`T=2\pi `$. In these simulations, we fix the time step $`\tau =T/100=2\pi /100`$ and $`\epsilon =0.05`$. Figure 4 shows the superposed phase-space plots of the $`N`$ trajectories (27) during the initial stage from which the bicluster structure emerges. All trajectories wind around the elliptic fixed points $`\theta =0,\pi `$. In figure 5 we plot some trajectories with initial value $`|\theta _0|<\pi /2`$, which show rapid oscillations of period $`T`$ superposed on a much slower oscillatory motion around the minimum in $`\theta =0`$ of the time averaged potential $`\mathrm{cos}2\mathrm{\Theta }`$ in Hamiltonian (24). The large time plots of figures 6 clearly reproduce the features of the full HMF $`N`$-degrees of freedom phase-space plots (e.g. figure 4 in reference ). The bicluster structure exhibits an overall oscillation on the period $`T`$ around the zero momenta fixed points $`\theta =0`$ and $`\theta =\pi `$.

It is also possible, using the effective Hamiltonian (24), to get at least a qualitative understanding of the behaviour of the position and momentum density functions, in particular of their singularities. Indeed, if we assume that all particles start from rest and are randomly distributed on the unit circle, we find for the density $`\rho (\theta )`$

$$\rho (\theta )=(2\pi Z)^1_0^{2\pi }𝑑\theta _0_{\mathrm{}}^{\mathrm{}}𝑑p\delta \left[\overline{H}(\theta ,p)\overline{H}(\theta _0,0)\right],$$
(28)
where $`Z`$ is the microcanonical phase space volume and where we have neglected the difference between the transformed variables and the original ones. The integral in (28) can be evaluated explicitly to yield

$$\rho (\theta )=𝒩K(\mathrm{cos}^2\theta ),$$
(29)
where $`K(m)`$ is the complete elliptic integral of the first kind as a function of the parameter $`m`$ and $`𝒩`$ is an appropriate normalization constant. The properties of the approximate solution (29) are very similar to the ones observed for the full $`N`$-body system, e.g. we can reproduce the logarithmic singularities at $`\theta =0,\pi `$. The general shape is also similar (see figure 7), but the agreement is not as satisfactory as that of the empirical formula (5), shown in figure 2.

A similar calculation for the momentum density yields

$$\rho (p)=𝒩K\left(1\frac{w^4}{16}\right)(|w|2),$$
(30)
where $`w`$ is the appropriately scaled momentum $`p/(\epsilon ^2/8\omega ^2)^{1/2}`$. Again, this has roughly the right form (see figure 8), but, quantitatively, it is not quite satisfactory.

Model (22) seems therefore to be a good representation of the dynamics of the antiferromagnetic HMF at vanishing energy, even if the assumption of a constant frequency $`\omega :=2\pi /T`$ is (slightly) violated due to self-consistency in the real system.

Let us now show that, although simple, this model has a non trivial dynamics. First note that KAM theorem, which guarantees the preservation of quasiperiodic motions under small perturbations, cannot be applied in its classical form to this system. Renaming the time $`t`$ as a phase variable $`\phi `$, the dynamics (22) derives from the two-degrees of freedom Hamiltonian

$$H(\theta ,p,\phi ,u)=\frac{1}{2}p^2+u\epsilon \mathrm{cos}\omega \phi \mathrm{cos}\theta $$
(31)
where the variable $`u`$ is conjugated to time $`\phi `$ and does not appear in (22). KAM theorem requires a nondegeneracy condition on the frequency vector $`𝝎:=(_pH,_uH)`$ to ensure that a large set of actions $`(p,u)`$ have “sufficiently irrational” frequencies. This condition is $`det\left(_{(p,u)}𝝎\right)0`$. This is trivially not satisfied here as $`𝝎=(p,1)`$. Yet, within the theory of averaging, KAM theory may be applied to an averaged Hamiltonian exponentially close to (24) with the result that, for $`\epsilon `$ small enough, the Poincaré section of $`\overline{H}`$ is filled up to a residue of exponentially small measure by invariant curves that are close to the level lines of $`\overline{H}`$ . This appears to be confirmed numerically, although in a rather unexpected fashion. Figure 9 displays the enlargement of a long-time Poincaré plot with initial conditions equally spaced on the $`\theta `$-axis and $`\epsilon =0.05`$. Over the time scales probed by the simulation ($`410^7`$ timesteps), the particle does not visit the whole phase space, but remains within a well-defined region in agreement with KAM predictions. In this region, however, it appears to lie on an extremely convoluted but smooth invariant surface. Indeed, the Poincaré surface shows a dense pattern of parallel lines on which the points lie. These presumably represent the intersection of a very high order torus with the Poincaré surface. Such behaviour is certainly peculiar, but may be well related to the degenerate nature of the unperturbed system. From this follows, in this particular case, that no short periodic orbits can exist in this part of the phase space.

If this picture is correct, the system is not ergodic. This then strongly suggests that, in this reduced one-particle model, the bicluster structure lasts forever. The relevance of this conclusion for the $`N`$ particle system, however, remains an open question. It is in principle conceivable that, among the class of initial conditions we consider, a non-vanishing measure lies on some invariant surface of the $`N`$ particle system. However, this certainly does not follow from our previous arguments concerning adiabatic invariance, which only holds over very long, but not over infinite times. Furthermore, the simulations discussed in Section 2 speak against such a possibility.

## 6 Conclusions

Summarizing, we have found a physically appealing approximate description of the low-temperature dynamics of the repulsive Hamiltonian mean-field model in terms of free particles on a ring performing torsional vibrations. The physical analog of the ring oscillation is the collective plasmon oscillation of the variables $`X_1`$ and $`X_2`$. The interaction between particles and ring arises solely from the influence of the particle positions on the moment of inertia of the ring. This provides a straightforward interpretation of the clustering phenomena observed in previous work: if the initial condition is such that the initial ring oscillation strongly dominates the motion of the individual particles, then a parametric instability sets in, driving the particles towards the part of the ring which is at rest. This creates, as shown in Section 4, an effective potential of the form $`\mathrm{cos}2\mathrm{\Theta }`$ in which the particles move. This same picture can also be used to obtain a description of the statistical mechanical equilibrium. In this case, it is clear that no clusters form, either in the canonical or in the microcanonical ensemble. Furthermore, this effective potential allows to form a qualitatively correct picture of both the particle density in position and in momentum space, though clear discrepancies remain, showing that the time-dependent nature of the problem is essential. Finally, we studied the chaotic properties of the time-dependent effective model. For a small, yet finite, modulus of the collective plasmon variable, this model was shown to present localization on very high order tori, and hence lack of ergodicity. This localization is presumably inherent to the small dimension of the system and might not occur in the finite $`N`$ particle model. Indeed, even if the tori present in the effective one-particle system actually survived in the $`N`$ particle system, Arnold diffusion might still act as a process that would eventually drive the system towards thermal equilibrium.

MCF thanks European Commission for support through a Marie Curie individual fellowship Contract No HPMFCT-2000-00596. This work was partially funded by grant IN112200 of DGAPA as well as CONACYT grant number 32173-E. Finally, we acknowledge financial support from INFN, the University of Florence and the MURST-COFIN00 project Chaos and localization in quantum and classical mechanics.

## Appendix A Derivation of the Approximate Lagrangian

We start from (10). We now introduce polar coordinates as follows:

$$x_k=\rho _k\mathrm{cos}\varphi _ky_k=\rho _k\mathrm{sin}\varphi _k.$$
(32)
Here the angles $`\varphi _k`$ are measured from the vector $`(X_1,X_2)`$. From this follows readily via the cosine theorem, that

$$\rho _k^2+X^22X\rho _k\mathrm{cos}\varphi _k=1,$$
(33)
where $`X`$ is defined as the norm and $`\mathrm{\Phi }`$ as the angle of the vector $`(X_1,X_2)`$. Within the approximations stated in (11), this yields

$$\rho _k=1+X\mathrm{cos}\varphi _k.$$
(34)
From this follows for the velocities

$$\dot{\rho }_k=\dot{X}\mathrm{cos}\varphi _kX\dot{\varphi }_k\mathrm{sin}\varphi _k=\dot{X}\mathrm{cos}\varphi _k,$$
(35)
where again the last equality uses (11) and $`\dot{X_1},\dot{X_2}\dot{\varphi }_k`$. The Lagrangian (10) now reads

$$L=\frac{1}{2}\underset{k=1}{\overset{N}{}}(\dot{\rho }_k^2+\rho _k^2\dot{\varphi }_k^2)+\frac{N}{2}\left[\dot{X}^2+X^2\dot{\mathrm{\Phi }}^2JX^2\right]$$
(36)
Substituting (34) and (35) into (36) and using systematically the approximations (11) yields the desired result.

## Appendix B Transformation to Action-angle Variables

In this appendix, we compute the canonical transformation that generates the local transformation to action-angle variables (17). We shall do this by computing its generating function as follows: we first determine a function $`S_0(I,X;\varphi _k)`$ such that the conditions

$$\frac{S_0}{I}=\psi \frac{S_0}{X}=P$$
(37)
are equivalent to (17). This is simply the generating function of the transformation to action-angle variables of the harmonic oscillator, where the $`\varphi _k`$ enter merely as parameters. If we now define

$$S(I,X;\varphi _k,P_k)=\underset{k=1}{\overset{N}{}}P_k\varphi _k+S_0(I,X;\varphi _k),$$
(38)
this determines in the usual way a canonical transformation on the whole space. To compute it explicitly, note that $`S_0`$ depends only on the combination $`\omega (\mathit{\varphi })I`$, so that

$`p_k`$ $`=`$ $`P_k+{\displaystyle \frac{I\psi }{\omega (\mathit{\varphi })}}{\displaystyle \frac{\omega (\mathit{\varphi })}{\varphi _k}}`$ (39)
$`=`$ $`P_k+{\displaystyle \frac{I\psi }{\omega (\mathit{\varphi })}}\mathrm{cos}\varphi _k\mathrm{sin}\varphi _k\sqrt{{\displaystyle \frac{NJ}{\left(1+_{l=1}^N\mathrm{cos}^2\varphi _l\right)^3}}}`$
Note that these correction terms are of order $`1/N`$.