text:
# A High-Resolution Rotation Curve of NGC 6822: A Test-case for Cold Dark Matter

## 1 Introduction

NGC 6822 is a dwarf irregular local group member. Located at a distance of $`490\pm 40`$ kpc (Mateo, 1998), NGC 6822 is the most nearby dwarf irregular apart from the LMC/SMC system. Due to the small distance the galaxy appears very extended on the sky: its optical angular diameter is over a quarter of a degree; the Hi disk measures close to a degree (de Blok & Walter, 2000). NGC 6822 is a member of an extended cloud of irregulars (Mateo, 1998) known as the “Local Group Cloud”. The galaxy has a total luminosity of $`M_B=15.8`$ (Hodge et al., 1991) and a total Hi mass of $`1.3\times 10^8M_{}`$ (de Blok & Walter, 2000), making it relatively gas-rich. It is a metal poor galaxy, with an ISM abundance of about 0.2 $`Z_{}`$ (Skillman et al., 1989) and has a star formation rate of $`0.06M_{}\mathrm{yr}^1`$ (based on H$`\alpha `$ and FIR fluxes) (Mateo, 1998; Israel et al., 1996). Hodge (1980) found evidence for increased star formation between 75 and 100 Myr ago, while Gallart et al. (1996) showed that the star formation in NGC 6822 increased by a factor of 2 to 6 between 100 and 200 Myr ago. This is consistent with the mostly constant but stochastic recent star formation histories often derived for dwarf and LSB galaxies (Grebel, 2001; Gerritsen & de Blok, 1999). NGC 6822 can be regarded as a rather average and quiescent dwarf irregular galaxy.

de Blok & Walter (2001) found that the outer Hi disk is dominated by what appears to be a super-giant Hi shell with extended, apparently tidal features in the outer disk. In contrast, the inner main Hi disk is remarkably uniform and circularly symmetric. There is no evidence for large kinematical disturbances, as we will show in this paper.

In this paper we address the detailed dynamics and kinematics of NGC 6822 and derive high-resolution rotation curves of this dwarf galaxy. The physical resolution of our highest resolution curve is $`20`$ pc or $`0.03h`$ where $`h`$ is the exponential scale length of the stellar disk (see Sect. 6.1), making it the highest resolution Hi rotation curve of any undisturbed dwarf galaxy.

In Section 2 we summarize the “small-scale crisis” in current cosmological simulations. Section 3 briefly describes the observations and data reduction, while Sect. 4 and 5 describe the derivation and analysis of the rotation curves in detail. In Section 6 we describe the models used in rotation curve fitting, while Section 7 presents the results. Section 8 discusses the results, and tries to fit them in the CDM picture. We summarize our conclusions in Section 9.

## 2 Dark Matter in Dwarfs

Cosmological numerical Cold Dark Matter (CDM) simulations predict a specific and universal shape for the dark matter mass density distributions (and therefore the rotation curves) of galaxies (Dubinski & Carlberg, 1991). This was investigated in detail by Navarro, Frenk & White (1996, 1997) who found that the mass-density distribution in the inner parts of simulated CDM halos could best be described by a $`r^1`$ power-law. This steeply increasing density cusp towards the center naturally translates into a steeply rising rotation curve. This kind of curve is however not observed in late-type disk galaxies (e.g. Swaters, 1999; Verheijen, 1997) . The rotation curves of dwarf and Low Surface Brightness (LSB) galaxies are shallow and rise linearly, more consistent with a dark matter distribution that is dominated by a constant-density core with a size of a few kpc (Moore, 1994; de Blok & McGaugh, 1997; de Blok et al., 2001b; Blais-Ouellette, Amram & Carignan, 2001). Especially in the case of LSB galaxies there has been much discussion about the reality of the observed shallow curves (Swaters, Madore & Trewhella, 2000). It was argued that systematic resolution effects (“beam-smearing”) could hide the steep CDM curves and lead to an erroneous conclusion that they would not be consistent with CDM.

Subsequent high-resolution follow-up studies, measuring the rotation curves of LSB galaxies in H$`\alpha `$ at a resolution of $`0.2`$ to $`1`$ kpc have now confirmed the conclusions derived from the early data (de Blok et al., 2001b; McGaugh, Rubin & de Blok, 2000; de Blok et al., 2001a; de Blok & Bosma, 2001; Blais-Ouellette, Amram & Carignan, 2001; Bolatto et al., 2002) Dwarfs and LSB galaxies seem to be dominated by a dark matter distribution that is best described in the inner parts by a powerlaw $`\rho r^\alpha `$ with $`\alpha =0.2\pm 0.2`$ (de Blok et al., 2001b).

With a spatial resolution of only $`20`$ pc, the current data set enables us to measure the rotation curve of this dwarf galaxy at a resolution which is an order of magnitude better than the observations described above. Among other things, the high resolution makes it possible to investigate the effect beam size has on mass models and test whether the inferred distribution of dark matter depends on resolution.

## 3 Observations and data reduction

NGC 6822 was observed with the Australia Telescope Compact Array for $`15\times 12`$ hours in its 375 (1 $`\times `$ 12h), 750D ($`2\times 12`$h), 1.5A (4 $`\times `$ 12h), 6A and 6D ($`8\times 12`$h) configurations over the period from June 1999 to March 2000. A total of 8 pointings was observed covering the entire Hi extent of the galaxy. We used a bandwidth of 4 MHz with a channel separation of 0.8 $`\text{km}\text{s}^1`$. Additionally, to obtain zero-spacing information, NGC 6822 was observed with the Parkes single dish radio telescope using the multibeam receiver in its narrowband mode in December 1998. The correlator configuration was identical to that used for the Compact Array observations.

The data were reduced and mosaicked together using the miriad data reduction package. Super-uniform weighting, reducing side lobes in individual pointings prior to mosaicing, was used. The resulting data cubes were cleaned with the miriad mossdi task. We combined the single dish and synthesis data to correct for the missing zero-spacings and produced data cubes at various spatial and velocity resolutions. For the rotation curve analysis presented here we use the data cubes with channel separations of 1.6 $`\text{km}\text{s}^1`$ (with an effective velocity resolution of 1.9 $`\text{km}\text{s}^1`$). These cubes gave the best compromise between velocity resolution and signal-to-noise.

We analyzed the 1.6 $`\text{km}\text{s}^1`$ data at five different spatial resolutions. The synthesized beam sizes used were $`96^{\prime \prime }\times 349.4^{\prime \prime }`$, $`48^{\prime \prime }\times 174.7^{\prime \prime }`$, $`24^{\prime \prime }\times 86.4^{\prime \prime }`$, $`12^{\prime \prime }\times 42.4^{\prime \prime }`$ and $`8^{\prime \prime }\times 28.3^{\prime \prime }`$. The position angle of the beam was 0 at all resolutions. We will refer to the 5 different resolution data sets as B96, B48, B24, B12 and B08, respectively. Table 1 gives some more information on the data sets. We retained the elliptical beamshape to obtain the highest possible resolution along the major axis. As the PA of the major axis of NGC 6822 is $`110^{}`$, as shown later, this means that the resolution along the major axis is to a few percent equal to the minor axis beam size. The large size of the galaxy compared to the beam rules out any beam smearing effects due to non-major axis information entering the beam.

## 4 Moment maps and velocity field

All subsequent analysis was done using the gipsy package. The low-resolution cube B96 was clipped at the $`2\sigma `$ level, and remaining noise peaks were removed by hand. Then for B48 and subsequent resolutions, we used the $`2\sigma `$-clipped cube of the previous resolution as a mask, after which remaining spurious noise peaks were removed by hand. This ensures a consistent selection of features at all resolutions. Furthermore using the lower resolution cube as a mask ensures that possible extended low-level structures are retained in the map.

The integrated Hi surface density maps were made in the usual manner by adding together all clipped channel maps. As a second step we then isolated the high signal-to-noise (S/N) regions of the maps as follows. For uniformly tapered maps in velocity $`\sigma _{tot}=\sqrt{N}\sigma _{ch}`$, where $`\sigma _{tot}`$ is the noise in a pixel in the integrated column density map, $`N`$ is the number of channels contributing to that pixel, and $`\sigma _{ch}`$ is the noise in one channel at that pixel. We constructed noise maps for each integrated column density map, and used these together with the column density maps to isolate those pixels in the column density maps where $`S/N>10`$. These high $`S/N`$ maps were used as masks for the velocity fields. Figure 1 shows the integrated column density map derived for the B12 data.

There are two standard ways to produce a velocity field. One commonly used procedure is to determine the intensity weighted mean of the velocities along each profile. It is known that this method can produce spurious velocities at low resolution and low S/N ratios. Given the high resolution and S/N of our data set, this is unlikely to be a problem here. An alternative method is to fit (a) Gaussian(s) to each profile. The latter method usually gives better results, as it is less affected by systematic effects due to skewed profiles etc. It is however computationally more expensive which becomes important for data cubes of the size we are dealing with here. The high quality of the current data puts us firmly in the regime where the intensity weighted mean gives accurate results.

To test this we have compared both methods for the B12 data set (Fig. 2) and find that for our data the differences are negligible. A histogram of the differences is very well described by a Gaussian with an average of $`0.3`$ $`\text{km}\text{s}^1`$ and an RMS of $`1.0`$ $`\text{km}\text{s}^1`$, i.e., both methods produce identical results to within better than a channel spacing. The number of pixels where the residuals reach $`5`$ $`\text{km}\text{s}^1`$ or higher is only $`0.1`$ percent of the total number of pixels.

Figure 3 shows a major axis position velocity diagram of the B12 data, taken along a position angle of 120. Over-plotted are the velocities as found in the intensity weighted and Gaussian velocity fields. The velocity fields produce the same results, without large-scale systematic deviations. We will therefore proceed to use intensity weighted velocity fields in subsequent analysis

## 5 Deriving Rotation Curves

The rotation curves were produced using the gipsy task rotcur. For a well-resolved and high S/N data set such as the current one, the tilted ring procedure is by far the best way to determine the rotation curve. Other methods which are often used, such as adjusting the tilted ring parameters by hand using position velocity diagrams as a guide (e.g. using the gipsy command inspector), are more subjective and results can depend on guesses on e.g. the magnitude and effects of beam smearing.

A look at Fig. 2 shows that the high-resolution velocity fields contain an extra-ordinary amount of detail. To get a good feel for the large scale structure in NGC 6822 we started by deriving the B96 rotation curve, and used the results for each resolution as the initial estimate for the next higher resolution, thus gradually refining the curve.

We adopted rings with a width of equal to the minor axis of the beam. The beam is elongated roughly along the minor axis of NGC 6822 with an effective radius of $`1.9`$ times the minor axis. A ring width equal to the minor axis is a good compromise between beam shape and complete sampling. Most information contributing to the rotation curve comes from near the major axis where the beam is narrow.

In general the procedure involved making a fit with all parameters (systemic velocity, centre position, position angle \[PA\], inclination and rotation velocity) free. The central position and systemic velocity were then fixed and several runs with either inclination or PA or both fixed were made to find the best model.

As the minor axis generally provides little information regarding the rotation curve we excluded an angle of 30 around the minor axis from the fits. One usually also down-weights the data around the minor axis by applying a $`|\mathrm{cos}\theta |`$ weighting, where $`\theta `$ is the angle with respect to the major axis in the plane of the galaxy. Other alternatives are to apply a uniform weighting or a $`\mathrm{cos}^2\theta `$ weighting. We have experimented with different weighting schemes as well as varying the exclusion angle around the minor axis between 15 and 45, but found no difference in the curves produced (this is due to the large number of independent points that do contribute at full weight around the major axis). We therefore adopt a uniform weighting scheme and a free angle around the minor axis of 30.

Once a satisfactory curve using both sides of the velocity field was produced, we also calculated two curves using only the approaching and receding sides of the velocity field. Comparing the three curves gives information about the symmetry of the system.

### 5.1 The rotation curves

For each resolution we derived the curve using the procedure described above. The results are shown in Table 2 and Fig. 4. The trends found in PA and inclination for the different resolutions agree well with each other. Only the B96 data does not show the PA trend as strongly. The variation of inclination with radius is only small; only the highest resolutions show weak evidence of a slight warp in the outer parts: choosing a constant inclination instead therefore does only marginally affect the curves. The trend in PA is pronounced and real; deriving curves with a constant PA results in velocities and inclinations that are inconsistent with the data, especially at the higher resolutions. The small-scale kinks in the inclination are not physical and mainly due to dispersion effects.

Beyond $`R1000^{\prime \prime }`$ (2.5 kpc) the difference between the approaching and receding sides becomes more pronounced with increasing resolution. This radius corresponds to the edge of the inner Hi disk. The receding curve shows no significant features there, but the approaching curve shows a dip in velocity. This coincides with the interface between the main Hi disk and the NW cloud (de Blok & Walter, 2001). However, at larger radii the receding curve converges with the approaching one, indicating that the galaxy and cloud (NW side) as well as the “tails” (SE side) are embedded in one symmetrical halo. The global dynamics of the system appear undisturbed despite the morphology of the Hi . The small kink at $`R150^{\prime \prime }`$ (0.38 kpc) is caused by a number of high velocity dispersion regions that over a small radial range happen to be aligned with the tilted rings. At $`R500^{\prime \prime }`$ a similar small kink corresponds with the inner edge of the large hole.

Figure 5 shows an overlay of the rotation curve on various position-velocity slices. The rotation curve is a good representation of the dynamics of the galaxy as a whole.

For all resolutions we have compared the velocity fields with model velocity fields constructed from the tilted ring fits, and find no systematic large scale residuals. Figure 6 compares the B12 velocity fields. There are only a few localised regions of slightly higher residuals in the higher resolution models, coinciding with regions of high velocity dispersion. Other resolutions give similar results.

### 5.2 Final curves

The error-bars derived from the tilted ring fits are an indication of the scatter in velocity around the best-fitting velocity at each tilted ring. As such they do not take into account large-scale asymmetries, differences between approaching and receding sides etc. To define more realistic errors, we regard the absolute difference between the rotation velocities of the approaching and receding sides as $`2\sigma `$ errors. We adopt the maximum of the tilted ring uncertainties and the asymmetry uncertainties as our error-bars, with one further modification. In some cases both sides of the galaxy are very symmetrical, and these usually are also the radii where the tilted ring uncertainties are very small. These combine to give very small errors $`(0.1\text{km}\text{s}^1)`$. We have therefore imposed a minimum error of half a channel width on the curves. After corrections for inclination this translates in a minimum error of $`0.8`$ $`\text{km}\text{s}^1`$.

In principle the Hi rotation curves also need to be corrected for the pressure gradients in the gas to derive the true rotation velocities. This correction for asymmetric drift is given by

$$v_c^2=v_\varphi ^2\sigma ^2\left(\frac{\mathrm{ln}\rho }{\mathrm{ln}R}+\frac{\mathrm{ln}\sigma ^2}{\mathrm{ln}R}\right),$$

where $`v_c`$ is the true rotation velocity, $`v_\varphi `$ the observed gas rotation velocity, $`\sigma `$ the velocity dispersion in the gas, and $`\rho `$ the volume density. The median (modal) velocity dispersion in the disk of NGC 6822 is $`5.8`$ ($`5.7`$) $`\text{km}\text{s}^1`$, as measured in the 1.6 $`\text{km}\text{s}^1`$ channel spacing data, and does not show any obvious radial trends. Assuming a constant scale height, we can derive the correction for asymmetric drift. We consider the B08 curves where, due to the large gradients in the Hi surface density, the corrections are largest. The corrected and uncorrected curves are shown in Fig. 7. Both curves are very similar, with the largest corrections of about $`3\text{km}\text{s}^1`$ occuring between 850′′ and 1050′′. It should be kept in mind that the correction itself is also uncertain and depends on assumptions of e.g. constant scale height, which may not hold in the outer parts of the gas disk. As the uncertainties due to asymmetric drift in the inner parts are negligible, ignoring them will not affect the results. Based on this analysis, we decided not to apply this correction.

### 5.3 Dynamical centres

The central position as derived above varies slightly from resolution to resolution. This is understandable given the different beam sizes. Fig. 8 shows the central positions over-plotted on a B12 column density map of NGC 6822, where the error-bars span the major and minor axis FWHM of the respective beams. Also plotted are $`K`$-band 2MASS (see below) and $`R`$-band isophotal centers. Due to the large size, low galactic latitude and the irregular morphology the optical center in NGC 6822 is difficult to determine. From different fits performed at various isophotal levels we estimate an uncertainty of $`0.5^{}`$. These positions are offset from the dynamical center by about 105′′ or $`0.25`$ kpc.

Though the dynamical center is *by definition* the zero-point of the rotation curve, it is instructive to re-derive the rotation curve under the assumption that the optical center is the true center of the galaxy. This is implicitly assumed in many emission line observations of rotation curves where the optical center of the galaxy is frequently used to line up the slit of the spectrograph. We re-derived the B12 rotation curve fixing the position of the center to that of the optical center. We kept the run of position angle and inclination identical to that in the original model. Fig. 9 over-plots the two curves. The overall shape of the curve is insensitive to the precise position of the center. The small-scale differences are of the same order of magnitude as the differences between the approaching and receding sides of the original curve. One should however not exaggerate the importance of this offset: had NGC 6822 been at a more typical distance of a few tens of Mpc, then the angular size of the offset would only be $`1^{\prime \prime }`$ and would have been unnoticed111 As an aside, we refer ahead to the discussion on central mass-density slopes in Sect. 8.1, where we determine the power-law slope $`\alpha `$ of the mass-density profile $`\rho r^\alpha `$ within a radius of 0.8 kpc. For the B12 case we there find a slope of $`\alpha =0.04\pm 0.09`$. For the curve presented here, determined using the optical center and also using an outer radius of 0.8 kpc, a similar analysis gives $`\alpha =0.23\pm 0.24`$. The uncertainties are larger, but the value is consistent with a soft core in NGC 6822. The exact position of the dynamical center does not significantly influence the results..

## 6 Mass components

### 6.1 Stellar component

We used a mosaic of full-resolution 2MASS $`K_s`$ images to derive the surface brightness profile of NGC 6822. The orientation parameters and the centers of the tilted ring fits to the velocity fields were used to measure the average surface brightness at each radius. Prior to the ellipse integration bright stars were masked out. The average surface brightness contribution of remaining foreground stars in the field away from NGC 6822 was determined and subtracted from the profile.

The $`K_s`$ surface brightness of NGC 6822 is low and at large radii the scatter is substantial. At radii beyond $`750^{\prime \prime }`$ the signal of the galaxy can no longer be distinguished from the sky background. Fig. 11 shows the azimuthally averaged surface brightness profile, corrected for inclination, determined using the B08 tilted ring parameters. The surface brightness was corrected for Galactic foreground extinction, though this was only a small correction of 0.088 mag in the $`K`$-band. We assumed NGC 6822 to be optically thin in $`K`$, and only applied a geometric inclination surface brightness correction.

The profile can be well described by an exponential disk. As the outer parts of the profile could be marginally affected by small uncertainties in the sky background we have fitted an exponential disk to the profile at $`R<400^{\prime \prime }`$ and find an exponential scale length of 286′′ or 0.68 kpc. The central surface brightness is $`\mu _0^i(K_s)=19.8`$ mag arcsec-2. The total absolute magnitude we derive assuming an exponential disk and extrapolating to infinity is $`M_K=17.9`$ or $`L_K=3.8\times 10^8L_{}`$.

For each resolution we determine the surface brightness profiles using the relevant tilted ring parameters determined. At $`R>240^{\prime \prime }`$ we replace the surface brightness profiles with that of the exponential disk fit as indicated in Fig. 11. The rotation curve of the stellar disk was then computed from the extended $`K`$-band profile following Casertano (1993) and Begeman (1987). The disk was assumed to have a vertical sech2 distribution with a scale height $`z_0=h/6`$ (van der Kruit & Searle, 1981). The rotation curves of the stellar component were sampled at the same radii as the rotation curves. We assume $`\mathrm{{\rm Y}}_{}`$ is constant with radius. While one expects some modest variation in $`\mathrm{{\rm Y}}_{}`$ with radius in optical bands (de Jong, 1996), the colour gradients in dwarf LSB galaxies in the $`K`$ band tend to be small, so this effect is not likely to be significant (see also Sect. 6.4).

### 6.2 Gas component

We used the orientation parameters of the tilted ring fits to determine the azimuthally average Hi surface density profiles using the total column density maps for each resolution. The Hi surface density profiles are presented in Fig. 10. In deriving the corresponding rotation curve with the gipsy task rotmod we scaled the surface density profiles by a factor of 1.4 to take the contribution of helium and metals into account.

We assume that the gas is distributed in an infinitely thin disk. This assumption is not crucial. Though one might expect dwarf galaxies to have a more puffed up disk than a typical spiral galaxy, the shape and amplitude of the gas rotation curve depends only very mildly on the thickness of the disk.

To use the B12 curve as an example, for the infinitely thin gas disk case the maximum rotation velocity of the gas component is $`17`$ $`\text{km}\text{s}^1`$. When increasing the thickness of the disk, this maximum rotation velocity drops only very slowly: for an exponential scale height $`z=0.5`$ kpc, it still measures $`14.5`$ $`\text{km}\text{s}^1`$. This value for $`z`$ is already approaching the value of the radial scale length of the stellar distribution and likely already too large. We thus adopt the thin disk approximation.

### 6.3 Halo models

The dark halo component differs from the previous two in that we are interested in parametrising this component assuming some fiducial model. The choice of this model is the crux of most of the DM analyses in the literature, and many models exist. These can all be broadly distinguished in two groups: halo models with a core, and halo models with a cusp. An example of the first category is the pseudo-isothermal halo, an example of the latter the CDM NFW halo.

#### 6.3.1 Pseudo-isothermal halo

The spherical pseudo-isothermal (ISO) halo has a density profile

$$\rho _{ISO}(R)=\frac{\rho _0}{1+(R/R_C)^2},$$
(1)
where $`\rho _0`$ is the central density of the halo, and $`R_C`$ the core radius of the halo. The corresponding rotation curve is given by

$$V(R)=\sqrt{4\pi G\rho _0R_C^2\left[1\frac{R_C}{R}\mathrm{arctan}\left(\frac{R}{R_C}\right)\right]}.$$
(2)
The asymptotic velocity of the halo, $`V_{\mathrm{}}`$, is given by

$$V_{\mathrm{}}=\sqrt{4\pi G\rho _0R_C^2}.$$
(3)
To characterise this halo only two of the three parameters $`(\rho _0,R_C,V_{\mathrm{}})`$ are needed, as equation (3) determines the value of the third parameter.

#### 6.3.2 NFW halo

The NFW mass density distribution takes the form (Navarro, Frenk & White, 1996)

$$\rho _{NFW}(R)=\frac{\rho _0}{(R/r_0)[1+(R/r_0)]^2}$$
(4)
where $`r_0`$ is the characteristic radius of the halo and $`\rho _0`$ the characteristic density (Navarro, Frenk & White, 1996; Blais-Ouellette, Amram & Carignan, 2001). This mass distribution gives rise to a halo rotation curve

$$V(R)=V_{200}\left[\frac{\mathrm{ln}(1+cx)cx/(1+cx)}{x[\mathrm{ln}(1+c)c/(1+c)]}\right]^{1/2},$$
(5)
where $`x=R/R_{200}`$. It is characterised by a concentration parameter $`c=R_{200}/R_s`$ and a radius $`R_{200}`$. These are directly related to $`R_s`$ and $`\rho _i`$, but are used instead as they are a convenient way to parametrise the rotation curve. The radius $`R_{200}`$ is the radius where the density contrast exceeds 200, roughly the virial radius (Navarro, Frenk & White, 1996). The characteristic velocity $`V_{200}`$ of the halo is defined in the same way as $`R_{200}`$. These parameters are not independent and are set by the cosmology.

### 6.4 Stellar Mass-to-light ratios and weighting

One of largest uncertainties in any mass model is the value of the stellar mass-to-light ratio $`\mathrm{{\rm Y}}_{}`$. Though broad trends in $`\mathrm{{\rm Y}}_{}`$ have been measured and modelled (e.g. Bottema, 1997; Bell & de Jong, 2000), the precise value for an individual galaxy is not well known, and depends on extinction, star formation history, Initial Mass Function, etc. The value of $`\mathrm{{\rm Y}}_{}`$ cannot be constrained using the rotation curve alone (van Albada & Sancisi, 1986; Lake & Feinswog, 1989) and some assumptions must be made. We therefore present disk-halo decompositions using four different assumptions.

Minimum disk. This model assumes that the observed rotation curve is due entirely to DM. This gives an upper limit on how concentrated the dark mass component can actually be. It is not a realistic model *per se*, because it ignores the gas disk which is obviously present. However, we present it here for ease of comparison with cosmological simulations, a large fraction of which tend to only model the dark matter component.

Minimum disk + gas. The contribution of the atomic gas (Hi and He) is taken into account, but $`\mathrm{{\rm Y}}_{}`$ is assumed to be zero.

Constant $`\mathrm{{\rm Y}}_{}`$. Here we choose a value for $`\mathrm{{\rm Y}}_{}`$ appropriate for NGC 6822 based on its stellar content, colours and likely star formation history. Unfortunately reliable integrated colours for NGC 6822 are rare due to its large angular size. We can however use indirect arguments to arrive at a value for $`\mathrm{{\rm Y}}_{}`$. Furthermore we are interested here only in the $`K`$-band value which is rather insensitive to effects of extinction and recent star formation.

Verheijen (1997) presents $`K^{}`$-band rotation curve fits of a large sample of HSB and LSB galaxies in the Ursa Major cluster. We consider the rotation curve fits derived using the Bottema disk prescription (Bottema, 1997). This recipe states that the maximum rotation velocity of the disk is 63% of the maximum rotation velocity observed (and is derived from measurements of stellar velocity dispersions). The distribution of $`\mathrm{{\rm Y}}_{}_K`$ ratios derived with the Bottema disk for the dwarf and LSB galaxies in the UMa sample peaks at $`\mathrm{{\rm Y}}_{}{}_{K}{}^{}=0.4\pm 0.05`$. If the colours of NGC 6822 are comparable to those found for other dwarfs and LSB galaxies (namely $`BV0.5`$, $`BR0.8`$ and $`VI0.8`$; de Blok, van der Hulst & Bothun 1995) we can use Table 3 in Bell & de Jong (2000) to derive an approximate value for $`\mathrm{{\rm Y}}_{}`$. For their formation epoch model with burst we find values for $`\mathrm{{\rm Y}}_{}_K`$ between 0.35 and 0.40. Similar values are found for the infall models and the hierarchical models. From detailed population synthesis modelling of NGC 6822, Gallart et al. (1996) find that the total mass of stars and stellar remnant ever formed in NGC 6822 must be $`910^7`$ $`M_{}`$. Combined with the $`K`$-band luminosity derived earlier, this gives a value $`\mathrm{{\rm Y}}_{}{}_{K}{}^{}=0.24`$.

Taking the above into account we adopt a value of $`\mathrm{{\rm Y}}_{}{}_{K}{}^{}=0.35`$ as representative of the stellar population in NGC 6822.

Maximum disk. The rotation curve of the stellar component is scaled to the maximum value allowed by the observed rotation curve, but with the restriction that the DM density gradient is required to be flat or negative at all radii (thus avoiding a “hollow halo”) (van Albada & Sancisi, 1986).

Each of the rotation curves was fitted using the gipsy task rotmas. The program determines the best-fitting combination of $`R_C`$ and $`V_{\mathrm{}}`$ (for the pseudo-isothermal halo) or $`c`$ and $`V_{200}`$ (for the NFW halo), using a least squares fitting routine. We assigned weights to the data points inversely proportional to the square of their uncertainty.

## 7 Results

Figures 12-16 show the rotation curves at all resolutions plotted for the four assumptions on $`\mathrm{{\rm Y}}_{}`$ for both NFW and ISO halo models. To rule out the possibility that the outer parts of NGC 6822, which might be affected by tidal effects, affect our conclusions we also present separate fits to just the inner part of the galaxy out to to $`R=1000^{\prime \prime }`$ (2.4 kpc) or the edge of the main inner disk. Tables 3 and 4 present the fit parameters for all models presented here.

An inspection of Tables 3 and 4 shows that the goodness-of-fit is in all cases better for ISO than for NFW. The only exceptions are some of the maximum disk fits where both ISO and NFW are difficult to fit, as NGC 6822 is not a maximum disk galaxy. A comparison between the $`\chi ^2`$ values is shown in Fig. 17.

The maximum disk values of $`\mathrm{{\rm Y}}_{}`$ generally have an uncertainty of $`\pm 0.1`$. Changing the values within this range has no discernible effect on the quality of the fits. One point of note is that for the B08 case the maximum disk value for $`\mathrm{{\rm Y}}_{}`$ is less than that for the constant $`\mathrm{{\rm Y}}_{}`$ case. The B08 constant $`\mathrm{{\rm Y}}_{}`$ fit thus results in a hollow halo.

Many of the NFW fits have fit parameters that do not make physical sense. The large majority of the fits prefer $`c0`$ and $`V_{200}\mathrm{}`$, which is another way of saying that the fit procedure is trying to fit a $`VR^{1/2}`$ curve to a $`VR`$ curve by stretching it to infinity, thus trying to take out the curvature in the model. Where this happened we fixed the $`c`$ parameter to $`c=0.1`$ with only $`V_{200}`$ as the free parameter. Such a value is still outside the range predicted by cosmological simulations where values of $`c10`$ would be expected. There is a general trend for the value of the maximum disk $`\mathrm{{\rm Y}}_{}`$ to decrease with increasing resolution. It is not clear what the cause of this effect is. The inner slopes of the rotation curve and stellar curve change subtly with resolution, but not enough to be the sole cause. It is likely that small differences between the tilted ring models also play a role.

We note that at all resolutions we find that (restricting ourselves to the ISO model for a moment) the min+gas model usually has the smallest value of $`\chi ^2`$, and that for all cases $`\chi _{m+g}^2<\chi _{con}^2`$. We have tried to constrain the value of $`\mathrm{{\rm Y}}_{}`$ by making a “best fit” to the rotation curves using ISO models: i.e. we also let $`\mathrm{{\rm Y}}_{}`$ be a free parameter in the fits. Unfortunately the different resolutions make it difficult to constrain $`\mathrm{{\rm Y}}_{}`$ in this manner. The B96 data is best fitted with a stellar component with $`\mathrm{{\rm Y}}_{}{}_{K}{}^{}=0.10\pm 0.13`$, while the B48 and B24 curves both demand $`\mathrm{{\rm Y}}_{}{}_{K}{}^{}=0.31\pm 0.04`$. For the B08 and B12 data such a “best fit” gives unphysical results (negative $`\mathrm{{\rm Y}}_{}`$). A constrained fit with $`\mathrm{{\rm Y}}_{}0`$ yields $`\mathrm{{\rm Y}}_{}=0`$ (i.e. the minimum disk+gas case) as best fit for these resolutions. Though it is clear that $`\mathrm{{\rm Y}}_{}`$ is difficult to constrain using only rotation curve information, the results do suggest that the stellar population is dynamically unimportant.

## 8 Discussion

### 8.1 Inner slopes

With the high-resolution rotation curves we can determine the inner slope of the mass-density distribution. The rationale behind this is described in de Blok et al. (2001b). In summary, the various halo models make specific predictions about the slope of the dark matter mass density distribution in the inner parts. If we approximate this distribution with a power-law $`\rho r^\alpha `$, then the ISO halo predicts $`\alpha =0`$, while the NFW halo predicts $`\alpha =1`$. In de Blok et al. (2001b) the inner mass density slopes of a large sample of LSB galaxies was derived, and the most representative value was $`\alpha =0.2\pm 0.2`$, i.e. closer to pseudo-isothermal than to CDM.

Here we derive the slopes of the mass density distribution of NGC 6822 at various resolutions using the inversion method described in de Blok et al. (2001b) where we compute the mass density profile and fit a power-law to the mass-density profile. Fig. 18 shows two examples. We show the B96 and B24 mass-density profiles, with the best fitting NFW and ISO (minimum disk) models over-plotted. The NFW models have central densities that are discrepant by a factor of 10 or more.

We measure the slope within $`R=0.8`$ kpc. This is an arbitrary choice, motived by the fact that we want to use at least 3 data points to determine the slope of the B96 curve. Other choices are obviously possible, and may result in slightly different values of the slope, but will not affect the conclusions. The resulting values are given in Table 5. For the B08 data a meaningful slope could not be derived due to the large gradients in the rotation velocity that are found at small scales. The values are at all resolutions significantly different from the NFW values. Note that there is no systematic steepening of the slope with resolution, which one would expect if the cusp was hidden by resolution effects. With a resolution of only tens of parsecs it is not clear that there is still room for a cusp in NGC 6822.

### 8.2 Dark matter, feedback and cosmology

In the previous sections we have made NFW fits to the rotation curves disregarding any of the predictions that CDM makes for $`c`$, $`V_{200}`$ and their relation. We will compare the observed rotation curve with cosmological predictions. If we assume that $`V_{200}V_{\mathrm{max}}55`$ $`\text{km}\text{s}^1`$, we find for $`\mathrm{\Lambda }`$CDM that $`c9.5`$ (Navarro, Frenk & White, 1997). The total mass $`M_{200}`$ of this halo is $`M_{200}=5.210^{10}`$ $`M_{}`$. This predicted halo is over-plotted on the observed curve in Fig. 19, and over-predicts the rotation velocity in the inner parts by a significant amount.

To look at the problem from a different perspective, let us assume $`\mathrm{{\rm Y}}_{}`$=0.35. This gives a stellar mass $`M_{}=1.3\times 10^8M_{}`$. The total (atomic) gas mass is $`1.5\times 10^8M_{}`$ (HI+He) and the total observed baryonic mass is $`2.8\times 10^8M_{}`$. A universal baryon fraction of $`0.09`$ (see e.g. White & Fabian 1995) implies a total mass $`M_{200}=3.1\times 10^9M_{}`$. For a NFW halo this implies $`V_{200}=21.5`$ $`\text{km}\text{s}^1`$ and $`R_{200}=28`$ kpc. These values are not consistent with the observed curve, and imply (apart from the inability of the model to fit the data) that either the baryon fraction in NGC 6822 differs significantly from the universal value or that large amounts of baryons have been expelled, as we will discuss below in more detail.

Let us work out the *observed* matter fractions in NGC 6822. We use the ISO models, and determine the amount of dark matter within a radius of 5 kpc, the outer edge of the Hi disk. This yields a total dark mass out to 5 kpc of $`3.2\times 10^9M_{}`$. With an observed baryonic mass of $`2.8\times 10^8M_{}`$ we find that NGC 6822 is heavily dark matter dominated with $`M_{vis}/M_{dark}(R<5\text{kpc})=0.09`$. This is comparable to values found in LSB galaxies and other dwarf galaxies (de Blok & McGaugh 1997). The corresponding baryon fraction is $`f_b=0.080`$, close to the universal baryon fraction.

The visible matter in NGC 6822 is thus a minor component of the total galaxy system. This makes it harder to explain the observed core-dominated dark matter distribution as resulting from a NFW halo modified by feedback. The term feedback is often used indiscriminately to indicate star formation and evolutionary processes that affect the mass distribution in a galaxy. In fact, there are two distinct forms of feedback. One of them is commonly observed in galaxies, also in NGC 6822, and usually shows itself as small-scale redistribution of mass in the disk.

The second form of feedback is more catastrophic, mostly theoretical, and was introduced to explain the discrepancy between observed and theoretical dark matter distributions. In essence it invokes large-scale and violent star formation, resulting in massive blow-outs that drag dark matter out due to gravitational interactions, thus destroying the primordial CDM cusp. The exact physics of feedback are not understood, and the models usually resorts to an empirical description that is fine-tuned to e.g. fit the Tully-Fisher relation or other observational constraints.

Nevertheless, let us explore the implications of feedback: using the amount of energy produced by supernovae and making some standard assumptions one can work out an expression for the amount of matter expelled by the effects of supernovae (see e.g. Eq. (A4) in van den Bosch et al. 2000). Feedback then requires the baryonic mass to drag along a similar amount of dark matter (modulo some efficiency factor that is unlikely to be much larger than unity).

Assuming CDM halos, we find for $`V_{200}55`$ $`\text{km}\text{s}^1`$ that $`M_{200}=5.2\times 10^{10}M_{}`$. Equation A4 in van den Bosch et al. (2000) then yields that the expelled baryonic mass is $`3.6\times 10^9M_{}`$. Assuming that the efficiency factor is unity (every unit mass of baryons drags along one unit mass of dark matter) we find that the amount of dark matter relocated is about equal to the amount of dark matter currently observed within the outermost radius of the observed disk. We are thus dealing with a major structural re-organisation of a galaxy.

Even if we assume the maximum possible maximum disk $`\mathrm{{\rm Y}}_{}`$=0.95 we still only find an observed baryonic mass of $`5.4\times 10^8M_{}`$ (though not physically motivated, this is a hard upper limit on the observed baryonic mass). Assuming an efficiency factor of unity, it seems NGC 6822 must have expelled over $`6`$ times its currently observed baryonic mass.

There has been some discussion in the past whether or not stellar feedback can remove a significant fraction of the mass of a galaxy to the Intergalactic Medium. Gnedin & Zhao (2002), as well as Mac Low & Ferrara (1998) showed that this process is unable to explain the observed cores, and is also inconsistent with other observational constraints. According to numerical simulation of stellar feedback by Mac Low & Ferrara (1998) the potential mass loss is a function of the total galaxy mass; for galaxies with gas masses $`<10^6M_{}`$ a galaxy can potentially destroy iself (‘blow-away’ scenario) whereas for gas masses between $`10^6`$ and $`10^9M_{}`$ some material may be removed (‘blow-out’, see their Fig. 1). According to these models we do not expect signifcant mass loss in NGC 6822 ($`M_{gas}=1.510^8`$) due to violent star formation.

From an observer’s perspective, there is relatively little evidence that mass loss does indeed occur in dwarf galaxies. Although strong outflow signatures are observed in some cases in dwarf galaxies there are only few cases where outflow material may permanently escape the galactic gravitational potential (e.g. Martin 1998).

A further problem with expelling a lot of mass are the global dynamics of the galaxy: if one were to have a huge mass loss of order a few times the total galaxy mass, one would certainly expect to find some clear evidence in the velocity field. The blown-out gas is such a case unlikely to end up as an undisturbed regularly rotating Hi disk. Furthermore, the stellar record shows no evidence for increases in the SFR with large factors at any time during the past $`9`$ Gyr (Gallart et al., 1996). It thus seems somewhat surprising that if feedback had occurred, it would have managed to hide all traces of a violent past.

In summary, we conclude that potential feedback of several times the total baryonic mass of NGC 6822 is unlikely. The observed rotation curve of NGC 6822, combined with the $`K`$-band and Hi data thus lead to a picture inconsistent with a combination of cuspy halo and feedback.

## 9 Conclusions

We obtained high resolution rotation curves for the local group dwarf irregular galaxy NGC 6822. By fitting both pseudo-isothermal and NFW halo models, it was found that the pseudo-isothermal model fits the data successfully. NGC 6822 is not maximum disk and is dominated by dark matter. The highest resolution rotation curve has some 250 independent points and beam-smearing effects are not an issue. The position of the dynamical center is unambiguously determined. There is no indication of a dark matter cusp down to scales of $`20`$ pc.

Trying to explain the observed curve with a combination of feedback and CDM halos leads to several inconsistencies. Feedback and CDM seem to imply that NGC 6822 has lost about 6 times the currently observed amount of baryons. There is however no evidence for such violent processes. We therefore conclude that (i) the dark matter distribution in NGC 6822 is unlikely to have been affected by feedback, and (ii) is best described by a model with a constant density core.

## Acknowledgements

We thank Steve Schneider for providing us with the 2MASS images. We thank Sébastien Blais-Ouellette for constructive comments that helped improve this paper. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.


format: text

text:
# 1 Lyapunov spectra

## 1 Lyapunov spectra

Recently, molecular dynamics simulations have been used to study many body systems representing simple fluids or solids from the point of view of dynamical systems theory. Due to the convex dispersive surface of the atoms, the phase trajectory of such systems is highly unstable and leads to an exponential growth, or decay, of small (infinitesimal) perturbations of an initial state along specified directions in phase space. This so-called Lyapunov instability is described by a set of rate constants, the Lyapunov exponents $`\{\lambda _l,l=1,\mathrm{},D\}`$, to which we refer as the Lyapunov spectrum. Conventionally, the exponents are taken to be ordered by size, $`\lambda _l\lambda _{l+1}`$. There are altogether $`D=2dN`$ exponents, where $`d`$ and $`D`$ denote the dimensions of space and of phase space, respectively, and $`N`$ is the number of particles. For fluids in nonequilibrium steady states close links between the Lyapunov spectrum and macroscopic dynamical properties, such as transport coefficients, irreversible entropy production, and the Second Law of thermodynamics, have been established . This important result provided the motivation for us to examine the spatial structure of the perturbed states associated with the various exponents. Here we present some of our results for two simple many-body systems representing dense two-dimensional fluids in thermodynamic equilibrium. The first model consists of $`N`$ hard disks (HD) interacting with hard elastic collisions, the second of $`N`$ soft disks interacting with a purely repulsive Weeks-Chandler-Anderson (WCA) potential.

The instantaneous state of a planar particle system is given by the $`4N`$-dimensional phase space vector $`𝚪=\{𝐫_i,𝐩_i,;i=1,\mathrm{},N\}`$, where $`𝐫_i`$ and $`𝐩_i`$ denote the respective position and linear momentum of molecule $`i`$. An infinitesimal perturbation $`\delta 𝚪=\{\delta 𝐫_i,\delta 𝐩_i;i=1,\mathrm{},N\}`$ evolves according to motion equations obtained by linearizing the equations of motion for $`𝚪(t)`$. For ergodic systems there exist $`D=4N`$ orthonormal initial vectors $`\{\delta 𝚪_l(0);l=1,\mathrm{},4N\}`$ in tangent space, such that the Lyapunov exponents

$$\lambda _l=\underset{t\mathrm{}}{lim}\frac{1}{t}\mathrm{ln}\frac{|\delta 𝚪_l(t)|}{|\delta 𝚪_l(0)|},l=1,\mathrm{},4N.$$
(1)
exist and are independent of the initial state . Geometrically, the Lyapunov spectrum describes the stretching and contraction along linearly-independent phase space directions of an infinitesimal hypersphere co-moving with the flow. For equilibrium systems the symplectic nature of the motion equations assures the conjugate pairing rule to hold : the exponents appear in pairs with a vanishing pair sum, $`\lambda _l+\lambda _{4N+1l}=0`$. Thus, only the positive half of the spectrum $`\{\lambda _{1l2N}\}`$ needs to be calculated. The sum of all Lyapunov exponents vanishes, which according to Liouville’s theorem is a manifestation of the conservation of phase volume by Hamiltonian systems. Six of the exponents, $`\{\lambda _{2N2l2N+3}\}`$, always vanish as a consequence of the conservation of energy, momentum, and center of mass, and of the non-exponential time evolution of a perturbation vector parallel to the phase flow.

For the computation of a complete spectrum a variant of the classical algorithm by Benettin et al. and Shimada et al. is commonly used . It requires following the time evolution of the reference trajectory and of an orthonormal set of tangent vectors $`\{\delta 𝚪_l(t);l=1,\mathrm{},4N\}`$, where the latter is periodically re-orthonormalized with a Gram-Schmidt (GS) procedure after consecutive time intervals $`\mathrm{\Delta }t_{GS}`$. The Lyapunov exponents are determined from the time-averaged renormalization factors. For the hard disk systems the free evolution of the particles is interrupted by hard elastic collisions and a linearized collision map needs to be calculated. An algorithm based on this idea has been developed by Dellago et al. . Although we make use of the conjugate pairing symmetry and compute only the positive branch of the spectrum, we are presently restricted to about 1000 particles by our available computer resources.

For our numerical work reduced units are used. In the case of the Weeks-Chandler-Anderson interaction potential,

$$\varphi (r)=\{\begin{array}{cc}4ϵ[(\sigma /r)^{12}(\sigma /r)^6]+ϵ,\hfill & r<2^{1/6}\sigma \hfill \\ 0,\hfill & r2^{1/6}\sigma .\hfill \end{array},$$
(2)
the particle mass $`m`$, the particle diameter $`\sigma `$, and the time $`(m\sigma ^2/ϵ)^{1/2}`$ are unity. In this paper we restrict our discussion to a thermodynamic state with a total energy per particle, $`E/N`$, also equal unity. For the hard-disk fluid, $`(Nm\sigma ^2/K)^{1/2}`$ is taken as the unit of time, where $`K`$ is the total kinetic energy, equal to the total energy $`E`$ of the system. There is no potential energy in this case. For the reduced temperature we have $`T=K/N`$, where Boltzmann’s constant is also unity. In the following, Lyapunov exponents for the two model fluids will be compared for equal temperatures (and not for equal total energy). This requires a rescaling of the hard-disk exponents by a factor of $`\sqrt{(K/N)_{WCA}/(K/N)_{HD}}`$ to account for the difference in temperature. All our simulations are for a reduced density $`\rho N/V=0.7`$, where the simulation box is a square with a volume $`V=L^2`$ and a side length $`L`$. Periodic boundaries are used throughout.

As an example, in Fig. 1 we compare the Lyapunov spectrum of a 100-particle WCA fluid to a spectrum for an analogous hard disk system at the same temperature ($`T=0.75`$) and density ($`\rho =0.7`$). A reduced index $`l/2N`$ is used on the abscissa. It is surprising that the two spectra differ so much in shape and size. The difference persists in the thermodynamic limit. The step-like structure displayed by the hard disk spectrum for $`l/2N`$ close to 1 is an indication of a coherent wave-like shape of the associated perturbations, to which we refer as Lyapunov modes . We defer the discussion of these modes to Section 4.

## 2 Time-dependent local exponents

We infer from Equ. (1) that the Lyapunov exponents are time averages over an (infinitely) long trajectory and, therefore, are global properties of the system. They can be rewritten as

$$\lambda _l=\underset{\tau \mathrm{}}{lim}_0^\tau \lambda _{}^{}{}_{l}{}^{}(𝚪(t))𝑑t\lambda _{}^{}{}_{l}{}^{},$$
(3)
where the (implicitly) time-dependent function $`\lambda _{}^{}{}_{l}{}^{}(𝚪)`$ depends on the state $`𝚪(t)`$ which the system occupies in phase space at time $`t`$. Thus, $`\lambda _{}^{}{}_{l}{}^{}(𝚪)`$ is called a local Lyapunov exponent. It may be estimated from

$$\lambda _{}^{}{}_{l}{}^{}(𝚪(t))=\frac{1}{\mathrm{\Delta }t_{GS}}\mathrm{ln}\frac{|\delta 𝚪_l(𝚪(t+\mathrm{\Delta }t_{GS})|}{|\delta 𝚪_l(𝚪(t)|},$$
(4)
where $`t`$ and $`t+\mathrm{\Delta }t_{GS}`$ refer to times immediately after consecutive Gram-Schmidt re-orthonormalization steps. Its time average along a trajectory is denoted by $`\mathrm{}`$ and gives the global exponent $`\lambda _l`$.

The local exponents fluctuate considerably along a trajectory. This is demonstrated in Fig. 2, where we have plotted the second momens $`\lambda _{}^{}{}_{l}{}^{2}`$ as a function of the Lyapunov index $`l,1l4N`$, for a system of 16 particles, both for the WCA and HD models. $`l=1`$ refers to the maximum, and $`l=64`$ to the most-negative exponent. The points for $`30l35`$ correspond to the 6 vanishing exponents and are not shown. We infer from this figure that for the soft WCA particles the fluctuations of the local exponents become very large for the Lyapunov exponents describing relatively-weak instabilities with near-zero growth rates, $`l2N`$. For the hard disk system, however, the relative importance of the fluctuations becomes minimal for the same exponents. We shall return to this point in Section 4.

We note that the computation of the second moments $`\lambda _{}^{}{}_{l}{}^{2}`$ for the hard disk system requires some care. Due to the hard core collisions they depend strongly on $`\mathrm{\Delta }t_{GS}`$. The variance of the fluctuating local exponents, $`\lambda _{}^{}{}_{l}{}^{2}\lambda _{}^{}{}_{l}{}^{}^2`$, varies with $`1/\mathrm{\Delta }t_{GS}`$ for small $`\mathrm{\Delta }t_{GS}`$, as is demonstrated in Fig. 3 for $`l=1`$. However, the shape of the fluctuation spectrum is hardly affected by the size of the renormalization interval.

## 3 The maximum exponent

The maximum Lyapunov exponent is the rate constant for the fastest growth of a phase space perturbation in a system. Thus, it is dominated by the fastest dynamical events in the fluid, binary collisions. There is strong numerical evidence for the existence of the thermodynamic limit $`\{N,V\mathrm{},N/V`$constant$`\}`$ for $`\lambda _1`$ and, hence, for the whole spectrum . Furthermore, the perturbation belonging to $`\lambda _1`$ is strongly localized in space . This may be demonstrated by projecting the tangent vector $`\delta 𝚪_1`$ onto the four-dimensional subspaces spanned by the perturbation components of the individual particles. The squared norm of this projection, $`\delta _i^2(t)(\delta 𝐫_i)_l^2+(\delta 𝐩_i)_l^2`$, indicates how active a particle $`i`$ is engaged in the growth process of the pertubation associated with $`\lambda _1`$ .

In Fig. 4 $`\delta _i^2(t)`$ is plotted along the vertical axis for all particles of a hard disk system at the respective positions $`(x_i,y_i)`$ of the disks in space, and the ensuing surface is interpolated over a periodic grid covering the simulation box. The figure refers to an instantaneous condition for a well-relaxed system, and no averaging is involved. A strong localization of the active particles is observed at any instant of time. Similar, albeit slightly broader peaks are observed for the WCA system or other soft disk potentials.

This localization is a consequence of two mechanisms: firstly, after a collision the delta-vector components of two colliding molecules are linear functions of their pre-collision values and have only a chance of further growth if their values before the collision were already far above average. Secondly, each renormalization step tends to reduce the (already small) components of the other non-colliding particles even further. Thus, the competition for maximum growth of tangent vector components favors the collision pair with the largest components. The active zone moves in space in a diffusive manner.

The localization also persists in the thermodynamic limit. To show this we follow Milanović et al.  and order all squared components $`[\delta 𝚪_1]_j^2;j=1,\mathrm{},4N`$ of the perturbation vector $`\delta 𝚪_1`$ according to size. By adding them up, starting with the largest, we determine the smallest number $`A`$ of terms, required for the sum to exceed a threshold $`\mathrm{\Theta }`$. Then, $`C_{1,\mathrm{\Theta }}A/4N`$ may be taken as a relative measure for the number of components actively contributing to $`\lambda _1`$:

$$\mathrm{\Theta }\underset{s=1}{\overset{4NC_{1,\mathrm{\Theta }}}{}}\left[\delta 𝚪_1\right]_s^2,[\delta 𝚪_1]_i^2[\delta 𝚪_1]_j^2\mathrm{for}i<j.$$
(5)
Obviously, $`C_{1,1}=1`$.

In Fig. 5 $`C_{1,\mathrm{\Theta }}`$ is shown for $`\mathrm{\Theta }=0.90`$ as a function of the particle number $`N`$, both for the WCA fluid and for the hard disk system. It converges to zero if our data are extrapolated to the thermodynamic limit, $`N\mathrm{}`$. This supports our assertion that in an infinite system only a vanishing part of the tangent-vector components (and, hence, of the particles) contributes significantly to the maximum Lyapunov exponent at any instant of time.

## 4 Lyapunov modes

We have already mentioned the appearance of a step-like structure in the Lyapunov sepctrum of the hard disk system for the positive exponents closest to zero, which was first observed in Ref. . They are a consequence of coherent wave-like spatial patterns generated by the perturbation vector components associated with the individual particles .

In Fig. 6 this is visualized by plotting the perturbations in the $`x`$ (bottom surface) and $`y`$ directions (top surface), $`\{\delta x_i,i=1,\mathrm{},N\}`$ and $`\{\delta y_i,i=1,\mathrm{},N\}`$, respectively, along the vertical axis at the instantaneous particle positions $`(x_i,y_i)`$ of all particles $`i`$. This figure depicts a transverse Lyapunov mode of type $`T_1`$, for which the perturbation components are perpendicular to the respective wave vectors with one wave length equal to $`L`$, the linear extension of the simulation box. The system consists of $`N=1024`$ hard disks, and the perturbation vector considered in the figure belongs to the smallest positive exponent $`\lambda _{2045}`$. An analogous plot of $`\delta p_x`$ and $`\delta p_y`$ for the same perturbation vector is identical in shape to that of $`\delta x`$ and $`\delta y`$ in Fig. 6, with the same phases of the waves. This is a consequence of the fact that the perturbations are solutions of linear first-order equations instead of second order equations such as the wave equation. Furthermore, the exponents for $`2042l2045`$ are equal. This four-fold degeneracy of non-propagating transversal modes, and an analogous eight-fold degeneracy of propagating longitudinal modes, are responsible for a complicated step structure for $`l`$ close to $`2N`$, which has been studied in detail in Refs. .

The wave length of the modes and the value of the corresponding exponents are determined by the size $`L`$ of the square simulation box. There is a kind of linear dispersion relation  indicating that the smallest positive exponent is proportional to $`1/L`$. This assures that for a square simulation box with aspect ratio 1 there is no positive lower bound for the positive exponents of a hard disk system in the thermodynamic limit .

So far, our discussion of modes applies only to hard disk fluids. In spite of a considerable computational effort we have not yet been able to indentify modes for two-dimensional fluid systems with a soft interaction potential such as WCA or similar potentials . This surprising result seems to be due to the very strong fluctuations of the local exponents discussed in Section 2. They tend to obscure any mode in the system in spite of considerable averaging and make a positive identification very difficult.

We are grateful to Christoph Dellago, Robin Hirschl, Bill Hoover, and Ljubomir Milanović for many illuminating discussions. This work was supported by the Austrian Fonds zur Förderung der wissenschaftlichen Forschung, Grants P11428-PHY and P15348-PHY.


text:
# Cooling non–linear lattices toward energy localisation

## I Introduction

Nonlinearity has revealed one of the key ingredients for describing many relevant features of different states of matter. In the realm of lattice dynamics, scattering processes among phonons, propagation of solitary waves and slow energy relaxation are typical examples. Recently, considerable efforts have been devoted to the study of periodic, localised, non–linear lattice excitations, named ”breathers” breath1  . They are quite peculiar (but generic) objects emerging from the interplay of nonlinearity and space discreteness breath2 . Their mathematical properties, like existence, stability, mobility etc. have been progressively unveiled, while they have been identified in many different scenarios.

Their role in non–equilibrium dynamics seems to be particularly fascinating. An example is the relaxation to energy equipartition of short–wavelength fluctuations chao\_breath . Another interesting scenario where breathers are found to emerge spontaneously is observed upon cooling the lattice at its boundaries Aubry1 ; Aubry2 ; noi , i.e. by considering a non–equilibrium process in which energy exchange with the environment is much faster at the surface than in the bulk. The numerical simulations neatly show that the energy dissipation rate may be significantly affected by the spontaneous excitation of breathers. As the latter exhibit a very weak interaction among themselves and with the boundaries, the energy release undergoes a sudden slowing down and is hardly detected on the time scales of a typical simulation. Thus, the lattice remains frozen in a pseudo–stationary, metastable configuration which is far from thermal equilibrium, which we shall refer to as residual state. Despite the absence of disorder, there is a close similarity with the glassy behaviour observed in disordered systems, as was already pointed out Aubry1 ; Aubry2 .

But what are the mechanisms leading to spontaneous localisation? As we have shown in a previous paper noi , the latter is intimately related to how dissipation acts on vibrational modes of different wavelengths. If long–wavelength phonons can be efficiently damped out, the modulation instability of short lattice waves Peyrard  becomes highly favoured and breathers can easily emerge from an interacting soliton–gas chao\_breath ; Kosevich . In this respect, this type of non–equilibrium condition is much more effective in exciting localised modes than an equilibrium one (see e.g. Bishop  for a related discussion).

Spontaneous localisation upon cooling has been observed for both on–site Aubry1 ; Aubry2  and pure nearest–neighbour (n.n.) interactions noi ; reigada1 . The two classes of models are known to behave differently, e.g. as to the mobility of breathers, which is much higher in the absence of local coupling. This is also confirmed by relaxations experiments, where differences in the energy decay laws have been observed. In particular, the approach to the residual state follows a simple exponential law in the Fermi-Pasta-Ulam (FPU) system noi , but turns to stretched exponential in the presence of local coupling Aubry2 .

Our aim in this paper is to provide an up–to–date overview of the phenomenon of breather localisation by cooling. Generalities about the energy release process, including models and indicators, will be presented in Section II, along with the discussion of the cooling process in an harmonic lattice and the results obtained for 1D non–linear systems. In particular we shall survey the different cooling pathways appearing in pure n.n. and on–site potentials. A specific subsection will be devoted to clarify the interpretation of some recent results reigada1 ; reigada2 . In Section III, we present the results for the n.n. 2D FPU model. We want to point out that the 2D is not a straightforward extension of the 1D case. Indeed, the existence of an energy activation threshold for breather solutions 2D\_tresh  leads to conjecture that a thermalised state may be cooled to the residual state only above some initial energy/temperature. Moreover, the residual state is a static multi–breather state, whereby in the 1D case it usually contains a single breather. Finally, we perform an analysis of the distribution of breather energies, based on a simple statistical model.

## II  Relaxation and localisation in 1D

In this section we will describe the lattice models and the generalities of a typical relaxation experiment as well as the relevant indicators used throughout the paper. Let us consider a chain of $`N`$ atoms of unit mass and denote with $`u_p`$ the displacement of the $`p`$–th particle from its equilibrium position $`pa`$ ($`a`$ is the lattice spacing). The atoms are labelled by the index $`p=0,1,\mathrm{}N1`$ and their dynamics is given by the equations of motion

$$\begin{array}{cc}\hfill \ddot{u}_p=& V^{}(u_{p+1}u_p)V^{}(u_pu_{p1})U^{}(u_p)\hfill \\ & \gamma \dot{u}_p\left[\delta _{p,0}+\delta _{p,N1}\right],\hfill \end{array}$$

where $`V(x)`$ and $`U(x)`$ are the interparticle and on–site potentials, respectively. We assume that $`V^{}(0)=U^{}(0)=0`$, denoting with a prime the derivative with respect to $`x`$. For convenience, we shall adopt non–dimensional units such that $`a`$ and $`V^{\prime \prime }(0)`$ are set to unity. Moreover, as we want to deal with systems in a finite volume, we will impose either free–end ($`u_1=u_0`$, $`u_{N1}=u_N`$) or fixed–end ($`u_1=u_N=0`$) boundary conditions (BC).

The last term in eq. (II) represents the interaction of the atoms with a “zero–temperature” heat bath in the form of a linear damping with rate $`\gamma `$. Following Aubry1 ; Aubry2 , we restricted to the case in which dissipation selectively acts only on the atoms located at the chain edges. As mentioned above, this choice is crucial and should model a physical situation in which energy exchange with the environment is much faster at the surface than in the bulk.

The general layout of a simulation can be summarised as follows . First, an equilibrium microstate is generated by, say, Nosè-Hoover (canonical) method Aubry1  or by letting the Hamiltonian (microcanonical) system ($`\gamma `$ = 0) evolve for a sufficiently long transient noi . In the following we will follow the second strategy. The initial condition for the transient is assigned by setting all displacements $`u_p`$ to zero and by drawing velocities at random from a Gaussian distribution. The velocities are then rescaled by a suitable factor to fix the desired value of the energy density (energy per particle) $`e_0`$. The resulting set of $`u_p`$ and $`\dot{u}_p`$ is then used as initial condition to integrate Eqs. (II) with $`\gamma >0`$ (see again Ref. noi  for details).

The result of the relaxation process is the residual state, i.e. a state with a finite fraction of the initial energy stored in the form of breathers. Remarkably, such residual state is characterised by a decay time scale which is several orders of magnitude longer than that of the localisation process.

### II.1  The harmonic lattice

A number of features which characterise the relaxation dynamics of non–linear lattices can be better understood by first studying the harmonic lattice, namely

$$V(x)=\frac{1}{2}x^2,$$

We anticipate that the presence of an harmonic on–site potential does not alter the main conclusions. We thus take $`U=0`$ throughout this section.

For small damping, an approximate analytical solution can be found by time–dependent perturbation theory. Let $`\omega _\alpha ^2=4\mathrm{sin}^2(q_\alpha /2)`$ and $`\eta ^\alpha `$ ($`\alpha =0,1,\mathrm{}N1`$) denote the eigenvalues and normalized eigenvectors of the unperturbed Hamiltonian problem, where the allowed wave–numbers are $`q_\alpha =\alpha \pi /N`$ and $`q_\alpha =(\alpha +1)\pi /(N+1)`$ for free–end and fixed–end BC, respectively. To the lowest order in $`\gamma `$ one obtains noi

$$u_p(t)=\underset{\alpha =0}{\overset{N1}{}}c_\alpha e^{\left({\displaystyle \frac{1}{\tau _\alpha }}+i\omega _\alpha \right)t}\eta _p^\alpha .$$
(1)
where $`\eta _p^\alpha `$ is the $`p`$-th component of the eigenvector $`\eta ^\alpha `$ and $`c_\alpha `$ are constant amplitudes fixed by the initial conditions. Notice that, since the number of the damped particles is fixed, the perturbative approximation is expected to improve upon increasing the system size $`N`$. The wavenumber-dependent damping rates are found to be

$$\frac{1}{\tau _\alpha }=\{\begin{array}{cc}\frac{1}{\tau _0}\mathrm{cos}^2\left(\frac{q_\alpha }{2}\right)\hfill & \text{for free–end BC}\hfill \\ \frac{1}{\tau _0}\mathrm{sin}^2(q_\alpha )\hfill & \text{for fixed–end BC}.\hfill \end{array}$$
(2)
with $`\tau _0=N/2\gamma `$. The free–end and fixed–end BC systems thus show considerably different behaviours. In the former case, the least damped modes are the short–wavelength ones ($`\alpha N`$), the largest lifetime being $`\tau _{N1}2N^3/\pi ^2\gamma `$, while the most damped modes are the ones in the vicinity of $`\alpha =0`$, with $`\tau _0`$ being the shortest decay time. On the contrary, for fixed ends the most damped modes are those around the band center ($`\alpha N/2`$) while short– and long–wavelength ones dissipate very weakly, being $`\tau _{N1}N(N+1)^2/2\pi ^2\gamma `$. As we shall illustrate in the following, such difference is crucial for the localisation in the non–linear system.

Using the solution (1), we evaluated the chain energy in the case of equipartition i.e. by replacing $`c_\alpha ^2\omega _\alpha ^2`$ with its average value $`2e_0`$. Finally, recalling eqs. (2) and approximating for large $`N`$ the sum over $`\alpha `$ with an integral ($`\tau _\alpha \tau (q)`$), we obtained noi

$$\frac{E(t)}{E(0)}=\frac{1}{\pi }_0^\pi e^{2t/\tau (q)}𝑑q=e^{t/\tau _0}I_0\left(\frac{t}{\tau _0}\right)=\{\begin{array}{cc}e^{t/\tau _0}\hfill & \mathrm{for}t\tau _0,\hfill \\ \frac{1}{\sqrt{2\pi \left(t/\tau _0\right)}}\hfill & \mathrm{for}t\tau _0.\hfill \end{array}$$
(3)
where $`I_0`$ is the modified zero–order Bessel function, whose asymptotic expansions have been used to obtain the last equality. We conclude that relaxation is exponential only at short times while a crossover to a power–law decay occurs at $`t\tau _0`$ . Moreover, despite the differences in the relaxation times (2), the decay law for $`E`$ turns out to be the same in the free and fixed–end systems.

The above formulas must be taken with caution when dealing with a large but finite system. In such a case, there is of course a lower cutoff at wavenumbers of order $`\pi /N`$. Therefore, after a time of the order of the lifetime of the longest–lived mode, $`\tau _{N1}N^3/\gamma `$, formula (3) no longer holds and a further crossover to the exponential decay law $`\mathrm{exp}(2t/\tau _{N1})`$ occurs. This has been also verified numerically noi .

### II.2 Non–linear lattices: the role of discreteness

Since the first results of relaxation experiments in non–linear lattices were reported Aubry1 ; Aubry2 , there has been some debate as to the nature of the decay law of the system energy and its relation with the spontaneously emerging localised modes. In particular, it has been claimed that the phenomenon of energy pinning in the form of discrete breathers would result in a glassy–like relaxation, i.e. stretched exponential decay law. Here we shall describe how the main role in determining the relaxation properties of a non–linear lattice is indeed played by the type of potential energy. In particular, we will highlight the importance of the relative strength of the inter–particle and on–site potentials, i.e. the degree of discreteness of the system. In order to do so, we study the general class of non–linear lattices obeying the equations of motion (II). In particular, we shall report here the results of numerical simulations performed with

$$V(x)=\frac{1}{2}x^2+\frac{1}{4}x^4U(x)=\frac{1}{4}\kappa x^4.$$
(4)
In this case, the relative strength of local and inter–particle non–linearities is accounted for by the parameter $`\kappa `$.

In Fig. 1 we show the space–time contour plots of the symmetrised site energies, defined as

$$h_p=\frac{1}{2}\dot{u}_p^2+\frac{1}{2}\left[V(u_{p+1}u_p)+V(u_pu_{p1})\right]+U(u_p).$$

The instantaneous total energy of the system is then given by $`E=_{p=0}^{N1}h_p`$. We see that a single localised excitation emerges from the relaxation process for both small and large values of $`\kappa `$. The pathway to localisation is the same as the one observed in the FPU model noi . In particular, the basic mechanism leading to localisation is modulational instability of short wavelength modes. The latter can only be effective if dissipation of long wavelength phonons is fast enough. As it shows from Eq. (2), this occurs only in the case of free–ends BC, whereas fixed–end BC strongly inhibit such process. On the other hand, breather mobility is strongly reduced in the highly discrete system (Fig. 1 (b)). Such difference is even more evident in the decay law of the normalized energy $`E(t)/E(0)`$, which turns from pure exponential to stretched–exponential behaviour upon increasing $`\kappa `$. This is better appreciated by plotting the indicator

$$𝒟(t)=\mathrm{ln}[E(t)/E(0)]$$
(5)
in a log–log scale, so that a stretched–exponential law of the form $`E(t)=E(0)\mathrm{exp}[(t/\tau )^\sigma ]`$ becomes a straight line with slope $`\sigma `$ that intercepts the $`y`$–axis at $`\sigma \mathrm{ln}\tau `$ (Fig. 2 (a)). The scaling regions correspond to the onset of energy localisation, as can be seen by plotting the localisation parameter

$$(t)=N\frac{{\displaystyle \underset{p}{}}h_p^2(t)}{\left[{\displaystyle \underset{p}{}}h_p(t)\right]^2}.$$
(6)
(Fig. 2 (b)). From its very definition, the fewer sites the energy is localised onto, the closer $``$ is to $`N`$ . On the other hand, the more evenly the energy is spread among all the particles, the closer $``$ is to a constant of order 1. It should be noted, that $``$ is a relative quantity, and bears no information on the amount of energy that is localised. As a matter of fact, the fraction of $`e_0`$ that gets trapped in the residual state is found to be greater the higher the degree of discreteness of the system.

To summarise, if the on–site force is weak enough, the behavior is basically the same as the FPU chain which, in turn, exhibits in the approach to the residual state the same exponential decay law of its linearised counterpart, i.e. $`E(t)/E(0)\mathrm{exp}(t/\tau _0)`$ (see Eqs. (3)). The reasons for such behaviour are the same discussed for the FPU model noi . First of all, an effective harmonic Hamiltonian with energy–dependent renormalized frequencies is known to account for several equilibrium properties of the FPU chain Alabiso ; Lepri . Second, the harmonic approximation becomes increasingly accurate as time elapses, simply because more and more energy is extracted from the system by the reservoir. On the contrary, if the system is highly discrete, the first relaxation stage is indeed described by a stretched exponential law. A quantitative explanation of it is still lacking, but it is clear that in this case an effective description of such genuine non–linear phenomenon in terms of quasi–harmonic modes breaks down. In other words, the resulting interactions among linear modes inhibit the energy flow from the bulk to the boundaries, thus slowing down the dissipation.

The onset of a pseudo–stationary state corresponds in Fig. 2 to the almost flat final portions of the curves. Actually, being the system globally dissipative, energy is still at that stage exponentially decreasing but with a huge time constant $`\tau _b`$. The latter is expected to be inversely proportional to the breather amplitude at the lattice edges, which, in turn, should be of order $`\mathrm{exp}(N/\mathrm{})`$, with $`\mathrm{}N`$ being the localisation length.

The above scenario is generic for large enough systems. Nonetheless, failure to spontaneously localise energy occurs the more frequently the smaller is the lattice. This is because the probability of occurrence of large enough energy fluctuations at a given energy density decreases with the number of particles. In other words, given a number of different initial equilibrium conditions, only a fraction $`n`$ of them will give rise to breathers. To give a quantitative idea, for the FPU model with $`e_0=1`$, $`n`$ is only about 50 % for a chain of $`N=20`$ particles but rapidly approaches 100 % already for $`N100`$. In the absence of spontaneous localisation, the non–linear systems show no difference with respect to the linear ones and $`E(t)`$ decays according to Eqs. (3).

### II.3 The effective–exponent analysis

We have repeatedly mentioned that for the FPU model no stretched exponential relaxation is observed, even in presence of a weak on–site force. In this section we wish to remark this statement by reporting a more detailed analysis of the energy decay. This can be accomplished by defining an effective exponent $`\zeta (t)`$ through the logarithmic derivative

$$\zeta (t)=\frac{d[\mathrm{ln}𝒟(t)]}{d[\mathrm{ln}t]}.$$
(7)
It is clear from definition (7) that if a portion of the energy decay curve goes as $`E(t)=E(0)\mathrm{exp}[(t/\tau )^\sigma ]`$ this would result in a plateau of height $`\sigma `$ in the corresponding $`\zeta (t)`$ curve. For the problem at hand, interpretation of numerical data with this indicator is however a delicate matter. Actually, we want to remark that a naive analysis can lead to misleading or even incorrect conclusions reigada1 ; reigada2 .

The most convincing way to illustrate this is to discuss the case of the harmonic chain where localisation and stretched-exponential behaviour are obviously excluded a priori. From Eq. (3), it is easily seen that $`\zeta `$ should start from 1 and monotonically vanish as $`1/\mathrm{ln}t`$ for $`t\tau _0`$. On the other hand, we already learnt that in a finite chain a further crossover to an exponential law $`\mathrm{exp}(2t/\tau _{N1})`$ occurs so that $`\zeta `$ must approach again 1 for $`t\tau _{N1}`$. The net results of those competing trends is that, for a finite $`N`$, $`\zeta `$ displays a broad minimum $`\zeta _{min}`$, at $`t\tau _{N1}`$ . The value $`\zeta _{min}`$ (which is independent of $`\gamma `$) can be estimated by noting that from Eqs. (7) and (3) one has $`\zeta (t)1/\mathrm{ln}[2\pi t/\tau _0]`$ for $`t\tau _0`$, and hence $`\zeta _{min}1/2\mathrm{ln}[4N/\pi ]`$. A simple plot of the curves shows how the the minimum $`\zeta _{min}`$ can be incorrectly interpreted as a plateau if observed on a too short time scale. Of course, the same scenario is likely to appear in relatively short FPU chains whenever the chosen initial conditions do not yield localisation (see Fig. 3).

When localisation does occur, the residual state decays exponentially noi . Therefore one expects $`\zeta `$ to converge to 1 as the residual state is approached. Actually, one observes first very small values of $`\zeta `$, while such a convergence does occur only on much longer time scales. To see why, let us consider a one–breather residual state. The total energy decays as $`E(t)=E_be^{t/\tau _b}`$, where $`E_b`$ is the initial energy of the breather and $`\tau _b`$ its characteristic decay time. The time constant $`\tau _b`$ is huge: in fact, it is roughly inversely proportional to the breather amplitude at the edge sites of the chain. From definitions (7) and (5) one has

$$\zeta =\frac{t/\tau _b}{t/\tau _b\mathrm{ln}[E_b/E(0)]}.$$
(8)
Since $`\mathrm{ln}[E_b/E(0)]`$ is a number of order 1, we see that, at times such that $`t\tau _b\mathrm{ln}[E_b/E(0)]`$, $`\zeta `$ is very small. The numerical results in this case are just a deceptive byproduct of the wrong normalization. The exponent $`\zeta `$ will start converging to 1 only when $`t\tau _b\mathrm{ln}[E_b/E(0)]`$, which is a huge time, being also $`E_bE(0)`$. In this case the $`\zeta `$ analysis may be misleading in identifying the nature of the decay law. The way around such problem would be to calculate the effective exponent by normalizing the total energy to $`E_b`$ . However, this solution is in turn complicated by the intrinsic uncertainty in locating the time origin for the decay of the breather state.

## III  2D lattices

Since a careful study of the 2D Klein-Gordon model has already been presented Aubry2 , we focus here on the 2D FPU model. In this respect, the most natural extension of model (II) should involve a two–component displacement vector. For the sake of simplicity, we rather consider here a scalar model with only one degree of freedom $`u_{i,j}`$ per lattice site. Actually, a preliminary series of simulations of the 2D vector version showed no appreciable differences with what reported here, and will be discussed elsewhere.

With the same choice of non–dimensional units introduced in section II, the model on a $`N\times N`$ lattice ($`i,j=0,\mathrm{},N1`$) with damping on all edges is defined by the equations of motion

$$\begin{array}{cc}\hfill \ddot{u}_{i,j}=& V^{}(u_{i+1,j}u_{i,j})V^{}(u_{i,j}u_{i1,j})+\hfill \\ & V^{}(u_{i,j+1}u_{i,j})V^{}(u_{i,j}u_{i,j1})\underset{p,q=0}{\overset{N1}{}}\mathrm{\Gamma }_{i,j}^{p,q}\dot{u}_{p,q},\hfill \end{array}$$

where $`\mathrm{\Gamma }_{i,j}^{p,q}=\gamma [g_{i,p}\delta _{j,q}+\delta _{i,p}g_{j,q}g_{i,p}g_{j,q}]`$, with $`g_{i,p}=\delta _{i,p}[\delta _{p,0}+\delta _{p,N1}]`$. Since localisation has been shown to be strongly inhibited by fixed–ends boundary conditions, we consider here free–ends BC.

Calculations identical to the one described in section II.1 can be easily extended to the case of a simple 2D harmonic lattice. It turns out that the contributions from the two spatial coordinates are identical for an $`N\times N`$ lattice, and therefore factorize, yielding

$$\frac{E(t)}{E(0)}=\left[e^{t/\tau _0}I_0\left(\frac{t}{\tau _0}\right)\right]^2\{\begin{array}{cc}e^{2t/\tau _0}\hfill & \text{for }t\tau _0\hfill \\ \frac{1}{2\pi \left(t/\tau _0\right)}\hfill & \text{for }t\tau _0\hfill \end{array},$$
(9)
with $`\tau _0=N/2\gamma `$. In particular, the law (9) applies in the absence of localisation, as well as during the first stage before the onset of localisation. In Fig. 4 we plot the energy decay curves of an FPU lattice with $`N=32`$ relaxing from different values of the energy density $`e_0`$. The first stage of the decay indeed follows the theoretical prediction Eq. (9). Moreover, it is clear that localisation does not occur below a finite value of $`e_0`$ (located between $`e_0=0.2`$ and $`e_0=0.3`$ in this case). In particular, in the absence of localisation the energy decay is again described by the two–crossover scenario. Notice that, also in this case, failure to spontaneously localise energy results in a minimum of the effective exponent $`\zeta `$, due to the competition of the power law decay $`1/t`$ and the exponential law $`\mathrm{exp}(2t/\tau N1)`$.

### III.1 Fluctuation–activated localisation

The relaxation process described in terms of the energy decay curves is identical to the two–crossover picture which applies to the 1D FPU case. Nonetheless, the scenario markedly differs from the 1D case in two respects. First of all, the mobility of the localised excitations which spontaneously emerge is drastically reduced. In fact, after a first stage in which strong interaction is observed, they arrange on a “random lattice”, which, on the time scale of our typical simulations, appears indeed to be frozen. One of such states is illustrated in Fig. 5 (a). This is presumably due to a smaller “scattering section” in 2D. Obviously, due to the presence of dissipation, this state decays on a much longer time scale, which (as discussed above) scales exponentially with the linear size of the lattice.

The second important difference with the 1D case stems from the appearance of a finite energy threshold $`\mathrm{\Delta }`$ for the existence of breathers (akin to envelope solitons in the small–amplitude limit) in 2D 2D\_tresh . As a consequence, we expect that localised modes are generated in the relaxation dynamics only by fluctuations that are large enough to overcome such threshold. The spontaneous excitation of breathers can thus be seen as an activated process. Accordingly, their number will be exponentially small in the ratio between $`\mathrm{\Delta }`$ and some quantity measuring the strength of fluctuations. Therefore, one expects the average density of breathers (i.e. the average number of breathers per lattice site) in the residual state to follow an Arrhenius law of the form

$$n_B\mathrm{exp}(\beta \mathrm{\Delta }).$$
(10)
From Eq. (10), it is tempting to identify $`\beta `$ with some inverse temperature. However, one has to keep in mind that we are dealing with a non–equilibrium process and this identification can only make sense if energy release is adiabatically slow. If this is true, $`1/\beta `$ should be proportional to (and smaller than) the initial temperature, which, in turn, is roughly proportional to the initial energy density.

In our numerical simulations we clearly observe that $`n_B`$ strongly depends on the initial specific energy. In Fig. 5 (b) we plot $`n_B`$ for two system sizes as a function of $`e_0`$, along with a fit performed with the law $`n_B=C\mathrm{exp}(\mathrm{\Delta }^{}/e_0)`$, with $`\mathrm{\Delta }^{}\mathrm{\Delta }`$. As it shows, the agreement is good. We conclude that the spontaneous localisation of energy in our system is indeed an activated process. In particular, the fit gives $`\mathrm{\Delta }^{}0.9`$. We note that a quantitative analysis of the time $`t_0`$ required to reach the residual state (localisation time) reveals that $`t_01/n_B\mathrm{exp}(\beta \mathrm{\Delta })`$ Esc\_proc . By fitting the numerical data for the localisation time, it is possible to get an independent confirmation of the thermal activation scenario. In particular, we get $`\mathrm{\Delta }^{}0.7`$, in good agreement with the value obtained above.

### III.2 Statistics of breathers’ energies

It is important to realize that the “pseudo–stationary” distribution of breathers in the residual state emerges as a result of two competing mechanisms, namely the birth out of fluctuations from the homogeneous state and the decay due to both the coupling with the external reservoir and to the various inelastic interactions (phonon–breather, breather–breather etc). Although the former can be safely neglected, it is very difficult to understand the latter in full detail. Indeed, to the best of our knowledge, only a few studies are available Johansson ; Bang  and a full description in terms of elementary interactions seems hardly feasible. Nonetheless, simple statistical arguments can be of help in understanding the energy distribution to a greater detail. Let us consider the fraction of breathers $`𝒫(ϵ,t)`$ having energy between $`ϵ`$ and $`ϵ+dϵ`$ with $`ϵ\mathrm{\Delta }`$ (in the FPU model breathers of arbitrarily large energy exist due to the unboundedness of the interaction potential) and let us define $`B(ϵ)`$ and $`D(ϵ)`$ to be some phenomenological birth and decay rates, respectively. Furthermore, let us assume that birth of a breather can only occur by a spontaneous fluctuation from the ground state ($`ϵ=0`$). Stationarity requires the flux–balance condition

$$B(ϵ)𝒫(0)=D(ϵ)𝒫(ϵ)$$
(11)
Since the birth is an activated process we expect $`B(ϵ)\mathrm{exp}[\beta ϵ]`$ <sup>1</sup><sup>1</sup>1In principle, one could have a further energy–dependent prefactor in front of the exponential. That would affect the form of the distribution. We neglect it for simplicity. . This assumption yields (for a constant $`𝒫(0)`$)

$$𝒫(ϵ)=\frac{𝒫(0)}{D(ϵ)}\mathrm{exp}(\beta ϵ)$$
(12)
In this simplified description, the prefactor of the exponential term is thus interpreted as a measure of the effective breather lifetime $`\tau (ϵ)`$ ($`𝒟(ϵ)1/\tau (ϵ)`$) and is a priori unknown. In all simulations we observe that small–amplitude breathers decay more easily. This observation agrees with the accepted existence of a preferred energy flow from small–amplitude breathers to large–amplitude ones in breather–breather collisions Bang . Moreover, Eq. (12) requires that $`\tau (\mathrm{\Delta })=0`$. Thus, we postulate

$$D(ϵ)(ϵ\mathrm{\Delta })^z,$$
(13)
where the exponent $`z>0`$ can be estimated by measuring the distribution of breather energies in the vicinity of the threshold $`\mathrm{\Delta }`$. We can estimate the average density of breathers in the residual state from Eqs. (12) and (13) as

$$n_B=_\mathrm{\Delta }^{\mathrm{}}𝒫(ϵ)𝑑ϵ\beta ^{(z+1)}e^{\beta \mathrm{\Delta }}.$$
(14)
This result is consistent with our initial hypothesis (Eq. (10)).

From Eqs. (12) and (13) it follows that the average breather energy $`ϵ`$ can be expressed as a function of $`\beta `$ and $`\mathrm{\Delta }`$. Hence, we can write the normalised distribution Eq. (12) as a function of the three parameters $`z`$, $`ϵ`$ and $`\mathrm{\Delta }`$ in the following fashion

$$𝒫(ϵ)=\frac{(z+1)^{z+1}}{\mathrm{\Gamma }(z+1)(ϵ\mathrm{\Delta })}\left[\frac{ϵ\mathrm{\Delta }}{ϵ\mathrm{\Delta }}\right]^ze^{(z+1)\left[\frac{ϵ\mathrm{\Delta }}{ϵ\mathrm{\Delta }}\right]},$$
(15)
where $`\mathrm{\Gamma }(z)`$ is the gamma function and $`ϵ=\mathrm{\Delta }+(z+1)/\beta `$. From the point of view of the numerics, it is more accurate to deal with the cumulative (integrated) distributions. These are easily obtained from the histograms without actually performing the integration, by noting that they are nothing but rank–size plots with the axes inverted. The cumulative distribution of the function (15) is

$$\begin{array}{cc}\hfill 𝒞(ϵ)=_ϵ^{\mathrm{}}𝒫(ϵ^{})𝑑ϵ^{}=& \frac{1}{\mathrm{\Gamma }(z+1)}\hfill \\ & \times \gamma (z+1,(z+1)\left[\frac{ϵ\mathrm{\Delta }}{ϵ\mathrm{\Delta }}\right])\hfill \end{array}$$
(16)
where $`\gamma (z,ϵ)`$ is the incomplete gamma function, defined as

$$\gamma (z,ϵ)=_ϵ^{\mathrm{}}y^{z1}e^y𝑑y.$$
(17)
In principle one could evaluate $`\mathrm{\Delta }`$ explicitly by calculating the energy of exact breather solutions of vanishing amplitude 2D\_tresh . However, we expect that $`\mathrm{\Delta }ϵ`$. Therefore, if the population of breathers used to calculate the experimental distribution $`𝒞(ϵ)`$ is large enough, we can assume that the smallest breather energy recorded $`ϵ_{min}`$ is close to the threshold $`\mathrm{\Delta }`$. Hence, we can perform a two–parameter fit with the theoretical prediction Eq. (16) with $`\mathrm{\Delta }=ϵ_{min}`$.

The quality of such fits is illustrated in Fig. 6 for a lattice with $`N=70`$. As it shows, the agreement of our simple model with the numerics is excellent. In particular, the analysis of the small–energy portion of the distributions yields $`z=2`$ (see inset in Fig. 6). We stress that the result of our fits were well reproduced irrespective of the particular realisation of $`ϵ_{min}`$. The fit of the experimental distributions over the whole range of available energies yields a slightly greater value, $`z2.6`$. This is because the large–energy tail of the distributions accounts for the rarest events, which are under–represented within reasonably large breather populations <sup>2</sup><sup>2</sup>2The main constraint comes from the time required to numerically integrate the equations of motion of large lattices for a large number of realisations of a given initial condition.. In order to convince oneself that this is the case, it is enough to introduce an energy–dependent weight function in the fit, which gauges the relative weights of the small–energy and large–energy portions of the curves. In this case, the best estimates of $`z`$ monotonically decrease towards the value $`z=2`$ upon lowering the relative weight assigned to the large–energy region.

The analysis of energy distributions also provides an independent confirmation that localisation is in the present case an activated process. To check this, we fitted the distributions obtained by letting systems of given size relax from different values of $`e_0`$. Then, we evaluated $`\beta `$ from the best–fit values of the fitting parameters $`z`$ and $`ϵ`$ according to the relation

$$\frac{1}{\beta }=\frac{ϵϵ_{min}}{z+1}.$$
(18)
As shown in Fig. 7, $`\beta `$ turns out to be inversely proportional to $`e_0`$ as expected.

## IV Conclusions

The mechanisms yielding the spontaneous formation of localised periodic excitations, i.e. breathers, in spatially discrete non–linear dynamical systems are of primary importance for the understanding of their physical interest. In this paper we have focused on the phenomenon of relaxation to breather states by cooling from the boundaries 1D and 2D lattices. We have provided a detailed description of the many facets of this phenomenon by combining numerical studies with theoretical arguments. In particular, we have shown that in models with sufficiently weak “substrate” forces the energy loss process can be explained on the basis of a simple perturbative analysis, that applies independently of the boundary conditions and of the nature of the nonlinearity. In this sense, we can claim that this is a general feature of this class of models, where an initial exponential decay is followed by a power–law behaviour. In the second stage of the decay process, boundary conditions play a crucial role in allowing for the formation of a long–living breather state, when they are taken free. On the contrary, fixed ends yield complete cooling of the lattice eventually ruled by the exponential decay rate of the longest–lived Fourier mode. We have also pointed out that the presence of a sufficiently strong “substrate” force can turn the initial exponential decay to a stretched–exponential law, while maintaining all the other features of the previously described scenario. A theoretical explanation of the origin of this stretched exponential behaviour and its dependence on the strength and on the nature of the local potential is a problem that deserves further and more refined investigations. We have also commented about the technical difficulties that can be encountered for identifying the correct time behaviour in these models, where the crossover between different regimes can be easily mistaken for an indication in favour of other scaling laws.

Another crucial aspect that we have widely analysed concerns the main differences between 1D and 2D systems. At variance with the former case, in 2D we have found numerical evidence of the existence of a finite energy threshold for breather formation. We have also proposed a statistical model that provides an effective agreement with numerical results concerning the breather residual state. This state eventually sets in as an almost static stationary configuration, where breathers of different amplitudes coexist. The statistical model is based on the hypothesis that this state emerges as a result of the competition between breather activation and mutual interaction. Moreover, it allows to determine an effective breather lifetime and its dependence on energy. In this framework the initial energy density plays the role of the control parameter regulating the strength of fluctuations.

In our opinion these results provide a satisfactory understanding of the phenomenon of spontaneous breather formation by cooling. As a final remark, it seems worth investigating experimentally the capability of real systems to store energy in the form of long–lived localised excitations.

###### Acknowledgements.

This work has been supported by the European Union under the RTN project LOCNET, Contract No. HPRN-CT-1999-00163. Francesco Piazza wishes to acknowledge numerous useful discussions with Katja Lindenberg.