text:
# Tracing the Cluster Internal Dynamics with Member Galaxies

## 1. Introduction

Already in 1931, Hubble & Humason (1931) provided rough estimates of the velocity dispersions of three galaxy clusters, ranging from 500 to 1200 km s-1. This is more or less the range of values found for rich clusters nowadays. A few years later, Zwicky (1933, 1937) and Smith (1936), determined the total masses of the Coma and the Virgo cluster, respectively, by application of the virial theorem to the cluster galaxies. The cluster masses they determined111Throughout this paper I use a Hubble constant H0=70 km s-1 Mpc-1. Quoted values are scaled accordingly., $`>3.3\times 10^{14}M_{}`$ (Coma), and $`7.7\times 10^{14}M_{}`$ (Virgo), are in reasonable good agreement with nowadays estimates. Of course, these early estimates used galaxies as unbiased tracers of the underlying mass distribution, but this was just a convenient assumption. Relaxing this assumption can have a large effect on the cluster mass estimate (see, e.g., Merritt 1987). In fact, galaxies are unlikely to be distributed like the mass, simply because different cluster galaxy populations have different distributions (e.g. emission-line galaxies, ELG hereafter, and non-ELG, see Biviano et al. 1997). That is why we need to determine the cluster mass distribution, if we want to determine its total mass. And if we use the cluster galaxies for this purpose, we must at the same time determine their orbital anisotropies.

In this contribution, I review our current understanding of the internal dynamics of clusters, as determined from the spatial and velocity distributions of cluster member galaxies. Due to space limitation, I focus my review on the determination of the cluster mass, and mass-to-light profiles, and on the orbital anisotropy of cluster galaxies. Readers interested in cluster substructures are referred to the recent review of Girardi & Biviano (2001), and those interested in a complete historical review of this topic are referred to Biviano (2000).

## 2. The Jeans approach

Here I briefly describe the traditional Jeans approach (see Binney & Tremaine 1987). For the sake of simplicity, I assume to deal with a non-rotating spherically symmetric system, in a steady-state. The Abel equations relate the projected density profile, $`I(R)`$, and projected velocity dispersion profile, $`\sigma _p(R)`$ (the observables), to the de-projected density profile, $`\nu (r)`$ and, respectively, the radial and tangential components of the velocity dispersion profile, $`\overline{v_r^2}`$, and $`\overline{v_t^2}`$. $`R`$ and $`r`$ are the projected and 3D clustercentric radial distances, respectively. When $`\nu (r)`$, $`\overline{v_r^2}`$ and $`\beta (r)`$ are determined, the mass profile follows from the Jeans equation.

In principle, one should determine the de-projected profiles from the observed ones, but this is not simple, unless the anisotropy term, $`\beta (r)1\overline{v_t^2}/\overline{v_r^2}`$ is zero. The general case requires solving an integro-diffential equation (eq.4-62 in Binney & Tremaine 1987; see Binney & Mamon 1982; Solanes & Salvador-Solé 1990). This “direct” approach is complicated when the data are noisy. Therefore, the simpler “inverse” route is usually chosen. One adopts suitable functional forms for the 3D profiles – $`\nu (r)`$, $`\overline{v_r^2}`$, $`\beta (r)`$, or the mass profile, $`\rho (r)`$ –, projects them in 2D, and determines their best-fit parameters by comparison with the observables (see, e.g., Carlberg et al. 1997a, and Sect. 3.2).

The Jeans approach has several potential problems. First of all, one generally assumes the system to be spherically symmetric, but clusters are not. However, the number of available galaxies in each single cluster is limited, and it is often necessary to build a “composite” cluster by joining together the data for several clusters, after proper normalization (see, e.g., Biviano et al. 1992, 1997; Carlberg et al. 1997a). This composite cluster is spherically symmetric by construction.

Evidence for cluster rotation has been claimed for a couple of clusters (see, e.g., Biviano et al. 1996; Dupke & Bregman 2001) but in general we can safely assume that rotational motions are not important.

The assumption of steady-state is not always valid. E.g. galaxy groups are probably still in the phase of collapse (Giuricin et al. 1988). On the other hand, in a low-density Universe, clusters evolve quite rapidly, and their formation process should be over by now. Following Ellingson et al. (2001), we estimate that the infall rate of field galaxies in nearby clusters is $`<1`$% of the total cluster population per Gyr. Time-dependent terms in the Jeans equation could be needed in more distant clusters, or when galaxies are subject to tidal disruption (e.g. spirals in cluster centres, see Katgert et al. 2002).

The Jeans equation is derived from the collisionless Boltzmann equation, and its applicability is limited to collisionless fluids. Clusters of galaxies are not collisionless in their inner regions (e.g. Biviano et al. 1992; den Hartog & Katgert 1996), and dynamical friction and mergers could be important for very massive custer galaxies. It is then wise to exclude the brightest cluster members from the sample (Katgert et al. 2002).

Interlopers are not a serious problem, when one deals with a composite sample of several hundred galaxies. Different methods of interloper rejection give similar results (Fadda et al. 1996). Subclustering is a more serious issue, since galaxies in subclusters are distributed in a different way from the smooth cluster population (Biviano et al. 2001; see also Katgert’s contribution in these proceedings). It is better to identify galaxies belonging to subclusters, and to exclude them from the sample (Katgert et al. 2002).

Finally, the basic problem of the “Jeans” approach is that the solution is intrinsically degenerate; we can either assume $`\beta (r)`$ and determine $`M(<r)`$, or viceversa, but we cannot simultaneously constrain $`M(<r)`$ and $`\beta (r)`$ from the observables (see, e.g., King 1972; Bailey 1982). A possible solution consists in modelling the whole velocity distribution of cluster galaxies, and not only its second moment (Merritt 1987; Merritt & Saha 1993). Another possibility is to break this degeneracy by the use of several different tracers of the same gravitational potential (Katgert et al. 2002; see Sect. 3.6).

## 3. Mass profiles and orbital anisotropies: results so far

Seven decades since its first application to Coma (Zwicky 1933) the virial theorem has not found simple alternatives yet (see Perea et al. 1990 and references therein). The main improvement over the classical virial method has been the introduction of the “surface-term” (The & White 1986; see also Sect. 3.4). The solutions to the direct Jeans approach provided by Binney & Mamon (1982) and Solanes & Salvador-Solé (1990), have never been applied to galaxy cluster data-sets so far, because these data are noisy. Most analyses therefore model the profiles instead of deriving them directly from the data. If the models are not general enough, this procedure can lead to biased results. That is why non-parametric approaches have been developed, mostly by Merritt and collaborators (Merritt & Saha 1993; Merritt & Tremblay 1994; Merritt & Gebhardt 1994; Gebhardt & Fisher 1995).

Despite these technical improvements, the question of whether or not the cluster mass is distributed like the galaxies has not received a definite answer yet. In clusters, as well as in groups, recent works provide evidence both in favour (Carlberg et al. 1997a; Mahdavi et al. 1999; van der Marel et al. 2000), and against (Koranyi et al. 1998; Rines et al. 2000; Carlberg et al. 2001) the hypothesis of a constant mass-to-light ratio with radius. Cluster mass density profiles are found to be in agreement with the ’universal’ profile of Navarro, Frenk & White (1997; NFW hereafter), but other models are acceptable as well (Carlberg et al. 1997b; van der Marel et al. 2000; Geller, Diaferio & Kurtz 1999). From a theoretical point of view, the situation is not much clearer. Both decreasing and increasing mass-to-light profiles are predicted, resulting from the processes of tidal stripping (Mamon 1999), or, respectively, of dynamical friction and merging (Fusco-Femiano & Menci 1998).

Numerical simulations show that dark matter particles are characterized by a moderate radial anisotropy, increasing out to $`2r_{200}`$ (e.g. Diaferio 1999). However, the idea of shells of collapsing material around clusters has been around since the work of Gunn & Gott (1972). This work stimulated several investigations, and the evidence for the infall of spirals into clusters has been accumulating over the years. Indirect evidence comes from the larger velocity dispersion of spirals as compared to that of early-type galaxies (Moss & Dickens 1977; Sodré et al. 1989). More direct evidence is provided by the analyses of the Tully-Fisher distance-velocity diagram in nearby clusters (Tully & Shaya 1984; Gavazzi et al. 1991).

Additional constraints on the orbital anisotropy of cluster galaxies come from the analysis of their velocity distributions. Most studies conclude that the orbits of early-type galaxies are quasi-isotropic, and those of late-type galaxies moderately radial (Mohr et al. 1996; Biviano et al. 1997; Adami et al. 1998; van der Marel et al. 2000, see also Sect. 3.3), with the exception of Ramírez & de Souza (1998, see Sect. 3.1). Similar results are also found by the application of dynamical models (e.g. Natarajan & Kneib 1996; Carlberg et al. 1997c; Mahdavi et al. 1999).

Other methods for constraining the galaxy orbits include Pryor & Geller’s (1984) and Merrifield’s (1998). Pryor & Geller try to constrain the orbits of HI-deficient galaxies by noting that cluster core-crossing is a necessary condition for gas stripping. Merrifield suggests using the radio or X-ray trails of galaxies in order to constrain their orbits.

### 3.1. Ramírez & de Souza

Ramírez & de Souza (1998) propose to use $`u<v_j>/\sigma _p`$ as a robust shape estimator of the velocity distribution of cluster members. The value of the anisotropy of galaxy orbits is obtained by comparing the observed values of $`u`$ with those derived from a kinematical model. The model is rather over-simplified, and not very realistic: the velocity distribution is Gaussian at any radial distance, and both $`\overline{v_r^2}`$ and $`\overline{v_t^2}`$ are constant.

Van der Marel et al. (2000) point out that the assumption of Gaussianity is not acceptable for anisotropic models. Another problem is that Ramírez & de Souza (1998) do not compute $`\sigma _p`$ and $`<v_j>`$ on the same sample. E.g., when they consider the subsample of ellipticals, $`\sigma _p`$ is still computed on all cluster galaxies. This effectively makes the $`u`$ parameter a scale rather than a shape estimator. They find that ellipticals have a lower $`u`$ than spirals, but this is just a restatement of the finding of Moss & Dickens (1977; see also Adami et al. 1998), that ellipticals have a lower velocity dispersion than spirals, and do not support Ramírez & de Souza’s (1998) claim of radial anisotropy.

### 3.2. Carlberg et al.

Carlberg et al. (1997a,b) analyse the internal dynamics of galaxy clusters from CNOC, by the “inverse” Jeans approach (see Sect. 2). They adopt ad hoc functional forms for $`\nu (r)`$, $`\overline{v_r^2}`$, and $`\beta (r)`$, project them in 2D, and constrain their parameters via a $`\chi ^2`$ comparison with the observables. The advantage of their approach is that the chosen functions are well behaved, and easy to invert, but the solution for $`M(<r)`$ is to some extent imposed a priori (Merritt & Gebhardt 1994).

Their results are based on $`1000`$ galaxy members in a composite of 14 “regular” clusters. They find that galaxies trace the mass, to within $`\pm 30`$%, and the anisotropy is null or at most mildly radial, $`0<\beta <0.5`$. The mass profile is found to be consistent with a NFW profile, but other profiles (Hernquist 1990) provide equally acceptable, if not better, fits.

Carlberg et al. (2001) have recently applied the same technique on a composite sample of $`700`$ galaxies in $`200`$ groups from the CNOC2. They find a strongly increasing mass-to-light ratio with radius, at variance with the findings of Mahdavi et al. (1999). The construction of a composite group sample is anyway problematic, due to uncertainties in the mean velocity and velocity dispersions of individual groups. Moreover, removing interlopers in small galaxy groups is a difficult task. Finally, groups may be far from virialization (Giuricin et al. 1988). Therefore, I feel that Carlberg et al.’s result needs confirmation by independent analyses.

### 3.3. van der Marel et al.

The CNOC sample has recently been re-analysed by van der Marel et al. (2000) using the method developed by van der Marel (1994). They assume a three-parameter family of models for the mass density profile, $`\rho (r)`$, and a set of values for the anisotropy $`\beta `$, taken to be a constant over radii, and determine the model parameters from the best-fit to $`\sigma _p(R)`$. The galaxy density profile is derived via Abel-inversion from a parametric representation of $`I(R)`$. Their approach improves over Carlberg et al.’s, since both the adopted model for $`\rho (r)`$, and the function representing $`I(R)`$, are quite general. On the other hand, a constant $`\beta `$ is unlikely to be a realistic model.

In order to break the intrinsic degeneracy between $`\rho (r)`$ and $`\beta `$, van der Marel et al. take into account the shape of the whole velocity distribution of cluster members (following the idea of Merritt & Saha 1993). They are thus able to conclude that mass follows light to within $`\pm 30`$%, $`\beta `$ ranges from $`0.8`$ to 0.1, and a NFW profile is an acceptable fit to $`\rho (r)`$ – but not unique.

### 3.4. Girardi et al.

Girardi et al. (1998) work out the masses of 170 clusters, the largest sample analysed so far, using the virial theorem with the surface-term correction. In order to work out this correction, they first determine $`\nu (r)`$ from the Abel inversion of a (simple) parametric representation of $`I(R)`$. Assuming $`\rho (r)\nu (r)`$, they then determine which of three $`\beta (r)`$-models provides the best-fit to $`\sigma _p(R)`$. From $`\beta (r)`$ and $`\rho (r)`$, $`\overline{v_r^2}(r)`$ and the surface-term correction are derived. In principle, it is then possible to compute the corrected virial mass at any radius, hence the mass profile, but the solution is certainly biased by the rather arbitrary assumptions.

They find that the surface-term correction is important, of order 20% at the virial radius (in agreement with The & White 1986, and Carlberg et al. 1997b). They find good agreement between their optically-determined cluster masses and those derived from the X-ray data. Using the same algorithm, Girardi & Mezzetti (2001) claim no evidence for the evolution of the internal cluster dynamics up to $`z0.3`$.

### 3.5. Diaferio & Geller

Based on the results of numerical simulations of structure formation, Diaferio & Geller (1997) and Diaferio (1999) have devised a new method for determining cluster mass profiles. They show that the velocity field of cluster surroundings is determined by the cluster mass distribution, via $`GM(<r)=_0^r𝒜^2(x)_\beta (x)𝑑x`$ where $`𝒜`$ is the amplitude of the caustics in the $`(R,v)`$ space, and $`_\beta `$ is a function of the gravitational potential and of the anisotropy profile, assumed to be a constant. The caustics are determined by an adaptive kernel method.

The “caustic method” does not rely upon the hypothesis of virial equilibrium, and thus it allows estimating $`M(<r)`$ to large radii, where other methods fail. On the other hand, the choice of the caustic in the $`(R,v)`$ plane is rather subjective, and the approximation $`_\beta `$ constant is valid only beyond 1–2 $`r_{200}`$.

Using this method Geller et al. (1999) estimate $`M(<r)`$ of the Coma cluster out to 14 Mpc, and find it to be well fitted by a NFW profile. Rines et al. (2000) determine the mass-to-light profile of Abell 576 out to 6 Mpc, and find it to be decreasing with radius. A similar indication for little mass outside the virialized cluster region is found by Reisenegger et al. (2000) in the Shapley supercluster. However, the method can constrain $`M(<r)`$ only to within a factor two (Diaferio 1999), and a constant mass-to-light profile is probably still consistent with the data.

### 3.6. Katgert & Biviano

Peter Katgert and myself, we are currently analysing the internal dynamics of a composite data-set of $`3000`$ galaxies in clusters from ENACS. We build $`M(<r)`$ directly from $`I(R)`$ and $`\sigma _p(R)`$, instead of modelling $`\rho (r)`$ and $`\beta (r)`$ with ad hoc chosen functions. In fact, “using ad hoc parametrized models severely bias the answer regardless of the quality of the data” (Merritt & Gebhardt 1994). Our approach is fully non-parametric; we use LOWESS (Gebhardt & Fisher 1995) to provide a smoothed representation of the observed profiles, which we invert to provide the mass profile in the isotropic case. Using several tracers of the gravitational potential (i.e. ellipticals and S0’s, early spirals, late spirals and ELG) we expect to break the intrinsic degeneracy of the problem, to constrain the orbits of cluster galaxies, and to test whether different galaxy populations are in equilibrium within the same gravitational potential.

The main problem of this method is numerical instability. This can be solved by smoothing the data, but excessive smoothing must be avoided, not to erase useful information in the observed profiles.

Preliminary results (Biviano et al. 1999; Mazure et al. 2000) indicate that different cluster galaxy populations are characterized by different anisotropy profiles.

## 4. Four new methods and one data-set

In this Section I compare the results obtained by four different methods on the same data-set. I consider 37 ENACS (Katgert et al. 1996, 1998) clusters with at least 45 galaxy redshifts, totalling 3365 galaxies. The observational data are shown in Fig. 1, where I also display the location of the caustics, determined with the method of Diaferio (1999). From these clusters I remove interlopers (following den Hartog & Katgert 1996), and galaxies with a significant probability of being in subclusters (following Biviano et al. 2001). I then build a composite cluster, in which galaxy velocities are scaled with the cluster $`\sigma _p`$’s, and galaxy clustercentric distances are scaled with the cluster $`r_{200}`$’s. There are 1557 non-ELG and 217 ELG in this “cleaned” composite cluster. I apply the methods of Carlberg et al., Girardi et al., and van der Marel et al. to the sample of 1557 non-ELG. The last method is also applied to the sample of 217 ELG. Diaferio & Geller’s metod is instead applied on the uncleaned sample of 3365 galaxies, since it does not require preliminary removal of interlopers, or subclusters. Results of the method of Katgert & Biviano will be presented elsewhere (Katgert et al. 2002).

The advantage of using a composite cluster is that one reduces asymmetries in the galaxy distributions, thus effectively imposing the circular symmetry to the problem. But of course, the composite cluster is an artificial construction, and as such it may bear no resemblance to any real cluster. On the other hand, we may reasonably expect that results obtained for this composite cluster are representative of the average cluster population. It is important to use a homogeneous sample (like ENACS, or CNOC) for properly dealing with observational selection functions. This is particulary important for the determination of $`I(R)`$ (for further details, see Biviano et al. 2001).

In Fig. 2a,b I plot the $`I(R)`$ and $`\sigma _p(R)`$ for the 1557 non-ELG. Also shown are the three fits to these profiles, determined using the three different methods labelled in the Figure, and assuming $`\beta (r)0`$. The most remarkable differences occur where the data are less constraining, i.e. at small and large radii. They are caused by different choices of the fitting functions.

The mass profiles determined with the four different methods (assuming $`\beta (r)0`$ where needed) are shown in Fig. 3. The mass scale is arbitrary and the different values of $`M(<r_{200})`$ indicate that different effective values of $`r_{200}`$ are implied by the different methods. Approximate 1-$`\sigma `$ confidence intervals are determined by splitting the sample in two, and independently computing $`M(<r)`$ on each subsample. The four different estimates of the mass profile are not in agreement. the mass at any given radius is uncertain by a factor two, and systematic errors dominate over random errors.

In Fig.s 4a,b I show the integrated mass-to-number and mass-to-luminosity ratios as a function of radius, computed with the four different methods. On average, light seems to trace mass to within $`\pm 30`$ %, at radii $`>0.4r_{200}`$. In the central region, the mass seems to be more concentrated than the galaxies, but less than their light, an indication of luminosity segregation (see, e.g., Biviano et al. 1992).

The determination of the mass for the non-ELG is based on the assumption of isotropy. This assumption is supported by the analysis of the non-ELG velocity distribution, shown in Fig. 5a, alongwith its Gauss-Hermite representation. The Gauss-Hermite moments of this distribution clearly indicate that isotropy (or, at most, mild tangential anisotropy, $`\overline{v_r^2}/\overline{v_t^2}0.8`$) is a valid assumption.

I then consider the sample of 217 ELG. I assume that the non-ELG mass profile, in the isotropic assumption, is the true cluster mass profile, and solve for $`\beta (r)`$ of the ELG, using the method of van der Marel et al. I show in Fig. 5b the best-fit $`\chi ^2`$ solution for a constant $`\overline{v_r^2}/\overline{v_t^2}`$. For consistency, the solutions for non-ELG are also shown, where deviations from the expected $`\overline{v_r^2}/\overline{v_t^2}1`$ solution are indicative of the accuracy of the method. The data indicate a moderate degree of radial anisotropy for the ELG, $`1.3\overline{v_r^2}/\overline{v_t^2}1.6`$, at 68% confidence level. Zero anisotropy models for the ELG are rejected at $`>99.99`$% probability.

## 5. Conclusions

The use of member galaxies to determine the internal cluster dynamics is now reaching an unprecedented level of accuracy, thanks to the increasing amount of data (ENACS, Katgert et al. 1996, 1998; CNOC, Yee, Ellingson, & Carlberg 1996; Coma, Geller et al. 1999; Abell 576, Rines et al. 2000), and the development of new algorithms. The size of the data-sets must still be increased if we want to separately consider different cluster galaxy populations (e.g. early-types and late-types).

Whatever the sample size, it is difficult to discriminate among different models for the cluster mass density profile. This is because models either differ at large radii, where the cluster fades into the field, or near the centre, where one samples the cD rather than the cluster potential. The determination of cluster mass profiles appears to suffer from the systematics of different methods. Even on samples of $`1000`$ galaxies, the total mass of a cluster can be uncertain by a factor two. Such an uncertainty also affects our current knowledge of the relative distributions of mass and light. Taking the data at face value, there are now some indications that mass is indeed roughly traced by light, at least outside the very central region (where there is evidence for luminosity segregation). Finally, independent analyses suggest that different cluster galaxy populations are characterized by different orbits. The currently best-buy model has (spectral or morphological) early-type galaxies on fairly isotropic orbits, and late-type galaxies on moderately radial orbits.

Ever since Zwicky (1937) we have been aware of the need for estimating cluster masses by different means. Different methods are affected by different biases, which are sometimes difficult to properly correct for. The perspective for improving upon the current situation is therefore clear. We need to constrain the cluster internal dynamics by the complementary analyses of the distribution of cluster members, as well as of X-ray and gravitational lensing data. For this, coordinated multi-wavelength surveys of galaxy clusters, using several instruments, are needed. Most promising, in this respect, is the upcoming XMM Large Scale Structure Survey222http://vela.astro.ulg.ac.be/themes/spatial/xmm/LSS/index\_e.html, but a much more intensive spectroscopic follow-up than currently planned would be required to properly address the issue of internal cluster dynamics.

#### Acknowledgments.

I wish to thank Stefano Borgani, Marino Mezzetti and Riccardo Valdarnini for organizing such an interesting and entertaining workshop. I acknowledge useful discussions with Peter Katgert, Marisa Girardi and Alain Mazure.

This paper is dedicated to the memory of my parents, Angelo and Anita.

## References

Adami, C., Biviano, A., & Mazure, A. 1998, A&A, 331, 439

Bailey, M.E. 1982, MNRAS, 201, 271

Binney, J., & Mamon, G. 1982, MNRAS, 200, 361

Binney, J., & Tremaine, S. 1987, “Galactic Dynamics” (Princeton: Princeton Univ. Press)

Biviano, A. 2000, in “Constructing the Universe with Clusters of Galaxies”, ed. F. Durret & D. Gerbal, http://www.iap.fr/Conferences/Colloque/coll2000/contributions/

Biviano, A., Durret, F., Gerbal, D. et al. 1996, A&A, 311, 95

Biviano, A., Girardi, M., Giuricin, G., Mardirossian, F., & Mezzetti, M. 1992, ApJ, 396, 35

Biviano, A., Katgert, P., Mazure, A. et al. 1997, A&A, 321, 84

Biviano, A., Katgert, P., Thomas, T., Mazure, A., & Adami, C. 1999, 1st workshop of the italian network “Formazione ed evoluzione delle galassie”, eds.: C. Chiosi, L. Portinari, & R. Tantalo, p.45,
http://www.brera.mi.astro.it/docB/galaxy/news.html

Biviano, A., Katgert, P., Thomas, T., Mazure, A., & Adami, C. 2001, A&A, submitted

Carlberg, R.G., Yee, H.K.C., Ellingson, E. 1997b, ApJ, 478, 462

Carlberg, R.G., Yee, H.K.C., Ellingson, E. et al. 1997c, ApJ, 476, L7

Carlberg, R.G., Yee, H.K.C., Ellingson, E. et al. 1997a, ApJ, 485, L13

Carlberg, R.G., Yee, H.K.C., Morris, S.L. et al. 2001, ApJ, 552, 427

den Hartog, R. & Katgert, P. 1996, MNRAS, 279, 349

Diaferio, A. 1999, MNRAS, 309, 610

Diaferio, A. & Geller, M.J. 1997, ApJ, 481, 633

Dupke, R.A., & Bregman, J.N. 2001, ApJ, 547, 705

Ellingson, E., Lin, H., Yee, H.K.C. & Carlberg, R.G. 2001, ApJ, 547, 609

Fadda, D., Girardi, M., Giuricin, G., Mardirossian, F., & Mezzetti, M. 1996, ApJ, 473, 670

Fusco-Femiano, R. & Menci, N. 1998, ApJ, 498, 95

Gavazzi, G., Scodeggio, M., Boselli, A. & Trinchieri, G. 1991, ApJ, 382, 19

Gebhardt, K. & Fischer, P. 1995, AJ, 109, 209

Geller, M.J., Diaferio, A., & Kurtz, M.J. 1999, ApJ, 517, L23

Girardi, M. & Biviano, A. 2002, to appear in “Merging Processes in Clusters of Galaxies”, ed. L. Feretti, I.M. Gioia, & G. Giovannini (Dordrecht: Kluwer)

Girardi, M., Giuricin, G., Mardirossian, F., Mezzetti, M., & Boschin, W. 1998, ApJ, 505, 74

Girardi, M. & Mezzetti, M. 2001, ApJ, 548, 79

Giuricin, G., Gondolo, P., Mardirossian, F., Mezzetti, M., & Ramella, M. 1988, A&A, 199, 85

Gunn, J.E. & Gott, J.R. III 1972, ApJ, 176, 1

Hernquist, L. 1990, ApJ, 356, 359

Hubble, E., & Humason, M.L. 1931, ApJ, 74, 43

Katgert, P., Biviano, A. et al. 2002, in preparation

Katgert, P., Mazure, A., den Hartog, R., Adami, C., Biviano, A., & Perea, J. 1998, A&AS, 129, 399

Katgert, P., Mazure, A., Perea, J. et al. 1996, A&A, 310, 8

King, I.R. 1972, ApJ, 174, L123

Koranyi, D.M., Geller, M.J., Mohr, J.J., & Wegner, G. 1998, AJ, 116, 2108

Mahdavi, A., Geller, M.J., Böhringer, H., Kurtz, M.J. & Ramella, M. 1999, ApJ, 518, 69

Mamon, G. A. 1999, in “Dynamics of Galaxies: from the Early Universe to the Present”, ASP Conf. Series 197, ed. F. Combes, G.A. Mamon, & V. Charmandaris, p.377

Mazure, A., Biviano, A., Katgert, P., Thomas, T. & Adami, C. 2000, in “Galaxy Dynamics: from the Early Universe to the Present”, ASP Conf. Series 197, ed. F. Combes, G.A. Mamon, & V. Charmandaris, p.373

Merrifield, M.R. 1998, MNRAS, 294, 347

Merritt, D. 1987, ApJ, 313, 121

Merritt, D. & Gebhardt, K. 1994, in “Clusters of Galaxies”, ed. F. Durret, A. Mazure & J. Tran Thanh Van (Gif-sur-Yvette: Editions Frontières)

Merritt, D. & Saha, P. 1993, ApJ, 409, 75

Merritt, D. & Tremblay, B. 1993, AJ, 106, 2229

Mohr, J.J., Geller, M.J., Fabricant, D.G., Wegner, G., & Thorstensen, J. 1996, ApJ, 470, 724

Moss, C. & Dickens, R.J. 1977, MNRAS, 178, 701

Natarajan, P. & Kneib, J.-P. 1996, MNRAS, 283, 1031

Navarro, J.F., Frenk, C.S. & White, S.D.M. 1997, ApJ, 490, 493 (NFW)

Perea, J., del Olmo, A., & Moles, M. 1990, A&A, 237, 319

Pryor, C., & Geller, M.J. 1984, ApJ, 278, 457

Ramírez, A.C. & de Souza, R.E. 1998, ApJ, 496, 693

Reisenegger, A., Quintana, H., Carrasco, E.R., & Maze, J. 2000, AJ, 120, 523

Rines, K., Geller, M.J., Diaferio, A., Mohr, J.J., & Wegner, G.A. 2000, AJ, 120, 2338

Smith, S. 1936, ApJ, 83, 23

Sodré, L.J., Capelato, H.V., Steiner, J.E. & Mazure, A. 1989, AJ, 97, 1279

Solanes, J.M. & Salvador-Solé 1990, A&A, 234, 93

The, L.S. & White, S.D.M. 1986, AJ, 92, 1248

Tully, R. B. & Shaya, E. J. 1984, ApJ, 281, 31

van der Marel, R. 1994, MNRAS, 270, 271

van der Marel, R., Magorrian, J., Carlberg, R.G., Yee, H.K.C., & Ellingson, E. 2000, AJ, 119, 2038

Yee, H.K.C., Ellingson, E. & Carlberg, R.G. 1996, ApJS, 102, 269

Zwicky, F. 1933, Helv. Phys. Acta, 6, 110

Zwicky, F. 1937, ApJ, 86, 217


format: text

text:
# A search for high redshift clusters associated with radio galaxies at 2 < z < 4

## 1. Introduction

Searching for clusters of galaxies or their progenitors (proto clusters) at redshifts $`>2`$ is important for constraining galaxy evolution, large scale structure formation and cosmological models (e.g. Giavalisco et al. 1998, Bahcall, Fan, & Cen 1997). Radio galaxies can act as beacons of these distant structures.

Luminous radio galaxies are the oldest and most massive known galaxies in the early Universe (De Breuck et al. 2000). In standard cosmological models massive galaxies form in the deepest of the first potential wells, in which - at either a later or earlier time depending on the galaxy formation model - other galaxies will form to assemble a cluster. There exists a diverse collection of evidence that high redshift radio galaxies (HzRGs) are actually residing in dense cluster environments, which includes: (i) the excess of companion galaxies around HzRGs (e.g. Röttgering et al. 1996), (ii) the extreme radio rotation measures of HzRGs indicating that they are embedded in hot magnetized gas (Carilli et al. 1997), and (iii) giant Ly$`\alpha `$ emitting gas halos surrounding HzRGs which can extent up to 200 kpc (e.g. Kurk et al. 2001).

In principle, using narrow band imaging in rest-frame Ly$`\alpha `$ star forming galaxies can be found effectively over a narrow redshift range. However, since there is only a small amount of dust needed to extinguish Ly$`\alpha `$ emission from recombining hydrogen ions, it was not clear whether Ly$`\alpha `$ photons could escape at all from star forming galaxies (Pritchet 1994). Now, using large telescopes several searches for Ly$`\alpha `$ emitters are underway and show good results (Rhoads et al. 2000, Stiavelli et al. 2001). With the presumption that powerful radio galaxies are located in overdense regions, we have started a search for high redshift clusters using the narrow band imaging technique to target Ly$`\alpha `$ emitting galaxies around radio galaxies at redshifts $`2<z<4`$.

## 2. The sample

From our collection of almost 150 radio galaxies at $`z>2`$, we have chosen eight objects with redshifts suitable for Ly$`\alpha `$ imaging with the available narrow band VLT/FORS filters: four at $`z2`$ and three at $`z3`$, while for two more galaxies at $`z=4.1`$ a custom filter was manufactured. These galaxies were also selected on their large radio, optical/IR continuum and Ly$`\alpha `$ luminosities, which are characteristics of massive galaxies. For each of the ten radio galaxy fields we (plan to) do spectroscopy of $`40`$ emitter candidates and aim to confirm $`>15`$ emitters. With this sample, we can study what fraction of HzRGs are located in galaxy overdensities and how properties of these overdensities (e.g. velocity dispersions, sizes) and its cluster-galaxies (e.g. rotation curves, star formation rates) change with redshift. Simultaneously, we can study the giant gas halos associated with the radio galaxies.

Currently, imaging of five fields has been carried out, resulting in the detection of fifteen to sixty emitters in each field with Ly$`\alpha `$ rest frame equivalent width $`>20`$ Å and significant signal-to-noise. Multi object spectroscopy has been completed for two fields: MRC 1138-262 and TN J1338-1942. The latter has the largest redshift of our sample and is described in the following section. The former was the subject of our pilot project and is described in detail by Pentericci et al. (this volume). Briefly, the existence of 14 Ly$`\alpha `$ emitters in this field was confirmed, with redshifts close to the radio galaxy but distributed in two groups (grey histogram in Fig. 1).

## 3. TN J1338-1942 at z = 4.11

About 25 candidates with above mentioned selection criteria were discovered in the field around TN J1338-1942. Subsequent spectroscopy of 22 of these led to the confirmation of 20 Ly$`\alpha `$ emitters at $`z4.1`$ (Fig. 2). Their spatial distribution is inhomogeneous: it is extended in the direction of the radio axis and Ly$`\alpha `$ halo of the radio galaxy (Fig. 2). Their velocity distribution (black histogram in Fig. 1) is narrower than that of MRC 1138-262 and has a dispersion ($`\sigma _\mathrm{v}`$) of $``$ 350 km s-1 which implies a total virial mass (Mvir) of $`2\times 10^{14}`$ M$`_{}`$, assuming spherical symmetry with a virial radius (R) of 1.5 Mpc, using Mvir = 5 R $`\sigma _\mathrm{v}^2`$ / G, where G is the gravitation constant. To determine whether we actually observe an overdensity here with respect to a field which does not contain a (proto) cluster, we have compared the number density of emitters around 1338-1942 with observations of blank fields by Cowie & Hu (1998) and Rhoads et al. (2000). The number density of emitters around 1338-1942 is five times higher than the field density and comparable to the density of emitters around 1138-262.

## 4. Conclusions and outlook

Using narrow band imaging and subsequent multi object spectroscopy on the VLT, we have confirmed the existence of concentrations of galaxies around radio galaxies at $`z=2.2`$ and 4.1. These observations are consistent with the idea that HzRGs are brightest cluster galaxies in forming clusters. Furthermore, it shows that narrow band imaging at the wavelength of redshifted Ly$`\alpha `$ is an effective method to find (proto) clusters at high redshift. We expect to uncover several more of these structures in the coming years. Obviously, these systems will be subject to further study at a multitude of wavelength bands, contributing important clues to outstanding questions in galaxy formation and cosmology.

## 5. References

Bahcall, N. A., Fan, X., Cen, R. 1997, ApJ, 485, L53
Carilli, C. L., Röttgering H. J. A., van Ojik, R., Miley, G. K., & van Breugel, W. J. M. 1997, ApJS, 109, 1
Cowie, L. L. & Hu, E. M. 1998, AJ, 115, 1319
De Breuck, C., van Breugel, W., Röttgering, H. J. A., & Miley, G. K. 2000, A&AS, 143, 303
Giavalisco, M., Steidel, C. C., Adelberger, K. L., Dickinson, M. E., Pettini, M., & Kellog, M. 1998, ApJ, 503, 543
Röttgering, H. J. A., West, M. J., Miley, G. K., & Chambers, K. C. 1996, A&A, 307, 376
Kurk, J. D. et al. 2000, A&A, 358, L1
Kurk, J. D., Röttgering, H. J. A., Miley, G. K., & Pentericci, L. 2001 in RevMexAA Conf. Proc. ’Emission Lines from Jet Flows’, ed. W. Henney, W. Steffen, L. Binette, & A. Raga, astro-ph/0102337
Pentericci, L. et al. 2000, A&A, 361, L25
Pritchet, C. 1994, PASP, 106, 1052
Rhoads, J. E., Malhotra, S., Dey, A., Stern, D., Spinrad, H., & Jannuzi, B. T. 2000, ApJ, 545, L85
Stiavelli, M., Scarlata, C., Panagia, N., Treu, T., Bertin, G., & Bertola, F. 2001, ApJ in press, astro-ph/0105503


text:
# Determining the location of Gamma-Ray Bursts through the evolution of their soft X-ray absorption

## 1 Introduction

An important way to unveil the nature of the Gamma-Ray Burst (GRB) progenitors is to investigate the properties of the ambient medium surrounding them. The most popular models for GRB formation involve either the collapse of a single massive star or the coalescence of two compact objects such as two neutron stars or a neutron star and a black hole (Eichler et al. 1989; Narayan et al. 1992; Woosley 1993) These models could be constrained by knowing the GRB environment. Massive stars have very short lives, thus they will explode in star-forming regions, which are typically characterized by very dense environments. On the other hand, most merging neutron stars would be very old and would have typically traveled far from their birthplace. The scenario of compact merger progenitors could thus be suggested by medium to low density environments. Recent observations (see e.g. Kulkarni et al. 1999) for the subclass of long bursts seem to suggest their association with star-forming regions.

Within the context of GRB progenitors involving the collapse of a massive star, several proposals have been made: the so-called “hypernova” scenario (Paczyński 1998; MacFadyen & Woosley 1998) predicts a simultaneous explosion of the burst with a Type Ic supernova, while the “Supranova” model by Vietri & Stella (1998) predicts a two-step explosion, in which the burst onset follows a supernova explosion by several months. While the rebrightening of the optical afterglows detected in several bursts $`30`$ days after the burst explosion (Bloom et al. 1999; Reichart 1999; Galama et al. 2000) is suggestive of the first scenario, the detection of iron emission and absorption features in the X-ray emission of the bursts and afterglows (Amati et al. 2000; Piro et al. 2000) can be better explained with a two-step explosion (Lazzati et al. 2001a). Independent tests of these different scenarios are therefore needed. A study of the density distribution in the close environment of the bursts can be a powerful diagnostic of the type of progenitors. In fact, if the burst explodes simultaneously with the supernova, the photons propagate through the pre-explosion stellar wind \[$`n(r)r^2`$\], while in the two-step scenario a high density metal enriched supernova remnant is expected to surround the burst explosion site (Lazzati et al. 1999).

Various methods have been proposed to investigate the properties of the GRB environment. It has been shown that the effect of the surrounding material on the propagation of photons can be a probe of the material itself through the detection of absorption features and their time evolution (Perna & Loeb 1998; Meszaros & Rees 1998; Böttcher et al. 1999; Ghisellini et al. 1999; Lazzati et al. 2001abc; Ramirez-Ruiz et al. 2001a). In X-rays, it is particularly interesting the possibility to detect iron $`K_\alpha `$ absorption edges (Ghisellini et al. 1999; Böttcher et al. 1999, Weth et al. 2000) and their time evolution (Lazzati et al. 2001c; Böttcher et al. 2001). In the optical, the evolution of the equivalent width of absorption features has been investigated by Perna & Loeb (1998). On the other hand, modelling the afterglow spectral and intensity evolution in the framework of the external shock synchrotron model can also constrain the type of environment (Wijers & Galama 1999; Berger et al. 2000; Panaitescu & Kumar 2001, Ramirez-Ruiz et al. 2001b). In the first case, however, the search for the signatures of the material is made difficult by the low signal-to-noise ratio of the lines, while the second procedure heavily relies on a particular model, and hence depends on its assumptions.

In this paper, rather than considering the feature produced by a particular transition, we consider the continuum soft X-ray absorption caused by the photoionization of low-intermediate atomic number ($`Z`$) atoms along the line of sight to the burst (Morrison & McCammon 1983). This absorption is usually quantified by means of the quantity $`N_H`$, the column density of neutral material with solar metallicity that would cause the observed absorption (in the simulations we use abundances from Anders & Ebihara 1982). The presence of $`N_H`$ is generally much easier to detect with respect to a Fe $`K_\alpha `$ absorption edge. In fact, the opacity of iron at $`7.1`$ keV is less than one hundredth of the opacity at 1 keV of a material with solar iron abundance (Morrison & McCammon 1983). On the other hand, the precise measurement of the absolute $`N_H`$ value can be model dependent. In fact, the $`N_H`$ measurement is performed by computing the difference of the observed soft X-ray flux with respect to an assumed model in the same band. A spectral break at $`12`$ keV may for this reason induce a high spurious $`N_H`$ measurement.

In the following sections we present the results of numerical simulations for the evolution of the measured $`N_H`$ as a function of the time elapsed after the onset of a source of ionizing photons in its centre. We consider in particular the possibility of detecting a variation of the column density during the $`100`$ seconds of duration of a typical GRB. By imposing that the $`N_H`$ should remain above the detection threshold for a reasonable time (1 s) but vary in less that the GRB duration, we can strongly constrain the geometrical properties of the material surrounding the GRB. We then compare our results with the $`N_H`$ measurements performed in the early emission of GRB 980329 (Frontera et al. 2000) and GRB 780506 (Connors & Hueter 1998).

## 2 Numerical Simulations

In the case of a single atomic transition it is possible to derive analytical equations to describe the time evolution of its opacity in the optically thin regime (Lazzati et al. 2001c). In the case of continuum absorption the calculation is complicated by the need of keeping track of all the ionization states of the various elements and relative cross sections. For this reason the problem can be faced only with the use of numerical simulations.

We have used the time-dependent photoionization code described in Lazzati et al. (2001c), which is based on routines developed by Raymond (1979) and specialized to the case of GRBs and their remnants by Perna, Raymond & Loeb (2000) and Perna & Raymond (2000). As described by Perna et al. (2000), this code computes in a fully time-dependent fashion the ionization state of the 12 most abundant astrophysical elements, while updating, at each time step, the temperature of the plasma111See also Schwarz 1973 for related work.. The ionizing continuum is simulated in the range \[0.01-10\] keV with a spectrum relevant to GRB emission222Note that the ionizing continuum has to be simulated in a spectral band wider than the one in which the absorption is measured since soft photons can ionize ions that contribute to the opacity at higher frequencies.. We adopted a photon index $`\mathrm{\Gamma }=1`$ and constant luminosity $`L_X=1.3\times 10^{49}`$ erg s-1 in the \[2–10\] keV band for $`tt_0`$. We considered three different radial distributions of the density of absorbers: a uniform cloud with $`R=3\times 10^{18}`$ cm, a shell with $`R=3\times 10^{18}`$ cm and $`\mathrm{\Delta }R/R=0.1`$ and a wind with $`n(r)=n_0(R/R_{\mathrm{min}})^2`$ for $`RR_{\mathrm{min}}=10^{13}`$ cm and $`n(r)=0`$ for smaller radii. For the uniform and shell geometries, we performed simulations with initial column densities ranging from $`3\times 10^{19}`$ cm-2 to $`3\times 10^{24}`$ cm-2, logarithmically spaced by a factor $`\sqrt{10}`$. In the case of the wind geometry, we considered only high values of the initial column density because of its very rapid decrease with time, and we performed simulations for the initial values $`N_H(0)=10^{22}`$, $`10^{23}`$ and $`10^{24}`$ cm-2.

Even though the simulations for the uniform and shell geometries are performed for a single value of the outer radius of the absorber distribution (see above), these can be applied to rescaled matter distributions. In fact, given a distribution of absorbers $`n(r)`$ and the ionizing luminosity $`L_X`$, the results of the simulations can be extended to any distribution $`n^{}(r)=a^1n(r/a)`$ (where $`a`$ is a scale parameter) and any luminosity $`L_X^{}=bL_X`$. The time of photoionization of an ion is proportional to the flux, which scales as $`FL/r^2`$, where $`r`$ is the distance of the ion from the photon source and $`L`$ the luminosity of the ionizing photons. For this reason

$$N_{H_{\{a,b\}}}(t)=N_H\left(\frac{bt}{a^2}\right)$$
(1)
where $`N_H(t)`$ is the column density as a function of time from the simulation and $`N_{H_{\{a,b\}}}(t)`$ is the column density as a function of time for the ionizing luminosity $`L_X^{}=bL_X`$ and the distribution of absorbers $`n^{}(r)=a^1n(r/a)`$.

### 2.1 Computation of the “observed” column density evolution

The above described numerical code produces frequency and time-dependent opacities $`\sigma (\nu ,t)`$ by fully accounting for the spectral variation with time of the opacity as the gas makes the transition from neutral to fully ionized; as such, however, the output of the code contains much more information than the observational data. This is particularly true for the GRB prompt X-ray emission, usually detected with proportional counter detectors, such as the BeppoSAX Wide Field Cameras (WFC). In particular, what is usually measured from the observed spectra is the amount of soft X-ray absorption, parameterized with the hydrogen column density in the assumption that all the material is cold, i.e. neutral (Morrison & McCammon 1983). We hence applied an analogous procedure to our synthetic dataset, averaging the opacity in the observed band in order to reduce the information contained in the simulated spectrum to a column density to be compared with the observations. This is equivalent to assume that the time and frequency-dependent opacity can be decomposed as $`\tau (\nu ,t)=N_H(t)\sigma (\nu )`$, where $`\sigma (\nu )`$ is the average cross section at frequency $`\nu `$, independent of time, weighted with the abundance of each element. Our procedure is also justified by the fact that the most widely used software for X-ray data analysis (XSPEC; Arnaud 1996) neglects the change of $`\sigma (\nu )`$ with ionization to compute the column density (model WABS) when the single ion features cannot be resolved333Note that there are cases, such as that of the warm absorbers, where the effects of the spectral variations of the opacity can actually be observed through the enhancement of the absorption edges of intermediate $`Z`$ elements (see, e.g., Fabian et al. 1994 for a more refined treatment of warm absorbers). In our case, however, the situation is made more difficult by the fact that the absorber is not in the same ionization state at a given time $`t`$. The inner regions will be completely ionized at a time in which the outer regions are still cold, while gas at intermediate distance is warm.. The consequence of this assumption is that the column density at a given time has a small dependence on the band in which it is measured. In the present paper we consider the energy range \[2–10\] keV, since our aim is to compare our results with the measurements performed with the BeppoSAX WFC.

In the above approximation, we compute the column density as:

$$N_H(t)=N_H(0)\frac{\tau (\nu ,t)}{\tau (\nu ,0)}_{[210]}$$
(2)
where the symbol $`_{[210]}`$ represents the average over the frequency range \[2–10\] keV. In Fig. 1 we show the results of a simulation for a uniform cloud with $`N_H(0)=3\times 10^{21}`$ cm-2 and $`R=3\times 10^{18}`$ cm. The opacities as a function of frequency are plotted for several times after the burst onset, while the inset shows the evolution of the column density in the first 100 s of the simulation.

In Fig. 2 we show the results of two simulations with identical initial conditions, but for different geometries of the absorbing material. The evolution of the absorbing column density in a uniform cloud and a shell are shown with a solid and a dashed line, respectively. The column density of a shell environment has a roughly linear decay, slower than that of the uniform cloud. This is due to the fact that in a shell geometry all the material is located at a large distance from the source of ionizing photons, and hence the ionization time is longer. With a good dataset, this difference can be succesfully used to investigate the radial structure of the absorbing medium.

## 3 Burst environment

In several bursts, detection of variable low-energy X-ray absorption has been claimed (GRB 780506, Connors & Hueter 1998; GRB 980329, Frontera et al. 2000; GRB 990705, Amati et al. 2000; GRB 010222, in’t Zand et al. 2001). In this section we examine under which conditions such an evolution of the measured column density can be detected during the prompt emission of the bursts.

In order to detect an $`N_H`$ variation it is sufficient (but also necessary) that the following two constraints be satisfied. First, the spectrum of the GRB must be integrated in order to perform a measurement of column density. The minimum integration time (see e.g. Frontera et al. 2000) is of the order of several seconds. Variations of the column density can then be detected only if

$$N_H(1)>N_{H_{\mathrm{Gal}}}10^{21}\mathrm{cm}^2,$$
(3)
where $`N_H(1)`$ is the measured $`N_H`$ one second after the burst onset, and $`N_{H_{\mathrm{Gal}}}`$ is the Galactic column density. With this condition we ensure that the column density in the surroundings of the GRB gives a detectable signal in the first integration bin.

On the other hand, in order to detect a variation, the column density must sensibly decrease during the GRB emission time. Our second condition is then that the observed column density decrease by a factor of two in less than 100 seconds, the average burst duration444The factor $`2`$ that we use in Eq. 4 has been chosen since a sizable decrease of the absorbing column is needed in order to have a significant measurement of its variation (see below). We tested different values (namely 1.1 and 10) finding that the allowed region changes slightly, but the conclusions are not significantly affected.:

$$N_H(100)N_H(1)/2.$$
(4)
We must also consider that the burst itself has to be detected. Since the soft $`\gamma `$-rays, in which GRBs are usually detected, are not attenuated by photoabsorption, we consider Thomson scattering as the most effective attenuation process. We then have $`N_H(0)1/\sigma _T1.5\times 10^{24}`$ cm-2. If $`N_H(0)\stackrel{>}{}10^{23}`$ cm-2, the \[2–10\] keV flux is completely absorbed and only a lower limit to the column density can be measured. In the extreme case of a large and dense cloud (line-shaded region in Fig 3, 4 and 6) the opacity in the \[2–10\] keV range is larger than unity for the whole duration of the bursts, and no flux can be detected in this band.

In Fig. 3, 4 and 5 we show that conditions (3) and (4) can be satisfied only in a limited range of cloud radii, geometries and initial column densities. In Fig. 3 we show the column density at time $`t=1`$ s and $`t=100`$ s as a function of the cloud radius and initial column density with solid and dashed contours, respectively. The shaded area underlines the region in which both conditions are satisfied. This shows that, in order to detect $`N_H`$ variations during the prompt emission of the burst, the absorbing region surrounding the GRB must be compact, with a radius $`R<5`$ pc and initial column density $`N_H(0)>10^{21}`$ cm-2. This translates to volume densities $`n\stackrel{>}{}10^3`$ cm-3 (dotted lines). Figure 4 shows the same results but for a shell with $`\mathrm{\Delta }R/R=0.1`$.

The possible contribution of a wind environment to the column density evolution is considered in Fig. 5. In the figure we plot the measured $`N_H`$ as a function of time for three different initial column densities: from a thin $`N_H(0)=10^{22}`$ cm-2 wind to a maximum value of $`N_H(0)=10^{24}`$ cm-2. Figure 5 shows how, no matter the initial value of the column density, it drops to an undetectable level in less then a fraction of a second. This is due to the fact that, in the case of a wind geometry, the column density is dominated by the region very close to $`R_{\mathrm{min}}`$, for which the ionization time is very short.

In order to compare the inferred cloud properties with the properties of Galactic molecular clouds, in Fig. 6 we compare the gray region in the parameter space with the properties of Galactic molecular clouds and their denser regions. We assume that the burst explodes in the centre of these regions. In Fig. 6, the circles show the properties of a sample of cloud radii and masses presented in Leisawitz et al. (1989), from which we derived average column densities. There are no circles that fall in the allowed region shaded in grey. Inside molecular clouds, however, there are overdense regions, where star formation is supposed to take place. The triangles show the position in the diagram of a sample of dense molecular cloud cores, selected for their presence of H2O masers and in which massive star formation is taking place (Plume et al. 1997). Their radii and column densities are consistent with the shaded region. In some cases the column density is even too high. A sample of Bok globules (Launhardt et al. 1998), thought to be directely related to star formation, is shown with asterisks. Their properties are fully consistent with the shaded region. Finally, we plot with diamonds the properties of dense cores in the Taurus molecular clouds, where low-mass star formation is though to happen (Zhou et al. 1994). These are too compact with respect to the shaded regions. It is however worth mentioning that the properties of these overdense regions at the time of the burst explosion may be somewhat different, due to the radiation pressure exerted by the luminosity of the massive star progenitor to the burst.

Should these $`N_H`$ variations be confirmed, they indicate that the burst progenitors (or at least some of them) are located inside overdense regions of molecular clouds. The properties of these regions are consistent with those of the places in which massive star formation takes place.

### 3.1 GRB 980329

To date, three GRBs detected by the WFC on board BeppoSAX showed evidence of variable soft X-ray absorption. Published data are available for GRB 980329 (Frontera et al. 2000). These data are shown, along with their uncertainties, in Fig. 7. The error bars on the last two measurements, non quoted in the original paper, have been assumed equal to the measurements themselves. The timescale has been reported to the beginning of the flaring activity of the burst, which starts $`23`$ s after the trigger (see Fig. 1 of Frontera et al. 2000)

We fitted our uniform-cloud and shell models for the evolution of the column density to the data. For the uniform model we obtained the best values $`N_H(0)=1.8\times 10^{23}`$ cm-2 and $`r=0.13`$ pc ($`n4.5\times 10^5`$ cm-3), while in the shell case we obtain $`N_H(0)=1.8\times 10^{23}`$ cm-2 and $`r=0.066`$ pc. In Fig. 7, the best fit models are shown with a solid line (uniform) and dashed line (shell). The shell model has a better $`\chi ^2`$ value, but due to the large uncertainties in the measurement, the two geometries cannot be significantly disentangled. The upper right inset shows the $`1\sigma `$, 90 per cent and 99 per cent confidence contours for the two parameters in both cases. Since the statistical significance of the variation of the column density is not compelling, the 99 per cent confidence contour is not closed. The different evolution of the two models shows, however, that with good quality data it is in principle possible to understand what is the geometry of the medium surrounding GRBs.

A problem in the fit of GRB 980329 data is the lack of a redshift measurement of the burst. In fact all our simulations are performed in the rest frame. The model should be however modified in two ways if the source of absorption lies at redshift $`z>0`$. First, the timescale is stretched by a factor $`(1+z)`$. Most importantly, if a column density $`N_{H_z}`$ surrounds a source at redshift $`z`$, the soft X-ray absorption measurement will yield a measured value $`N_H=(1+z)^{2.5}N_{H_z}`$, where the exponent $`2.5`$ holds in the redshift range $`0<z<4`$. This implies that for a source at redshift $`z=1`$ the measured column density is underestimated by a factor $`5`$. The modelling of the data with a known redshift may hence give a best fit which is an even more compact and dense cloud.

### 3.2 GRB 780506

Besides the three bursts detected by the BeppoSAX WFC, a measurement of column density was possible in a burst serendipitously detected in 1987 by the HEAO-1 satellite with the A-4 and A-2 instruments (GRB 780506, Connors & Hueter 1998). The burst spectrum was fitted with a three component model: a power-law, a thermal component and an absorption component, with the column density free to vary. The lightcurve was integrated in time bins of 10.24 s each and the six spectra were fitted indipendently. The measured $`N_H`$ values are shown in Fig. 8. In the first time bin a very significant ($`6\sigma `$) soft X-ray absorption was detected, while the fits in the remaining 5 time intervals yielded only upper limits. Fitting our uniform cloud model to the data, again without the possibility to correct for the redshift of the burster, we find a best fit cloud with $`R=2\times 10^{17}`$ cm and $`N_H(0)=2\times 10^{23}`$ cm-2. Since only one positive measurement was possible, as shown by the inset in Fig. 8, the radius and initial column density of the cloud cannot be precisely determined simultaneously. We can safely conclude, however, that the density of the cloud was $`n\stackrel{>}{}10^6`$ cm-3. Since for this burst only one positive measurement was obtained, it is not possible to distinguish between a uniform and a shell geometry.

## 4 Discussion and Conclusions

We have performed numerical simulations to compute the evolution of the soft X-ray absorption of a cloud under the ionization flux of a Gamma-Ray Burst exploded at its centre. We have considered a uniform cloud, a shell and a wind geometry, and we have shown that variations of the column density during the duration of a GRB can be detected only under well-defined geometrical conditions of the absorbing material.

The column density of a wind, in which the bulk of the material responsible for the absorption is close to the GRB, cannot survive the first milliseconds of the illumination and is undetectable with the present instrumentation. On the other hand, if the material is spread uniformly over a region with size $`0.05\stackrel{<}{}R\stackrel{<}{}5`$ pc, the soft X-ray absorption can be observed to evolve in time. Present data do not allow to firmly detect the evolution of the column density, nor to tell if the material is uniformly spread around the source or rather concentrated in a shell at distance $`R`$ from the burster. By comparing the allowed properties of molecular clouds with the results presented here, we conclude that the weak evidence of $`N_H`$ variation during the prompt X-ray emission of some GRBs can be accounted for if they are associated to compact overdense regions similar to the cocoon of star formation within molecular clouds or to Bok globules.

A number of approximations have been applied to this analysis and are worth a discussion. First, we consider an ionizing flux with fixed spectrum and luminosity. Both these assumptions are coarse approximations of a GRB lightcurve and spectral evolution. Even though GRB are frantically variable, the results of our simulations are exact if the time $`t`$ is substituted with the quantity $`t^{}=_0^tL(t)𝑑t/_0^{t_{\mathrm{GRB}}}L(t)𝑑t`$, where $`t_{\mathrm{GRB}}`$ is the total burst duration. Since in this work the detailed time variability of the $`N_H`$ lightcurve is not addressed, the approximation of constant luminosity is reasonable.

A more important issue is the variation of the spectral slope with time. More complete simulations of time-dependent column densities will be presented in a forthcoming paper (Perna et al. in preparation), in which the dependence of the observed column density on the ionizing spectrum will be discussed in detail. However, $`\mathrm{\Gamma }=1`$ is a good value for the average time integrated low-energy slope of GRBs (Preece et al. 2000).

We have assumed that the medium is at rest during the photoionization process. While Lazzati et al. (2001c) showed that the fireball interaction is not important, radiation pressure may accelerate the medium. The ratio of the timescale of photoionization $`t_{\mathrm{pi}}`$ over the timescale of acceleration $`t_{\mathrm{acc}}`$ is given, for a $`\mathrm{\Gamma }=1`$ spectrum, by:

$$\frac{t_{\mathrm{pi}}}{t_{\mathrm{acc}}}=\frac{h\nu _\mathrm{b}\sigma _T}{m_pc^2\sigma _{\mathrm{pi}}}2\times 10^6\left(\frac{h\nu _b}{511\mathrm{keV}}\right)\left(\frac{h\nu _{\mathrm{th}}}{1\mathrm{keV}}\right)^{2.3}$$
(5)
where $`\nu _b`$ is the break energy of the GRB spectrum and $`\sigma _{\mathrm{pi}}2\times 10^{22}(h\nu _{\mathrm{th}}/1\mathrm{keV})^{2.3}`$ is the average photoionization cross section at the threshold energy $`h\nu _{\mathrm{th}}`$ weighted with the element solar abundance (Anders & Ebihara 1982). Even for iron, the relevant element with the highest threshold energy, the acceleration timescale is several orders of magnitude longer than the photoionization one.

Finally, as discussed in §4.1, the redshift of the burst introduces time dilation and decreases the observed column density at a given time. All these factors have not been taken into account in Fig. 3 and Fig. 4, so that the gray shaded regions should be considered indicative, with uncertainties up to a factor of order unity. Nevertheless we can firmly conclude that if a column density variation is observed in the first $`100`$ seconds of the prompt emission of a GRB, then it must be located within a cloud with properties roughly consistent with those of Fig. 3 or inside a shell (Fig. 4). We have not considered here the possibility of a beamed fireball. In the case of a conical jet, as discussed to date for GRB fireballs (Sari et al. 1999, Ghisellini & Lazzati 1999, Panaitescu & Kumar 2001), beaming does not affect our calculations, since all processes take place along the line of sight. In the case of a cylindrical jet, the radial dilution of photons would be different, and our computations are not appplicable.

Variations in the soft X-ray absorption have indeed been claimed in several bursts. Only one of them, however, shows this evidence quite clearly. GRB 780506, serendipitously observed by the HEAO-1 satellite (Connors & Hueter 1998) has a highly significant soft X-ray absorption in the first time integration bin, which vanishes in the rest of the burst lightcurve (Fig. 8). GRB 980329 (Frontera et al. 2000), shows a very well defined trend of the $`N_H`$ to decrease with time, which can be fit in our model (see Fig. 7), but with large statistical uncertainties. GRB 990705 (Amati et al. 2000) shows evidence of time dependent iron absorption features. An evolution of the column density is also claimed, but the data are not published. Finally, GRB 010222 (in ’t Zand et al. 2001) shows clear evidence of $`N_H`$ variations, but the column density appears to rise and fall, a behavior that cannot be explained in our model and is likely due to the presence of a spectral break that moves inside and outside the observing window.

As a consequence of the high density ($`n10^3`$) described above, the afterglow of the bursts should evolve more rapidly than in a very low density environment. In particular, the time at which the relativistic to non-relativistic transition takes place is $`t_{\mathrm{nr}}10(E_{53}/n_6)^{1/3}`$ d (see, e.g. Wijers, Rees & Meszaros 1997), where $`E_{53}`$ is the isotropic fireball energy in units of $`10^{53}`$ erg and $`n_6`$ the density of the surrounding medium in units of $`10^6`$ cm-3. Unfortunately, there are not accurate afterglow lightcurves for any of the bursts described in this paper. Indeed, a very high density of the surrounding medium is required in GRB 990705 (Lazzati et al. 2001a). In the case of GRB 010222, which displays $`N_H`$ variations but cannot be straightforwardly included in our model, a very high density of the surrounding medium is claimed (in ’t Zand et al. 2001) on the base of spectral and lightcurve considerations.

The HETE2 satellite has on-board two experiment, the Soft X-ray Camera (SXC) and the Wide Field X-ray Monitor (WXM) potentially sensitive to $`N_H`$ variations. The SXC, has an effective area of $`10`$ cm2, is sensitive in the \[0.5–10\] keV range and has a field of view of 0.91 sr. The WXM has a larger effective area (175 cm2) and field of view ($`1.6`$ sr) but is sensitive in the harder \[2–25\] keV band, where continuum absorption is less severe. The combination of these two instruments, however, will hopefully clarify the properties and time variability of the soft X-ray absorption in the early prompt X-ray spectra of GRBs.

## Acknowledgments

We wish to thank the anonymous referee for his/her careful reading of the manuscript and useful comments. We thank Gabriele Ghisellini, Enrico Ramirez-Ruiz, Martin Rees, Luigi Stella and Mario Vietri for very useful discussions and suggestions. We thank John Raymond for his help in the process of updating the numerical code used in this work.


text:
# Pairing Instability and Mechanical Collapse of a Bose Gas with an Attractive Interaction

## Abstract

We study the pairing instability and mechanical collapse of a dilute homogeneous bose gas with an attractive interaction. The pairing phase is found to be a saddle point and unstable against pairing fluctuations. This pairing saddle point exists above a critical temperature. Below this critical temperature, the system is totally unstable in the pairing channel. Thus the system could collapse in the pairing channel in addition to mechanical collapse. The critical temperatures of pairing instability and mechanical collapse are higher than the BEC temperature of an ideal bose gas with the same density. When fluctuations are taken into account, we find that the critical temperature of mechanical collapse is even higher. The difference between the collapse temperature and the BEC temperature is proportional to $`(n|a_s|^3)^{2/9}`$, where $`n`$ is the density and $`a_s`$ is the scattering length.

The system of Bose gas with attractive interaction has received a lot of attention recently due to the success in cooling 7Li and 85Rb systems. In these trapped systems, Bose Einstein condensation can be realized at zero temperature if the number of atoms in the condensate is below a critical value, because of a finite size effect. In a 7Li system, the condensate experiences repetitive growth and collapse as the system is cooled down. In a 85Rb system, the interaction can be varied from being repulsive to being attractive by tuning a Feshbach resonance. During this process, the system loses particles until the critical number is reached.

A homogeneous Bose gas with attractive interaction is generally unstable at zero temperature because the system can always lower energy by increasing density. At high temperatures entropy stabilizes the gaseous phase. There are several scenarios as to how the system becomes unstable when the temperature is lowered. The simplest picture is that the system collapses within the normal gaseous phase, which is supported by several theoretical studies. This collapse is mechanical collapse, where the compressibility becomes negative and the system can lower energy by separating into a low density phase and a high density phase. Another possibility is that the system goes into the Bose-condensed phase and then collapses in the condensed phase. However the condensation transition temperature is found to be lower than that of mechanical collapse. In addition to these scenarios, Evans and Imry proposed that there is a pairing phase similar to the BCS phase in superconductors at low temperature. However it was also concluded in the previous theoretical studies that this pairing phase transition temperature is lower than the mechanical collapse temperature of the normal phase.

It was not fully addressed in the previous studies whether the pairing phase could exist as a metastable phase or not. In this paper, we examine its stability in detail. In addition, we go beyond the Hartree-Fock approximation to see how the temperature of mechanical collapse is affected by fluctuations.

The interaction between atoms is singular at short distances. However, as for other cases in condensed matter and elementary particle physics, at low energies the dilute Bose gas system can be well described by the Hamiltonian with an effective contact interaction with the coupling constant proportional to the scattering length $`a_s`$, up to a momentum transfer which is determined by the effective range of the interaction. The unphysical ultraviolet divergences in this approach can be subtracted by the pseudopotential method. This renormalization scheme is justified in the case of a dilute Bose gas with weak interaction at zero temperature, where it can be shown that the divergences come from double counting of certain diagrams and the dominant perturbations come from the repetitive scatterings of two particles.

At finite temperature when there is no condensate, the medium effect becomes important and the coupling constant could depend on other quantities in addition to $`a_s`$. To avoid this complexity, here we use a simple model with a contact interaction which has coupling constant $`g`$ and an explicit cutoff $`\mathrm{\Lambda }`$ in the wave-vector space. The Hamiltonian of the system is given by

$$H=\frac{\mathrm{}^2}{2m}\psi ^{}\psi +g\psi _{}^{}{}_{}{}^{2}\psi ^2.$$
(1)
A length scale $`a`$ can be defined from the coupling constant, $`am/(4\pi \mathrm{}^2g)`$, but $`a`$ is not the scattering length $`a_s`$. By solving the two-body scattering problem, we obtain the relation between $`a`$ and $`a_s`$, $`a_s=a/(1+2a\mathrm{\Lambda }/\pi )`$. The cutoff $`\mathrm{\Lambda }`$ is proportional to $`1/\sqrt{|a_s|r_s}`$, where $`r_s`$ is the effective range. Here we consider the case of a weakly attractive interaction with $`a<0`$, $`2a\mathrm{\Lambda }/\pi <1`$, and $`|na^3|1`$, where $`n`$ is the density.

To study this system, we use the Peierls variational method which gives an upper bound for the free energy or grand thermodynamic potential of a quantum system. To discuss both the normal phase and the pairing phase, we assume the quasiparticles are superpositions of particles and holes, and they have certain occupation numbers at finite temperature. In momentum space, the density matrix is given by $`<\psi _𝐤^{}\psi _𝐤>=(f_k+1/2)\mathrm{cosh}\theta _k1/2`$, and $`<\psi _𝐤\psi _𝐤>=(f_k+1/2)\mathrm{sinh}\theta _k`$, where $`\mathrm{cosh}\theta _k`$ and $`\mathrm{sinh}\theta _k`$ are coherence factors and are chosen to be real. The occupation number is given by $`f_k`$.

Within this approach, the expectation value of the grand potential is given by

$`\mathrm{\Omega }(\{\theta _k\},\{f_k\})={\displaystyle \underset{𝐤}{}}(ϵ_k\mu )[(f_k+1/2)\mathrm{cosh}\theta _k1/2]+{\displaystyle \frac{g}{2V}}\{{\displaystyle \underset{𝐤}{}}(f_k+1/2)\mathrm{sinh}\theta _k\}^2`$
$`+{\displaystyle \frac{g}{V}}\{{\displaystyle \underset{𝐤}{}}[(f_k+1/2)\mathrm{cosh}\theta _k1/2]\}^2T{\displaystyle \underset{𝐤}{}}[(f_k+1)\mathrm{ln}(f_k+1)f_k\mathrm{ln}f_k].`$
The normal phase and the pairing phase are given by the saddle points of the variational space, $`\mathrm{\Omega }/f_k=0`$, $`\mathrm{\Omega }/\theta _k=0`$. The saddle point equations can be further written as

$`\mathrm{tanh}\theta _k`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }}{ϵ_kA\mu }},`$ (2)
$`f_k`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{coth}{\displaystyle \frac{E_k}{2T}}{\displaystyle \frac{1}{2}},`$ (3)
where

$$\mathrm{\Delta }=\frac{g}{V}\underset{𝐤}{}(f_k+1/2)\mathrm{sinh}\theta _k,$$

$$A=\frac{2g}{V}\underset{𝐤}{}[(f_k+1/2)\mathrm{cosh}\theta _k1/2],$$

$`E_k=\sqrt{\eta _k^2\mathrm{\Delta }^2}`$, and $`\eta _k=ϵ_kA\mu `$. The pairing phase corresponds to the solution with finite $`\mathrm{\Delta }`$ and the normal phase corresponds to the solution with $`\mathrm{\Delta }=0`$. The transition from normal phase to pairing phase occurs when $`\mathrm{\Delta }`$ approaches zero.

To simplify the discussion without compromising the saddle point structure, we look at a more restricted variational space defined by Eq.(2) and Eq.(3) with $`\mathrm{\Delta }`$ and $`A`$ as variational parameters. In terms of these new variational parameters, the grand potential is given by

$$\mathrm{\Omega }(\mathrm{\Delta },A)=\underset{𝐤}{}(ϵ_k\mu )\left[\frac{\eta _k}{2E_k}\mathrm{coth}\frac{E_k}{2T}\frac{1}{2}\right]+\frac{V\stackrel{~}{\mathrm{\Delta }}^2}{2g}+\frac{V\stackrel{~}{A}^2}{4g}\underset{𝐤}{}\frac{E_k}{2}\mathrm{coth}\frac{E_k}{2T}+T\underset{𝐤}{}\mathrm{ln}\left(2\mathrm{sinh}\frac{E_k}{2T}\right),$$
(4)
where the parameters $`\stackrel{~}{\mathrm{\Delta }},\stackrel{~}{A}`$ are defined by

$$\stackrel{~}{\mathrm{\Delta }}(\mathrm{\Delta },A)=\frac{g}{V}\underset{𝐤}{}(f_k+\frac{1}{2})\mathrm{sinh}\theta _k=\frac{g}{V}\underset{𝐤}{}\frac{\mathrm{\Delta }}{2E_k}\mathrm{coth}\frac{E_k}{2T},$$

and

$$\stackrel{~}{A}(\mathrm{\Delta },A)=2\frac{g}{V}\underset{𝐤}{}\left[(f_k+\frac{1}{2})\mathrm{cosh}\theta _k\frac{1}{2}\right]=\frac{g}{V}\underset{𝐤}{}\left[\frac{\eta _k}{E_k}\mathrm{coth}\frac{E_k}{2T}1\right].$$

The new saddle point equations are given by $`\stackrel{~}{A}=A`$, $`\stackrel{~}{\mathrm{\Delta }}=\mathrm{\Delta }`$.

The stable phases correspond to the local minimum points of the grand potential. If the pairing phase is stable, one of the necessary conditions is that it should be stable against fluctuations in pairing, which is given by

$$\frac{^2\mathrm{\Omega }}{\mathrm{\Delta }^2}|_{\stackrel{~}{A}=A,\stackrel{~}{\mathrm{\Delta }}=\mathrm{\Delta }}>0.$$
(5)
Since the transition into the pairing phase would be a second-order phase transition, at the critical temperature $`T_P`$,

$$\frac{^2\mathrm{\Omega }}{\mathrm{\Delta }^2}|_{\stackrel{~}{A}=A,\stackrel{~}{\mathrm{\Delta }}=\mathrm{\Delta }}=0,$$
(6)
and the stability condition Eq.(5) becomes

$$\frac{^4\mathrm{\Omega }}{\mathrm{\Delta }^4}|_{\stackrel{~}{A}=A,\mathrm{\Delta }=0}>0.$$
(7)
However, we find that the opposite inequality holds

$$\frac{^4\mathrm{\Omega }}{\mathrm{\Delta }^4}|_{\stackrel{~}{A}=A,\mathrm{\Delta }=0}=\frac{3V}{2g}\left[\frac{^2\stackrel{~}{A}}{\mathrm{\Delta }^2}\right]^2\underset{𝐤}{}\frac{3}{2\eta _k^3}(1+2f_k^2)<0.$$
(8)
It clearly shows that the pairing phase is generally unstable. The instability of the pairing phase can also be illustrated by plotting the profile of the grand potential as shown in Fig.(1). The pairing phase is located at a saddle point and the normal phase is located at a minimum point. The pairing saddle point is a maximum in the direction of the normal phase and a minimum in the perpendicular direction.

This instability can be studied in more detail by constructing an effective Lagrangian near the critical temperature $`T_P`$ based on Landau’s theory of phase transitions

$$(\mathrm{\Delta })\mathrm{\Omega }(\mathrm{\Delta },\stackrel{~}{A})\mathrm{\Omega }(0,\stackrel{~}{A})\alpha \mathrm{\Delta }^2+\beta \mathrm{\Delta }^4,$$
(9)
where $`\alpha =\frac{^2\mathrm{\Omega }}{2\mathrm{\Delta }^2}|_{\stackrel{~}{A}=A,\mathrm{\Delta }=0}`$ and $`\beta =\frac{^4\mathrm{\Omega }}{24\mathrm{\Delta }^4}|_{\stackrel{~}{A}=A,\mathrm{\Delta }=0}`$, $`\beta <0`$ from Eq.(8). When $`T>T_P`$, $`\alpha >0`$, there is a pairing saddle-point solution $`\mathrm{\Delta }^2=\alpha /(2\beta )`$. When $`T<T_P`$, $`\alpha <0`$, the normal phase becomes unstable in the pairing channel and there is no pairing saddle-point solution. It is clear that the system is unstable in the pairing channel when the temperature is below $`T_P`$. Even above $`T_P`$, the normal phase can become unstable in the pairing channel by tunneling into the region beyond the pairing saddle point if any place in that region has lower energy than the normal phase. The pairing saddle point serves as an energy barrier in this case. However, in our limited numerical calculations, we have not been able to determine the condition for the normal phase to become metastable. When the temperature is far above $`T_P`$, the pairing saddle point disappears. In this case the effective Lagrangian description breaks down and a better method is needed to find the general expression of the critical temperature below which the pairing saddle point appears.

The critical temperature of pairing instability $`T_P`$ is given by the following equation

$$1=g\frac{d^3k}{(2\pi )^3}\frac{1}{2\eta _k}\mathrm{coth}\frac{\eta _k}{2T_P},$$
(10)
where $`\eta _k=ϵ_k+2gn\mu `$. This equation has solution only when the temperature is close to the ideal BEC temperature $`T_0`$. In this case it can be further simplified,

$$1+\frac{2a\mathrm{\Lambda }}{\pi }g\frac{d^3k}{(2\pi )^3}\frac{T_P}{\eta _k^2}$$
(11)
and

$$1+\frac{2a\mathrm{\Lambda }}{\pi }\frac{T_Pa}{\mathrm{}}\sqrt{\frac{2m}{\eta _0}}.$$
(12)
From the density equation, we obtain

$$\eta _0\frac{9}{16\pi }\zeta ^2(\frac{3}{2})T_0(\frac{T}{T_0}1)^2,$$
(13)
where

$$T_0=\frac{2\pi \mathrm{}^2}{m}(\frac{n}{\zeta (\frac{3}{2})})^{\frac{2}{3}}$$

and $`\zeta (x)`$ is the Riemann-Zeta function. Therefore the solution of Eq.(10) is given by

$$\frac{T_P}{T_0}1\frac{2\pi }{3\zeta ^{\frac{4}{3}}(\frac{3}{2})(1+\frac{2a\mathrm{\Lambda }}{\pi })}n^{\frac{1}{3}}|a|=\frac{8\pi }{3\zeta ^{\frac{4}{3}}(\frac{3}{2})}n^{\frac{1}{3}}|a_s|.$$
(14)
In addition to the pairing instability, mechanical collapse can also occur in the normal phase as discussed in several references. It is important to find out which instability takes place first. The mechanical collapse happens when the compressibility goes to zero $`\mu /n=0`$, which can be more explicitly written as

$$\frac{g}{2T_C}\frac{d^3k}{(2\pi )^3}\mathrm{csch}^2\frac{\eta _k}{2T_C}=1,$$
(15)
where $`T_C`$ is the temperature of mechanical collapse. Similar to Eq.(10), Eq.(15) can only be satisfied when $`T_C`$ is close to $`T_0`$, and can be simplified as

$$1g\frac{d^3k}{(2\pi )^3}\frac{2T_C}{\eta _k^2}$$
(16)
and

$$1\frac{2T_Ca}{\mathrm{}}\sqrt{\frac{2m}{\eta _0}}.$$
(17)
The solution is given by

$$\frac{T_C}{T_0}1\frac{16\pi }{3\zeta ^{\frac{4}{3}}(\frac{3}{2})}n^{\frac{1}{3}}|a|.$$
(18)
Similar results have been obtained in the previous studies.

By comparing the two critical temperatures Eq.(14) and Eq.(18), we find that when $`4a\mathrm{\Lambda }/\pi >1`$, $`T_P>T_C`$; when $`4a\mathrm{\Lambda }/\pi <1`$, $`T_P<T_C`$. However, this result is obtained within the variational approach, which is equivalent to the Hartree-Fock approximation. To obtain a better quantitative result, fluctuations need to be considered.

There are infrared divergences in all orders of perturbation at $`T_0`$. To obtain the correct renormalization near $`T_0`$, summing the leading divergent terms at each order is necessary, which can be effectively performed in a self-consistent approach. Here we consider the simplest one-loop renormalization with self consistency. In the Hartree-Fock approximation, the self-energy in the normal phase is simply given by $`2gn`$. When fluctuations are taken into account, the coupling constant $`g`$ is renormalized. In one-loop order, $`g`$ is re normalized by particle-particle scattering. It is well known that at zero temperature, such scatterings renormalize $`g`$ to $`4\pi \mathrm{}^2a_s/m`$. In this approach, the renormalized coupling constant is given by the $`t`$-matrix

$$t(𝐪,\mathrm{\Omega })=gT\underset{\omega }{}\frac{d^3k}{(2\pi )^3}gt(𝐪,\mathrm{\Omega })G(𝐤,\omega )G(𝐪𝐤,\mathrm{\Omega }\omega ),$$
(19)
where $`G(𝐤,\omega )=1/(i\omega \mathrm{\Sigma }(𝐤,\omega )ϵ_k+\mu )`$ is the finite-temperature Green’s function. Under this approximation, the new self-energy is given by

$$\mathrm{\Sigma }(𝐤,\omega )=2T\underset{\mathrm{\Omega }}{}\frac{d^3q}{(2\pi )^3}t(𝐪,\mathrm{\Omega })G(𝐤𝐪,\omega \mathrm{\Omega }).$$
(20)
The infrared behavior of the integrals in Eq.(19) and Eq.(20) is dominated by $`\mathrm{\Sigma }(0,0)\mu `$, which is small and finite at the collapse temperature. The $`k`$-dependence and $`\omega `$-dependence of $`\mathrm{\Sigma }`$ just provides a renormalization factor to the propagator in the leading order. So here we only consider only the constant part of self-energy,

$$\mathrm{\Sigma }\mathrm{\Sigma }(0,0)2t(0,0)n,$$

and get

$$\frac{g}{t(0,0)}1+g\frac{d^3k}{(2\pi )^3}\frac{1}{2\eta _k^{}}\mathrm{coth}\frac{\eta _k^{}}{2T},$$
(21)
where $`\eta _k^{}=ϵ_k+\mathrm{\Sigma }\mu `$. Using the same technique in solving Eq.(10), we obtain $`t(0,0)`$ in terms of $`\mathrm{\Sigma }`$

$$t(0,0)\frac{g}{1+\frac{2a\mathrm{\Lambda }}{\pi }+\frac{Ta}{\mathrm{}}\sqrt{\frac{2m}{(\mathrm{\Sigma }\mu )}}}.$$
(22)
Therefore the self-energy $`\mathrm{\Sigma }`$ is given by the self-consistent equation

$$\mathrm{\Sigma }\frac{2gn}{1+\frac{2a\mathrm{\Lambda }}{\pi }+\frac{Ta}{\mathrm{}}\sqrt{\frac{2m}{(\mathrm{\Sigma }\mu )}}}.$$
(23)
The condition of mechanical collapse $`\mu /n=0`$ now leads to the following equation

$`1`$ $`={\displaystyle \frac{d^3k}{(2\pi )^3}\frac{1}{4T}\frac{\mathrm{\Sigma }}{n}\mathrm{csch}^2\frac{ϵ_k+\mathrm{\Sigma }\mu }{2T}}`$ (25)
$`{\displaystyle \frac{Tm}{4\pi \mathrm{}^3}}{\displaystyle \frac{\mathrm{\Sigma }}{n}}\sqrt{{\displaystyle \frac{2m}{\mathrm{\Sigma }\mu }}},`$
where from Eq.(23) we obtain

$$\frac{\mathrm{\Sigma }}{n}\frac{2g[1+\frac{2a\mathrm{\Lambda }}{\pi }+\frac{Ta}{\mathrm{}}\sqrt{\frac{2m}{(\mathrm{\Sigma }\mu )}}]}{[1+\frac{2a\mathrm{\Lambda }}{\pi }+\frac{Ta}{\mathrm{}}\sqrt{\frac{2m}{(\mathrm{\Sigma }\mu )}}]^2\frac{Ta}{\mathrm{}}gn\sqrt{\frac{2m}{(\mathrm{\Sigma }\mu )^3}}}.$$
(26)
The solution of Eq.(25) is approximately given by

$$(1+\frac{2a\mathrm{\Lambda }}{\pi })^2\frac{8\pi \mathrm{}a^2nT}{\sqrt{2m(\mathrm{\Sigma }\mu )^3}}.$$
(27)
Close to $`T_0`$, the density equation yields

$$\mathrm{\Sigma }\mu \frac{9\zeta ^2(\frac{3}{2})\delta T^2}{16\pi T_0},$$
(28)
where $`\delta T=TT_0`$. Combining Eq.(27) with Eq.(28), we obtain the new critical temperature of mechanical collapse $`T_C^{}`$

$$\frac{T_C^{}}{T_0}14\frac{(2\pi |a_s|n^{\frac{1}{3}})^{\frac{2}{3}}}{3\zeta ^{\frac{8}{9}}(\frac{3}{2})}.$$
(29)
The difference between $`T_C^{}`$ and $`T_0`$ is now proportional to $`(n|a_s|^3)^{2/9}`$. In the weakly-interacting limit, it is much bigger than $`n^{1/3}|a_s|`$ which is proportional to the temperature difference in Hartree-Fock approximation. The strong renormalization of the collapse temperature is a result of the enhancement of the Hartree-Fock attractive interaction due to fluctuations. In contrast, the critical temperature of pairing instability is not renormalized by the $`t`$-matrix. We have used the one-loop and constant self-energy approximations which may limit the accuracy of the numerical coefficient. Future works with better approximations may improve the numerical prefactor and provide an estimate of the next order term.

Our result shows that the collapse temperature is higher than the temperature of pairing instability. However the pairing saddle point exists above the temperature of pairing instability. The saddle point may still be important to collapse dynamics if it appears above the collapse temperature. In this case, if the normal phase is only a local minimum and not a global minimum of the grand potential, there is a finite probability for the system to pass the pairing saddle point and tunnel into the unstable region. Discussing the dynamical process of collapse is beyond the scope of this paper. However, within this variational framework, the collapse process can probably be studied by the general method of analyzing the stability of metastable states in the future.

We would like to thank the Aspen Center for Physics where some parts of this work were produced. This work is supported by the National Science Foundation under Grant No. DMR 98-15932.


text:
# Magnetotransport properties of untwinned YBa2Cu3Oy single crystals: Novel 60-K-phase anomalies in the charge transport

## I Introduction

In YBa2Cu3Oy (YBCO) compound, an increase in the oxygen content $`y`$ from 6 to 7 causes the hole doping into the CuO2 planes (above $`y6.2`$) and leads to superconductivity (above $`y6.4`$). The dependence of $`T_c`$ on $`y`$ is non-trivial and there is a well-known plateau with $`T_c`$ 60 K for $`y`$ of around 6.7 (“60-K plateau” or “60-K phase”). Both an oxygen-ordering scenario  and an 1/8-anomaly scenario  have been discussed for the origin of the 60-K plateau, but the case remains controversial. To understand the true nature of the 60-K phase in YBCO, we have conducted systematic measurements of the transport properties across the 60-K phase, which reveal clear electronic anomalies in the 60-K phase and thus indicate an electronic origin of the 60-K plateau. In this paper, we summarize the novel 60-K-phase anomalies we found in the charge transport in YBCO and discuss their possible origin in conjunction with the self-organization of the holes into stripes.

## II Experimental

It should be noted that the oxygen arrangement in the Cu-O chain layers is essentially random, causing complications to the study of YBCO; for example, for a given $`y`$ the actual hole doping can differ depending on the arrangement of the O atoms, and the O atoms in the Cu-O chains can rather easily rearrange at room temperature, which causes the “room-temperature (RT) annealing effect” . For this work, the crystals are always quenched at the end of the high-temperature annealing (which tunes the oxygen content) and detwinning is performed at temperatures below 220C after the annealing. The samples are left at room temperature for at least a week for the oxygen arrangement to equilibrate before the measurements; therefore, the oxygen atoms on the chain sites are expected to be locally ordered (because of the RT annealing) but macroscopically uniform (because of the quenching) . (Note that a long-time annealing at relatively low temperature ($``$100C) causes a macroscopic phase separation in heavily-underdoped samples  and messes up the transport properties.) Our procedure ensures very good reproducibility of the transport properties, as has been demonstrated in Ref. . In particular, we have determined the absolute values of the resistivity and the Hall coefficient for a given $`y`$ with the accuracy of 5%; this gives us confidence in discussing the systematics of the transport properties across the 60-K phase. Details of the measurement techniques are described in Refs. .

## III Results and Discussions

Let us start with the resistivity behavior. Figure 1(a) shows the temperature dependences of $`\rho _a`$ for a wide range of oxygen contents, $`y=6.457.0`$. Remember that the Cu-O chains run along the $`b`$-axis and thus $`\rho _a`$ is not complicated with the conductivity of the chains . We emphasize that at least 3 samples are measured for each $`y`$ and the data are reproducible within 5%; in fact, Fig. 1(a) is a summary of the measurements of more than 30 samples. One can immediately notice in Fig. 1(a) that the $`\rho _a(T)`$ data for $`y=6.656.80`$ show clear overlap below $``$130 K. Note that in the underdoped YBCO the pseudogap opening can be inferred from a downward deviation from the high-temperature $`T`$-linear dependence , and thus the overlapping of $`\rho _a(T)`$ is observed in the pseudogap state. Figure 1(b) shows the corresponding evolution of $`\rho _b(T)`$ with $`y`$; $`\rho _b`$ is generally smaller than $`\rho _a`$ for the same $`y`$, which is believed to be due to the finite chain conductivity . The $`\rho _b(T)`$ data do not show as clear overlap in the 60-K phase as the $`\rho _a(T)`$ data.

The overlap of the $`\rho _a(T)`$ data is very unusual. Unless the effective mass $`m^{}`$ of the charge carrier in YBCO is anomalously changing with $`y`$ (which is very unlikely), the $`y`$-independence of $`\rho _a`$ can have only two possible origins: (i) both the carrier concentration $`n`$ and the scattering time $`\tau `$ remain unchanged with $`y`$, or (ii) a change in $`n`$ is compensated by a change in $`\tau `$. To clarify which of the two is actually the case, we measured the room-temperature thermopower $`S(290\mathrm{K})`$, which is generally believed to reflect the change in the hole concentration and thus may be used as a guide to estimate $`n`$ . Figure 2 shows the $`y`$ dependence of the thermopower measured along the $`a`$-axis at 290 K, $`S_a(290\mathrm{K})`$, which is expected to be free from the contribution of the Cu-O chain transport. The $`S_a(290\mathrm{K})`$ data show a continuous change across the 60-K phase, which strongly suggests that $`n`$ is continuously changing with $`y`$ in our samples. Therefore, among the two possibilities listed above, it is more likely that a change in $`n`$ is somehow compensated by a change in $`\tau `$ in the 60-K phase.

Examination of the Hall channel in the in-plane transport gives us a clue to the (phenomenological) origin of the anomaly in $`\rho _a(T)`$. In YBCO, it is expected that $`\sigma _{xy}`$ is governed by the properties of the planes (since the Cu-O chains contribute little to $`\sigma _{xy}`$ because of their one-dimensionality), while the Hall resistivity $`\rho _{xy}`$ \[which is expressed as $`\rho _{xy}\sigma _{xy}/(\sigma _{xx}\sigma _{yy})`$ with $`\sigma _{xx}1/\rho _a`$ and $`\sigma _{yy}1/\rho _b`$\] is affected by the properties of the chains through $`\sigma _{yy}`$. Therefore, $`\sigma _{xy}`$ is a better indicator of the properties of the planes compared to $`R_H`$. Perhaps reflecting this situation, $`\sigma _{xy}`$ at 125 K shows a clearly anomalous $`y`$ dependence while such anomalies are smeared out in the $`y`$ dependence of $`R_H`$ \[Figs. 3(a) and 3(b)\]. The nature of this anomaly in the Hall channel is best understood by the plot of the Hall mobility in the planes, $`\mu _H=\sigma _{xy}/(B\sigma _{xx})`$ (Fig. 4), which reflects $`\tau /m^{}`$ and does not include $`n`$ in the Drude picture. One can clearly see that $`\mu _H`$ is anomalously enhanced near $`y`$=6.65, particularly at 125 K \[where the overlap of $`\rho _a(T)`$ is observed\]. Therefore, it appears that the scattering time $`\tau `$ gets enhanced upon reducing $`y`$ from 6.80 to 6.65 and this change in $`\tau `$ cancels the change in $`n`$, causing the $`y`$-independent resistivity in the CuO2 planes for $`y=6.656.80`$. It should be noted that this anomalous enhancement of $`\tau `$ takes place only when the pseudogap opens, and thus the overlap of $`\rho _a(T)`$ is observed only below $``$130 K.

The Hall effect reveals yet another aspect of the anomalies in the 60-K phase. Figure 5(a) shows the temperature dependences of $`R_H`$ (measured with the current along the $`a`$-axis) for a wide range of $`y`$; the data points in Fig. 3(a) are extracted from this series of data. One may notice in Fig. 5(a) that only the data for $`y=6.656.85`$ show noticeable drop in $`R_H(T)`$ with lowering temperature upon approaching $`T_c`$; this situation becomes obvious in Fig. 5(b), where the $`R_H(T)`$ data for $`y=7.0`$, 6.75, and 6.50 are compared, with clearly indicated $`T_c`$. It is apparent that only the 60-K-phase sample shows a marked drop in $`R_H(T)`$ from well above $`T_c`$ and the data seem to be heading towards zero. To demonstrate that this anomaly is not extrinsically caused by the effect of the Cu-O chains on $`R_H`$, we show similar comparison of $`\sigma _{xy}(T)`$ for the three compositions \[Fig. 5(c)\], which again shows that the anomalous $`T`$ dependence is observable (though a bit weakened) only in the 60-K phase sample. This remarkable drop in $`R_H`$ is somewhat reminiscent of the $`R_H(T)`$ behavior in Nd-doped La2-xSrxCuO4 , where it has been discussed that some peculiarity of the charge transport  imposed by the stripes  is responsible for the drop in $`R_H`$. Thus, the anomaly observed in the $`R_H(T)`$ behavior is suggestive of a role of the stripes in the 60-K phase of YBCO.

We also measured the magnetoresistance in $`\rho _a`$ in both the transverse and the longitudinal geometries, from which we extract orbital magnetoconductivity (MC). The orbital MC at high temperatures, whose temperature dependence is very well fitted with $`[bT^2+c]^2`$ , is due to the normal-state contribution; the data show additional contribution at lower temperatures, which can be identified as the Aslamazov-Larkin orbital (ALO) component of the superconducting fluctuation MC . The full details of the data and their analysis will be published elsewhere , and here we show only the fluctuation MC ($`\mathrm{\Delta }\sigma _{\mathrm{FL}}/\sigma `$) obtained after the analysis. Figure 6(a) shows how the fluctuation MC evolves in the $`(TT_c)/T_c`$ vs. $`y`$ plane. One can see that there is a general growth of $`\mathrm{\Delta }\sigma _{\mathrm{FL}}/\sigma `$ with decreasing $`y`$, but on top of this trend there is a marked enhancement of $`\mathrm{\Delta }\sigma _{\mathrm{FL}}/\sigma `$ near $`y6.7`$. An enhancement of $`\mathrm{\Delta }\sigma _{\mathrm{FL}}/\sigma `$ means that the characteristic magnetic-field scale to suppress the superconducting fluctuations is smaller, which suggests that $`H_{c2}`$ is reduced in the 60-K phase. In fact, when we extract the in-plane coherence length $`\xi _{ab}`$ by fitting the $`\mathrm{\Delta }\sigma _{\mathrm{FL}}/\sigma `$ data to the ALO formula, obtained $`\xi _{ab}`$ \[Fig. 6(b)\] shows an anomalous enhancement in its $`y`$ dependence near $`y6.7`$; this gives rise to a reduction of the mean-field upper critical field $`H_{c2}^{\mathrm{MF}}`$ \[$`\mathrm{\Phi }_0/(2\pi \xi _{ab}^2)`$\] as shown in Fig. 6(c).

## IV Stripes?

All the above results indicate that the electronic state in the CuO2 planes near $`y6.7`$ is somewhat anomalous compared to other composition and that the 60-K-phase anomalies are definitely not simply due to oxygen ordering. The origin of the electronic anomaly in the 60-K phase is not clear at this stage, but an intriguing possibility is the charge stripes, because the drop of $`R_H`$ with decreasing temperature observed in the 60-K-phase samples is most easily associated with the stripe physics . Below we propose a highly speculative scenario, which just describes one possibility: Suppose that the physics in the whole phase diagram is governed by the fluctuating stripes or “electronic liquid crystals” (as is proposed by Kivelson et al. ) and that the superconductivity is caused by the fluctuating stripes. If there is a quantum critical point (QCP) associated with the stripes near $`y6.7`$ , the superconductivity may be weakened at such QCP (which leads to a reduction in $`H_{c2}`$ and a plateau in $`T_c`$). Since the fluctuations are always enhanced near QCP, the charge mobility may also be enhanced. Thus, though highly speculative, it is possible to qualitatively understand the observed anomalies in the 60-K phase in terms of the electronic liquid crystal picture. In the above scenario, the drop in $`R_H(T)`$ should be due to the particle-hole symmetry  rather than the suppression of the transverse motion of charges .

## V Summary

We have systematically measured the transport properties of untwinned YBCO single crystals in a wide range of doping and found that the behaviors of $`\rho _a(T)`$, $`\mu _H(y)`$, $`R_H(T)`$, $`\mathrm{\Delta }\sigma _{\mathrm{FL}}(y)`$, and $`H_{c2}^{\mathrm{MF}}(y)`$ all show novel anomalies in the 60-K phase. These anomalies are clearly of electronic origin, and possibly related to the stripe physics in this compound.

## VI Acknowledgments

We thank S. A. Kivelson and A. N. Lavrov for helpful discussions.