text:
# The phase-diagram of the IGM and the entropy floor of groups and clusters: are clusters born warm?

## 1 Introduction

A key test of cosmological scenarios is to check whether they can reproduce the large-scale structure of the present universe. Many studies have already shown that the usual hierarchical scenarios (such as the standard CDM model) agree reasonably well with observations of galaxies (Valageas & Schaeffer 1999, Kauffmann et al. 1993, Cole et al. 1994), Lyman-$`\alpha `$ clouds (Valageas et al. 1999, Petitjean et al. 1992, Miralda-Escude et al. 1996), and quasars, as well as with reionization constraints (Valageas & Silk 1999a, Gnedin & Ostriker 1997, Haiman & Loeb 1998).

However some puzzles still remain unresolved. For instance, the mass of baryons observed at $`z=0`$ (stars, intracluster hot gas) yields $`\mathrm{\Omega }_\mathrm{b}0.01`$ (Fukugita et al. 1998), while the Lyman alpha forest accounts for about as much (Penton et al. 2000). This falls short of the CMB determination by Archeops (Benoit et al. 2003), $`\mathrm{\Omega }_\mathrm{b}0.052\pm 0.007`$, or the prediction of standard nucleosynthesis $`\mathrm{\Omega }_\mathrm{b}0.045\pm 0.006`$ (Tytler et al. 2000) 111Our choice of cosmological parameters is given in Section 2.. As argued in Cen & Ostriker (1999) and Dave et al. (2001), up to $`50\%`$ of the baryons could be in the form of a “warm” phase of the intergalactic medium (IGM), with temperatures in the range $`10^5<T<10^7`$ K, and provide the explanation for the “missing” baryonic matter.

A second problem is the temperature - X-ray luminosity relation of clusters. Indeed, the simplest models yield $`L_\mathrm{X}T^2`$ (Kaiser 1986) which disagrees with observations (Ponman et al. 1996). It has been suggested (e.g., Evrard & Henry 1991, Cavaliere et al. 1997, Ponman et al. 1999) that this discrepancy could be due to a “preheating” of the IGM which would raise its entropy before clusters form. Such an entropy “floor” breaks the standard scaling laws and yields a steeper $`TL_\mathrm{X}`$ relation, and may possibly play a role in the solution of the “overcooling” problem (Blanchard et al. 1992). Two sources of energy which could provide this preheating appear naturally in the usual cosmological scenarios: supernovae/stellar winds and quasars. However, realistic models (e.g., Valageas & Silk 1999b, Wu et al. 2000, Bower et al. 2001) show that the energy released by supernovae is unlikely to be sufficient (unless the energy transfer is highly efficient and targeted to the gas which will build clusters and groups). On the other hand, the energy provided by quasars is largely sufficient (Valageas & Silk 1999b), but it is not easy to estimate with a good accuracy the efficiency of the energy transfer. Although it is commonly agreed that an increase in entropy is required, its actual origin remains an open issue.

Some numerical simulations (Muanwong et al. 2001), on the other hand, show that inhomogeneous cooling of the baryonic gas, which removes the cooler gas as it condenses into galaxies while leaving only the hotter component available as intracluster gas, can play the same role as an increase in entropy. A related proposal (Voit & Bryan 2001) is that the entropy along the edge of the cooling curve matches, at a given virial temperature, the observed entropy of the intracluster gas. However, to be efficient this may imply an excessive amount of cooling (Muanwong et al. 2001).

In this work, we investigate another source of entropy: the shock-heating associated with the gravitational collapse of very large scales which are just turning non-linear. This is partly motivated by the recent finding in numerical simulations of a warm IGM phase generated by such a transfer of energy from collapsing very large scales to much lower scales, see Dave et al. (2001). This mechanism, driven by the non-linear length scale ($`9`$ Mpc at $`z=0`$), heats the IGM which is clustered within much smaller filaments (i.e. with a similar length but a smaller thickness of order of a few hundred kpc at $`z=0`$) and occupies a small fraction of the volume of the universe (e.g., Cen & Ostriker 1999, Valageas et al. 2002). Since clusters are located at the crossing of such filaments, their properties could also be influenced by this shock-heating. Therefore, in this study we investigate whether the entropy floor observed within clusters could be generated by this shock-heating from scales much above the cluster scale. This could occur through two mechanisms. First, the warm IGM embedded within neighbouring filaments, where it underwent this shock-heating, could be efficiently mixed within the cluster as the accretion which builds the halo proceeds. Then, the additional entropy brought with this flow of matter could increase the entropy of the cluster. In order to be relevant, this process requires a sufficiently efficient mixing of the gas so that at least the small clusters with a temperature below $`1`$ keV which display the anomaly are affected in their bulk gas distribution. We shall evaluate the level of mixing required by such a scenario. A second mechanism would be for the gas to be shock-heated in situ, that is the shocks generated by the large scales which are just turning non-linear may run through the cluster and heat the ICM itself. We shall not study this alternative scenario in detail here, but we note that it requires stronger shocks in order to significantly increase the entropy of this higher-density gas.

## 2 Phase-diagram of cosmological baryons

For consistency with our earlier work, we use a Hubble constant $`H_0=65`$ km s-1, and consider a standard $`\mathrm{\Lambda }`$CDM universe with $`\mathrm{\Omega }_\mathrm{m}=0.4`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.6`$, with a baryonic density parameter of $`\mathrm{\Omega }_\mathrm{b}=0.047`$. The $`\mathrm{\Lambda }`$CDM power-spectrum of the linear density fluctuations is normalized by $`\sigma _8=0.8`$.

First, we review in Fig.1 the properties of the IGM at redshift $`z=0`$, as obtained in Valageas et al. (2002). These can be described by the distribution of the baryons within the $`(\rho ,T)`$ parameter space, or equivalently the $`(1+\delta ,T)`$ space, where $`\delta =\rho /\overline{\rho }1`$ is the density contrast. The warm IGM is constrained to lie within the white central region in Fig.1 which is surrounded by exclusion domains. In counter-clockwise order, starting from the right side, these correspond to 1) a fast cooling region where the cooling time $`t_{\mathrm{cool}}`$ is shorter than the Hubble time $`t_\mathrm{H}`$ (so that the gas cannot remain in this region), 2) a high-temperature and high-density domain made of very rare and massive potential wells which are beyond the cut-off of the probability distribution function (pdf) $`𝒫(\rho )`$ (hence they may occur, but with a vanishingly small probability as one moves away from the boundary), 3) linear scales where the gravitational dynamics has not heated the gas yet through shocks, 4) rare voids below the low-density cutoff of the pdf $`𝒫(\rho )`$, and 5) a low-temperature region excluded by radiative heating since the UV background flux heats the gas (here the Lyman-$`\alpha `$ forest) up to $`T10^3`$ K.

The IGM can then be split into two phases. Firstly, a “cool” component described by a well-defined equation of state lies along the lower bound set by radiative heating in the range $`10^11+\delta 6`$. This corresponds to the Lyman-$`\alpha `$ forest where the temperature is set by the interplay of radiative heating by the UV background with the expansion of the universe. Secondly, a warm component is broadly distributed within the central region, around the curve labeled “warm IGM” in Fig.1. This gas is shock-heated by the gravitational dynamics, up to the temperature of the large-scale structures which are just turning non-linear. Finally, some of the gas which enters the cooling region forms stars and galactic disks (note that feedback from supernovae and stellar winds only allow a small fraction of this matter to turn into stars, e.g. Valageas & Schaeffer 1999). The phase-diagram displayed in Fig.1 agrees with the results of numerical simulations (e.g., Dave et al. 2001). We refer the reader to Valageas et al. (2002) for a detailed discussion of these processes and of their redshift evolution.

## 3 Entropy sources

As recalled in the Introduction, in order to break the $`L_\mathrm{X}T^2`$ rule, which holds (Kaiser 1986) when the cluster gravitational potential is the sole source of energy, we may look for additional heating mechanisms. In this paper, we investigate the gravitational shock-heating induced by the larger scales which are just turning non-linear. More precisely, we consider the first scenario described in the Introduction, where this transfer of energy from such large scales down to the ICM proceeds via the warm IGM. Thus, the gas is first shock-heated within surrounding filaments and next falls into the cluster potential well. Therefore, we compute below in Sect.3.1 the entropy $`S_\mathrm{w}`$ of the warm IGM. However, some of the gas within the ICM does not correspond to this late-falling matter. Indeed, part of the ICM observed at $`z=0`$ was already embedded within virialized structures at $`z=1`$ (within loose groups or cooled galaxies) and was later on ejected (through ram pressure, supernovae winds, tidal stripping,…) within the diffuse ICM. Hence we compute in Sect.3.2 the properties of this second ICM component, ejected from substructures within the cluster. Hereafter, following Ponman et al. (1999), we define222The thermodynamical specific entropy is thus, within an additional -irrelevant- constant, equal to $`3/2\mathrm{ln}S`$. the specific entropy of the gas as:

$$S=\frac{kTn_\mathrm{b}^{2/3}}{1\mathrm{keV}\mathrm{cm}^2},$$
(1)
where $`n_\mathrm{b}`$ is the baryon number density.

### 3.1 Entropy associated with the larger scales: the warm IGM contribution

From the phase-diagram shown in Fig.1, we compute the mean333To average the entropy over patches of different entropy $`S_i`$ but also different densities $`n_i`$ it is more relevant to bring them adiabatically to the same effective density, say $`n_{\mathrm{eff}}`$ and then take the mass-weighted average $`T=S_in_{\mathrm{eff}}^{2/3}`$ of their rescaled temperature $`S_in_{\mathrm{eff}}^{2/3}`$ and next to define the average entropy $`S`$ by eq.(1): $`S=T/n_{\mathrm{eff}}^{2/3}`$. This procedure is independent of the choice of $`n_{\mathrm{eff}}`$ and leads simply as the average specific entropy to the mass weighted average (2) of $`S`$ defined in eq.(1). A straightforward way to convince oneself that this is the correct average is to see that this procedure provides the correct final entropy in case the patches are physically mixed at the fixed density $`n_{\mathrm{eff}}`$. This choice is the one which is compatible with our final use of $`T_{\mathrm{ad}}`$ in Sect.4. specific entropy $`S_\mathrm{w}(z)`$ at redshift $`z`$ of the warm IGM, created by the just-collapsing large scales, as:

$$S_\mathrm{w}(z)=\frac{1}{F_\mathrm{w}}dxxh_\mathrm{w}(x)S(x),$$
(2)
with:

$$F_\mathrm{w}(z)=dxxh_\mathrm{w}(x).$$
(3)
The entropy $`S`$ entering the average (2) is the one defined in eq.(1) and the integral runs along the upward rising curve in Fig.1 labeled “Warm IGM”. The variable $`x(1+\delta )/\overline{\xi }`$ (where $`\overline{\xi }(R)=\delta _R^2`$ is the mean two-point correlation over a sphere of radius $`R`$) is a coordinate which, in eq.(2), is taken at a given redshift $`z`$ along the mean equation of state of the warm IGM. The quantity $`xh_\mathrm{w}(x)\mathrm{d}x`$ is the mass fraction enclosed within the interval $`[x,x+\mathrm{d}x]`$ and the mass function $`h_\mathrm{w}(x)`$ is proportional to the pdf $`𝒫(\rho )`$ as described in Valageas et al. (2002), normalized to the total mass within the warm IGM component. Thus, $`F_\mathrm{w}`$ is the fraction of matter enclosed within the warm IGM phase.

We display our result for $`S_\mathrm{w}`$ as the solid line in Fig.2. We see that the entropy $`S_\mathrm{w}`$ of the warm IGM increases with time, since the mean density of the universe declines while the characteristic temperature grows as larger gravitational structures form. This is due to the increasing depth of the potential well of the typical object which turns non-linear (i.e. the virial temperature of the late forming clusters is higher than that attached to the early forming galaxies). This is also apparent in Fig.1 of Valageas et al. (2002) which gives the time evolution of the phase diagram of Fig.1.

### 3.2 Entropy associated with the smaller scales: the substructure contribution

Here we turn to the second component of the ICM: the gas which was already embedded within virialized objects at earlier redshifts. To this order, we first consider the properties of the virialized halos which form at a given redshift $`z`$.

#### 3.2.1 Entropy within objects newly formed at redshift $`z`$

We first define as in eq.(2) the mean entropy $`S_{\mathrm{vir}}(z^{})`$ attached to the virialized halos newly formed at redshift $`z^{}`$:

$$S_{\mathrm{vir}}(z^{})=\frac{1}{F_{\mathrm{vir}}}dxxh_{\mathrm{vir}}(x)S(x),$$
(4)
with:

$$F_{\mathrm{vir}}(z^{})=dxxh_{\mathrm{vir}}(x).$$
(5)
The quantity $`xh_{\mathrm{vir}}(x)\mathrm{d}x`$ is again the mass fraction enclosed within the interval $`[x,x+\mathrm{d}x]`$, but the integration now runs along the line $`\delta =\mathrm{\Delta }_{\mathrm{vir}}`$, above the lower branch of the cooling curve. This is shown by the vertical dotted line in Fig.1, see also Fig.1 in Valageas et al. (2002). The density contrast threshold $`\mathrm{\Delta }_{\mathrm{vir}}200`$ is obtained from the usual spherical model. In particular, $`F_{\mathrm{vir}}(z^{})`$ is the fraction of matter embedded within these virialized halos at a redshift $`z^{}`$. In eq.(5) the lower bound of the integration is set by the lower branch of the cooling curve. This quantity was already computed in Valageas et al. (2002) (solid line in Fig.3 there) and checked against numerical simulations. To obtain the entropy $`S(x)`$ of a given halo we use the definition (1) at the adequate density $`n_b`$, setting $`T=T_{\mathrm{vir}}`$, where $`T_{\mathrm{vir}}`$ is the virial temperature:

$$T_{\mathrm{vir}}=\frac{𝒢\mu m_pM}{2kR},$$
(6)
$`M`$ is the mass and $`R`$ the radius of the halo. Therefore, the mean $`S_{\mathrm{vir}}`$ describes the mean entropy of just virialized halos formed at a given redshift $`z^{}`$, generated during their own gravitational collapse (whence the use of $`T_{\mathrm{vir}}`$).

Although we integrate eq.(5) along the density threshold $`\mathrm{\Delta }_{\mathrm{vir}}`$ at redshift $`z^{}`$, it allows us to describe the high-density substructures present within virialized clusters at later redshifts $`z<z^{}`$. Indeed, as noticed in Navarro et al. (1996) the core density within dark matter halos roughly scales as the density of the universe at the redshift when this scale turned non-linear. This scaling is also consistent with the stable-clustering Ansatz (Peebles 1980). Therefore, the matter described by eqs.(4)-(5) at redshift $`z^{}`$ actually corresponds to substructures with a density contrast $`\mathrm{\Delta }_{\mathrm{vir}}(z^{})(1+z^{})^3`$ at redshift $`z=0`$ within larger just-virialized cluster, see Sect. 3.2.2 below.

We estimate the entropy available through ejection of mass from these virialized objects by this specific entropy $`S_{\mathrm{vir}}(z^{})`$. Indeed, the gas which can escape from smaller-scale halos (e.g., galaxies) must have a kinetic or thermal energy above the potential well of these objects (e.g., the velocity of outflows must be larger than the local escape velocity). We assume that the additional energy poured into the ICM scales as the potential energy of this subunit. Here, the advantage of using the entropy as a basic variable is that shock-heating which may occur after virialization (when small-scale objects feed the intermediate scale structures) can only increase the entropy of the system. Moreover, it takes into account the expansion of the gas as it spreads within the cluster.

#### 3.2.2 Entropy of a collection of condensed objects formed prior a given epoch $`z`$

The substructures located within a new larger-scale cluster were formed in the past so that there is still an average to be made over their respective epoch of formation. This mean is expected to be dominated by contributions $`S_{\mathrm{vir}}(z^{})`$, at definitely earlier times $`z^{}>z`$, and $`S_{\mathrm{vir}}(z)`$ enters this sum only marginally.

We now proceed evaluating this entropy, $`S_{\mathrm{sub}}`$, within substructures existing at redshift $`z`$, averaged over their previous redshift of formation and weighted by the mass flux entering into virialized objects. More precisely, we write for the average entropy $`S_{\mathrm{sub}}`$ of these structures:

$$S_{\mathrm{sub}}(z)=\frac{1}{F_{\mathrm{vir}}(z)}_z^{\mathrm{}}dz^{}\frac{\mathrm{d}F_{\mathrm{vir}}}{\mathrm{d}z}(z^{})S_{\mathrm{vir}}(z^{}).$$
(7)
The meaning of $`S_{\mathrm{sub}}`$ is the following. A new cluster builds up at redshift $`z`$ as a large scale density fluctuation turns non-linear and older smaller-scale structures gather within this deeper potential well. Then, these substructures, already virialized within the just-forming cluster of interest, may also be the source of entropy for this new cluster. Indeed, the gas embedded within such smaller objects (e.g., galaxies) may escape from these subunits and become part of the intra-cluster medium (ICM), for any reason (supernovae feedback, tidal stripping, mergers), so it can bring some entropy to the more diffuse ICM.

We display $`S_{\mathrm{sub}}(z)`$ as the dotted line in Fig.2. As expected, $`S_{\mathrm{sub}}(z)`$ declines at higher redshift, as $`S_\mathrm{w}`$, because at earlier times typical potential wells were shallower. The average over redshift which enters the definition (7) also implies that $`S_{\mathrm{sub}}(z)`$ is significantly smaller than $`S_{\mathrm{vir}}`$. Thus, we obtain numerically $`S_{\mathrm{sub}}(z=0)S_{\mathrm{vir}}(z=1)`$ which means that an important fraction of the matter which enters virialized halos at $`z=0`$ was already embedded within collapsed structures at $`z1`$ (note that this is consistent with numerical simulations, e.g. Borgani et al. 2002). So, the major contribution to $`S_{\mathrm{sub}}`$ at the present epoch may be attributed to objects which form at $`z1`$, that somewhat loosely one may call galaxies. But it should be clear that in our average specific entropy all substructures are accounted for.

### 3.3 Relative efficiency of external energy sources

The entropy $`S_{\mathrm{sub}}`$ may provide an additional contribution to the cluster entropy floor, independently of $`S_\mathrm{w}`$. However, this entropy is initially “trapped” within the substructures, and we now discuss to what extent it is available to the intracluster gas. Some feedback from the gas trapped within galaxies, for instance, is expected, since the cluster gas is known to contain metals. But it is also expected that galaxies retain a non-negligible fraction of their gas. As an example to start with, we shall adopt in Sect.4 the simple rule that this is the source of about $`1/3`$ of the cluster gas, an estimate which is consistent with the typical metallicity, $`1/3`$ solar, of a cluster. Thus, we may estimate the external entropy level to be

$$S_{\mathrm{ext}}=\left(\frac{1}{3}S_{\mathrm{sub}}+\frac{2}{3}S_\mathrm{w}\right).$$
(8)
This is the essence of hierarchical clustering, with a contribution from the larger just-collapsing structures as well as from the smaller, already collapsed objects, the dark matter as well as the attached baryons being incorporated into the newly forming structures.

Note that both processes, leading to either $`S_\mathrm{w}`$ or $`S_{\mathrm{sub}}`$, break the simple scaling laws which led to $`L_\mathrm{X}T^2`$ and, if sufficiently strong, might explain the observed temperature-luminosity relation of clusters (as with other physical processes based on supernovae or quasar-driven energy injection which also yield an entropy floor). First, the formation of large-scale structures through gravitational instability builds a network of filaments and shock-heated regions which provide a background entropy and temperature for the gas which eventually gathers to build a cluster. Second, a collapsing group has inner substructures which may be relevant to its subsequent evolution.

These entropy floors however have a negligible impact on the formation of very rare and massive clusters, characterized by high temperatures. This may also be tied to the fact that rare events are known to be almost spherically symmetric (e.g., Bernardeau 1994, Valageas 2002), while the ones under the influence of external heating by the larger scale filaments are precisely those which are expected not to bear the spherical symmetry. Most important, let us note that our new heating process which we have assumed to contribute to structures ranging from rich clusters (where it turns out to be subdominant) down to small groups (where it modifies significantly the gas distribution of the whole cluster) needs only to be present in the latter for our explanation for the breaking of the $`L_\mathrm{X}T^2`$ rule to hold.

## 4 The temperature-luminosity relation of clusters

We compute the temperature - bolometric luminosity relation of clusters as in Valageas & Schaeffer (2000) (see Sect.5). Hence, we assume approximate hydrostatic equilibrium for the gas within a dark matter halo described by the density profile:

$$\rho (r)=\frac{\rho _c}{(r/r_s)^{1.5}\left[1+(r/r_s)^{1.5}\right]}\text{with}r_s=\frac{R}{c},$$
(9)
where $`c`$ is the concentration parameter and $`R`$ the virial radius. The latter is defined as usual as the radius where the overall density contrast reaches the threshold $`\mathrm{\Delta }_{\mathrm{vir}}200`$ given by the standard spherical model. We use $`c=4`$ as in Moore et al. (1999), although this parameter may depend slightly on the mass of the cluster. Besides, at each radius $`r`$ we can associate an effective temperature $`T_\mathrm{s}(r)`$ which measures the depth of the potential well, defined by:

$$T_\mathrm{s}(r)=\frac{𝒢\mu m_pM(<r)}{2kr}=T\frac{\mathrm{ln}\left(1+(r/r_s)^{1.5}\right)}{\mathrm{ln}\left(1+c^{1.5}\right)}\frac{R}{r}.$$
(10)
Thus, $`T_\mathrm{s}(r)`$ yields the temperature which can be reached by the gas through non-adiabatic shock-heating during the formation of the cluster. On the other hand, since the initial entropy (or temperature) of the gas is non-zero we can also define the temperature $`T_{\mathrm{ad}}(r)`$ associated with an adiabatic compression during the infall, from the typical density of the warm IGM (i.e. $`\delta 10`$) up to the density $`n_\mathrm{g}(r)`$ reached within the cluster (i.e. $`\delta >\mathrm{\Delta }_{\mathrm{vir}}`$). It is set by the entropy floor $`S_{\mathrm{ext}}=S_\mathrm{w}`$ through the relation (see eq.(1)):

$$T_{\mathrm{ad}}(r)=\frac{n_\mathrm{g}(r)^{2/3}}{1\mathrm{cm}^2}S_{\mathrm{ext}}\mathrm{keV},$$
(11)
since it corresponds to a constant-entropy (i.e. adiabatic) dynamics. Thus, $`T_{\mathrm{ad}}`$ is the temperature reached by some adiabatic gas falling from the outer warm IGM into the cluster, with an initial entropy $`S_{\mathrm{ext}}`$ acquired in a larger-scale filament. In order to take into account the fact that distinct regions of the warm IGM can have different entropies, within the range allowed by the phase-diagram shown in Fig.1, we consider that the gas thermalizes within the cluster at a typical density $`\delta \mathrm{\Delta }_{\mathrm{vir}}`$. From eq.(11), we see that this amounts to using as a unique initial entropy level the mean $`S_\mathrm{w}`$ defined in eq.(2). Moreover, through eq.(8) this also allows us to take into account the gas which did not come lately from the warm IGM but was ejected from substructures. Next, in order to take into account both these possible external sources and the shock-heating associated with the cluster itself, assuming all sources of entropy correspond to matter which is well-mixed throughout the whole cluster 444We indeed expect the small, somewhat irregular, groups where this contribution is most important to have undergone substantial mixing., we write:

$$T_\mathrm{g}(r)=T_\mathrm{s}(r)+T_{\mathrm{ad}}(r).$$
(12)
Thus, if $`T_\mathrm{s}>T_{\mathrm{ad}}`$, the physics of the gas is governed by the local potential well and $`T_\mathrm{g}T_\mathrm{s}`$. This corresponds to massive high-temperature clusters or to the outer shells. On the other hand, if $`T_\mathrm{s}<T_{\mathrm{ad}}`$ the behaviour of the gas is governed by the external entropy source. In fact, the density of the gas $`\rho _\mathrm{g}`$ then shows a plateau within a core radius $`R_{\mathrm{core}}`$ where both terms in eq.(12) are equal. Indeed, when $`T_\mathrm{g}T_\mathrm{s}`$ the gas no longer falls into the potential well and a core with a nearly constant baryon density appears. On the other hand, at larger radii the gas follows the dark matter and we take $`\rho _\mathrm{g}=\mathrm{\Omega }_\mathrm{m}/\mathrm{\Omega }_\mathrm{b}\rho `$. Thus, for smaller halos the preheating becomes more important and the gas density profile gets flatter, decreasing the X-ray luminosity. We refer the reader to Valageas & Schaeffer (2000) for a more detailed discussion of this model (and the role of radiative cooling).

Finally, we display in Fig.3 the $`L_\mathrm{X}T`$ relation of clusters obtained in this way at $`z=0`$. The solid line is $`S_{\mathrm{ext}}=S_\mathrm{w}`$, the dotted line $`S_{\mathrm{ext}}=S_{\mathrm{sub}}`$ while the lower dashed line is the mean defined in eq.(8). We recover the knee at low luminosities that is found in observations. This feature is common to all models which include some preheating of the gas which provides an entropy floor of order $`100`$ keV cm2. However, we see that the case $`S_{\mathrm{ext}}=S_{\mathrm{sub}}`$ provides an entropy level which is too small. Therefore, gas escaping from older substructures cannot explain the extra entropy of small groups. In particular, if the only source of entropy were galactic outflows (driven by supernovae or any other process), this would not suffice to account for the observed entropy floor. This is consistent with detailed studies of preheating by supernovae and stellar winds (e.g., Valageas & Silk 1999b, Wu et al. 2000). The warm IGM, on the other hand, provides an important contribution. In particular, we see that even in the case defined by eq.(8), where the energy brought by the warm IGM is diluted by a factor $`2/3`$ (i.e. only $`2/3`$ of the current X-ray emitting gas came from the outer IGM over a redshift range $`\mathrm{\Delta }z1`$), the entropy floor is sufficient to explain the knee of the $`L_\mathrm{X}T`$ relation. An interesting exercise is to consider the lowest possible entropy injection required by the data (left dashed-line in Fig.3). This may be achieved even with no entropy contribution from galactic outflows provided $`38\%`$ of the cluster gas originates from the recently formed warm IGM (the other $`62\%`$ having a negligible entropy). If we include our (modest, but probably quite realistic) estimate of the entropy associated with galactic outflows the requirement is that at least $`28\%`$ of the ICM bears the entropy of the warm IGM phase (the other $`72\%`$ bearing $`S_{\mathrm{sub}}`$). These modest requirements are quite plausible, indeed, and strongly back our suggestion.

This entropy floor agrees with the plateau associated with small groups (e.g., Lloyd-Davies et al. 2000). Thus, in this case the properties of the warm IGM would explain in a natural fashion the observed properties of groups and clusters. From the $`L_\mathrm{X}T`$ relation we can also compute the luminosity function of clusters and its evolution with redshift. Such a study is described in Valageas & Schaeffer (2000).

We should emphasize that the underlying picture of our model is that the entropy floor is connected with the warm IGM, and to its equation of state which involves regions of density contrast $`\delta 10`$ and scales much larger than the cluster scale. This is quite different from the alternative suggestion of Voit & Bryan (2001). In the latter model, for instance, it is seen from their Fig.1 that their entropy threshold is associated with the internal cooling of the cluster, in very high density regions of $`\delta 10^4`$, namely those actually observed in the centres of clusters which have completed their formation process. Their suggestion actually relies on a very specific feature of the cooling curve, which, drawn in the $`S(T)`$ variables, happens to exhibit a constant-entropy plateau in a very specific -and extremely narrow- range of temperatures. Of course, in our scenario the cooling curve $`t_{\mathrm{cool}}=t_\mathrm{H}`$ also comes into play, but as a whole as it defines the high-density boundary of the “warm” IGM phase in Fig.1. However, the entropy of the gas is not at all defined by the edge of this exclusion region (and hence does not depend on its local behaviour in a small region of phase-space) since the warm IGM component is actually spread over a broad band within the white central region drawn in Fig.1, as shown by numerical simulations (e.g., Dave et al. 2001). Besides, the gas is truly preheated before the cluster forms while it is embedded within the complex network of filaments at $`\delta 10`$.

Here we may also note that Muanwong et al. (2001) find that their numerical simulations which involve radiative cooling indeed recover a strong increase of the entropy of the X-ray gas within clusters which allows them to match the observed $`L_\mathrm{X}T`$ relation. However, they remark that the mass of cold gas they obtain is too large. Here we did not explicitly model this process but it simply means that some of the gas which has fallen into the cooling region has been reheated by some process (supernovae being one of the possibilities). This gas may then either remain associated with galactic halos (in which case it is irrelevant for our purpose here) or be ejected into the ICM. This corresponds to the contribution $`S_{\mathrm{sub}}`$ discussed in Sect.3.2.

## 5 Shock-heating of the ICM itself from larger scales

Finally, we briefly discuss here the second scenario we presented in the Introduction. This corresponds to the case where the gravitational shock-heating associated with the large scales which are just turning non-linear directly heats the ICM in situ, without involving the transfer of matter from the warm IGM. The characteristic temperature involved by this scenario is again the temperature (or shock velocity) associated with the non-linear length scale, which is of order $`T1`$ keV at $`z=0`$ (e.g., Cen & Ostriker 1999, Valageas et al. 2002). Therefore, such shock-heating may indeed influence the properties of groups with $`T<1`$ keV. In fact, this is the same “coincidence” which underlies the scenario we described in the previous sections, which relies on the warm IGM as an intermediate tool to realize this energy transfer from the non-linear length scale down to the ICM scale. However, in order to increase the entropy of the ICM by the same level as the first scenario this model requires a stronger shock-heating because of the higher density of the gas within the cluster as compared within filaments, see eq.(1). Indeed, we can see from Fig.1 that by heating the gas within filaments we do not need to systematically reach high temperatures $`T1`$ keV in order to achieve high entropies, because of the smaller density of the warm IGM (the density contrast $`\delta `$ typically runs from 6 to 200). On the other hand, although large-scale shocks (which may extend over $`15h^1`$ Mpc) are indeed seen in numerical simulations to envelop and penetrate deep inside clusters (e.g., Miniati et al. 2000) it is difficult to estimate their efficiency. We shall not investigate this scenario further in this work but we expect it can provide an alternative way of injecting energy into clusters from larger scales.

## 6 Conclusion

We have shown that the entropy of the warm IGM can be sufficient to explain the observed entropy of small clusters. This scenario relies on the assumption that the gas which has recently falled into the cluster from the surrounding filaments has been efficiently mixed within the halo so as to form at least $`28\%`$ of the ICM at the radii where the extra-entropy is observed. This point is clearly a matter of debate since it is not obvious to follow the behaviour of the gas within such complex objects as clusters. However, we feel this requirement is quite plausible. A natural scale appears in our scenario, $`1`$ keV, the potential energy of the scales just turning non-linear. Clusters below $`1`$ keV where this external entropy source can play a role are specifically the ones under debate. These rather small clusters should be rather irregular so that their accretion is not expected to be well-ordered or spherically symmetric. Then, the late gas coming from the warm IGM may indeed spread over the whole group, taking advantage of the general mixing (e.g., turbulence) brought by the chaotic dynamics of the various components (e.g., galaxies) of the group. On the other hand, we briefly notice that one could also imagine a second scenario. There, the gravitational energy associated with the large scales which are just turning non-linear could directly heat the ICM itself. This assumes that shocks generated by those large scales can spread within the cluster and heat at least $`28\%`$ of the ICM. However, in order to reach the same entropies this scenario requires a stronger shock-heating than in the first model because of the need to compensate for the higher density of the gas within the ICM as compared with filaments.

Altough the efficiency of these mechanisms remains to be evaluated with a better accuracy before one can safely determine whether they indeed play an important role in the physics of clusters, we note that such scenarios present the advantage to dispense with the need of a new external preheating source (e.g., supernovae or quasars). Indeed, they only rely on gravitational shock-heating which has recently been seen to explain the properties of the IGM in numerical simulations (e.g., Dave et al. 2001). Therefore, we merely propose here that the mechanisms which build the warm IGM may also govern the properties of small clusters. So far, numerical simulations have not really substantiated this claim since they usually recover the standard scaling $`L_XT^2`$, unless one adds some ad-hoc preheating or takes into account cooling processes (which lead to other problems like overcooling). However, they may not have been able to pin down all salient effects of shock-heating mechanisms yet.

Similar processes may be expected to be at work at the epoch of galaxy formation. They may for instance modify the hydrodynamical equilibrium of the smaller galaxies. In particular, as already noticed in Valageas & Silk (1999b) preheating models are strongly constrained by the contradictory requirements arising from the observed properties of clusters and galaxies. This will be studied in a forthcoming paper for the process presented here which exhibits a specific evolution with redshift.

In summary, there are three hypothesis for the entropy floor in a given cluster. It may arise 1) from external preheating by supernovae or quasars; 2) from a very specific feature of the cooling curve, as implied by the suggestion of Voit & Bryan (2001); 3) as we suggest here, from the pre-collapse phase of structures much larger than clusters which transfer their energy to the smaller cluster scales through shock-heating, which also builds the warm IGM phase.

We have stressed the difficulties encountered by the two first possibilities. The advantage of the latter explanation is that it is a natural outcome of hierarchical scenarios of large-scale structure formation. Numerical simulations and observations may be used to check the entropy evolution implied by our suggestion.

## acknowledgements

We are very grateful to the referee, Trevor Ponman, for his stimulating comments, helping us to improve this paper.


format: text

text:
# The ultimate state of thermal convection

## Abstract

The ultimate regime of thermal convection, the so called Kraichnan regime (R. H. Kraichnan, Phys. Fluids 5, 1374 (1962)), hitherto has been elusive. Here, numerical evidence for that regime is presented by performing simulations of the bulk of turbulence only, eliminating the thermal and kinetic boundary layers and replacing them by periodic boundary conditions.

Thermally driven turbulence is of tremendous importance in oceanography, geophysics, meteorology, astrophysics, or process technology. If the flow is confined to a box heated from below and cooled from above, the setup is called Rayleigh-Bénard convection. It is one of the classical problems in fluid dynamics . The control parameter are the dimensionless temperature difference between bottom and top (called Rayleigh number $`Ra`$), the ratio between kinematic viscosity $`\nu `$ and heat diffusivity $`\kappa `$ (called Prandtl number $`Pr`$), and the aspect ratio of the cell. The system answers with some heat transfer from the bottom to the top (called Nusselt number $`Nu`$ in dimensionless form) and the degree of turbulence in the cell (expressed in terms of the Reynolds number $`Re`$).

All classical experiments on this system show (effective) power law dependences with $`NuRa^{0.25}`$ up to $`NuRa^{0.33}`$. However, for very large $`Ra`$, Kraichnan  has predicted an ultimate scaling law $`NuRa^{1/2}`$. In this ultimate state the heat flux and the turbulent velocity are expected to be independent of the viscosity and heat diffusivity. Indeed, for a system with infinite aspect ratio a mathematically strict upper bound $`Nu0.167Ra^{1/2}1`$ could be found .

Hitherto, this ultimate state of thermal convection remained elusive , inspite of tremendous efforts to find it. Niemela et al.  measured up to $`Ra10^{18}`$ in helium gas close to the critical point ($`Pr1`$), still finding an effective power law $`NuRa^{0.31}`$. Glazier et al.  measured up to $`Ra10^{11}`$ in mercury ($`Pr=0.025`$), finding $`NuRa^{0.28}`$. The only evidence for a transition (at $`Ra10^{11}`$) towards the ultimate regime has been claimed by Chavanne et al. , but Niemela and Sreenivasan  argue that those data would be consistent with $`NuRa^{1/3}`$. In summary, it is therefore not clear

* whether in today’s experiments the Rayleigh number is only still too small for the ultimate state, or
* whether it does not exist at all.

However, it is crucial to know which of these two alternative is the correct one:

* for practical reasons in view of above mentioned applications of thermal convection in geophysics, industry etc, as upscaling from lab-scale models does not work if there is a transition to a new state of turbulence;
* for fundamental reasons, as is a “central dogma in turbulence”  that its effects become independent of viscosity and diffusivity for large enough Reynolds (or here Rayleigh) numbers .

Indeed, also within Grossmann and Lohse’s unifying theory on thermal convection  the Kraichnan or ultimate regime follows, once the thermal and kinetic boundary layers either break down due to their expected instability at large $`Ra`$, or hardly contribute to the global kinetic and thermal dissipation . According to ref.  the breakdown of the laminar kinetic boundary layer should happen at $`Ra10^{14}`$ for $`Pr=1`$ and at $`Ra10^{12}`$ for $`Pr=0.025`$ for an aspect ratio $`\mathrm{\Gamma }=1`$ cell; Niemela and Sreenivasan  find similar values. Then the total kinetic dissipation rate $`ϵ_u`$ equals the bulk kinetic dissipation $`ϵ_{u,bulk}`$ and the total thermal dissipation rate $`ϵ_\theta `$ equals the bulk thermal dissipation $`ϵ_{\theta ,bulk}`$ and the Kraichnan or ultimate regime with

$$NuRa^{1/2}Pr^{1/2}$$
(1)
and

$$ReRa^{1/2}Pr^{1/2}$$
(2)
results .

From the above it follows that an artificial destruction of the boundary layers should enhance the ultimate regime. Indeed, by performing Rayleigh-Bénard convection experiments in a cell with rough top and bottom walls Roche et al.  find evidence for the onset of the Kraichnan regime.

Here, we follow another approach: We perform numerical simulation of the bulk of turbulence only, eliminating the thermal and kinetic boundary layers and replacing them by periodic boundary conditions . Then the Kraichnan regime should follow immediately. The thermal bulk turbulence is forced by a mean temperature profile $`z\mathrm{\Delta }/L`$, where $`L`$ is the periodicity length and $`\mathrm{\Delta }`$ the temperature difference on that length scale. The temperature fluctuations $`\theta `$ then obey the advection equation

$$_t\theta +u_i_i\theta =\kappa _i_i\theta +u_3\frac{\mathrm{\Delta }}{L}.$$
(3)
The velocity field $`u_i(𝒙,t)`$ obeys the standard Boussinesq equation,

$$_tu_i+u_j_ju_i=_ip+\nu _j^2u_i+\beta g\delta _{i3}\theta .$$
(4)
$`\beta `$ is the thermal expansion coefficient, $`g`$ gravity, and $`p`$ the pressure. $`Ra`$ and $`Pr`$ are defined as usual as $`Ra=\beta g\mathrm{\Delta }L^3/(\nu \kappa )`$ and $`Pr=\nu /\kappa `$, respectively. The Prandtl number is 1 in our simulations presented here. Also $`Nu`$ is defined as usual as $`Nu=u_3\theta L/(\kappa \mathrm{\Delta })_3\theta _AL/\mathrm{\Delta }`$. The Reynolds number is defined with the rms mean velocity fluctuation $`u^{}`$, i.e., $`Re=u^{}L/\nu `$. A similar simulation has earlier been performed by Borue and Orszag . The focus of that work was on scaling properties and spectra and thus hyperviscosity has been employed.

From eqs. (3) and (4) one can derive by volume averaging and assuming statistically stationarity the exact  relations $`ϵ_\theta =\kappa L^2\mathrm{\Delta }^2Nu`$ and $`ϵ_u=\nu ^3L^4NuRaPr^2`$. As in refs.  we estimate $`ϵ_{u,bulk}u^3/L`$ and $`ϵ_{\theta ,bulk}u^{}\mathrm{\Delta }^2/L`$ and the Kraichnan regime (1) – (2) results. Note that $`L`$ is the relevant length scale for these estimates as for the Bolgiano length scale it holds $`l_B\text{ }\stackrel{>}{}\text{ }L`$.

This Kraichnan scaling is indeed seen in our numerical simulations of eqs. (3) and (4). We run four different $`Ra`$ numbers $`Ra=8.6410^5`$ (for 160 large eddy turnovers), $`Ra=2.2110^6`$ (for 60), $`Ra=4.5110^6`$ (for 180), and $`Ra=1.3810^7`$ (for 334), and average over space and time to obtain $`Nu`$ and $`Re`$ . The $`Ra`$-scaling is consistent with eqs. (1) – (2), see figure 1. To our knowledge it is the first realization of the ultimate regime of thermal convection. It is remarkable that it is realized in spite of the relatively low $`Ra`$ numbers of the simulations. The reason is that the boundary layers have been eliminated and the simulations focus on the bulk. Due to the lack of boundary layers and the more efficient driving, the Reynolds numbers achieved here are much larger than they would be in the standard Rayleigh-Bénard case. Vice versa, the required Rayleigh numbers to achieve the same degree of turbulence for the standard Rayleigh-Bénard case as here are much higher. E.g., to achieve $`Re=4996`$ (here for $`Ra=1.3810^7`$) one needs a Rayleigh number of $`Ra=7.210^8`$ in the standard case .

In conclusion, from above alternative the first one is favorable: Once $`Ra`$ is so large that the laminar kinetic boundary layers breaks down, the ultimate regime of thermal convection, which has been so elusive in experiment, should finally show up. Obviously, the ultimate proof for the existence of this regime can only come from experiment.

Acknowledgement: The work is part of the research program of the Stichting voor Fundamenteel Onderzoek der Materie (FOM), which is financially supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO). We acknowledge support from the European Union (EU) through the European Research Network on “Nonideal Turbulence” (contract HPRN–CT–200000162). We also acknowledge INFN, Sezione di Pisa for the APEmille computer resources.


text:
# Absorption losses in periodic arrays of thin metallic wires

## Abstract

We analyze the transmission and reflection of the electromagnetic wave calculated from transfer matrix simulations of periodic arrangements of thin metallic wires. The effective permittivity and the absorption is determined. Their dependence on the wire thickness and the conductance of the metallic wires is studied. The cutoff frequency or effective plasma frequency is obtained and is compared with analytical predictions. It is shown that the periodic arrangement of wires exhibits a frequency region in which the real part of the permittivity is negative while its imaginary part is very small. This behavior is seen for wires with thickness as small as 17 $`\mu `$m with a lattice constant of 3.33 mm.

Rapidly increasing interest in left-handed metamaterials (LHM) (for a recent review, see MSnato ) raised some interesting questions about the electromagnetic (EM) properties of composites which contain very thin metallic components. The most simple example of such a composite is a periodic array of thin metallic wires. Pendry el al. Pendry  predicted that such system behaves as a high pass filter with an effective permittivity

$$ϵ_{\mathrm{eff}}=1\frac{f_p^2}{f^2+2i\gamma f}.$$
(1)
In Eq. 1, $`f_p`$ is the effective plasma frequency, or cutoff frequency Sigalas , and $`\gamma `$ is the damping factor. Various theoretical formulas were derived for the dependence of the plasma frequency on the lattice period $`a`$ and the wire radius $`r`$. Pendry el al. Pendry  obtained that

$$f_p^2=\frac{c_{\mathrm{light}}^2}{2\pi a^2\mathrm{ln}(a/r)},$$
(2)
Shalaev and Sarychev Shalaev  obtained that

$$f_p^2=\frac{c_{\mathrm{light}}^2}{2\pi a^2[\mathrm{ln}(a/\sqrt{2}r)+\pi /23]},$$
(3)
while Maslovski et al. Maslovski  that

$$f_p^2=\frac{c_{\mathrm{light}}^2}{2\pi a^2[\mathrm{ln}a^2/4r(ar)]}.$$
(4)
In Eqs. 2-4, $`c_{\mathrm{light}}`$ is the velocity of light in vacuum.

Periodic arrangements of thin metallic wires are used as a negative-$`ϵ`$ medium Brown ; Rotman  in the left-handed structures Smith ; Shelby ; Claudio . It is therefore important to understand how the electromagnetic response - not only the effective plasma frequency, but also the factor $`\gamma `$ \- depends on the structural parameters of the wire system. Recently, Ponizhovskaya el al. PGN  claimed that for small wire radius, the absorption in the wire system is so large that the transmission losses do not allow any propagation of EM wave in left-handed structure. Very low transmission, measured in the original experiments on the LHM, Smith ; Shelby  seemed to agree with their pessimistic conclusion. However, it is not clear why the transmission was so low in the original experiments. Recent experimental measurements Claudio ; Ekmel  established that the transmission of LHM could be as good as in the right-handed systems.

Our aim in this Letter is to study numerically how the effective permittivity of the periodic arrangement of metallic wires depends on the wire radius and on the conductance of the wires. We present results for the real ($`ϵ_{\mathrm{eff}}^{}`$) and imaginary ($`ϵ_{\mathrm{eff}}^{\prime \prime }`$) part of the effective permittivity of the wire medium, estimate the transmission losses and the plasma frequency and compare our results with the analytical formulas given in Eqns. 1-4.

In our numerical simulations we use the transfer matrix method (TMM). Details of the method are given elsewhere MS-1 . Here we only point out the main advantage of the TMM, namely that it gives directly the transmission $`t`$, reflection $`r`$ and absorption $`A=1|t|^2|r|^2`$ of the EM plane wave passing through the system. Contrary to this, in the FDTD method, which was used in Ref. PGN , one obtains $`t`$ and $`r`$ form the time development of the wave packet, which is much more complicated and probably also less accurate. To be able to obtain $`ϵ_{\mathrm{eff}}`$, one also needs the phase of $`r`$ and $`t`$, in addition to their amplitudes. This is also easily achieved in the TMM.

From the obtained data for the transmission and reflection, we calculate the effective permittivity of the system. MS-2  The refraction index is given by the equation

$$\mathrm{cos}(nkL)=\frac{1}{2t}\left[1r^2+t^2\right].$$
(5)
Since we do not expect any magnetic response, we fixed the value of the permeability to be $`\mu 1`$. The permittivity is then found as $`ϵ_{\mathrm{eff}}=n^2`$.

The way we discrete the space in the TMM might control the accuracy of our results. To test how discretization influences our results, we repeated the numerical simulation for different discretizations. The wire is represented as a rectangular with square cross section $`2r\times 2r`$, $`r`$ being the “wire radius”. The obtained results for the effective permittivity $`ϵ_{\mathrm{eff}}`$ (both real and imaginary part) are almost independent of the discretization procedure.

Figure 1 shows how the effective permittivity depends on the wire radius. We analyzed four different wire arrays with period $`a=3.33`$ mm. For all of them the real part of the effective permittivity is negative and could be fitted by Eq. (1). This enables us to obtain easily the plasma frequency. Only for the smallest wire thickness studied ($`17\times 17\mu `$m) one gets a relatively large $`ϵ_{\mathrm{eff}}^{\prime \prime }`$ for small frequencies. In this region relation (1) is not valid. Nevertheless, for frequencies larger than 5 GHz, $`ϵ_{\mathrm{eff}}^{\prime \prime }`$ is small and $`ϵ_{\mathrm{eff}}^{}`$ negative.

The right lower panel of fig. 1 shows data for wires with the cross section of $`17\times 300\mu `$m. These parameters were used in the experiment of Shelby et al. Shelby . We again see that the formula given by Eq. 1 qualitatively agree with our data. Thus, there is no doubt that this array of wires really produces a medium with negative $`ϵ_{\mathrm{eff}}^{}`$, which then can be used in the creation of the left handed systems.

Figure 2 compares our data for plasma frequency with the analytical formulas given by Eqns. 2-4. Accepting some uncertainty in the estimation of the plasma frequency from the numerical data, we can conclude that for thin wires our data agree with the theoretical formula (3) of Shalaev and Sarychev. Shalaev  For thicker wires, our results are in agreement with the formula (4) of Maslovski et al. Maslovski

We also study the role of the conductance of the metallic wires. In the simulations shown in fig. 1, we consider the metallic permittivity to be $`ϵ_\mathrm{m}=(3+588i)\times 10^3`$. We are aware that this value of $`ϵ_\mathrm{m}`$ is smaller than the permittivity of realistic metallic wires: for instance for copper $`ϵ_\mathrm{m}5\times 10^7i`$, as follows from the relation between the permittivity and the conductance Jackson  (the conductivity of copper is $`\sigma 5.9\times 10^7(\mathrm{\Omega }\mathrm{m})^1`$). Our data in fig. 1 therefore underestimate losses, because transmission losses are smaller for higher values of the metallic permittivity Pendry-2 . This is clearly shown in fig. 3 where we present the ratio $`\kappa =|ϵ_{\mathrm{eff}}^{\prime \prime }/ϵ_{\mathrm{eff}}^{}|`$ vs frequency for two systems which differ only in the value of the imaginary part of the metallic permittivity. Fig. 3 shows also the frequency dependence of the absorption as obtained from the numerical simulations. The absorption also exhibits a maximum in the neighbor of the plasma frequency.

Fig. 1 also confirms that $`ϵ_{\mathrm{eff}}^{\prime \prime }`$ increases when the wire radius decreases. For instance, the parameter $`\gamma `$ is only 0.003 GHz for $`r=100\mu `$m. As we will see below, (fig. 4) $`\gamma `$ increases up to 1.2 GHz when the wire radius decreases to 15$`\mu `$m. Nevertheless, even for very thin wires, losses are much less than what was claimed in Ref. PGN . As it is shown in fig 1, also an array of wires with thickness $`17\times 17\mu `$m creates a negative-$`ϵ`$ medium. As this is in strong contrast with the results of Ref. PGN , we decided to study exactly the same system as that of PGN . Results of our simulations are shown in fig. 4. Although such systems are not used in experimental arrangements of the LHM, our results give a comparison between two different numerical treatments. Our data again clearly show that $`ϵ_{\mathrm{eff}}^{}`$ is negative for $`f<f_p`$. As it is shown in the inset, the transmission losses are also small.

We believe that the present data are more accurate, than those published previously in PGN , not only because they agree with the theoretical analysis but also because the TMM gives the reflection and its phase straightforward. This is important because the main difference between our results and those of PGN  seems to be in the estimation of reflection $`R=|r|^2`$. When comparing our data for absorption, given in the inset of fig. 4 with those given in fig 2b of Ref. PGN , we see that our absorption is much less than that estimated in Ref. PGN .

In conclusion, we analyzed numerically the transmission properties of a periodic arrangements of thin metallic wires. From the transmission and reflection data we calculate the effective permittivity and plasma frequency, which agrees qualitatively with theoretical predictions. Both the effective permittivity and the absorption data confirm that the array of thin metallic wires, used in the recent experiments on the left-handed metamaterials, indeed behaves as a negative permittivity medium with low losses.

Acknowledgments: This work was supported by Ames Laboratory (Contract. n. W-7405-Eng-82). Financial support of DARPA, NATO (Grant No. PST.CLG.978088), APVT (Project No. APVT-51-021602) and EU project DALHM are also acknowledged.


text:
# Nonlinear Dynamics of a Bose-Einstein Condensate in a Magnetic Waveguide

## Abstract

We have studied the internal and external dynamics of a Bose-Einstein condensate in an anharmonic magnetic waveguide. An oscillating condensate experiences a strong coupling between the center of mass motion and the internal collective modes. Due to the anharmonicity of the magnetic potential, not only the center of mass motion shows harmonic frequency generation, but also the internal dynamics exhibit nonlinear frequency mixing. Thereby, the condensate shows shape oscillations with an extremely large change in the aspect ratio of up to a factor of 10. We describe the data with a theoretical model to high accuracy. For strong excitations we test the experimental data for indications of a chaotic behavior.

Since the first realization of Bose-Einstein condensates in 1995 ideas have been discussed how to exploit this new form of coherent matter for atom optics. Atom interferometry with single atoms has been extensively investigated in the past decade Berman1997a  and finds important applications for high precision measurements of gravitational or rotational forces Riehle1991a ; Kasevich1991a . It was also possible to realize and demonstrate a number of atom optical elements such as mirrors Johnson1998a , gratings Rasel1995a , beamsplitters Cassettari2000a  and waveguides Key2000a ; Muller1999a . It is a fascinating vision to integrate these elements on a microfabricated atom chip. Extremely sensitive chip based sensors for forces and accelerations as well as for applications in quantum computing are conceivable. Since Bose-Einstein condensates have been loaded into microfabricated traps Ott2001a ; Hansel2001b  this vision has changed into a very real research topic. Integrated atom optics with condensates combine the nonlinear interaction of a Bose-Einstein condensates with trapping potentials, whose geometry is in general anharmonic. This is in contrast to both single atom experiments, where the interaction plays no important role, and conventional experimental arrangements for condensates, where the trapping potential is typically harmonic and consequently the center of mass (CM) motion is decoupled from the internal dynamics Dobson1994a . In microtraps one therefore expects not only a coupling between the external and internal dynamics but also pronounced nonlinear effects including even chaotic dynamics for large excitation amplitudes. The study of the dynamics of a condensate in micropotentials is an extension of previous work on collective excitations and provides essential information for future applications for integrated atom optics.

We present a detailed experimental study of the condensate dynamics in the most fundamental atom optical element — a magnetic waveguide. We regard waveguides which are closed at both ends such that the generic motion of the condensate is a periodic oscillation. The linear propagation of Bose-Einstein condensates along a non terminated guide was recently investigated in Refs. Leanhardt2002a  and Fortagh2002c . We have measured the time evolution of the CM motion and of the aspect ratio i. e. the ratio between the radial and the axial condensate radii. In contrast to previous experimental work on collective modes of a condensate, we realize a closed system with conserved energy and a self-sustained excitation scheme. This system allows us to realize high amplitude oscillations of the aspect ratio which exceed previously reported values by far. To describe the experimental data we have developed a theoretical model which is described in Ref. Ott2002b . The model is based on the separation of the CM motion and reduces the internal dynamics to that of a condensate in a fictitious time-dependent harmonic potential. Such a scenario is known to be exactly solvable in Thomas-Fermi approximation Kagan1997b ; Castin1996a . The model also allows us to identify an experimentally accessible regime where the dynamics become chaotic.

The waveguide potential is generated with microfabricated copper conductors at the surface of a chip. The chip is mounted upside down in a vacuum chamber at a background pressure of $`10^{11}`$ mbar. The 87Rb atoms are collected in a six-beam magneto-optical trap, pumped into the $`(F=2,m_\mathrm{F}=2)`$ ground state and trapped in a spherical magnetic quadrupole trap. The atoms are then adiabatically transferred into the waveguide potential at the chip surface by gradually changing the magnetic field geometry Ott2001a ; Fortagh2002b . By forced evaporation, Bose-Einstein condensates with 100.000 atoms are produced. For the presented experiments, the trapping frequencies are adjusted to $`110`$ Hz in the radial and $`8`$ Hz in the axial direction. The shape of the axial potential is sketched in Fig. 1a. The waveguide potential is characterized by an almost constant radial oscillation frequency (the maximal deviations are smaller than 10 percent) and a strongly varying axial curvature. With a numerical calculation of the magnetic field, we determine the axial potential $`U(x)`$ to

$$\frac{U(x)}{\mu }\left[G\right]=0.219x^20.013x^30.033x^40.007x^5,$$
(1)
where $`x`$ has to be quoted in millimeter and $`\mu `$ is the magnetic moment of the atom. To induce the oscillation we displace the condensate within $`400`$ ms with a magnetic field gradient and switch off the field gradient non-adiabatically ($`<60\mu \mathrm{s}`$). The experiment was performed for three different initial displacements $`d`$: $`0.45`$ mm (series A, duration $`t_\mathrm{A}=1.015`$ s), $`0.63`$ mm (B, $`t_\mathrm{B}=2.7`$ s) and $`0.79`$ mm (C, $`t_\mathrm{C}=0.765`$ s). In Fig. 1b the first $`1.5`$ oscillation periods of series C are shown. To determine the CM position and the aspect ratio, the absorption images are fitted to a Thomas-Fermi profile. The aspect ratio bears the advantage of being independent from the atom number. To exclude influences of the nearby copper conductors Fortagh2002d ; Kraft2002a , the experiment was carried out in a distance to the surface of $`220`$$`\mu \mathrm{m}`$.

Our theoretical model describes the condensate with a set of two coupled differential equation systems. The first system are the classical equations of motion for the CM position of the condensate

$$m\ddot{𝐑}=U(𝐑),$$
(2)
where $`m`$ is the mass of the atom, $`𝐑`$ is the CM position and $`U`$ is the confining potential. In Thomas-Fermi approximation, the density distribution of a condensate in a time-dependent harmonic potential can be written as follows Castin1996a :

$$n(𝐱,𝐭)=\frac{\mu }{g}\frac{1}{\lambda _1(t)\lambda _2(t)\lambda _3(t)}\left(1\underset{i=1}{\overset{3}{}}\left(\frac{x_i}{r_{i0}\lambda _i(t)}\right)^2\right),$$
(3)
where $`\mu `$ is the chemical potential and $`g`$ is the coupling constant. The $`\lambda _i(t)`$ contain the complete time evolution of the Thomas-Fermi radii $`r_i(t)=r_{i0}\lambda _i(t)`$. The $`\lambda _i`$ are governed by three ordinary differential equations

$$\ddot{\lambda }_i=\frac{\omega _{i0}^2}{\lambda _i\lambda _1\lambda _2\lambda _3}\omega _i^2(t)\lambda _i,i=1,2,3.$$
(4)
Here, $`\omega _i(t)`$ are the time-dependent oscillation frequencies of the harmonic potential. If the modulation of the potential starts at $`t=0`$, the initial conditions for the parameters are $`\omega _i(0)=\omega _{i0}`$ and $`r_i(0)=r_{i0}`$ with $`\lambda _i(0)=1`$. The anharmonicity of the potential now provides a coupling of the CM motion and the time evolution of the Thomas-Fermi radii, which is given by

$$\omega _i^2(t)=\frac{1}{m}\frac{^2}{x_i^2}U(𝐑),i=1,2,3.$$
(5)
To compare the experimental data with the model, (2), (4) and (5) are numerically solved with the initial conditions $`\lambda _i(0)=1`$ and $`𝐑(0)=(d,0,0)`$.

The experimental data for the CM motion of series B are presented in Fig. 1c . We cannot extract any damping and the data are consistent with a ($`1/e`$) – lifetime of at least 20 minutes. Due to the anharmonicity of the potential (1), the frequency spectrum of the CM motion shows second and third harmonic generation. For series C with the largest initial displacement, the nonlinear effects are most pronounced and lead to a relative amplitude of 5 percent for the second harmonic and 1 percent for the third harmonic frequency with respect to the amplitude for the fundamental frequency $`\nu _0`$. For small amplitudes $`\nu _0`$ amounts to $`8.5`$ Hz but is reduced with increasing displacement (A, B, and C) by 4, 6, and 11 percent respectively. Because of the extremely low damping, the CM motion is suited to detect small forces, which may shift the trapping potential. An estimate for our experimental parameters indicates a resulting sensitivity of $`F/m=10^4g`$, with the gravitational acceleration $`g`$.

In Fig. 2a-2d the measured aspect ratio is compared with the theoretical model. While the theory contains no free parameter, a constant multiplication factor of $`1.2`$ is needed to describe the data accurately. This factor arises from the fact, that the condensate is radially compressed during the switch-off of the magnetic trap which — due to the mean field repulsion — leads to a larger radial size in the time-of-flight image switchoff . The agreement between experiment and theory is remarkable and proves true the picture that the external motion in an anharmonic potential leads to a time-dependent harmonic potential in the rest frame. In Fig. 1b, the condensates marked with an asterisk cannot be spatially resolved by our absorption imaging system, indicating an axial extension of the condensate of less than $`15\mu \mathrm{m}`$. Therefore, the condensate may in fact follow the strong dynamics, which is predicted by the theory (Fig 2d). For the present data, we measure a change in the aspect ratio by more than a factor of 10, the largest value reported so far. The data of series B which cover $`2.7`$ s are promising to find effects arising from loss of atoms (Fig 2e). At the end of the series, the atom number has dropped by more than a factor of 15 with only 3000 atoms left in the condensate. The systematic increase of the aspect ratio is due to the crossover to the 1D regime as defined in Ref. Goerlitz2001a . The Thomas-Fermi approximation breaks down because the chemical potential of the condensate becomes smaller than the radial kinetic energy, which then dominates the radial expansion of the condensate.

We have performed a frequency analysis of the aspect ratio with an algorithm suitable for unequally spaced or missing data Lomb1976a . The results are shown in Fig. 3. The spectra show three major frequencies: the fundamental frequency $`\nu _0`$ (peak (i) in Fig. 3), the second harmonic (ii) and the resonance frequency $`\nu _q`$ of the lowest lying collective mode (iii). The latter does not deviate more than 4 percent from its small amplitude value for strongly anisotropic traps which is given by $`\sqrt{5/2}\nu _0`$. Furthermore, the spectrum for series C (Fig. 3) shows several mixed frequencies. Two of them, the sum and the difference between the second harmonic frequency and the lowest resonance frequency $`2\nu _0\nu _q`$ and $`2\nu _0+\nu _q`$ have been indentified ((iv) and (v)). This gives evidence for the nonlinear coupling of the collective excitations. Other than the coupled modes observed by Hechenblaikner et al. Hechenblaikner2000a , the mode coupling in our experiment is non-resonant, and the mixed frequencies vanish in the limit of small amplitudes. Because the condensate performs a free oscillation in an external potential the total energy of the system is conserved. Summing up only the contributions for the internal energy of the condensate, one finds that the sum is not constant. Adding the kinetic and potential energy of the CM motion, the energy is conserved, which can be verified by using a Hamiltonian formalism Ott2002b . Therefore, not only the internal dynamics is affected by the external motion, but there is also an energy exchange between the internal dynamics and the CM motion. In our experiment this effect leads to a theoretically calculated amplitude variation for the CM motion of $`1/4000`$, which is too small to be detected.

Fig. 4a shows absorption images for an initial displacement of $`0.94`$ mm. The data show a strong axial excitation of the condensate, which raises the question of the onset of chaos. A signature of a chaotic behavior in experimental data is an enhanced noise in the frequency spectrum of the data. The spectrum of the aspect ratio for series C in Fig. 3 clearly shows an enhanced background, compared to the two other series. However, this should not be regarded as an unchallengeable experimental proof. Such a proof would require a longer data set and an enhanced resolution of the imaging system. Nevertheless, the theoretical model can be analyzed for the experimental parameters. One of the strongest evidences for a chaotic dynamic is a positive Lyapunov exponent, which indicates an exponential sensitivity to the initial conditions. For this purpose neighboring trajectories for slightly different initial displacements have been simulated over 200 s. Fig. 4 shows the evolution of the Lyapunov exponent with increasing displacement. For a displacement of more than $`0.63`$ mm, the calculated Lyapunov exponent becomes positive and grows linearly. As the total energy of the system is conserved, one can easily test the numerical stability and exclude any numerical artefacts. The corresponding Poincaré maps and the frequency spectra of the trajectories indicate the same results. This gives strong evidence, that series C was performed in a regime where the longtime behavior of the trajectory is chaotic.

In conclusion, we have studied the dynamics of a Bose-Einstein condensate oscillating in an anharmonic trapping potential. We have identified several nonlinear effects including second and third harmonic generation of the center of mass motion, a coupling of the internal and external dynamics and nonlinear mode mixing. The theoretical model describes these effects with high accuracy. The model predicts a chaotic dynamic above a critical oscillation amplitude and the experimental data give indication for a chaotic motion too. Considering applications of magnetic microtraps and integrated atom optics, we have shown that the motion of a condensate in a magnetic waveguide necessarily leads to excitations of the inner degrees of freedom. A thorough understanding of these effects plays a key role to exploit the nonlinearity to build high precision sensors. It also shows fundamental limits for the use of Bose-Einstein condensates for technical applications.

We gratefully acknowledge financial support from the Deutsche Forschungsgemeinschaft under Grant No. Zi/419-5.