text:
# How different Fermi surface maps emerge in photoemission from Bi2212

## Abstract

We report angle-resolved photoemission spectra (ARPES) from the Fermi energy ($`E_F`$) over a large area of the ($`k_x,k_y`$) plane using 21.2 eV and 32 eV photons in two distinct polarizations from an optimally doped single crystal of Bi2Sr2CaCu2O8+δ (Bi2212), together with extensive first-principles simulations of the ARPES intensities. The results display a wide-ranging level of accord between theory and experiment and clarify how myriad Fermi surface (FS) maps emerge in ARPES under various experimental conditions. The energy and polarization dependences of the ARPES matrix element help disentangle primary contributions to the spectrum due to the pristine lattice from those arising from modulations of the underlying tetragonal symmetry and provide a route for separating closely placed FS sheets in low dimensional materials.

\]

Establishing the nature of the Fermi surface (FS) in the metallic state–the shape, size, connectivity and the number of FS sheets and the character of states at and near the Fermi energy ($`E_F`$), is crucially important for identifying the correct superconducting scenario in the cuprates. In this connection, high-resolution angle-resolved photoemission spectroscopy (ARPES) provides a unique window in low dimensional materials. However, even under identical conditions of doping, temperature and other physical parameters, the imprint of the FS in ARPES often varies dramatically with the energy and polarization of the incident photons. This article reports fundamental progress in understanding this remarkable behavior by considering the example of Bi2212 which has played the role of a touchstone and a workhorse in ARPES studies of the high-Tc’s. By comparing ARPES measurements over a large area of the FS with corresponding first-principles simulations, we show how myriad tapestries of FS maps originate in photoemission, and how the ARPES matrix element blends and highlights various parts of the electronic spectrum –pristine states, effects of modulations, different FS sheets and the like, depending upon the specifics of the photon energy and polarization.

High quality single crystals of optimally doped Bi2212 ($`T_c`$=90 K) were grown by the traveling solvent floating zone method. The FS maps were produced by measuring the photointensity within a narrow energy window at $`E_F`$ over the full 360 angular range in the ($`k_x,k_y`$) plane at the SU8 high-resolution beamline of the Super-Aco synchrotron ring at LURE.  The data were collected on a regular mesh in polar ($`\theta `$) and azimuthal angles ($`\varphi `$) of steps $`\mathrm{\Delta }\theta =1^o`$, $`\mathrm{\Delta }\varphi =1^o`$. Some measurements employed He-I light (unpolarized) from a He discharge lamp. The overall energy resolution varied from 12 meV to 60 meV depending on measurement conditions. Procedures to assure accurate sample alignment included a complete set of photoelectron diffraction azimuthal and polar scans of the Bi core level in the high-energy regime which were recorded insitu in order to define the $`\overline{\mathrm{\Gamma }}`$ point and the main high symmetry directions.  Two different detection geometries – referred to as ‘even’ and ‘odd’ – have been employed to obtain the ARPES spectra. In the even case, the detector is constrained to move in the horizontal plane, which is the same as the plane of the incident light, whereas for the odd case, the detector is chosen to move in the vertical plane. The polarization is defined as measured with respect to radial lines through $`\mathrm{\Gamma }`$: even (odd) polarization implies the polarization vector lies along (transverse to) this radial direction. It can then be shown that when emission from a mirror plane in the crystal lattice is involved (e.g. along the $`\mathrm{\Gamma }M`$ line), the even (or odd) experiment only probes initial states of even (or odd) symmetry. In the case of emission from the Cu-O planes in Bi2212, the relevant initial states are predominantly Cu $`3d_{x^2y^2}`$ hybridized with O $`2p_{x,y}`$ and states along the $`\overline{\mathrm{\Gamma }}\overline{M}`$ line (along the Cu-O-Cu bonds) are expected to be intense in the even experiment, while states along $`\overline{\mathrm{\Gamma }}\overline{X}(\overline{Y})`$ will be strong in the odd measurements (since deviations from tetragonal order are small).

All computations in this article are based on the one-step model of photoemission extended to treat arbitrarily complex unit cell materials . Effects of multiple scattering and the ARPES matrix element  are thus included realistically in the presence of a specific surface termination, taken here to be the Bi-O layer. The finite lifetimes of the initial and final states are incorporated by giving suitable imaginary parts to the self-energies of the associated initial and final state propagators. The crystal potential was obtained via a selfconsistent KKR procedure  assuming a perfect tetragonal lattice which yielded the well-known LDA-based band structure and FS of Bi2212. The actual potential used in this work however has been modified slightly such that the Bi-O pockets around the $`\overline{M}`$-point are lifted above $`E_F`$ to account for their absence in the experimental ARPES spectra. A number of simulations have also been carried out where repulsive barriers were placed between the two CuO2 planes in the Ca layers in order to mimic correlation effects beyond the LDA which are expected to reduce the bilayer splitting. Finally, in order to describe the effects of the orthorhombic modulation and the $`(1\times 5)`$ superstructure present in Bi2212, we have taken the spectrum computed for the tetragonal case and superposed on this the same spectrum with appropriate weights after suitable translations in the reciprocal space. The weights estimated from relative intensities in the experimental FS maps are summarized in Table 1 and described in the text below.

Figure 1 considers two sets of theoretical spectra at 21.2 eV in the even and odd detection geometries for the pristine tetragonal lattice, and helps set the stage for our discussion. The right hand Figs. 1(b,d) refer to simulations where the bilayer splitting has been artificially reduced to a near zero value and both FS sheets are hole-like around the $`\overline{X}(\overline{Y})`$ symmetry point. The case with a bilayer splitting of $``$ 200 meV where the first FS sheet is hole-like and the 2nd sheet is slightly electron like is depicted in the left hand Figs. 1(a,c). We see by comparing Figs. 1(a) and 1(b) that there is a striking increase in the intensity around the $`\overline{M}`$-point for even detection due to the presence of the aforementioned 2nd underlying FS sheet. In fact, we have carried out numerous other simulations in which the size of the bilayer splitting and/or the position of $`E_F`$ around the $`\overline{M}`$-point was varied; these results show clearly that about the only way in which one can obtain a large intensity around the $`\overline{M}`$-point in the even detection geometry for 21.2 eV photons is when there is an additional band present very close to $`E_F`$ around the $`\overline{M}`$-point. This is a crucially important observation which we will recall below in analyzing the experimental data. Fig. 1 also highlights the differences between the even and odd maps more generally. The 2nd FS sheet so prominent in the even map of Fig. 1(a) is essentially ‘invisible’ in the corresponding odd map of Fig. 1(c). Also, the odd map is complementary to the even map in the sense that regions of high intensity in one map are often replaced by those of low intensity in the other and vice versa.

The determination of the FS topology of the underlying pristine phase in Bi2212 has been complicated by the presence of secondary features tied to superlattice and shadow bands. The results of Table 1 offer insight into how a headway can be made in disentangling primary and secondary features by exploiting the energy and polarization dependences of the ARPES matrix element. We see that in the 21.2 eV even geometry, no secondary feature displays a weight greater than 20 $`\%`$ of the primary features, and therefore, this energy and polarization is well-suited for delineating the primary FS. In sharp contrast, the weights of secondary features are significantly greater at 32 eV in both the even and odd spectra as well as in the 21.2 eV odd spectrum. Interestingly, the shadow feature is most prominent in the 32 eV odd case, and the second superlattice image is visible in 32 eV even map, albeit weakly.

Although absolute intensities are not our focus, one aspect of these is particularly germane to this discussion. The theoretically computed maximum absolute ARPES intensities of the primary features are in the ratios: 21.2 eV (even): 21.2 eV (odd): 32 eV (even): 32 eV (odd) = 5:2:1.2:0.7, in reasonable accord with the corresponding experimental values of: 5:1.7:0.8:0.6, uncertainty in determining absolute experimental intensity notwithstanding. Applying these ratios to the data of Table 1 \[which gives the relative weights of the secondary features\], it is found that the absolute intensities of various secondary features are not all that different at the two photon energies considered for either the odd or the even maps. In other words, the primary emissions become highlighted in the 21.2 eV even measurements not because the secondaries are weak, but due to the primaries becoming much stronger.

Figures 2 and 3 compare experimental and theoretical spectra at 21.2 eV and 32eV respectively. While changes in the ARPES maps of the FS with photon energy and polarization are striking, an excellent overall agreement with theory is seen in all cases. The unpolarized (He I) data (Fig. 2(a)) and theory (Fig. 2(b)) both show a clear imprint of the large hole sheet centered around $`\overline{X}`$ or $`\overline{Y}`$, but relatively little trace of the 2nd FS sheet around $`\overline{M}`$. However, under even polarized light (Figs. 2(c,d)) this 2nd sheet comes alive giving intense emission around the $`\overline{M}`$-points. As emphasized already in connection with Fig. 1 above, the appearance of this large intensity is an unmistakable signature of the presence of a second band lying close to EF at $`\overline{M}`$. As noted above, under these experimental conditions secondary features are quite weak (Table 1) and can be ignored in the analysis of Figs. 2(c,d). We see thus how a systematic comparison between experimental and theoretical ARPES maps shows that the 21 eV even detection geometry is well suited for establishing unambiguously the presence of the second FS sheet with negligible interference from umklapps and shadow bands.

Under odd polarization (Figs. 2(e,f)) secondary features develop greater relative intensity, so the even and odd maps look quite different. As expected, the emission around $`\overline{M}`$ from the 2nd sheet is suppressed in this polarization. On the other hand, the trace of the FS sheets is reinforced along the lattice modulation direction via superposition with its Umklapp images, thereby distorting the image of the hole sheet in the odd map (Fig. 2(e)), so that caution must be exercised in deducing physical parameters (size and shape) from this spectrum. The mix of primary and secondary features is quite different in the 32eV data, Fig. 3, even though some of the characteristics are basically similar. For example, the 2nd FS sheet is emphasized in the even polarization in Figs. 3(a,b), but suppressed in the odd case of Figs. 3(c,d). Compared to 21.2 eV, there is generally a greater spectral weight at higher momenta; for example, the image of the 2nd sheet is more or less equally intense in the 1st and 2nd BZ’s in Figs. 3(a,b), but in the 21.2 eV Figs. 2(c,d), the 2nd BZ imprint is quite weak in relation to that in the 1st BZ. The combination of the two FS sheets with appropriately weighted secondary images yields two nearly parallel bands of intense emission which are oriented along the ($`\pi ,\pi `$) direction in Fig. 3 and appear quite striking. Interestingly, in even polarization (Figs. 2(c,d), or 3(a,b)) the spectra show a remarkable line of essentially zero intensity along the $`\overline{\mathrm{\Gamma }}\overline{Y}`$ direction, as expected for a mirror plane. Along the perpendicular line $`\overline{\mathrm{\Gamma }}\overline{X}`$, there is clear residual intensity – a signature of superlattice effects which break the symmetry selection rules.

In conclusion, we demonstrate that dramatic differences in the observed FS maps under various experimental conditions arise as a consequence of the ARPES matrix element which gives the primary FS features (arising from the pristine lattice) a weight which depends upon the energy and polarization of light in a very different manner from that of the secondary features associated with the superlattice and orthorhombic modulations. In this way the ARPES matrix element helps disentangle various aspects of the electronic structure and fermiology of this complex system. In particular, we confirm recent reports of two FS sheets in Bi2212, extending the result to optimal doping. There is one large hole-like standard sheet and a second electron-like sheet which arises from a band which lies very close the Fermi energy at the $`\overline{M}`$-point. The wide ranging accord between theory and experiment suggests that the band theory framework implicit in our computations captures the essential underlying physics. If so, the two FS pieces should be viewed as the result of conventional bilayer splitting. However, it would be wise to keep the door open for admitting other explanations, e.g. some form of a nano-scale phase separation or stripe order.

###### Acknowledgements.

We acknowledge support by the Spanish agency DGES under grant PB-97-1199, the U.S.D.O.E. Contracts W-31-109-ENG-38 and AC02-98CH10886, the large scale Installation Program of the EU, and the allocation of supercomputer time at the NERSC and the Northeastern University (NU-ASCC) Computation Centers.

Figure Captions

Fig. 1

Theoretical maps at 21.2 eV for emission from $`E_F`$ in pristine (tetragonal) Bi2212 for even and odd polarizations of light described in the text. The left hand simulations (a,c) refer to the band theory-based crystal potential with a bilayer splitting of $``$ 200 meV at the $`\overline{M}`$ point. In the right hand side simulations (b,d), the bilayer splitting has been artificially reduced to a nearly zero value. Intensities are normalized to the same value in all maps and plotted in hot colors (whites and yellows are highs) on a linear scale. The square in solid white lines marks the 2D Brillouin zone.

Fig. 2

Theory and experiment are compared directly at 21.2 eV for even and odd polarizations, as well as for the unpolarized He I light. Effects of superlattice and orthorhombic modulations are included in theory as discussed in the text. See caption to Fig. 1 for other pertinent details.

Fig. 3

Same as Fig. 2, except that this figure refers to 32 eV synchrotron light.


format: text

text:
# Comment on “Effects of Disorder on Ferromagnetism in Diluted Magnetic Semiconductors”

\[

\]

In a recent Letter, Berciu and Bhatt  have presented a mean-field theory of ferromagnetism in III-V semiconductors doped with Mn . The approach starts from uncorrelated Mn impurity states, which overlap to form an impurity band (IB). We argue here that this approach is not appropriate for (Ga,Mn)As containing 1–5% Mn. There is also a sign error in Ref. . After correcting this error, mean-field theory no longer yields stable ferromagnetism for the system sizes used.

For the overlap matrix elements between hole impurity states, Ref.  gives $`t(r)=2E_b(1+r/a_B)e^{r/a_B}`$ with hole binding energy $`E_b`$ and Bohr radius $`a_B`$. The calculations are done with $`t(r)>0`$. We have reproduced the numerical results of Ref. , for example, the dotted curves in the inset in Fig. 1 give the Mn and hole spin polarizations for three impurity configurations with a Mn concentration of $`x=0.05`$ and $`p=0.1`$ holes per Mn, using $`200`$ Mn spins. These are evidently the same parameters as used for Fig. 3 of Ref. .

However, the overlap matrix elements $`t(r)`$ should be negative, as we show now, using the sign convention of Ref. . For a parabolic valence band, the Hamiltonian for a hole and two impurities at $`𝐑_A`$ and $`𝐑_B`$ reads $`H=p^2/2m^{}+V_A(𝐫)+V_B(𝐫)`$ with $`V_i(𝐫)=e^2/ϵ|𝐫𝐑_i|`$ where $`i=A,B`$. $`m^{}`$ is the effective mass and $`ϵ`$ is the dielectric constant. The hydrogen-like Hamiltonian for a single impurity is $`H_i=p^2/2m^{}+V_i(𝐫)`$. Let $`\psi _i`$ be the ground state of $`H_i`$ with energy $`E_b<0`$. The overlap matrix element between $`\psi _A`$ and $`\psi _B`$ is $`t\psi _B|H|\psi _A`$, which we can also write as $`t=\psi _B|H_A+V_B(𝐫)|\psi _A=E_b\psi _B|\psi _A\psi _B|e^2/ϵ|𝐫𝐑_B||\psi _A`$. Since the wave functions $`\psi _i`$ are real and positive, we obtain $`t<0`$.

The sign change in $`t(r)`$ inverts the highly asymmetric IB. With the correct sign the IB has a long tail to negative hole energies (positive electron energies), resulting from strongly bound states. Fig. 1 shows the density of states (DOS) for $`x=0.05`$ with $`200`$ Mn impurities. The subband structure is due to the finite system size. We find that the IB width (independent of the sign of $`t`$) is about $`18\mathrm{eV}`$, and thus much larger than the gap of GaAs. Furthermore, the IB strongly overlaps with the valence band. These results show that the IB picture is not appropriate. The reason is that there are typically many Mn impurities within the range over which $`t(r)`$ falls off, i.e., the system is not weakly doped .

With the corrected $`t(r)`$ the Fermi energy no longer lies in a region of very large DOS as in Ref. , but in the band tail close to the conduction band edge, as shown in Fig. 1 for $`p=0.1`$. Average Mn and hole spin polarizations for $`x=0.05`$ and $`p=0.1`$ are shown in the inset. The maximum Curie temperature is only about $`20\mathrm{K}`$ and the stability of the ferromagnetic state strongly depends on the impurity configuration. If ferromagnetism is found at all, the hole polarization is at most $`0.05`$. For the numerical parameters this means that only a single hole spin is flipped, suggesting large finite-size effects.

To conclude, the impurity-band approach is not applicable to (Ga,Mn)As containing a few percent of Mn. However, it should be valuable for weak doping. We would like to thank J. König, M. Raikh, and J. Schliemann for interesting discussions.


text:
# Deterministic diffusion in one-dimensional maps – Calculation of diffusion constants by harmonic inversion of periodic orbit sums

## 1 Introduction

Deterministic diffusion plays an important rôle in a variety of physical phenomena and applications, e.g. particle confinement in magnetic fields , conductivity in metals , etc. A general theory of diffusion for systems with continuous time is presented in . In systems with periodic phase-space structure, diffusion may be induced by chaotic dynamics inside the elementary cells, causing an irregular jumping between the cells. The investigation of deterministic diffusion is an interesting task both from the theoretical point of view, concerning the question of how the typical stochastic properties arise from a purely deterministic dynamics, as well as with regard to its various applications in different physical contexts.

It has been observed that even the simplest model systems – one-dimensional maps – which exhibit deterministic diffusion are capable of explaining the behaviour of relevant physical systems . For such maps, it has been demonstrated that the diffusive properties are related to the periodic orbits of the map restricted to an elementary cell . In particular, the diffusion constant characterizing the diffusive spreading can be obtained from the leading zero of a zeta function consisting of a product over contributions from all periodic orbits of the reduced map. Similar to the dynamical zeta functions which appear in semiclassical quantization problems, this product usually does not converge. In former work , the convergence problem was handled by the application of cycle expansion techniques. However, cycle expansion strongly depends on the existence of a symbolic dynamics and the shadowing of long orbits by combinations of short ones.

On the other hand, it has been demonstrated  that the convergence problems of the zeta functions and corresponding response functions arising in semiclassical quantization problems can be circumvented by the application of harmonic inversion, which is a universal method in the sense that it does not rely on a symbolic code or other special properties of the system. The method has recently been extended from Hamiltonian systems to maps . In this paper we demonstrate that harmonic inversion can also be used for the calculation of diffusion constants of one-dimensional maps. In analogy with the procedure for semiclassical quantization, we rewrite the zeta function as a response function, from which we construct a periodic orbit signal including a finite set of periodic orbits. The analysis of the signal by harmonic inversion yields the zeros of the underlying zeta function, from which the diffusion constant is calculated. We apply this technique to an example map for which the periodic orbits can easily be calculated, and which has already been examined with the help of cycle expansion techniques , rendering a comparison with previous results possible.

We start by describing the general procedure. We consider one-dimensional maps $`\widehat{f}:`$ defined by an anti-symmetric function $`\widehat{f}`$,

$$\widehat{x}_{n+1}=\widehat{f}(\widehat{x}_n),\widehat{f}(\widehat{x})=\widehat{f}(\widehat{x})$$
(1)
for points $`\widehat{x}`$ inside the elementary cell $`[1/2,1/2)`$, and a discrete translational symmetry outside the elementary cell,

$$\widehat{f}(\widehat{x}+n)=\widehat{f}(\widehat{x})+n,n.$$
(2)
From this map, a reduced map $`f:[1/2,1/2)[1/2,1/2)`$ can be constructed by restricting the dynamics to the elementary cell via the definition

$`x`$ $`=`$ $`\widehat{x}[\widehat{x}+1/2]`$ (3)
$`f(x)`$ $`=`$ $`\widehat{f}(x)[\widehat{f}(x)+1/2],`$ (4)
where $`[]`$ denotes the integer part of the quantity in the brackets. Each periodic orbit of the reduced map $`f(x)`$, consisting of a sequence of fixed points $`\{x_1,x_2,\mathrm{},x_n\}`$, corresponds to an orbit of the full map $`\widehat{f}(\widehat{x})`$ which returns to an equivalent point,

$$\widehat{f}^n(x_j)=x_j+\sigma _p,\sigma _p.$$
(5)
The diffusion constant associated with the map is determined by the lattice translation $`\sigma `$ after $`n`$ applications of the map, averaged over initial conditions in the elementary cell

$$D=\underset{n\mathrm{}}{lim}\frac{1}{2n}\sigma ^2_n.$$
(6)
As was shown in Ref. , for hyperbolic or almost hyperbolic maps, the diffusion constant can be determined from the periodic orbits of the reduced map (4) by considering the zeta function

$$Z(z,\beta )=\underset{m=0}{\overset{\mathrm{}}{}}\underset{p}{}\left(1\frac{z^{n_p}e^{\sigma _p\beta }}{|\mathrm{\Lambda }_p|\mathrm{\Lambda }_p^m}\right),$$
(7)
where the second product runs over all primitive periodic orbits of the reduced map; $`n_p`$ and $`\mathrm{\Lambda }_p=_{i=1}^nf^{}(x_i)`$ are the topological length and the stability of the orbit, respectively, $`\sigma _p`$ denotes the lattice translation of the corresponding orbit of the full map, and $`\beta `$ is a free parameter. According to Ref. , the diffusion constant $`D`$ can be expressed in terms of the leading zero $`z_0(\beta )`$ of the zeta function as a function of $`\beta `$,

$$D\frac{1}{2}\left[\frac{^2}{\beta ^2}z_0(\beta )\right]_{\beta =0}.$$
(8)
Additionally, from the symmetry of the map, it follows that $`[z_0/\beta ]_{\beta =0}=0`$. Using these relations, an expansion of $`\mathrm{ln}z_0(\beta )`$ for small $`\beta `$ yields

$$\mathrm{ln}z_0(\beta )\beta ^2D,$$
(9)
which can directly be used to calculate the diffusion constant from the leading zero for small $`\beta `$.

The problem that remains is how to actually calculate the leading zero, since the expression (7) for the zeta function usually diverges. In past works, this problem has been tackled by application of cycle expansion techniques. As an alternative, we demonstrate in the following how the zeta function can be expressed in terms of a response function, and how the leading zero can then be obtained by harmonic inversion of a periodic orbit signal constructed from this response function.

## 2 Harmonic inversion for maps

We start by defining the quantity $`w=i\mathrm{ln}z`$ and rewriting the zeta function (7) as a function of $`w`$,

$$Z(w,\beta )=\underset{m=0}{\overset{\mathrm{}}{}}\underset{p}{}\left(1\frac{e^{iwn_p}e^{\sigma _p\beta }}{|\mathrm{\Lambda }_p|\mathrm{\Lambda }_p^m}\right).$$
(10)
The density of zeros of the zeta function (10) can be written in the form

$$\rho (w,\beta )=\frac{1}{\pi }\mathrm{Im}g(w,\beta ),$$
(11)
where the response function $`g(w,\beta )`$ is given by

$$g(w,\beta )=g_0(w,\beta )+\frac{}{w}\mathrm{ln}Z(w,\beta )$$
(12)
with $`g_0(w,\beta )`$ a smooth, slowly varying function of $`w`$. Using the definition

$$t_p^{(m)}=\frac{e^{iwn_p}e^{\sigma _p\beta }}{|\mathrm{\Lambda }_p|\mathrm{\Lambda }_p^m},$$
(13)
we rewrite the second term of the response function as

$`{\displaystyle \frac{}{w}}\mathrm{ln}Z(w,\beta )`$ $`=`$ $`{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{p}{}}(in_p){\displaystyle \frac{t_p^{(m)}}{1t_p^{(m)}}}`$ (14)
$`=`$ $`{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{p}{}}{\displaystyle \underset{r=1}{\overset{\mathrm{}}{}}}(in_p)\left(t_p^{(m)}\right)^r`$ (15)
$`=`$ $`{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}{\displaystyle \underset{p}{}}{\displaystyle \underset{r=1}{\overset{\mathrm{}}{}}}(in_p){\displaystyle \frac{e^{iwrn_p}e^{r\sigma _p\beta }}{|\mathrm{\Lambda }_p^r|\mathrm{\Lambda }_p^{rm}}}`$ (16)
$`=`$ $`{\displaystyle \underset{\mathrm{po}}{}}(in_0){\displaystyle \frac{e^{iwn_{\mathrm{po}}}e^{\sigma _{\mathrm{po}}\beta }}{|\mathrm{\Lambda }_{\mathrm{po}}|}}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}\mathrm{\Lambda }_{\mathrm{po}}^m,`$ (17)
where the first sum in (17) now runs over all periodic orbits (p.o.) including multiple traversals of primitive orbits, and $`n_0`$ denotes the topological length of the underlying primitive orbit. After carrying out the sum over $`m`$, the response function finally reads

$$g(w,\beta )=g_0(w,\beta )+\underset{\mathrm{po}}{}(in_0)\frac{e^{\sigma _{\mathrm{po}}\beta }}{|1\mathrm{\Lambda }_{\mathrm{po}}|}e^{iwn_{\mathrm{po}}}.$$
(18)
On the other hand, the response function can be written in terms of the zeros $`w_k`$ of the zeta function (10) and their multiplicities $`d_k`$,

$$g(w,\beta )=\underset{k}{}\frac{d_k}{ww_k(\beta )}.$$
(19)
Like the zeta function (10), the expression (18) for the response function will usually not converge in those regions of the complex variable $`w`$, where the poles of $`g(w,\beta )`$ are located, and the zeros of the zeta function cannot directly be calculated from (18). Instead, we apply harmonic inversion to obtain an analytical continuation . The poles $`w_k(\beta )`$ of the analytical continuation of $`g(w,\beta )`$ in (18) are the zeros of the analytical continuation of the zeta function $`Z(w,\beta )`$ in (10). To apply harmonic inversion, we now construct a periodic orbit signal by Fourier transformation of the response function (18),

$$C(s,\beta )=\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}g(w,\beta )e^{isw}𝑑w.$$
(20)
Considering only the second, oscillating part of $`g(w,\beta )`$, this results in

$$C(s,\beta )=\underset{\mathrm{po}}{}\frac{in_0e^{\sigma _{\mathrm{po}}\beta }}{|1\mathrm{\Lambda }_{\mathrm{po}}|}\delta (sn_{\mathrm{po}}).$$
(21)
The corresponding signal resulting from Fourier transformation of Eq. (19) reads

$$C(s,\beta )=i\underset{k}{}d_ke^{isw_k(\beta )},$$
(22)
which contains the zeros $`w_k(\beta )`$ of the zeta function (10) as frequencies. The zeros $`w_k(\beta )`$ can now be determined by adjusting the periodic orbit signal (21) – including all orbits up to a maximum topological length $`n_{\mathrm{max}}`$ – to the form of the corresponding signal (22) by harmonic inversion. According to Eq. (9) and with the relation $`w=i\mathrm{ln}z`$, the diffusion constant can be obtained from the leading frequency $`w_0`$ using the relation

$$D\frac{i}{\beta ^2}w_0(\beta ),$$
(23)
which for small $`\beta `$ should be independent of $`\beta `$.

## 3 Results

As an example we consider the following map, which has also been examined in Ref. : For points $`\widehat{x}`$ inside the interval $`[1/2,1/2]`$ the map is defined by

$$\widehat{f}(\widehat{x})=\widehat{x}(1+2|2\widehat{x}|^\alpha ),$$
(24)
where $`\alpha `$ is a free parameter with the restriction $`\alpha >1`$. Outside the elementary cell, the map is continued according to Eq. (2). Using Eqs. (3) and (4), a corresponding reduced map $`f(x)`$ can be constructed.

As was discussed in Ref. , the periodic orbits of the reduced map $`f(x)`$ can be described by a complete symbolic dynamics with a three-letter alphabet. For $`1<\alpha 0`$ the dynamics is completely hyperbolic. The 0 orbit (corresponding to the fixed point $`x=0`$) is infinitely unstable (i.e., $`\mathrm{\Lambda }_p=\mathrm{}`$) and does not contribute to the response function (18). For $`\alpha >0`$, the 0 orbit becomes marginally stable with $`\mathrm{\Lambda }_p=1`$, and the dynamics becomes intermittent. Furthermore, it was shown in  that at $`\alpha =1`$ the system undergoes a phase transition from normal to anomalous diffusion, and the diffusion constant as defined by Eq. (6) equals 0 for $`\alpha 1`$.

For $`\alpha >0`$, the marginally stable 0 orbit would give a diverging contribution to the zeta function (7) as well as to the periodic orbit signal (21). As was discussed in Ref. , it is appropriate to simply omit this orbit from the zeta function, as was also done in Ref. . In Ref. , the diffusion constant of the map (24) as a function of $`\alpha `$ was determined using explicit cycle expansion expressions for the average in Eq. (6) rather than making use of Eq. (8). In the following, we will determine the diffusion constant by harmonic inversion as outlined above.

We calculated the diffusion constant for different values of the parameter $`\alpha `$ in the region $`0.5\alpha 2`$ using the periodic orbits up to topological length $`N=13`$ (but omitting the 0 orbit). For each value of $`\alpha `$, we performed calculations for different small values of $`\beta `$ and checked the convergence of the results for the diffusion constant for $`\beta 0`$. The results for the diffusion constant are presented in Figure 1. In the region $`\alpha 0.5`$, our results are in excellent agreement with those from Ref.  obtained with the periodic orbits up to topological length $`N=10`$. For $`\alpha >0.5`$, the different techniques applied in Ref.  yielded quantitatively different values for the diffusion constant; here, our results show a good qualitative agreement with the different cycle expansion results. This is even true when we used the reduced set of orbits up to topological length $`N=10`$ in our calculations: The function $`D(\alpha )`$ in Fig. 1 becomes slightly more flat, but still shows a good qualitative agreement with the results of Ref. . As in , it was not possible to reproduce the theoretical value $`D=0`$ for $`\alpha 1`$, but one only obtains an asymptotical convergence to this value.

In our calculations, we observed that the condition that for given $`\alpha `$ the expression (23) for the diffusion constant should become independent of $`\beta `$ for $`\beta 0`$ was well fulfilled for $`\alpha <0`$. However, with increasing $`\alpha >0`$, the results became numerically unstable if $`\beta `$ was too close to zero, so that we had to extrapolate from the results for $`\beta `$ around $`0.5`$. The convergence of the frequencies and amplitudes in the harmonic inversion procedure itself was also very good for $`\alpha <0`$, but convergence was increasingly hard to obtain for $`\alpha >0`$. A criterion for the reliability of the results for a given value of $`\alpha `$ is the closeness of the leading frequency for $`\beta =0`$ to its theoretical value $`w_0=0`$, as well as the closeness of the corresponding amplitude $`d_0`$ (see Eq. (22)) obtained by harmonic inversion to its theoretical value 1 (the multiplicity of the leading zero of the zeta function). The absolute value of the leading frequency for $`\beta =0`$ and the corresponding amplitude as a function of $`\alpha `$ are presented in Figure 2. The absolute value of the leading frequency was smallest for $`\alpha <0`$ but stayed below $`|w_0|10^3`$ for all values of $`\alpha `$. The corresponding amplitudes were nearly exactly equal to their theoretical value 1 for $`\alpha <0`$, but began to deviate from their theoretical value for $`\alpha 0.3`$, indicating that the harmonic inversion results became less well converged. A possible reason for the deviations may be that the density of frequencies becomes too large for the signal length chosen. It can be expected that the accuracy of the results will be improved by including longer orbits.

In Ref.  it was proposed to truncate the periodic orbit sum by introducing a stability cutoff and to use the stability as an ordering parameter for the cycle expansion. It is worth noting that the harmonic inversion procedure has the advantage that no discrete ordering parameter is needed at all. However, a stability cutoff can also be used with the harmonic inversion technique because very unstable orbits give only a small contribution to the signal which can be neglected.

## 4 Conclusion

In conclusion, we have demonstrated how the diffusion constant of one-dimensional diffusive maps can be determined from a finite set of periodic orbits using harmonic inversion techniques. Starting from a zeta function whose leading zero possesses a known relation to the diffusion constant, we have constructed a periodic orbit signal which was analysed by harmonic inversion. The diffusion constant could then be determined from the leading frequency of the signal. Since, in contrast to other methods such as cycle expansion, harmonic inversion does not depend on any special properties of the map, the general procedure should work in the same way for all chaotic maps. We have tested the method for a simple map whose periodic orbits can easily be determined and for which results for the diffusion constant obtained by cycle expansion techniques are available in the literature. Our results are in excellent agreement with the ones from cycle expansion.

## Acknowledgement

We are grateful to R. Grauer for turning our attention to the subject. This work was supported by the Deutsche Forschungsgemeinschaft and Deutscher Akademischer Austauschdienst.


text:
# Influence of correlations on the velocity statistics of scalar granular gases

## I Introduction

Some twenty years ago, S.Ulam proposed an elegant and simple model  aimed to show how the Maxwell distribution for the velocity statistics in a gas may be derived directly from the conservation laws encoded in the particle dynamics. He considered the evolution of an assembly of $`N`$ particles, whose initial velocity distribution is arbitrary. At each instant a random pair of particles is selected and their velocities are updated as if they had collided elastically. Such a process converges spontaneously to a macro-state whose velocity distribution is a Gaussian, with the variance given by the constant energy of the system.

Relaxing the energy conservation law, i.e. considering inelastic collisions between particles, leads to different statistical properties, which in fact are related to a well defined physical problem, i.e. the dynamics of Fluidized Granular Materials . These are assemblies of moving macroscopic particles, whose evolution is controlled by the inelasticity of their mutual collisions. They present unusual structural and statistical properties which stimulated many studies not only because of their industrial and technological relevance, but also due to the fundamental questions they raise.

Roughly speaking, the inelasticity causes two phenomena on the velocity statistics: a) in the absence of energy injection the fluid cools down since the variance of the velocity distribution, the so called granular temperature, decreases; b) the form of the distribution may even cease to be Gaussian and approach a stationary shape under a suitable rescaling.

A large amount of recent studies  indicates that the evolution of a freely evolving granular system consists of three stages: a homogeneous cooling during which both the velocity and the density are spatially uniform; a second stage where the velocity field forms vortices and a third stage where the density becomes highly inhomogeneous and clusters appear. Unfortunately, many of the existing approaches have intrinsic limitations: MD numerical simulations are time demanding and hardly access to the correlated regimes; kinetic theories, based on the molecular chaos assumption, seem to be limited to the quasi elastic situation; hydrodynamic equations suffer from their phenomenological derivation. For these reason minimalistic models, which either lend themselves to analytical solutions or greatly reduce the computational effort, can provide useful hints.

Hereafter we shall focus on the physics of rapid granular flows employing a one dimensional collision rule in order to render the analysis as simple as possible. This choice, in spite of its apparent simplicity shares many features with the physics of higher dimensions. The paper is organised as follows: first, we shall revisit the inelastic hard rod model moving on a ring by means of event driven dynamics. Next, we shall compare these results with those of the inelastic Ulam model (which corresponds to a gas of pseudo-Maxwell molecules ), for which we find the exact asymptotic velocity statistics, showing how they disagree even in the homogeneous regime. Thirdly, we shall extend Ulam’s model including spatial fluctuations, and show that many features of the real one-dimensional gas, both in the homogeneous and in the inhomogeneous regime, are recovered.

## II Inelastic hard rods on a ring

We define a one dimensional granular gas as a system of hard rods, confined to a ring and colliding inelastically. After a binary collision the scalar velocities of the particles change according to:

$`v_i^{}`$ $`=`$ $`v_i{\displaystyle \frac{1+r}{2}}(v_iv_j)`$ (1)
$`v_j^{}`$ $`=`$ $`v_j+{\displaystyle \frac{1+r}{2}}(v_iv_j)`$ (2)
where $`r`$ is a restitution coefficient, which takes on the value $`1`$ for perfectly elastic systems and $`0`$ for completely inelastic particles. Few years ago McNamara and Young , and Sela and Goldhirsh , simulated such a model and observed a universal algebraic decay of the kinetic energy, $`E(t)t^2`$ (Haff’s Law), together with an anomalous behavior for the global velocity distribution, even in the early homogeneous regime. Furthermore, the appearance of strong inhomogeneities is the precursor of a numerically catastrophic event, named inelastic collapse: i.e. particles perform an infinite number of collisions in finite time interval. A renewal of interest on one-dimensional granular flows, has been generated by the recent work of Ben-Naim et al. . They eliminated the inconvenience of the inelastic collapse, by means of a physically motivated ’regularization’: they assumed the restitution coefficient to be elastic for small momentum transfer. This allows the system to enter into a new dynamical regime, independent of the choice of the regularization. During such a regime, the kinetic energy decays as $`E(t)t^{2/3}`$, for any $`r<1`$. A direct inspection of the hydrodinamic profiles, shows that such a regime is highly inhomogeneous, with density clusters and shocks in the velocity field. They suggest that inelastic systems behave asymptotically as a sticky ($`r=0`$) gas , which is known to be described by the Burgers equation in the inviscid limit . This should justify the presence of shocks and clustering and, moreover, predicts that the tails of the velocity distribution decay as $`\mathrm{exp}(v^3)`$ . On the other hand, an accurate numerical test of this prediction  seems problematic, because the bulk of the distribution does not deviate appreciably from a Gaussian (see Fig. 2). Such a behavior reflects the fact that the asymptotics is dominated by the dynamics of clusters of particles, which move through the system and coalesce, like sticky objects. Although the sticky nature of the inelastic granular dynamics seems a very appealing description, it is not clear if it applies to higher dimensions .

To the best of our knowledge, exact analytic treatments of the hard rod model are limited to the homogeneous regime, by means of solutions of the Boltzmann equation. Caglioti et al. solved asymptotically the Boltzmann equation for a 1d inelastic hard rod gas in the limit $`r1^{}`$ and $`\rho \mathrm{}`$ (where $`\rho `$ is the density). They obtained the asymptotic velocity distribution, which, if rescaled to have a constant variance, is the superposition of two delta functions (see also ). For MD simulations, the observation of two distinct peaks in the rescaled velocity distribution was reported in  (see Fig. 1): it represents a precursor of that singular asymptotic distribution. However, the appearance of inhomogeneities and, eventually, the inelastic collapse prevent the approach to the limiting distribution.

## III Ulam’s model

Alternatively, one may attack the Boltzmann equation, by assuming a simpler form for the scattering cross section in the collision integral, i.e. taking the relative velocity of the colliding pair to be proportional to the average thermal velocity: $`|vv^{}|\sqrt{E}`$. This is the so called pseudo-Maxwell model . The energy factor $`\sqrt{E}`$ can be eliminated via a time reparametrization, and one obtains a simpler equation:

$$_\tau P(v,\tau )+P(v,\tau )=\beta 𝑑uP(u,\tau )P(\beta v+(1\beta )u,\tau )$$
(3)
where $`\beta =2/(1+r)`$ and the $`\tau `$ counts the number of collisions per particle. Eq. (3) is the master equation of the inelastic version of Ulam’s scalar model: at each step an arbitrary pair is selected and the scalar velocities are transformed according to the rule of Eq. (2). In a stimulating paper, Ben-Naim and Krapivsky  considered such scalar model, and obtained the evolution of the moments of the velocity distributions. Since at large times, $`v^n\mathrm{exp}(\tau a_n)`$, and the decay rates $`a_nna_2/2`$ (they depend non-linearly on $`n`$), they argued that such a multiscaling behavior prevents the existence of a rescaled asymptotic distribution $`f`$ such that $`P(v,\tau )f(v/v_0(\tau ))/v_0(\tau )`$, for large $`\tau `$, where $`v_0^2(\tau )=v^2P(v,\tau )𝑑v=E(\tau )`$. On the contrary, we believe that such “multiscaling” behavior only indicates the fact that the moments of the rescaled distribution $`x^nf(x)𝑑x=v^n/v_0^n`$ diverge asymptotically for $`n3`$, and does not rule out the possibility of the existence of an asymptotic distribution with power law tails. In fact, the Fourier transform of Eq. (3)

$$_\tau \widehat{P}(k,\tau )+\widehat{P}(k,\tau )=\widehat{P}[k/(1\beta ),\tau ]\widehat{P}[k/\beta ,\tau ]$$
(4)
possesses several self-similar solutions of the kind $`\widehat{P}(k,\tau )=\widehat{f}(kv_0(\tau ))`$, which correspond to the asymptotic rescaled distribution $`P(v,\tau )=f(v/v_0(\tau ))/v_0(\tau )`$. Many of them do not correspond to physically acceptable velocity distributions . The divergence of the higher moments implies a non analytic structure of $`\widehat{f}`$ in $`k=0`$, since $`v^n/v_0^n=(i)^n\frac{d^n}{dk^n}\widehat{f}(k)|_{k=0}`$, and represents a guide in the selection of the physical solution. As shown in Fig. 3, our data collapse on the function

$$f(v/v_0(\tau ))=\frac{2}{\pi \left[1+(v/v_0(\tau ))^2\right]^2}.$$
(5)
corresponding to the self-similar solution $`\widehat{f}(k)=(1+|k|)\mathrm{exp}(|k|)`$ . Notice that (5) is a solution of Eq.(4) for every $`r<1`$, i.e. the asymptotic velocity distribution does not depend on the value of $`r<1`$, as shown in Fig. 3 (left frame). However, such solution (5) does not resemble the distribution obtained in the simulation of the inelastic hard rod model, suggesting that the approach á la Ulam fails to describe even the homogeneous cooling phase (see Fig. 1). We believe that such failure is due to the absence of velocity correlations.

## IV One-dimensional lattice model

To reinstate spatial correlations we considered a 1d lattice model where the $`N`$ sites have a “velocity” $`v_i`$. At every step of the evolution a pair of neighboring sites is chosen randomly and undergoes an inelastic collision according to Eq.(2) if $`(v_iv_j)\times (ij)<0`$. The latter condition, which avoids collisions between particles “moving” away from each other, is to be be referred below as the kinematic constraint and plays a key role in the formation of structures during the inhomogeneous phase. A unit of time $`\tau `$ correspond to $`N`$ of such steps (i.e. to $`1`$ collision per particle on average). Note that, since there is no real particle movement, the number of collisions per particle, $`\tau `$, remains the only measure of time in our model. This represents a problem in the direct comparison with MD, because $`\tau (t)`$ is a well defined “universal” function for the homogeneous regime only: in the inhomogeneous regime $`\tau (t)`$ diverges due to the inelastic collapse, or depends on the “regularization” of the restitution coefficient. We observe that initially the energy $`E(\tau )`$ decreases exponentially, as in the Haff regime for MD simulation, but after the formation of long range velocity correlations, the decay slows down, and a power law decay of the form $`E\tau ^{1/2}`$, appears. As stated above, this law cannot be directly compared to the $`t^{2/3}`$ behavior, observed in MD, since, in the non-homogeneous regime the function $`\tau (t)`$ depends on the regularization employed to avoid the inelastic collapse. However the problem can be circumvented if one assumes the energy, instead of $`\tau `$, as the physical clock of both dynamics.

A surprising agreement between the MD and the lattice gas is found in the comparison of the velocity distribution taken at instants corresponding to the same energy. Both the early two humps structure and the late Gaussian shape seen in MD simulations agree with those of the lattice gas, as shown in Figs. 1 and 2. The Gaussianity of the asymptotic field distribution, and the $`\tau ^{1/2}`$ asymptotic decay suggest that, in the lattice model, the scalar velocity field $`v(i,\tau )`$ is governed by a discretized version of the diffusion equation of the form: $`_\tau v(i,\tau )=\nu \mathrm{\Delta }v(i,\tau )`$, where $`\mathrm{\Delta }`$ is a lattice laplacian. We checked this possibility studying the structure factor $`S(k,\tau )=\widehat{v}(k,\tau )\widehat{v}(k,\tau )`$ where $`\widehat{v}(k,\tau )`$ is the Fourier transform of the field $`v(i,\tau )`$. In fact, at late times, $`S(k,\tau )=f_1(k\tau ^{1/2})`$ scales similarly to a diffusive process as shown by the data collapse in Fig. 4 (right frame). However the form of the scaling function differs from the Gaussian shape predicted by the diffusion equation. In same figure (left) we report the structure factor of the MD simulation (fourier transform of $`C(r,t)=v(i,t)v(i+r,t)`$, where $`v(i,t)`$ is the velocity of the $`i`$-th particle, see below). In this case a good collapse, $`S(k,t)=f_2(kt^{2/3})`$, has been obtained . The two collapses can be directly compared, because, in both models, the structure factor is a function of $`k/E`$, consistently with the previous discussion on $`t`$ and $`\tau `$. The algebraic tails of the structure factors observed in both models carry important information about the nature of the growing velocity correlations. A comparison of velocity profiles at instants having the same energy in the two models is presented in Fig. 5. The top frame displays $`v(x,t)`$ of the MD, with the characteristic shock and preshock structures. In the middle frame, it is shown the profile $`v(i,t)v(x_i(t))`$ of the MD, where $`i`$ is the particle label. The difference between $`v(x,t)`$ and $`v(i,t)`$ in MD is in the sign of large gradients: $`v(x,t)`$ displays large negative gradients, while $`v(i,t)`$ has large positive gradients. Finally, the bottom frame displays the $`v(i,\tau )`$ profile of the lattice gas which fairly compares with the analogous MD quantity. The power law behavior $`S(k)k^2`$ observed in both models are due to the presence of short wavelength defects, viz. shocks, as predicted by Porod’s law . Note that the presence of shocks in the lattice model is due to the kinematic constraint, since they disappear, together with the Porod’s tails, relaxing the constraint.

## V Conclusions

In summary, we have investigated the velocity statistics of the homogeneous and inhomogeneous cooling regimes of the 1d granular gas. The complex observed behavior, featuring two distinct dynamical regimes, the onset of long-ranged velocity correlations and shocks in the velocity profiles can be reproduced even disregarding the real kinematics of the particles. In fact we introduced a lattice model of immobile particles which reproduces quantitavely the velocity statistics both at short and long times. The structure factors and the velocity profiles can also be qualitatively compared. The crucial ingredient of the model is the onset of velocity correlations in the system. Such correlations are relevant for the velocity statistics even when short-ranged. We corroborated such a conclusion through the study of the mean-field limit of the lattice model, which corresponds to a pseudo-Maxwell model of an inelastic gas. In this limit the velocity distributions change drammatically, displaying large tails, at odds with the real gas. We found the exact expression for this power law tailed distribution, which represent a self-similar solution of a non-trivial integro-differential Boltzmann like equation . Although the one-dimensional character of our results may seems restrictive, they can be interesting in several contests, as kinetic theory and microscopic models for turbulence (Burgers equation) or for surface growth (KPZ equation).

We wish to thank E. Caglioti and A. Vulpiani for many fruitful discussions.


text:
# The Impact of Massive Stars on the Formation of Young Stellar Clusters

## 1. Introduction

Stars are formed not in isolation, but in associations, groups and clusters. Associations span tens of parsecs and are defined by the molecular cloud complexes from which their constituent stars form; some examples are the T association in the Taurus molecular cloud, and the extensive association of OB stars interacting with the Orion molecular clouds. Within these molecular clouds, star formation is typically not widespread, but concentrated into distinct, parsec–sized groups and clusters with memberships ranging from ten stars to a thousand stars111Since systems of more than 100 young stars may persist as gravitationally bound open clusters, Adam & Myers (2001) define clusters as systems of more than 100 stars, and groups as systems containing less than 100 stars.. The stars in associations and clusters are not formed simultaneously in a single star forming event, but over a span of 1 Myr in clusters (Palla & Stahler 2000), and a span of several to 10 Myr in associations (Briceño et al. 2001). The immediate implication of both the non-isolated and non-coeval nature of star formation is that the formation of stars can be affected by winds, radiation, and outflows from previous generations of stars.

OB stars have the most dramatic influence on star formation. These stars both ionize and photodissociate molecular clouds with their intense UV radiation and sweep up atomic and molecular gas with winds, ionization fronts, and – at the end of their lives – supernovae. In particular, stars with spectral types of B2 or earlier produce enough UV radiation to effectively ionize cluster–forming molecular cores. Since OB stars are relatively rare, it is important to assess their importance by examining unbiased near–IR surveys of molecular clouds. In an analysis of the 2MASS 2nd incremental release data covering the Perseus, MonR2, Orion A and Orion B clouds, Carpenter (2000) found a total of 1200 isolated stars, 200 stars in eight groups, and 2900 stars in six clusters. Hence, the majority of the stars in these clouds formed in clusters containing 100 or more stars. In the sample, four of the clusters, GGD12-15, Mon R2, NGC 2024 & the Orion Nebula Cluster, contain at least one star with a spectral type of B2 or earlier. One of the groups, NGC 2068222Our classification of NGC 2068 as a group is due to the partial coverage of this region in the 2nd incremental release of the 2MASS data. While Carpenter finds only 45 stars in the incomplete 2MASS data, Lada (1989) reports 192 stars in the region, indicating that NGC 2068 should be classified as a cluster., contains a B1.5 star. In total, 57% of the stars detected in the studied regions are in close proximity to ionizing OB stars. This number is a lower limit since it only counts OB stars which form within clusters. As we will demonstrate, external OB stars may also influence the formation of clusters.

In this paper, we discuss how massive stars influence the process of cluster formation on sizescales ranging from 0.1 pc to 300 pc, with a particular emphasis on our own study of the NGC 281 region. We start by examining the disruption of cluster–forming molecular cores by nascent OB stars within the cluster. As we demonstrated above, the majority of observed young clusters appear to contain OB stars. We then study the disruption of cluster forming cores by neighboring OB stars, which are situated a few parsecs away from the star forming region. This case occurs during the formation of an OB association, when OB stars which have dispersed their own parental cores begin to ionize neighboring star forming cores. Finally, we discuss a possible case of supernova driven star formation, where the influence of OB stars extends beyond the boundaries of molecular cloud complexes to distances of a few hundred parsecs.

## 2. Core Destruction by Internal OB Stars ($`<1`$ pc)

The birth sites of massive stars are 1000 M dense molecular cores with densities of $`10^5`$ cm-3 and sizes of 1 pc (Plume et al. 1997). Within these cores, entire clusters of lower–mass stars form. OB stars still embedded within their parental cores are thought to form ultracompact HII regions (UCHIIs), HII regions with sizes of $`<0.1`$ pc and electron densities $`>10^4`$ cm-3. This picture is supported by the growing number of observations showing UCHIIs surrounded by embedded clusters of lower–mass stars (Megeath et al. 1996; Tapia, Persi & Roth 1996; Pratap, Megeath & Bergin 1999). After OB stars dissipate their parental dense cores, the HII regions expand and are no longer classified as UCHIIs (Tieftrunk et al. 1997). Hence, the lifetime of UCHIIs is the time required for young OB stars, once they have formed, to destroy their parental cores and terminate the process of cluster formation.

The best current estimates of the lifetimes of UCHIIs are based on empirical arguments. Wood & Churchwell (1989) found that UCHIIs can be identified on the basis of their IRAS colors. By comparing the number of IRAS identified UCHIIs in the Galaxy with the observed number of OB stars in the solar neighborhood, they estimated a lifetime of 0.5 Myr. By eliminating lower–mass stars contaminating the IRAS sample, as well as taking into account variations in the density of UCHIIs with galactic radius, this lifetime has been revised by Comeron & Torra (1996) and Casassus et al. (2000) to $`0.050.1`$ Myr.

The estimated UCHII lifetimes indicate that once an OB star is formed in a cluster forming core, the core destruction occurs on a timescale of 0.1 Myr. Since the timescale for cluster formation is ten times longer, 1 Myr (Palla & Stahler 2000), only $`1/10`$ of the forming clusters should contain ultracompact HII regions. Interestingly, in the combined surveys of W3 (Tieftrunk et al. 1998), Perseus, Orion A, Mon R2 (Carpenter 2000) & Orion B (Lada 1992), $`4/10`$ of the clusters contain UCHIIs: Mon R2, W3 OH, W3 Main, and GGD12-15. Although the ratio is four times that expected, the current sample of molecular clouds is small. A larger sample is needed to firmly establish the fraction of star–forming clusters with UCHIIs.

## 3. Core Destruction by External OB Stars (1-10 pc)

Optical images of HII regions commonly show bright rimmed and cometary globules: dense molecular gas that survives in the HII region once the surrounding, lower density gas is swept away. These globules are directly exposed to the UV radiation from OB stars within the HII region. In a process called photoevaporation, the UV radiation produces a flow of ionized gas off the surfaces of the clump; the optical line emission producing the bright rims arises in this ionized flow. Despite the continual erosion of the molecular gas, a theoretical analysis by Bertoldi & McKee (1990) has demonstrated that these globules can survive for millions of years in the HII region, depending on their mass and distance from the OB stars. This relatively long timescale for the core destruction allows for the formation of groups and clusters in the exposed cores. Several examples of young stellar groups forming in bright rimmed globules have been reported by Sugitani, Tamura & Ogura (1995).

The formation of groups and clusters in externally ionized cores may not result solely from the long destruction timescales, but also from the triggering of star formation by the ionization. A likely example of triggered cluster formation is the NGC 281 infrared cluster, the richest cluster yet observed towards a bright rimmed globule. It is the most active region of ongoing star formation in NGC 281, containing at least one massive star, as evidenced by a $`10^4`$ L IRAS source (Henning 1994), and a cluster of 240 lower–mass stars, as detected in near–IR images (Megeath & Wilson 1997). The infrared cluster, located in the southwestern quadrant of the NGC 281 nebula, is not associated with the optical bright rimmed globules apparent in the nebula (Fig. 1).

The evidence for triggered star formation comes from combined centimeter, millimeter, and near–IR data. Since optical observations of this region are complicated by a spatially varying extinction, we obtained C18O maps, tracing the dense cores, and a VLA 20 cm map, tracing the ionized gas, to disentangle the structure of the H2/HII interface with angular resolutions $`<20^{\prime \prime }`$. The maps show that the molecular gas is concentrated into three main cores and several smaller clumps (Fig. 2). To the northern edges of each core, 20 cm emission tracing the ionized gas is detected, with particularly strong 20 cm emission adjacent to the two northernmost cores. These cores are radio bright rimmed globules. In support of the radio data, recent Hubble space telescopes images taken with the WFPC2 camera show the optical emission from H$`\alpha `$ and \[SII\] lines extending to the cores. These data demonstrate that the cores are exposed to the UV radiation from the ionizing O star, and are photoevaporating. Megeath & Wilson (1997) estimate that at the current rate of photoevaporation, the observed cores will be fully evaporated in 2.5 Myr.

A comparison of the pressure from the ionized, photoevaporating gas at the clump surfaces to the turbulent pressure of the molecular gas in the clumps, indicates that the pressure of the ionized gas is sufficient to drive shock waves into the molecular cores. The shocks would have a velocity of 1.5 km s-1 and a crossing time of 0.7 Myr. The C18O data show a ridge of blueshifted gas along the northern edge of the cores: this may be the detection of gas swept up by the shock (see Megeath & Wilson 1997).

Near–IR images of the two northermost cores show a steep gradient in the surface density of stars with most of the embedded stars concentrated near the photoevaporating surfaces of the cores (Fig. 3). The lack of a corresponding gradient in the extinction through the cores, as demonstrated by the C18O maps, implies that the observed gradient in stellar density is due to an asymmetric distribution of stars. Megeath & Wilson (1997) have proposed that the asymmetry is the result of triggered star formation. In this picture, the shock fronts have traversed roughly halfway through the cores. As the shocks compress the gas, they precipitate the formation of stars which then fill the post-shock gas. Near–IR imaging of groups in optical bright rimmed globules have detected similar asymmetries in the distribution of stars (Sugitani, Tamura & Ogura 1995). These data suggest that shock–triggered formation may be a common process in cores exposed to the UV radiation from neighboring OB stars.

## 4. Supernovae Driven Formation of Clusters ($`>100`$ pc)

The NGC 281 nebula is part of a complex of molecular clouds in the Perseus arm of the galaxy (Fig. 4). The measured distance to the nebula (2.9 kpc; Guetter & Turner 1997) implies that the complex is located 300 pc below the galactic plane; in contrast, the rest of the Perseus arm clouds are found within 150 pc of the plane. Another unusual property of the NGC 281 complex is the high dispersion in cloud velocities: the longitude-velocity diagram in Fig. 5 shows a dispersion of 45 km s-1. The high velocity dispersion suggests that the complex may only be a chance alignment of clouds in the sky. However, the lack of other Perseus arm clouds at these galactic latitudes, the linking of the CO clouds by the HI emission (Fig. 4 and 5), and the elliptical distribution of CO and HI emission in the longitude-velocity diagram, all support the notion that these clouds form a bona fide complex of molecular and atomic clouds.

The elliptical pattern evident in the longitude-velocity diagram is indicative of an expanding ring. The diameter of the ring is 270 pc, the expansion velocity in the line of sight is 22 km s-1, and the resulting dynamical time is 6 Myr. The mass in molecular gas, using a standard conversion from the velocity–integrated CO temperatures to $`H_2`$ column density, is $`10^5`$ M. The HI mass integrated over latitudes of -6o to -11o and velocities of -15 km s-1 to -60 km s-1 is $`3.5\times 10^5`$ M. The resulting kinetic energy is $`4.5\times 10^{51}`$ ergs, comparable to the energy of a supernova.

We are currently exploring the possibility that the ring has formed in gas lifted out of the plane by the expansion of a superbubble. The HI emission in the NGC 281 complex forms part of a loop of HI extending out of the plane; this loop may trace the shell of a superbubble driven by supernovae and winds from OB associations near the plane of the Galaxy (Fig. 4). Models of superbubbles expanding into a stratified galactic “atmosphere” predict that at a distance of a few hundred parsecs above the plane, the upward expansion of the bubble accelerates resulting in a blowout (MacLow & McCray 1988). Numerical simulations of blowouts predict that the fragmenting shell can have velocities as large as 50 km s-1 tangential to the galactic plane (MacLow, McCray & Norman 1989). It is then plausible that the observed expanding ring has been formed in a blowout, with the clouds located at the distance from the plane where the blowout started. This interpretation suggests that the formation of the entire NGC 281 complex has been triggered by previous generations of OB stars. This mode would operate on a totally different timescale (6 Myr vs 1 Myr) and sizescale (300 pc vs 4 pc) than the triggered star formation discussed in the previous section.

## References

Bertoldi, F. & McKee, C. F. 1990 ApJ 354, 529.

Briceño, C., Vivas, A. K., Calvet, N. et al. 2001 Sci, 291, 93.

Carpenter, J. M. 2000 AJ 120, 3139.

Casassus, S., Bronfman, L, & May, J. et al. 2000 A&A 358, 514.

Comeron, F. & Torra, J. 1996, A&A 314, 776.

Dame, T. M., Hartmann, D., & Thaddeus, P. 2001 ApJ 547, 792.

Guetter, H. H. & Turner, D. G. 1997 AJ 113, 2116.

Hartmann, D. & Burton, W. B. 1997, Atlas of Galactic Neutral Hydrogen (Cambridge: Cambridge Univ. Press).

Henning, Th., Martin, K., Reimann, H.-G., et al. 1994 A&A 288, 282.

Lada, E. A. 1992 ApJ 393, L25.

MacLow, M.-M., McCray, R. 1988 ApJ 324, 776.

MacLow, M.-M., McCray, R. & Norman, M. L. 1989 ApJ 337, 141.

Megeath, S. T., Herter, T., Beichman, C., et al. 1996 A&A 307, 775.

Megeath & Wilson 1997 AJ 114, 1106.

Palla, F., Stahler, S. W. 2000 ApJ 540, 255.

Plume, R., Jaffe, D. T., Evans, N. J., II. et al. 1997, ApJ 476, 730.

Pratap, P., Megeath, S. T., & Bergin, E. A. 1999 ApJ 517, 799.

Sugitani, K., Tamura, M., & Ogura, K. 1995 ApJ 455, L31.

Tapia M., Persi, P., & Roth, M. 1996 A&A 316, 102.

Tieftrunk, A. R., Megeath, S. T., Wilson, T. L., et al. 1998 A&A 336, 991.

Tieftrunk, A. R., Gaume, R. A., Claussen, M. J., et al. 1997 A&A 318, 931.

Wood, D. O. S. & Churchwell, E. 1989 ApJ 340, 265.


text:
# MULTIPLE PRODUCTION, TRANSPORT IN ATMOSPHERE AND DETECTION OF HIGH ENERGY COSMIC RAYS

## 1 Monte Carlo Collision generators for Cosmic Rays and extrapolation at ultra high energy

In order to simulate in a reasonable time the 4-dimensional development of Extensive Air Showers at very high energy, it is important to elaborate collision generators reproducing rapidly the detailed features of multiple production observed in accelerators and colliders. Among the variety of models implemented in CORSIKA, there are microscopic and phenomenological models. The first ones include all the steps of the parton momenta generation from the parton distribution functions related to valence quarks and diquarks, sea quarks and gluons, convoluted with the fragmentation of the respective strings into secondary hadrons. The phenomenological models provide fast and direct hadron sampling, but take into account some global characteristics of the Gribov-Regge theories with parameters adjusted to the experimental data. One common feature is the separation between the non single diffractive component and the diffractive component (single and double) with the inclusion of soft and hard mechanisms.
As an example (fig.1), the reproduction of the inclusive data with the model HDPM2 (2nd version of the hybrid dual parton model) with parameters tuned to fit the experimental data of FermiLab at $`\sqrt{s}=630GeV`$ gives a total average inelasticity for p-p collisions of $`0.7`$ instead of $`0.5`$ when adjusted to the previous measurements of UA5. The forward trajectory on Fig. 1 is not constrained above 5.5 units of pseudo rapidity and several models taken in option for CORSIKA can have different trajectories in the forward region with a good agreement to the observed data, but rather different inelasticities. A cosmic ray cascade is initiated by the unique interaction of the primary particle and the fluctuations of the elementary act have to be reproduced carefully taking into account the semi-inclusive data as well as the correlation between charged and neutral secondaries. The extrapolation at ultra high energy can be carried with the model HDPM2 taking into account recent features of collider physics such as $`p_t`$ versus central rapidity density (UA1-MIMI exp.) and recent results of FermiLab  for pseudo-rapidity up to 5.5. The pseudo-rapidity distributions obtained  with HDPM2 for 2000 collisions NSD are shown on Fig. 1. The distributions on Fig. 1 (histograms for HDPM2, full line for QSJET model) give for HDPM2 different extrapolations at $`10^{10}`$ GeV following the PDF assumed ($`B0`$ and $`B_{}`$).
The correlation between the central rapidity density and the average transverse momentum $`<p_t>`$ plays also one role in the muon and hadron radial distributions: the dependence of the transverse momentum on energy turns to a permanent increase on Fig. 1.

The p-p collision is transformed to p-A collision following Glauber’s considerations and Nuclei-Air collisions are treated by the abrasion evaporation procedure.

## 2 Simulation of giant EAS

The electromagnetic cascade is treated by EGS4 involving cross sections, taking into acount the LPM effect, for bremsstrahlung and pair production. The nuclear photo production has been implemented in the simulation and its contribution is compared versus the primary photon energy to the previous processes on the Fig. 2 .

The fluctuations of multiplicity are governed by the negative binomial distribution and at $`10^{20}eV`$, the average charged multiplicity is close to 300 for HDPM (1000 for QGSJET); this corresponds to energy densities of $`30GeV/fm^3`$ in average… and we have not considered here the important modifications from a possible phase transition to QGP.
The longitudinal developments obtained for proton, photon and iron induced extensive air showers are compared on Fig. 2. We notice that the showers (here averaged on 10 individual cascades simulated with CORSIKA) have their maximum near 1000–1500 m altitude for proton and iron. The $`\gamma `$ shower initiated at $`50g/cm^2`$ exhibits a different behaviour due to the LPM effect; a rule of thumb for hadronic showers is to divide by 2 the total energy in $`GeV`$ to get $`N_{max}`$. This maximum depth corresponds to a minimum of fluctuations, suggested by analytic cascade theory and confirmed by the Monte Carlo simulation, indicating the well adapted localisation of AGASA and AUGER experiments. The depth of the maximum depends on the logarithm of the primary energy, when $`N_{max}`$ remains proportional to energy with a similar factor for a large variety of interaction models.
A more delicate problem is the estimation of $`N_{max}`$; it can be read directly from the cascade curve derived from fluorescence measurements, as performed by the Fly’s Eye or Hires, but the approach with a small number of detectors hit at large distances (1 to 1.5 km) turns to a hopeless topological problem.

## 3 Topological aspects of radial distribution

The profile of the lateral distribution assumed and the method of core localization are especially important. We propose instead of NKG and other Euler Beta functions, the employment of the gaussian hypergeometric formalism giving also normalization and better skewness, under the form:

$`f(x)`$ $`=g(s)x^{sa}(x+1)^{sb}(1+dx)^c`$ (1)
which has the advantage (for values of parameters respecting the conditions of convergence $`sa+2>0`$ and $`c2s+b2>0`$) to be exactly normalized in terms of Gaussian Hypergeometric function $`F_{HG}`$ $`=`$ $`F(c,sa+2,c+bs;1d)`$ by:

$$g(s)=\frac{\mathrm{\Gamma }(c+bs)}{2\pi \mathrm{\Gamma }(sa+2)\mathrm{\Gamma }(c2s+b+a2)F_{HG}}$$
(2)
The empirical distributions, such as AGASA function , as underlined by Vishwanath  enters in the category of Hypergeometric Gaussian functions , under the general form of structure function

$$f(x)=C_ex^\alpha (1+x)^{(\eta \alpha )}(1+dx)^\beta $$
(3)
with the conditions $`2\alpha >0`$ and $`\beta +\eta 2>0`$. The value used in AGASA function for the coefficient $`C_e`$ is just an approximation; the exact value is

$$C_e=\frac{\mathrm{\Gamma }(\beta +\eta \alpha )}{2\pi \mathrm{\Gamma }(2\alpha )\mathrm{\Gamma }(\beta +\eta 2)}\frac{1}{F_{HG}}$$

where $`F_{HG}=F_{HG}(\beta ,2\alpha ,\beta +\eta \alpha ;1d)`$. The Hypergeometric Gaussian function can be easily calculated from the hypergeometric series:

$`F_{HG}(a,b,c;z)={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(a)_n(b)_n}{(c)_nn!}}z^n,\text{ }c0,1,2,\mathrm{},`$
$`(a)_n=\mathrm{\Gamma }(a+n)/\mathrm{\Gamma }(a),\text{ }(a)_0=1`$
This equation is equivalent to our version $`(3)`$ containing the age parameter s with the relations between respective coefficients:   $`x=\frac{r}{r_\mathrm{M}}`$,   $`d=\frac{r_\mathrm{M}}{r_0}`$,   $`s=1.03`$,   $`\alpha =as`$,   $`\eta =bs+\alpha `$ (the value of s is taken from the longitudinal development simulated).

We have adjusted with MINUIT the parameters of our hypergeometric function (Table 1) to the average lateral distributions (set JNC01 for charged, JNC02 for electrons) of groups of 10 showers at $`10^{20}`$ eV simulated with CORSIKA (QGSJET model), as shown on Fig. 3.

In each case, the adjustment has been performed with 50 points from the simulation distributed from 0.1 m up to 10 km from axis position for charged particles (muons and electrons). The advantages of JNC01 formula can be seen on Fig. 3 and on Table 2.

The major part of the particles is contained inside 200 m from the axis and only the skewness of the hypergeometric function allows a reliable relation between size and density at 600 m. This function has been applied to localizations of the showers contained in the catalogues of Volcano Ranch and Yakutsk. The core position has been obtained by minimization with Minuit program between different formulas available for lateral densities written versus the coordinates $`X`$, $`Y`$ as

$$\varrho (r)=\varrho (\sqrt{(XX_c)^2+(YY_c)^2})$$
(4)
where the core coordinates $`X_c`$ and $`Y_c`$ are taken as two additive parameters in the minimization. The adjustments are generally improved when compared to the original treatments, turning to lower sizes (in the case of Yakutsk formula) and better approximation of the density at 600 m. The situation of the most energetic event of AGASA  is given in Table 3.

## 4 Fluorescence, Cerenkov and Radio emissions

From the longitudinal development, we have developed a method of fast simulation based on the structural stability of the subshowers (in the sense of catastroph theory). This is an efficient alternative to the thinning technique. The lateral extension of the total amount of light received at a fixed distance from the axis is compared for both Cerenkov and fluorescence component and shown on Fig. 4.
The opportunity of CORSIKA which separates positive and negative particles taught us that a regular negative excess close to $`25\%`$ is present in the e.m. component of EAS, allowing the calculation of radio emission, following Askarian’s effect. In Yakutsk experiment, it is possible to estimate the primary energy from the density at $`600m`$ as well as from the Cerenkov component; in AUGER array, the direct calibration of showers detected simultaneously by the ground array and the fluorescence detector looks promising.

## 5 Diquark breaking mechanism and GZK cut off

The diquark breaking mechanism disturbs strongly the leading particle effect present in the different models used in cosmic rays. In the classical form of the dual parton model, the 3 valence quarks of the proton projectile are separated in a fast diquark and another valence quark slowed down. The diquark is recombined with one quark of the sea to produce, the most commonly, an outgoing leader baryon propagating the energy deeper in the cascade. The 3 valence quarks separated will be recombined in various meson structures in pairs $`|u\overline{d}>`$, $`|d\overline{u}>`$,… or neutral mesons as $`1/\sqrt{2}(d\overline{d}u\overline{u})`$. The configuration with the simultaneous final state for the valence quarks of $`3\pi ^0`$’s could be especially interesting with a probability of emergence that we can evaluate from the quark content and the quark additive model as $`1/27`$. Such configuration (with intermediate final states of higher probabilities, one pair of charged pions and one neutral, one pair of neutral and one charged…) will transfer a large part of energy to the electromagnetic component and this energy will be definitely missing for both hadronic and penetrating components. Remembering that for the same primary energy, the cascade theory shows that one primary photon produces at maximum, approximately, two times more electrons, we can expect a large electron excess for some cascades initiated with diquark breaking. The longitudinal development calculated for protons of the same energy of $`310^8`$ GeV is compared to a classical development, here the model HDPM2 with D–option, on Fig. 4. The electron size at maximum is doubled in the assumption of diquark breaking and relatively rare recombination simultaneously in $`3`$ neutral pions . We have simulated 104 showers, between $`10EeV`$ and $`50EeV`$, where the generation, following the primary energy differential spectrum, is cut by brute force. The showers simulated with axis distributed randomly are treated following the method of AGASA and the primary spectrum reconstructed is superimposed on the original one. It appears an extra component above $`50EeV`$, due to the showers including the diquark mechanism, considered previously, and such artefact could appear as a GZK violation. With AUGER array, such showers could be recognized very easily by an energy 3 or 4 times larger from the fluorescence detector than estimated by the water tank detectors. We observed also that such artificial component could appear from spacial topological situations at $`50EeV`$; for instance, with an axis falling near one detector, all the others are hit by densities near the low density threshold and the answer very unstable of those surrounding detectors can turn to an important overestimation of the primary energy.

## 6 Conclusions

The hypergeometric approach gives a better accuracy in the interpolation of densities at 600–1000 m, a more reliable estimation of the shower size (when the axis is in the array), a better axis localization and finally a more correct constraint of the primary energy. The position of the axis is crucial as the density near $`600m`$ varies as $`r^4`$ and an uncertainty of $`50m`$ on axis turns to an error of $`40\%`$ at least on primary energy. Topological difficulties and also non standard aspect of multiple production needs probably more attention to be able to rule out definitely the cut off of GZK.

## Acknowledgments

This work has been partly supported (JNC, FC and CLG) by INTAS contract 1339.


text:
# The New Emerging Model for the Structure of Cooling Cores in Clusters of Galaxies

## 1. Introduction

X-ray imaging observations have shown that the X-ray emitting, hot gas in a large fraction of all galaxy clusters reaches high enough densities in the cluster centers that the cooling time of the gas falls below the Hubble time, and gas may cool and condense in the absence of a suitable fine-tuned heating source (e.g. Silk 1976, Fabian & Nulsen 1977). From the detailed analysis of surface brightness profiles of X-ray images of clusters obtained with the Einstein, EXOSAT, and ROSAT observatories, the detailed, self-consistent scenario of inhomogeneous, comoving cooling flows emerged (e.g. Fabian et al. 1984, Nulsen 1986, Thomas, Fabian, & Nulsen 1987, Fabian 1994). The main assumptions on which the cooling flow model is based and some important implications are: (i) Each radial zone in the cooling flow region comprises different plasma phases covering a wide range of temperatures. The consequence of this temperature distribution is that gas will cool to low temperature and condense over a wide range of radii. (ii) The gas features an inflow in which all phases with different temperature move with the same flow speed. (iii) There is no energy exchange between the different phases, between material at different radii, and no heating.

Now the first analysis of high resolution X-ray spectra and imaging spectroscopy obtained with XMM-Newton has shown to our surprise that the spectra show no signatures of cooler phases of the cooling flow gas below an intermediate temperature which constitutes a problem for the interpretation of the results in the conventional cooling flow picture (e.g. Peterson et al. 2001, Tamura et al. 2001). Another result is that the spectroscopic data are better explained with local isothermality in the cooling flow region (e.g. Böhringer et al. 2001a, Matsushita et al. 2001, Molendi & Pizzolato 2001) also in conflict with the inhomogeneous cooling flow model. Here, we discuss these new spectroscopic results and their implications and point out the way to a new possible model for this phenomenon. The results are mostly based on the detailed observations of the M87 X-ray halo. A detailed description of this study is provided by Böhringer et al. (2001b).

## 2. Spectroscopic Diagnostics of Cluster Cooling Cores

XMM Reflection Grating Spectrometer (RGS) observations of several cooling core regions show signatures of different temperature phases ranging approximately from the hot virial temperature of the cluster to a lower limiting temperature, $`T_{low}`$. Clearly observable spectroscopic features of even lower temperature gas expected for a cooling flow model are not observed. A1835 with a bulk temperature of about 8.3 keV has $`T_{low}`$ around 2.7 keV (Peterson et al. 2001) and similar results have been derived for A1795 (Tamura et al. 2001). These results are very well confirmed by XMM observations with the energy sensitive imaging devices, EPN and EMOS, providing spectral information across the entire cooling core region, yielding the result that (for M87, A1795, and A1835) single temperature models provide a better representation of the data than cooling flow models (Böhringer et al. 2001a, Molendi & Pizzolato 2001) also implying the lack of low temperature components. The very detailed analysis of M87 by Matsushita et al. (2001 and the contribution to this workshop) has shown that the temperature structure is well described locally by a single temperature over most of the cooling core region, except for the regions of the radio lobes and the very center ($`r1`$ arcmin, $`5`$ kpc).

Among the spectroscopic signatures which are sensitive to the plasma temperature in the relevant temperature range, the complex of iron L-shell lines is most important. Fig. 1 shows simulated X-ray spectra as predicted for the XMM EPN instrument in the spectral region around the Fe L-shell lines for a single-temperature plasma at various temperatures from 0.4 to 2.0 keV and 0.7 solar metallicity. There is a very obvious shift in the location of the peak making this feature an excellent thermometer.

For a cooling flow with a broad range of temperatures one expects a composite of several of the relatively narrow line blend features, resulting in a quite broad peak. Fig. 2a shows for example the deprojected spectrum of the M87 halo plasma for the radial range 1 - 2 arcmin (outside the inner radio lobes) and a fit of a cooling flow model with a mass deposition rate slightly less than 1 M yr-1 as expected for this radial range from the analysis of the surface brightness profile (e.g. Stewart et al. 1984, Matsushita et al. 2001). It is evident that the peak in the cooling flow model is much broader than the observed spectral feature. For comparison Fig. 2b shows the same spectrum fitted by a cooling flow model where a temperature of 2 keV was chosen for the maximum temperature and a suitable lower temperature cut-off (1.44 keV) was determined by the fit. The very narrow temperature interval (almost isothermality) is well consistent with the narrow peak. A similar result is obtained for other clusters, e.g. A1795 as shown in Fig. 2d.

Since this diagnostics of the temperature structure is essentially based on the observation of metal lines, an inhomogeneous distribution of the metal abundances in the cluster ICM and a resulting suppression of line emission at low temperatures was suggested as a possible way to reconcile the above findings with the standard cooling flow model by Fabian et al. (2001a and contribution in these proceedings). As shown by Böhringer et al. (2001b) such a scenario will still result in a relatively broad Fe L-line feature and does not solve the problem in this case of M87.

## 3. Internal absorption

Another possible attempt to obtain consistency is to allow the absorption parameter in the fit to adjust freely. This is demonstrated in Fig. 2c with the same observed spectrum where the best fitting absorption column density is selected in such a way by the fit that the absorption edge limits the extent of the Fe-L line feature towards lower energies. This is actually the general finding with ASCA observations which has shown two possible options for the interpretation of the spectra of cluster core regions: (1) an interpretation of the results in form of an inhomogeneous cooling flow model which than necessarily includes an internal absorption component (e.g. Allen 2000, Allen et al. 2001), or (2) an explanation of the spectra in terms of a two-temperature component model (e.g. Ikebe et al. 1997, 1999, Makishima 2001) where the hot component is roughly equivalent to the hot bulk temperature of the clusters and the cool component corresponds approximately to $`T_{low}`$. Thus for the cooling flow interpretation to work and to produce a sharp Fe-L line feature as observed, the absorption edge has to appear at the right energy and therefore values for the absorption column of typically around $`310^{21}`$ cm-2 are needed (e.g. Allen 2000 and Allen et al. 2001 who find values in the range $`1.5510^{21}`$ cm-2).

It is therefore important to perform an independent test on the presence of absorbing material in the cluster cores. Thanks to CHANDRA and XMM-Newton we can now use central cluster AGN as independent light sources for probing. Using the nucleus and jet of M87 (with XMM, see Fig. 3) and the nucleus of NGC1275 (with CHANDRA) we find no signature of internal absorption. Thus at least for these two cases it is difficult to argue for internal absorption to obtain consistency of the observations with the cooling flow model.

## 4. Heating Model

In view of these difficulties of interpreting the observations with the standard cooling flow model, we may consider the possibility that the cooling and mass deposition rates are much smaller than previously thought, that is reduced by at least one order of magnitude. To decrease the mass condensation under energy conservation some form of heating is clearly necessary. Three forms of heat input into the cooling flow region have been discussed: (i) heating by the energy output of the central AGN (e.g. Pedlar et al. 1990, Binney & Tabor 1993, McNamara et al. 2000,

(ii) heating by heat conduction from the hotter gas outside the cooling flow (e.g. Tucker & Rosner 1983, Bertschinger & Meiksin 1986), and (iii) heating by magnetic fields, basically through some form of reconnection (e.g. Soker & Sarazin 1990, Makishima et al. 2001). The latter two processes depend on poorly known plasma physical conditions and are thus more speculative. The energy output of the central AGN, however, can be determined as shown below.

A heating scenario can only successfully explain the observations if among others the two most important requirements are met: (i) The energy input has to provide sufficient heating to balance the cooling flow losses, that is about $`10^{60}`$ to $`10^{61}`$ erg in 10 Gyr or on average about $`310^{43}310^{44}`$ erg s-1, and (ii) The energy input has to be fine-tuned. Too much heating would result in an outflow from the central region and the central regions would be less dense than observed. Too little heat will not reduce the cooling flow by a large factor. Therefore the heating process has to be self-regulated: mass deposition triggers the heating process and the heating process reduces the mass deposition.

Further constraints are discussed by Böhringer et al. (2001b). The total energy input into the ICM by the relativistic jets of the central AGN can be estimated by the interaction effect of the jets with the ICM by means of the scenario described in Churazov et al. (2000). It relies on a comparison of the inflation and buoyant rise time of the bubbles of relativistic plasma which are observed e.g. in the case of NGC 1275 (Böhringer et al. 1993, Fabian et al. 2001b) and requires as observational input parameters the bubble size, $`r_B`$, the ambient pressure, $`P_{th}`$, and the Keplerian velocity at the bubble radius in the cluster, $`v_K`$. The parameters and the estimated total energy output is given in Table 1 for three examples, M87, Perseus, and Hydra A. These values for the energy input have to be compared with the energy loss in the cooling flow, which is of the order of $`10^{43}`$ erg s-1 for M87 and about $`10^{44}`$ erg s-1 for Perseus. Thus in these cases the energy input is larger than the radiation losses in the cooling flow for at least about the last $`10^8`$ yr. We have, however, evidence that this energy input continued for a longer time with evidence given by the outer radio halo around M87 with an outer radius of 35 - 40 kpc (e.g. Kassim et al. 1993, Rottmann et al. 1996). Owen et al. (2000) give a detailed physical account of the halo and model the energy input into it. They estimate the total current energy content in the halo in form of relativistic plasma to $`310^{59}`$ erg and the power input for a lifetime of about $`10^8`$ years, which is also close to the lifetime of the synchrotron emitting electrons, to the order of $`10^{44}`$ erg s-1, consistent with our estimate. The very characteristic sharp outer boundary of the outer radio halo of M87, noted by Owen et al. (2000), has the important implications, that this could not have been produced by magnetic field advection in a cooling flow.

Thus, we find a radio structure providing evidence for a power input from the central AGN into the halo region of the order of about ten times the radiative energy loss rate over at least about $`10^8`$ years (for this representative example of M87). The energy input could therefore balance the heating for at least about $`10^9`$ years. The observation of active AGN in the centers of cooling flows is a very common phenomenon. E.g. Ball et al. (1993) find in a systematic VLA study of the radio properties of cD galaxies in cluster centers, that 71% of the cooling flow clusters have radio loud cDs compared to 23% of the non-cooling flow cluster cDs. Therefore we can safely assume that the current episode of activity was not the only one in the life of M87 and its cooling flow.

The mechanism for a fine-tuned heating of the cooling flow region should most probably be searched for in a feeding mechanism of the AGN by the cooling flow gas. The most simple physical situation would be given if simple Bondi type of accretion from the inner cooling core region would roughly provide the order of magnitude of the power output that is observed and required. Using the classical formula for spherical accretion from a hot gas by Bondi (1952) we can obtain a very rough estimate for this number. For the proton density near the M87 nucleus ($`r15`$ arcsec) of about 0.1 cm-3, a temperature of about $`10^7`$ K (e.g. Matsushita et al. 2001), and a black hole mass of $`310^9`$ M (e.g. Ford et al. 1994) we find a mass accretion rate of about 0.01 M yr-1 and an energy output of about $`710^{43}`$ erg s-1, where we have assumed the canonical value of 0.1 for the ratio of the rest mass accretion rate to the energy output. The corresponding accretion radius is about 50 pc ($`0.6`$ arcsec). This accretion rate is more than a factor of 1000 below the Eddington value and thus no reduction effects of the spherical accretion rate by radiation pressure has to be expected. Small changes in the temperature and density structure in the inner cooling core region will directly have an effect on the accretion rate. Therefore we have all the best prospects for building a successful self-regulated AGN-feeding and cooling flow-heating model.

## 5. Conclusions

Several observational constraints have let us to the conclusion that the mass deposition rates in galaxy cluster cooling cores are not as high as previously predicted. The new X-ray spectroscopic observations with a lack of spectral signatures for the coolest gas phases expected for cooling flows and the lower mass deposition rates indicated at other wavelength bands than X-rays are more consistent with mass deposition rates reduced by one or two orders of magnitude below the previously derived values. This can, however, only be achieved if the gas in the cooling flow region is heated. The most promising heating model is a self-regulated heating model powered by the large energy output of the central AGN in most cooling flows.

Most of the guidance and the support of the heating model proposed here (based on concepts developed in Churazov et al. 2000, 2001) is taken from the detailed observations of a cooling core region in the halo of M87 and to a smaller part from the observations in the Perseus cluster. These observations show that the central AGN produces sufficient heat for the energy balance of the cooling flow, that the most fundamental and classical accretion process originally proposed by Bondi (1952) provides an elegant way of devising a self-regulated model of AGN heating of the cooling flow, and that most of the further requirements that have to be met by a heating model to be consistent with the observations can most probably be fulfilled. Since these ideas are mostly developed to match the conditions in M87, it is important to extent such detailed studies to most other nearby cooling flow clusters.

In this new perspective the cooling cores of galaxy clusters become the sites where most of the energy output of the central cluster AGN is finally dissipated. Strong cooling flows should therefore be the locations of AGN with the largest mass accretion rates. While in the case of M87 with a possible current mass accretion rate of about 0.01 M y-1 the mass addition to the black hole (with an estimated mass of about $`310^9`$ M) is a smaller fraction of the total mass, the mass build-up may become very important for the formation of massive black holes in the most massive cooling flows, where mass accretion rates above 0.1 M y-1 become important over cosmological times.

## References

Allen, S.W., 2000, MNRAS, 315, 269

Allen, S.W., Fabian, A.C., Johnstone, R.M., et al., 2001, MNRAS, 322, 589

Ball, R., Burns, J.O., Loken, C., 1993, AJ, 105, 53

Bertschinger, E. & Meiksin, A., 1986, ApJ, 306, L1

Böhringer, H., Nulsen, P.E.J., Braun, R., Fabian, A.C., 1995, MNRAS, 274, L67

Böhringer, H., Belsole, E., Kennea, J., et al., 2001a, A&A, 365, L181

Böhringer, H., Matsushita, K., Ikebe, Y., et al. 2001b, A&A, submitted

Bondi, H., 1952, MNRAS, 112, 195

Churazov, E., Forman, W., Jones, C., Böhringer, H., 2000, A&A, 356, 788

Churazov, E., Brüggen, M., Kaiser, C.R., Böhringer, H., & Forman, W., 2001, ApJ, 554, 261

Fabian, A.C., 1994, ARA&A, 32, 277

Fabian, A.C., Nulsen, P.E.J., 1977, MNRAS, 180, 479

Fabian, A.C., Nulsen, P.E.J., & Canizares, C.R., 1984 Nature, 310, 733

Fabian, A. C., Sanders, J. S., Ettori, S., et al. 2000, MNRAS 318, 65

Fabian, A.C., Mushotzky, R.F., Nulsen, P.E.J., Peterson, J.R., 2001, MNRAS, 321, 20

Ford, H.C., et al. 1994, ApJ, 435, L27

Ikebe, Y., Makishima, K., Ezawa, H., et al., 1997, ApJ, 481, 660

Ikebe, Y., Makishima, K., Fukazawa, Y., et al., 1999, ApJ, 525, 58

Kassim, N., Perley, R.A., Erickson, W.C., Dwarakanath, K.S., 1993, AJ, 106, 2218

Makishima, K., Ezawa, H., Fukazawa, Y., et al., 2001, PASJ, 53, 401

Matsushita, K., Belsole, E., Finoguenov, A., Böhringer, H., 2001, A&A, submitted

McNamara, B., Wise, M., Nulsen, P.E.J., et al. 2000, ApJ, 534, L135

Molendi, S. & Pizzolato, F., 2001, ApJ (in press), astro-ph/0106552

Nulsen, P.E.J., 1986, MNRAS, 221, 377

Owen, F.N., Eilek, J.A., Kassim, N.E., 2000, ApJ, 543, 611

Pedlar, A., Ghataure, H. S., Davies, R. D., et al., 1990, MNRAS, 246, 477

Peterson, J.R., Paerels, F.B.S., Kaastra, J.S., et al., 2001, A&A, 365, L104

Rottmann, H., Mack, K.-H., Klein, U., Wielebinski, R., 1996, A&A, 309, L9

Silk, J., 1976, ApJ., 208, 646

Soker, N. & Sarazin, C.L., 1990, ApJ, 348, 73

Stewart, G.C., Canizares, C.R., Fabian, A.C., & Nulsen, P.E.J., 1984, ApJ, 278, 536

Tabor, G. & Binney, J., 1993, MNRAS, 263, 323

Tamura, T., Kaastra, J.S., Peterson, J.R., et al., 2001, A&A, 365, L87

Thomas, P.A., Fabian, A.C., & Nulsen, P.E.J., 1987, MNRAS, 228, 973

Tucker, W.H. & Rosner, R., 1983, ApJ, 267, 547