text:
# The mass function of the Las Campanas loose groups of galaxies

## 1 Introduction

One important constraint of cosmological models is provided by the mass function of clusters and groups of galaxies. In the case of popular $`\mathrm{\Lambda }`$CDM models the amplitude and shape of the mass function depend primarily on the mean density of matter in the Universe, $`\mathrm{\Omega }_m=\mathrm{\Omega }_c+\mathrm{\Omega }_b`$, where $`\mathrm{\Omega }_c`$ and $`\mathrm{\Omega }_b`$ are the mean densities of the cold dark matter and baryonic matter in units of the critical cosmological density, respectively. The amplitude of the mass function also depends on the amplitude of the power spectrum of density fluctuations, which can be characterised by the $`\sigma _8`$ parameter (the linearly extrapolated rms density fluctuations in a sphere of of 8 $`h^1`$ Mpc radius). The local abundance of rich clusters of galaxies can be used to estimate $`\sigma _8`$ (Bahcall, Fan & Cen bac:fan (1997), de Theije, van Kampen & Slijkhuis the:the (1998), Cen cen:cen (1998)).

Determination of the cluster mass function consists of two tasks: calculation of cluster masses, and estimation of the spatial density of clusters. Masses of clusters can be derived using either X-ray data and the mass-temperature relation or data from optical surveys using the velocity dispersion of galaxies in clusters (virial masses). In a pioneering study by Bahcall and Cen (bac:cen (1993); hereafter BC93), both mass determination methods were used. Masses were derived for Abell clusters of all richness classes and for groups of galaxies. Biviano et al. (biv:biv (1993)) and Girardi et al. (gir:bor (1998), hereafter G98) calculated virial masses of nearby clusters. Reiprich & Böhringer (rei:rei (2002)) used a X-ray flux-limited sample of galaxy clusters (HIFLUGCS, the HIghest X-ray FLUx Galaxy Cluster Sample) to obtain a mass function. The spatial density of massive clusters of galaxies, according to G98, exceeds the density found by BC93 by a factor of almost ten. As shown by G98, this difference is mainly due to the fact that BC93 assumed a one-to-one correspondence between the richness of a cluster and its mass, whereas in reality this relationship has a large scatter. A new determination of the mass function of loose groups of galaxies has been provided by Girardi and Giuricin (gir:giu (2000)).

The sample of loose groups of galaxies identified in the Las Campanas Redshift Survey (LCRS) by Tucker et al. (tuc:tuc (2000), hereafter TUC) provides a possibility to derive a new independent estimate of the mass function of clusters and groups of galaxies. Loose groups represent the most numerous component of galaxy clustering, and so we should be able to determine their statistical properties much better than those for rare rich clusters of galaxies. The virial masses of loose groups of galaxies in the LCRS range from $`M10^{12}\mathrm{M}_{\mathrm{}}`$ to $`10^{15}\mathrm{M}_{\mathrm{}}`$; these masses were estimated by TUC using velocity dispersions of galaxies and harmonic radii of groups.

In this paper we shall estimate the mass function of LCRS loose groups. Our main task is to derive the spatial density of groups of various masses. The volume of the LCRS sample is well determined, thus we hope to get an unbiased estimate of the spatial density of groups. We shall compare the empirical mass function with theoretical mass functions found using numerical simulations of structure evolution. In our simulations we assign to a halo all the particles identified as members of the halo by halo finder. In observations this means that the velocity dispersion of galaxies is assumed to be identical to that of the dark matter. This comparison allows us to check the consistency of structure evolution models with empirical data and to find the set of cosmological parameters which brings models into agreement with data. We also determine the mass interval where the sample of groups is not influenced by selection effects.

## 2 The Data

### 2.1 Observations

The LCRS (Shectman et al. she:she (1996)) is an optically selected galaxy redshift survey where a multi–object spectrograph was used to measure simultaneously redshifts of 50 or 112 galaxies. Extending to a redshift of $`z0.2`$, the catalogue covers 6 slices of size on average $`1.5^{}\times 80^{}`$, containing 23,967 galaxies with measured redshifts within the official survey photometric and geometric limits. Three slices are located in the northern Galactic hemisphere, centred at declinations $`\delta =3^{},6^{},12^{}`$, and the other three slices are located in the southern Galactic hemisphere, centred at declinations $`\delta =39^{},42^{},45^{}`$. TUC applied a friends-of-friends (FoF) percolation algorithm to extract the catalogue of Las Campanas loose groups of galaxies (hereafter LCLG). The linking parameter was chosen to get a density enhancement limit of $`\delta n/n80`$. The minimum group membership was chosen to be three.

Because the spectroscopy was carried out for each field either via a 50 or a 112 fibre multiobject spectrograph, the selection criteria varied from field to field. The nominal apparent magnitude limits for 50 fibre fields were $`16.0R17.3`$ and for 112 fibre fields the limits were $`15.0R17.7`$. According to TUC great effort was put forth in accounting for these field-to-field sampling variations. The general properties of the 50 fibre and the 112 fibre groups agree well with properties of groups found from other surveys, thus we shall use the whole TUC group sample.

Only groups with redshifts $`10,000cz45,000`$ km s-1 were included into the LCLG sample. TUC concluded that the LCLG is a useful sample for a variety of studies requiring an unbiased collection of loose groups. It is based on the LCRS galaxy sample, which is the first redshift survey that can claim to enclose a reasonably fair sample of the nearby universe.

The complete LCLG list includes 1495 groups. TUC also introduced a “clean sample”, where groups with four potential bias factors are excluded: 1) groups which are too close to a slice edge, 2) groups which have the crossing time greater than the Hubble time, 3) groups with the corrected velocity dispersion zero what can happen since TUC subtracted an redshift measurement error of 67 km/s in quadrature from each group velocity dispersion, and 4) groups containing a 55 arc-sec orphan galaxy, i.e. a galaxy with no measured redshift. The last bias was caused by technical reasons (the fibre separation limit, which prevents the observation of neighbouring galaxies within 55 arc-seconds of each other). In the full sample this effect was compensated for by reintroducing lost galaxies, assigning to them a redshift equal to that of its nearest neighbour, convolved with a Gaussian of width $`\sigma =200`$ km s-1.

### 2.2 Simulations

A low density CDM universe with a cosmological constant ($`\mathrm{\Lambda }`$CDM) is widely regarded as the best model compatible with most of the currently available data; e.g. with the microwave background anisotropy measured by BOOMERANG (de Bernardis et al. ber:ber (2000)) and MAXIMA I (Hanany et al. han:han (2000)) experiments, and with data on the large scale structure of the universe: two-point correlation functions and power spectra, and the spatial density of mass-limited samples of galaxy clusters (Governato et al. gov:gov (1999), Colberg et al. col:col (2000), Pierpaoli et al. pie:sco (2001)).

For the present study we employ a flat cosmological model ($`\mathrm{\Omega }_m+\mathrm{\Omega }_\mathrm{\Lambda }=1`$) with the following parameters: the matter density $`\mathrm{\Omega }_m=0.3`$, the baryonic density $`\mathrm{\Omega }_b=0.04`$, the vacuum energy density (cosmological constant) $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, and the Hubble constant $`h=0.7`$ (here and throughout this paper $`h`$ is the present-day Hubble constant in units of 100 km s-1 Mpc-1). The simulations were performed using a P3M code (Couchman couch91 (1991)) for a cube of 200 $`h^1`$ Mpc size and a $`256^3`$ mesh and for the same number of particles; each particle has a mass of $`4.0\times 10^{10}h^1\mathrm{M}_{\mathrm{}}`$. The transfer function was computed using the CMBFAST code by Seljak and Zaldarriaga (sel:zal (1996)). The rms mass density fluctuation parameter of this model is $`\sigma _8=0.87`$.

We also calculated an additional model using the same code with $`128^3`$ particles in a cube of size $`L=100`$ $`h^1`$ Mpc, and the cosmological parameters $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`\mathrm{\Omega }_b=0.05`$, and $`h=0.65`$; this COBE normalised model has the density fluctuation parameter $`\sigma _8=0.78`$.

We used the FoF algorithm to identify CDM halos. The only free parameter in the FoF method is the linking length, which is defined as the maximum separation between particles which are still joined into groups. In our case the linking length was chosen as 0.23 in units of the mean particle separation, which approximately selects the matter density contrast $`\delta n/n=80`$. This density contrast was the one chosen by TUC to extract the group catalogue from LCRS, and it is typical of the values used in extracting groups from observational galaxy catalogues. While in the Einstein-de Sitter universe the overdensity within the virial radius of a cluster is $`\delta n/n=178`$, in a low density spatially flat universe with a cosmological constant it is $`\delta n/n178\mathrm{\Omega }_m^{0.6}`$ (White et al. whi:efs (1993)). The corresponding overdensity for our simulation is $`\delta n/n366`$. This means that the group sample determined using a low density contrast ($`\delta n/n=80`$) contains groups which can be in an uncertain dynamical state and do not have to be virialised. When particles outside the virialised core are included in the groups, the masses of simulated groups may be overestimated. For comparison we also calculated the mass function of publicly available numerical simulations of the Virgo consortium111http://www.MPA-Garching.MPG.DE/Virgo/, where Jenkins et al. jen:fre (1998) selected FoF-groups with a linking length of 0.2 of the mean inter-particle separation.

## 3 Results

### 3.1 Selection effects in the LCLG sample

At first we discuss selection effects and show how we took them into account when calculating the mass function. There are two main selection effects in the LCRS: 1) observations are performed in a fixed apparent magnitude interval which transforms to a distance dependent absolute magnitude interval, galaxies fainter or brighter than this interval are not included in the survey; 2) depending on the field 50 or 112 galaxies were measured for redshifts, the actual number of galaxies in the magnitude window could be larger, thus the sample is diluted.

The first selection effect makes it impossible for galaxies outside the window to enter the survey. Consequently groups consisting of faint galaxies occur only in the nearest region of the survey; with increasing distance fainter groups gradually disappear from the sample. This effect is seen in Figure 1 where the total luminosity of LCRS loose groups is shown as a function of distance. This selection effect can be statistically corrected using individual weights in calculation of the group mass function, following the procedure by Moore, Frenk & White (moo:moo (1993)) and Girardi & Guiricin (gir:giu (2000)). Namely, we weight each group by $`w=1/\mathrm{\Gamma }`$, where $`\mathrm{\Gamma }`$ is the volume, which corresponds to the absolute magnitude limit of the third brightest galaxy of a group; if $`\mathrm{\Gamma }>\mathrm{\Gamma }_0`$ (the full volume of a slice), we take $`\mathrm{\Gamma }=\mathrm{\Gamma }_0`$. The data to calculate these limits have been tabulated by TUC separately for each of the 327 spectroscopic fields of the LCRS. Figure 2 shows the distribution of the normalised weights $`w_0`$ (in units of the full volume of a slice) as a function of group mass. We see that there is a weak relation between cluster mass and the relative weight: massive groups contain as a rule sufficiently bright galaxies and can be observed in the whole volume of the survey, thus for massive clusters the relative weight is closer to unity.

Now we consider the second selection effect. To estimate the degree of dilution of the LCRS we compare the surface density of galaxies of the LCRS with the surface density of the Sloan Digital Sky Survey. Here the goal was to measure all galaxies in the magnitude window. The number of galaxies with measured redshifts in the Early Data Release of SDSS is 70 per square degree, while the LCRS has only 20 galaxies per square degree in the $`6^{}`$ slice, and $`3038`$ per square degree in the rest of slices. The faint end limit in most slices of the LCRS is almost same as for the SDSS ($`17.7`$), thus the difference in number density is due to dilution of the LCRS. Most slices of the LCRS are diluted by a factor of about 2, and the slice $`6^{}`$ by a factor of 3.5. Dilution decreases not only the number of galaxies but also the number of clusters and groups, as the number of galaxies in a group may fall below 3 and the group will not be included in the catalogue. The number of groups in the $`6^{}`$ slice is actually lower than in other slices by a factor of about 2.

However, uniform dilution does not influence strongly the mass function of groups. Firstly, virial masses of groups practically do not change. When determining virial masses, group members play the role of test particles, which move in a common gravitational potential. Using less test particles gives the same estimate for the potential and for the virial mass as before, only with a larger variance. We checked this conclusion by diluting simulated group catalogues.

Dilution may affect the group catalogue only by reducing the group richness of smaller groups below the catalogue limit. These groups will drop out of the catalogue and the number density of smaller (less massive) groups will fall. We estimated this correction, using group catalogues obtained in our numerical simulations and randomly diluting the groups by 50%, as estimated above. We found that the conditional probability distribution $`P(\mathrm{log}N|M)`$ ($`N`$ is the richness of a diluted group, $`M`$ is the virial mass of the group, and we use decimal logarithms here) is close to a Gaussian with a mean

$$\overline{\mathrm{log}N(M)}\mathrm{log}(M/\mathrm{M}_{\mathrm{}})11.5$$

and a rms error $`\sigma `$ that is also a function of the virial mass,

$$\sigma (M)0.3(14\mathrm{log}(M/\mathrm{M}_{\mathrm{}})).$$

The fraction $`\alpha `$ of the groups, which remain in the catalogue after dilution, is given then by the error integral:

$$\alpha (M)=0.5\text{erfc}\left(\frac{\mathrm{log}N_L\overline{\mathrm{log}N(M)}}{\sqrt{2}\sigma (M)}\right),$$

where $`N_L`$ is the catalogue limit (3 and 5, in our case) and $`\overline{\mathrm{log}N(M)}`$ and $`\sigma (M)`$ are given above. In order to restore the original differential mass function, we have to multiply it by the factor $`\alpha ^1`$. This dilution correction is larger than 1 only in the mass interval $`10^{12}`$$`10^{13}\mathrm{M}_{\mathrm{}}`$, reaching 1.95 at its lowest end (for $`N_L=3`$) and 2.67 for $`N_L=5`$. As the differential mass function is small in this mass interval, this correction does not change appreciably the final integral mass function. Also, the dilution correction does not change the high mass end of the mass function at all.

### 3.2 Mass function of LCLGs

The mass function (MF) of galaxy clusters/groups is defined as the number density of clusters above a given mass $`M`$, $`n(>M)`$. To construct the group mass function from a group sample one needs accurate group masses and well defined volume and selection functions of the sample. We have used the masses estimated by TUC, who assumed that the groups were virialised and calculated virial masses of groups as:

$$\mathrm{M}_{vir}=\frac{d\sigma _{los}^2R_h}{G},$$
(1)
where $`G`$ is the gravitational constant, $`R_h`$ is the harmonic radius of the group, $`\sigma _{los}`$ is the group line-of-sight velocity dispersion, and $`d=6`$ in the case of a spherically symmetric velocity distribution of the group.

We calculated the volumes of slices as follows:

$$V=\frac{\mathrm{cos}\delta _m\mathrm{\Delta }\alpha \mathrm{\Delta }\delta }{3}(r_2^3r_1^3);$$
(2)
here $`\delta _m`$ is the mean declination of the slice, $`\mathrm{\Delta }\delta `$ and $`\mathrm{\Delta }\alpha `$ are the widths of the slice in declination and right ascension (in radians), and $`r_1`$ and $`r_2`$ are the lower and upper distance used in calculations. This formula is valid for spatially flat (k=0) cosmologies.

The comoving cosmological distance $`r`$ is a function of measured redshift and depends on the cosmological model as (e.g. Peebles peebles (1993)):

$$r=\frac{c}{H_0}_0^z\frac{dz^{}}{E(z^{})},$$
(3)
where the function $`E(z)`$ is given by

$$E^2(z)=\mathrm{\Omega }_m(1+z)^3+(1\mathrm{\Omega }_m\mathrm{\Omega }_\mathrm{\Lambda })(1+z)^2+\mathrm{\Omega }_\mathrm{\Lambda }.$$
(4)
The LCLG survey extends up to the redshift $`z0.15`$. Tucker et al. (tuc:tuc (2000)) used the Einstein-de Sitter cosmology (assuming $`\mathrm{\Omega }_m=1,\mathrm{\Omega }_\mathrm{\Lambda }=0`$) to determine distances in the LCLG survey, while modern data prefer the parameters $`\mathrm{\Omega }_m`$=0.3 and $`\mathrm{\Omega }_\mathrm{\Lambda }`$=0.7. Thus for the present study we recalculated comoving distances for all galaxies of the LCLG survey. For the particular redshift range $`10000cz45000`$ km s-1 the ratio between the comoving volumes in these two models is about 0.8. The numbers of groups for all slices, for the clean sample and the volume of each slice in the LCLG catalogue are shown in Table 1. Only groups with defined mass are included (for some groups the velocity dispersion, corrected for random errors, is zero and mass determination is impossible).

The MF was constructed by sorting the masses in logarithmic bins of width 0.2. Errors in group masses (an estimate of the rms error in virial mass has been found by TUC) were accounted for by convolving group masses with a log-normal distribution of the same mode and variance. Figure 3 shows the mass function of LCLG for all six slices. To have a general picture we added slices together and found the average mass function, denoted by a thick solid line. We found that the average mass function and the mass function of the clean sample distributions are very similar, which indicates that the differences between the clean sample and the full sample are unimportant for the MF. We also tried excluding the slice $`\delta =6^{}`$ and found that it has only a weak influence on the total MF. In order to better show the differences between various slices we did not correct these mass functions for dilution. The dilution-corrected final mass function is given in Table 2 and in Fig. 7.

There is a clear difference in the MF between slices. The two northern slices at $`\delta =6^{}`$ and $`\delta =12^{}`$ differ from others, having the lowest and the highest amplitudes at the faint end of the MF. We suppose that this discrepancy between mass functions are due to differences in dilution, since the slice $`\delta =6^{}`$ was observed by the 50 fibre spectrograph only and the $`\delta =12^{}`$ slice by the 112 fibre spectrograph only. This discrepancy cannot be removed by a luminosity-weighting procedure. The remaining differences are probably due to cosmic variance. The slices of the LCRS survey intersect large voids and superclusters. In such slices, loose groups form a variety of structures that may influence the observed mass function at its massive end, since the most massive groups (clusters) are associated with supercluster structures.

Figure 3 shows clearly the flattening of mass functions at the low mass end. Such flattening is usually interpreted as the result of incompleteness of a sample (Girardi et al. gir:bor (1998)), although it could also be real, describing a diminishing number of groups of smaller masses. Thus, we have to study the completeness problem in more detail.

To estimate the mass completeness limit of the LCLG we can study the volume-density of the groups (within different mass intervals) as a function of distance. In Figure 4 lines denote mean values of the density in various mass ranges. Any rapid fall of the lines indicates incompleteness for that mass interval. We find from the Figure 4 that below masses of $`10^{13.5}10^{14}\mathrm{M}_{\mathrm{}}`$ our sample is probably incomplete.

Another way to estimate the mass completeness limit of the LCLG is to use the velocity dispersion distribution function, VDF (Fadda et al. fed:fed (1996)). The VDF is defined similarly to the MF, but for the velocity dispersion of groups. Assuming that the VDF can be described by a simple power law function (note that this assumption does not have a clear physical justification), Fadda et al. (fed:fed (1996)) set the completeness limit for their sample at the point where the power-law exponent of the VDF starts to change. For the LCLG this point is reached at the line-of-sight velocity dispersion $`\sigma _{los}250`$ km/s.

Figure 5 shows the relation between $`\sigma _{los}`$ and mass $`M`$ for the LCLG catalogue. Equation (1) defines the relation between these quantities via the harmonic radius of the group. The dotted line illustrates a power-law fit:

$$\sigma _{los}=a\left(\frac{hM}{10^{15}\mathrm{M}_{\mathrm{}}}\right)^c\mathrm{km}/\mathrm{s},$$
(5)
where $`a=819`$ and $`c=0.35`$. The result is in quite good agreement with previous results by del Popolo & Gambera (pop:popgam (2000)) and del Popolo et al. (pop:popetal (2000)), $`a=842\mathrm{km}/\mathrm{s}`$ and $`c=0.33`$, based on N-body simulations and X-ray observations. A value 1/3 for the exponent $`b`$ means that the system is virialised. The fact that the LCLG follow rather closely the virial relation shows that the possible contamination is small and the groups are physical. Setting the completeness limit to $`\sigma _{los}=250`$ km/s, based on the VDF, in equation (5), we get an estimate for the completeness limit in mass between $`10^{13}M10^{14}h^1\mathrm{M}_{\mathrm{}}`$, with the mean value about $`M5\times 10^{13}h^1\mathrm{M}_{\mathrm{}}`$. This result agrees with the result of the first analysis.

### 3.3 Reliability of groups

Ramella et al. (ram (1997)) and Girardi & Giuricin (gir:giu (2000)) noted that groups with a few (less than 5 members) may constitute not real, but pseudo-groups. On average, in pseudo-groups the velocity dispersion is considerably larger than in real groups (Ramella et al. ram (1997)). Indeed in the LCLG sample, groups with less than 5 galaxies have a very large scatter of velocity dispersions, and some velocity dispersions are extremely large. Analysis of a subsample, where poor groups with only 3 or 4 members are excluded, shows that in the $`wM`$ plane (as in Figure 2) the scatter for this subsample is somewhat smaller than for the whole sample. However, exclusion of all poor groups decreases the mass function considerably in the whole mass range (see Figure 7, where the mass function is given for two cases, for all groups and for groups with the number of members ($`N_m5`$). This indicates that at least part of small groups with a large velocity dispersion (and thus a large mass) are probably pseudo-groups (studies by Einasto et al. ein02 (2002) indicate that such groups are located mostly in low density regions). On the other hand, most of small groups have velocity dispersions in the same range as richer groups, and are probably real. In other words, exclusion of all small groups introduces an additional selection effect, which is difficult to quantify.

## 4 Discussion

### 4.1 Comparison with numerical simulations

Now we compare our mass function with mass functions found for groups in numerical simulations. Figure 6 shows the velocity dispersion of dark matter particles in groups as a function of the halo (group) mass in our numerical simulations. We see that in small groups with mass less than $`10^{13}h^1\mathrm{M}_{\mathrm{}}`$ (which contain less than $`250`$ particles) the scatter of the velocity dispersion is very large. This phenomenon can be explained if we assume that small groups are not yet virialised and are contaminated by non-virialised particles, or alternatively are composed of several small subgroups moving fast relative to each other. Another explanation for this effect is that it is a result of the chosen $`\delta n/n`$ — the FoF algorithm collects particles from outskirts of small clusters although they do not belong dynamically to the group. The fraction of interlopers in simulated loose groups can become rather high, about 20–30%, as found by Diaferio et al. (diaferio (1999)).

To test the second assumption we built groups, using a lower linking length $`b=0.17`$ (in units of the mean particle separation). This did not change considerably the scatter of the velocity dispersion among the groups with less than 250 particles. Thus the first explanation is more likely.

To eliminate the effect of unbound groups in simulations we calculated both the potential and kinetic energy for all groups and excluded unbound particles, using the condition $`|E_{pot}|<0.8E_{kin}`$ ($`E_{pot}`$ is the potential energy and $`E_{kin}`$ the kinetic energy of a group) for an unbound particle. We split groups with a high kinetic energy as compared to the potential energy to separate interacting or merging groups. To this aim we employed a 6D-group finder with a linking measure $`\mathrm{\Delta }r^2/r^2+\mathrm{\Delta }v^2/\sigma _v^2`$, where $`r`$ and $`\sigma _v`$ are the group effective (half-mass) radius and the velocity dispersion. This procedure affected only about 10% of the halos, with many of the halos having a small number of DM particles. We show the results of our clean-up procedure in the right panel of Figure 6, which is similar to the left panel, only non-virialised particles and bypassing groups are excluded. We see that our procedure to eliminate non-virialised particles works well — the low-mass groups follow the same trend as observed for the high-mass groups and clusters.

We have calculated the mass functions for the sample of all groups and for the sample of virialised groups in our simulations. These functions are very similar, only for groups of very small mass ($`M<10^{12}h^1\mathrm{M}_{\mathrm{}}`$) the MF of virialised groups differs from that of all groups by less than 1%. Thus the contribution of unbound groups to the total mass function is practically negligible.

Figure 7 shows the mass function of dark matter halos in our simulation together with the MF of the LCLG. The Poisson error bars (square root of counts per bin multiplied with the mean weight) for the full LCLG mass function are shown by the dash-dotted lines. In order not to overcrowd the figure, we have not shown the error bars for the $`N>=5`$ sample; they are very close to those of the full sample. For small masses ($`M<510^{13}\mathrm{M}_{\mathrm{}}`$) the error bars are smaller than the size of the plotting symbols. For comparison, Figure 7 shows also the mass function of the simulations by the Virgo Consortium, $`\mathrm{\Lambda }`$CDM-gif (Jenkins et al. jen:fre (1998), Kauffmann et al. 1999a ). In the Virgo simulation the side of the box was 141.3 $`h^1`$ Mpc, the matter density $`\mathrm{\Omega }_m=0.3`$, the cosmological term $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, the particle mass $`1.4\times 10^{10}\mathrm{M}_{\mathrm{}}/h`$, and the number of the particles $`256^3`$. Figure 7 shows that our $`128^3`$ simulation with the amplitude parameter $`\sigma _8=0.78`$ gives the best fit among the models we used to simulate the LCLG subsample at the massive end ($`210^{14}\mathrm{M}_{\mathrm{}}10^{15}\mathrm{M}_{\mathrm{}}`$) of the MF. For smaller masses, down to $`410^{13}\mathrm{M}_{\mathrm{}}`$ the model with a higher $`\sigma _8=0.87`$ gives a better fit. The cluster MFs in the Virgo Consortium simulation and in our $`256^3`$ simulation are quite similar. For yet smaller masses, where we suspect that the LCLG sample is incomplete, the LCLG mass function flattens out, in contrast to the rising simulated mass functions.

Most of earlier results give $`\sigma _80.91`$ for the $`\mathrm{\Omega }_m=0.3`$ scenario (Eke et al. eke (1996), Viana and Liddle viena (1996) etc.). This value agrees with the recent result by Pierpaoli et al. (pie:sco (2001)), $`\sigma _8=1.02`$. Lower values are obtained based on X-ray cluster samples. Reiprich & Böhringer (rei:rei (2002)) derived for $`\mathrm{\Omega }_m=0.3`$, $`\sigma _80.68`$. From the ROSAT Deep Cluster Survey Borgani et al. (borgani (2001)) derived $`\mathrm{\Omega }_m=0.35`$ and $`\sigma _80.66`$. Recently Lahav et al. (lah (2002)) obtained $`\sigma _8=0.73`$ using the 2dFGRS+CMB data. Finally, Bahcall et al. (bah (2002)) used the SDSS data and got $`\sigma _8=0.72`$ for $`\mathrm{\Omega }_m=0.3`$. In principle, $`\sigma _8`$ can be determined on the basis of the CMB data only; future experiments as MAP and PLANCK will help to pinpoint this parameter.

### 4.2 Comparison with earlier work

The mass function determined by Girardi & Giuricin (gir:giu (2000)) is shown in Figure 7. This group catalogue was constructed using a subsample of galaxies of the Lyon-Meudon Extragalactic Database. Groups were identified using two different methods: the FoF algorithm (P-groups, 453 groups) and the hierarchical method (H-groups, 498 groups). These group samples extend until $`cz5500`$ km s-1.

Here we used the P-data from figure 3 in Girardi & Giuricin (gir:giu (2000)). In that case they chose the subsample of the nearby groups with $`cz2000`$ km s-1. , which they assumed to be a complete sample and a good representation of the total population. However, this sample covers only the local supercluster, while the LCLG catalogue extends up to distances $``$ 450 $`h^1`$ Mpc. This volume contains already a variety of structures.

We see that in the high mass range ($`M10^{14}h^1\mathrm{M}_{\mathrm{}}`$) the MFs are rather similar. Around masses $`M510^{13}`$ the mass functions of the simulations and the LCLG start to differ from each other. Also, at these lower masses the mass function by Girardi & Giuricin (gir:giu (2000)) is considerably higher than the LCLG mass function.

This is an intriguing result that certainly needs further analysis. We see that although loose groups promise to provide better statistics, there are a number of unsolved problems that do not permit a realisation of this promise at the moment. There are two problems in the LCRS group sample: reliable determination of masses of small groups, and selection effects. Virial mass estimates can be in error; Girardi and Giuricin (gir:giu (2000)) try to correct for this. Their correction is rather complicated and could induce its own errors.

The simplest explanation of the difference between the mass functions found in the present paper and that of Girardi and Giuricin (gir:giu (2000)) is a mass-dependent selection effect. The LCRS is a strongly diluted sample of galaxies, and it is biased against low surface brightness galaxies. However, as shown above (Sec. 3.2), dilution practically does not affect the integral mass function.

The second effect, the absence of low-brightness galaxies is probably much less present in the Lyon-Meudon galaxy catalogue, from which Garcia (garcia (1998)) extracted the catalogue of local groups that served as a source for the Girardi and Giuricin (gir:giu (2000)) mass function. If low surface brightness galaxies were more prevalent in low mass groups – which is plausible – the LCLG mass function would flatten out as seen in Fig. 7.

In a recent paper Martínez et al. (martinez (2002)) found a similar flattening of the mass function of groups, extracted from the 2dF catalogue. This is also a deep sample, similar to the LCLG catalogue.

We should also look for possible problems with our simulated mass functions.

One problem is that we are comparing dark matter halos with galaxy-populated halos. It is probable that halos of smaller masses ($`10^{12}10^{13}`$) can not easily host three or more galaxies. This would explain the deficiency of small-mass galaxy groups. Thus, better simulations of group catalogues will have to include galaxy formation.

There could also be a problem with mass determination. We devoted much attention to accurately evaluate the mass estimates in simulations; in particular, we succeeded in filtering out unbound groups. Still, the mass-velocity dispersion relation of virialised groups widens toward the small mass end, showing that virial masses of small groups may have large intrinsic errors. This widening could be caused, first, by the predominance of radial orbits in small groups, but this should be a rather small effect. Another possibility is that small groups might have a higher intrinsic velocity spread than large groups due to different formation histories, or in other words, due to a mixture of young and old groups. This is a question that could be only answered by a detailed dynamical study of observed groups.

Another numerical effect could be an enhanced group production in simulations, caused by gravitational two-body collisions, which were discussed by Suisalu & Saar suisalu (1995) and Splinter et al. splinter (1998)). More recent analysis concerning the convergence of dissipationless dark matter $`N`$-body codes by Knebe et al. (knebe (2000)) and Power et al. (power (2002)) demonstrates that these effects are small if the gravitational softening and the time steps are chosen carefully.

## 5 Conclusion

In this work we have presented the mass function of the Las Campanas Loose Groups selected by Tucker et al. (tuc:tuc (2000)). We have also studied the mass-velocity dispersion relation that shows the expected scaling relation for virialised systems. We conclude that the loose galaxy groups are basically physical systems. The mass function was compared with results of dark matter simulations.

Our main conclusions are:

1) The completeness interval of the LCLG sample is found to be around $`510^{13}\mathrm{M}_{\mathrm{}}/hM8\times 10^{14}\mathrm{M}_{\mathrm{}}/h`$.

2) The high end of the mass function of the LCLG sample lies between the simulated mass functions corresponding to $`\sigma _8=0.78`$ and $`\sigma _8=0.87`$.

3) At small masses $`M<510^{13}\mathrm{M}_{\mathrm{}}/h`$ the mass function of the LCRS sample of loose groups flattens out.

###### Acknowledgements.

The present study was supported by the Estonian Science Foundation grant 4695. P.H. was supported by the Academy of Finland (grant 46733). D.L.T. was supported by the US Department of Energy under contract No. DE-AC02-76CH03000. We thank Mirt Gramann, Sahar Allam and Gert Hütsi for comments. We acknowledge the use of the AP3M code made public by H. Couchman. Finally, we thank our referee for a number of useful and constructive comments.


format: text

text:
# Force-induced unfolding of a homopolymer on fractal lattice: exact results vs. mean field predictions

## Abstract

We study the force-induced unfolding of a homopolymer on the three dimensional Sierpinski gasket. The polymer is subject to a contact energy between nearest neighbour sites not consecutive along the chain and to a stretching force. The hierarchical nature of the lattice we consider allows for an exact treatment which yields the phase diagram and the critical behaviour. We show that for this model mean field predictions are not correct, in particular in the exact phase diagram there is not a low temperature reentrance and we find that the force induced unfolding transition below the theta temperature is second order.

64.60.Cn, 87.15-v, 87.15.He, 05.10.Cc

The recent development of single molecule techniques have given experimentalists the opportunity to grab with suitable handles and mechanically manipulate, by means of optical tweezers or cantilevers such as atomic force microscopes, proteins, molecular motors or DNA molecules. In this way it has become possible to measure or exert on these molecules forces in the piconewton range. In particular the effect of a stress on the giant molecule titin has been studied in where it was found that a force induced unfolding transition takes place between a compact native like state and an extended state. The presence of strong hysteresis together with the rather abrupt jumps observed in the force vs. elongation curves suggest that the unfolding is a first order phase transition. On the theoretical side, many simple models have been proposed to make contact with these experiments (see e.g.  for the stretching of proteins and homopolymers,  for DNA unzipping and  for RNA unzipping). In particular, in Refs. the authors study the force-induced unfolding of a homopolymer and a heteropolymer in the mean field approximation and find that in both cases the critical line separating the globule from the coil is re-entrant at low temperatures. This is analogous to what had been found for the phase boundary in the phase diagram valid for the DNA unzipping, where it has been proved exactly that in the presence of a pulling force mean field is correct.

When self-avoidance is incorporated in the models, exact results are rare and mean field treatments become popular. In this work, we study exactly the force-induced unfolding of a self-avoiding walk (SAW) on a fractal lattice, the three-dimensional Sierpinski Gasket (SG) . This gasket has a fractal dimension $`D_f=2`$, and this renders the system interesting from the theoretical point of view because we are below the upper critical dimension for theta collapse. Mean field theory would predict a reentrance here because the hamiltonian compact walks in this lattice have nonzero entropy as compared to the zero entropy of the completely stretched coil (same argument as in ). However, the exact critical line that we find in the SG shows no reentrance. Similarly, by comparing the free energies of the globule and that of the coil one would naively expect a first order transition, whereas the exact calculation yields a second order transition. Though we cannot conclude from this calculation that also on Euclidean lattices the situation will be analogous, we feel that this calculation should give a warning that the mean field prediction need not be correct. The presence or not of reentrance on Euclidean lattices together with the nature of transition is thus an intriguing question which deserves further work.

We consider a SAW on the three-dimensional SG (3DSG), a hierarchical lattice with ramification number $`4`$ and $`D_f=2`$ (Fig. 1). We study the combined effect of a force $`\stackrel{}{f}`$ which stretches the polymer along one edge of the SG (Fig. 1) together with a compacting self-attractive term, obtaining by assigning a weight $`\mathrm{exp}(\beta ϵ)`$ ($`ϵ>0`$) every time two non consecutive sites of the SAW are nearest neighbour in the 3DSG. We call $`\beta T^1`$ the inverse temperature. To describe the effect of the stretching force, we give each step of the walk an orientation and a weight $`\mathrm{exp}(\beta f\mathrm{\Delta }a)`$, where $`f|\stackrel{}{f}|`$ is the modulus of $`\stackrel{}{f}`$ and $`\mathrm{\Delta }a`$ is the projection of the oriented step along $`\stackrel{}{f}`$.

The calculation of the partition function and other thermodynamic quantities involves the evaluation of twenty-five generating functions (Fig. 2). Twelve generating functions involve contributions of a SAW which enters in the $`n`$-th order SG at one vertex and goes out at one other vertex. Other twelve generating functions arise when an (oriented) SAW starts from a vertex, goes out from another one and then re-enters the $`n`$-th order SG at a later stage. The last generating function is the void. The one-leg generating functions are labelled $`A_i`$, $`i=1,\mathrm{},12`$, and the two-leg contributions $`B_i`$, $`i=1,\mathrm{},12`$. One can write recursion relations for the generating functions by exactly enumerating all walks in the SG on the computer: one has to sum up the contributions of SAWs at the $`n`$th order SG in order to generate SAWs in the $`n+1`$-th order SG.

The initial conditions for the generating functions are:

$`A_1`$ $`=`$ $`zy^1A_{2,3,4,8}=zy^{1/2}`$ (1)
$`A_{5,9,10,11}`$ $`=`$ $`zy^{1/2}A_{6,7}=zA_{12}=zy`$ (2)
$`B_{1,2,3,6}`$ $`=`$ $`z^2w^4y^1B_{4,5,8,9}=z^2w^4B_{7,10,11,12}=z^2w^4y`$ (3)
where $`z`$ is the SAW step fugacity, $`y\mathrm{exp}\left(\beta f\right)`$, and $`w\mathrm{exp}\left(\beta ϵ\right)`$, is the weight responsible for the theta collapse in absence of force.

Note that, as is apparent from the above equations, we have adopted the convention that interactions are restricted to sites within the first order SG and moreover that after one SAW has touched one vertex in the SG it is obliged to exit from the SG at that order. It has been proved that this approach is equivalent to the more general approach as regards universality of the phase transition .

When there is no pulling force the recursion relations simplify because:

$`A_iAB_iBi=1,\mathrm{},12,`$ (4)
so one gets:

$`A^{}`$ $`=`$ $`A^2+2A^3+2A^4+4A^3B+6A^2B^2,`$ (5)
$`B^{}`$ $`=`$ $`A^4+4A^3B+22B^4.`$ (6)
At non-zero $`\stackrel{}{f}`$, the equations retain the same structure but every $`B`$ and every $`A`$ has to be labelled by the appropriate number as found in the exact enumeration. To solve the model, we proceed as follows. First, we must find the fixed points of the recursion relations: in general, for every phase or critical point (line) there is a fixed point. Then, for every $`f`$ and $`T`$ fixed, we have to find the critical step fugacity, $`z_c(\beta ϵ,\beta f)`$ such that the flux defined by the twenty-five recursion relations and initial conditions will bring the system in the fixed point corresponding to those values of $`f`$ and $`T`$. This critical step fugacity $`z_c(\beta ϵ,\beta f)`$ allows us to have a numerical expression for all the quantities we are interested in, namely, the free energy $`F(\beta ϵ,\beta f)\frac{1}{\beta }\mathrm{log}(z_c(\beta ϵ,\beta f))`$, the average elongation along the direction of the pulling force, $`x(\beta ϵ,\beta f)`$ and the average number of contacts, $`n(\beta ϵ,\beta f)`$, found by calculating the appropriate derivative of $`F(\beta ϵ,\beta f)`$.

From the linearized recursion equations , it is possible to get the critical exponents of the phase transition. Every eigenvalue $`\lambda `$ (there are at most as many such eigenvalues as there are recursion relations) , of the linearized flux at given $`f`$ and $`T`$ defines two critical exponents, $`Y`$ and $`\nu `$, through:

$$\lambda 2^Y2^{\frac{1}{\nu }}$$
(7)
where $`\nu \frac{1}{Y}`$ is defined in terms of the average squared elongation $`x^2`$ as $`x^2N^{2\nu }`$ for large values of the number of steps in the SAW, $`N`$.

Performing this analysis for our model, we find six different fixed points. At zero force, we re-obtain the fixed points given in . In particular, for $`T<T_\theta \frac{2}{\mathrm{log}(3)}`$, the theta collapse temperature in the absence of force, there is a zero force compact fixed point for the recursion relations:

$`A_i=0B_i=22^{1/3}B^{}i=1,\mathrm{},12.`$ (8)
For $`T>T_\theta `$, i.e. in the swollen phase, the fixed point that is approached is:

$`A_i0.4294\mathrm{}B_i0.04998\mathrm{}i=1,\mathrm{},12.`$ (9)
Just at criticality at the theta temperature, $`T=T_\theta `$, the fixed point is:

$`A_i=1/3B_i=1/3i=1,\mathrm{},12.`$ (10)
The corresponding critical indeces $`\nu `$ (corresponding to the largest eigenvelues, $`\lambda _1`$, in every one of the three regimes) are $`1/2`$, $`0.5294\mathrm{}`$ and $`0.7294\mathrm{}`$ respectively in the collapsed phase, at the theta point and in the swollen phase.

At $`T<T_\theta `$, there exists a critical line $`f=f_c(T)`$ separating a compact from an open phase. When $`f<f_c(T)`$ the recursion relations display another compact fixed point:

$`A_i`$ $`=`$ $`0i=1,\mathrm{},12`$ (11)
$`B_{4,5,8,9}`$ $`=`$ $`22^{1/3}B_{1,2,3,6}=0B_{7,10,11,12}=+\mathrm{}`$ (12)
The divergence arises from the pulling force and one can convince that the two leg generating functions are diverging (vanishing) for $`L\mathrm{}`$ ($`L`$ is the system size, and one has $`L=2^{n1}`$ at the $`n`$-th level of iteration in the construction of the 3DSG) as $`y^L`$ ($`y^L`$). Consequently, we can eliminate the divergences if e.g. we multiply the diverging functions for $`B_1`$, and the vanishing ones for $`B_7`$, and then take the square root of the result: if we do this all the two leg generating functions converge to the fixed point $`22^{1/3}`$ as in the zero force compact phase. When the SAW is in the open phase ($`f>f_c(T)`$), the fixed point is:

$`A_i`$ $`=`$ $`0i=1,\mathrm{},11;A_{12}=1`$ (13)
$`B_i`$ $`=`$ $`0i=1,\mathrm{},12`$ (14)
Finally, when the force is exacly tuned at the critical line, the fixed point which is approached is:

$`A_i`$ $`=`$ $`0i=1,\mathrm{},11;A_{12}={\displaystyle \frac{1}{1+622^{2/3}}}0.5668\mathrm{}A^{}`$ (15)
$`B_{4,5,8,9}`$ $`=`$ $`22^{1/3}B_{1,2,3,6}=0B_{7,10,11,12}=+\mathrm{}`$ (16)
again the square root of the products of one diverging two leg function times a vanishing two leg function is $`22^{1/3}`$. We will see later that it is possible to justify the value found numerically also for the nonzero force system. The exponent $`\nu `$ is equal to $`1/2`$ for $`ff_c(T)`$ and is $`1`$ in the open phase.

In practice, by taking advantage of the knowledge of the fixed points, one can devise a convenient way to calculate numerically the critical force as a function of temperature. Once the critical step fugacity $`z_c(\beta ϵ,\beta f=0)`$, of the compact phase at a fixed temperature $`T`$, is known one can further fix $`z`$ to this value, and then tune $`f`$ in order to pass from the compact fixed point to the critical line fixed point (within the precision allowed by the computer). In order to have a sufficient accuracy in these calculations, it proved necessary to use a computer with quadruple precision.

The phase diagram obtained numerically is shown in Fig. 3. Two remarks are in order. First, we notice that there is no re-entrance in the critical line, the slope at $`T=0`$ being approximately $`0.05`$. This is at variance with the prediction doable on the basis of mean field-like treatments (as those done e.g. in Refs., see also below). Second, the behaviour of the critical line near $`T_\theta `$ is:

$$f_c(T)(T_\theta T)^{a=0.87\pm 0.01},$$
(17)
in agreement with the prediction $`a=\frac{\nu _\theta }{\varphi }=0.868`$, where $`\nu _\theta `$ is the critical exponent of the end-to-end distance at the theta temperature and $`\varphi `$ is the theta transition crossover exponent (see below for an argument, leading to the behaviour in Eq. 17).

We report in Fig. 4 the plot of $`x`$ vs. $`f`$ at $`T=0.35`$. This figure supports the hypotesis that the unfolding transition for the homopolymeric SAW on the 3DSG is second order. This is in agreement with the argument based on the recursion relations given below. Our belief supported by the exact numerics is that the transition is second order for any nonzero $`T`$ and is first order only at $`T=0`$.

We now present some arguments to interpret our results (based on the renormalization group (RG) flux). Even though to write explicitly the equations and to cope numerically with them it was necessary to put the force in step by step, it is equivalent to evolve the $`f=0`$ generating functions (Eq. 5) and then put in the dependence on $`\stackrel{}{f}`$ at every iteration by multiplying the generating functions times suitable powers of $`y`$. If we do this, we can exploit the symmetries of the problem which make so that before multiplying the generating functions times $`y`$ dependent term, all one leg and two leg diagrams are separately equivalent. By noting the structure of the nonzero force fixed points, one can argue that the one leg generating functions must vanish in the compact phase like $`ACy_c^L`$, where $`y_c\mathrm{exp}\left(\beta f_c(T)\right)`$ and $`C`$ is a constant, for now undetermined. By matching the exponentials in Eq.5, we obtain that the two equations, when $`TT_\theta `$ ($`T0`$), can be approximated by:

$`A^{}=A^2+6A^2B^2,B^{}=22B^4.`$ (18)
The flux corresponding to these RG equations is shown in Fig. 5. The non-trivial fixed point, $`(A^{},B^{})`$ with both $`A^{}`$ and $`B^{}`$ non-zero, is obtained with $`B^{}=22^{1/3}`$ and $`A^{}=\frac{1}{1+6\left(B^{}\right)^2}`$. Thus $`ACy_c^L`$ as $`L\mathrm{}`$, with $`C=A^{}`$. When $`T0`$, we notice that the fixed points can be approached at once, i.e. since the initial condition, also in the original (not simplified) recursion relations, because the terms that should evolve to $`0`$ have an initial value which is zero for $`T0`$. As a consequence, as $`T0`$, the phase boundary of Fig. 3 is found by matching the initial conditions with the fixed point values for $`A_{12}`$ and any $`B`$. Thus the critical line is found by solving this system:

$`z_c`$ $`=`$ $`\mathrm{exp}\left(2\beta ϵ\right)\left(B^{}\right)^{1/2},`$ (19)
$`z_c`$ $`=`$ $`\mathrm{exp}\left(\beta f\right)A^{},`$ (20)
where $`B^{}`$ and $`A^{}`$ are the nonzero fixed points of Eq. 18, corresponding to the critical force (such that $`y^LA1`$ for large $`L`$, see Eq.15). The critical line for $`T0`$ is:

$$f_c(T)2ϵ+T\mathrm{log}\left(\frac{A^{}}{(B^{})^{1/2}}\right)2ϵ0.0525\mathrm{}T,$$
(21)
so that $`f_c(T)`$ starts with negative slope as found numerically. For $`f>f_c(T)`$, the fixed point to be reached in the $`(A,B)`$ plane of Fig. 5 is $`(1,0)`$. We observe that for a small deviation from $`A^{}`$ the RG flux takes the SAW to the open fixed point by moving along the line $`B(A)`$, which we can find in the neighbourhood of $`(A^{},B^{})`$ by requiring that it be a fixed line under the flux defined by the recursion equations. We require that $`B(A)B^{}+c(\delta A)^\alpha `$ in the neighbourhood of $`A^{}`$. We thus need to solve the system:

$`A^{}`$ $`=`$ $`\left[A^{}+(\delta A)\right]^2\left[1+6\left[B^{}+c(\delta A)^\alpha \right]^2\right]`$ (22)
$`B(A^{})`$ $`=`$ $`22\left[B^{}+c(\delta A)^\alpha \right]^4`$ (23)
in the unknown quantities $`c`$ and $`\alpha `$. We obtain $`c=1/12A_{}^{}{}_{}{}^{3}B^{}1.2826\mathrm{}`$ and $`\alpha =2`$ so that the fixed line $`B(A)`$ approaches smoothly the fixed point at critical force. By inserting these values in the expression $`B(A)B^{}+c(\delta A)^\alpha `$ together with the ansatz:

$$z_c(\beta f_c,\beta ϵ)z_c(\beta f,\beta ϵ)(ff_c)^\gamma ,$$
(24)
one gets $`\gamma =2`$ which implies $`x(ff_c)`$ for $`f\stackrel{>}{}f_c`$ and the transition is second order. This argument is strictly valid for low $`T`$. However it is unlikely that the order of the transition could change along the phase boundary and this is also confirmed by our numerics. The situation at the point $`T=0`$ is somewhat special: the entropy vanishes and balancing the energetic terms gives a first order transition.

One should notice the importance of the term $`6A^2B^2`$. Let us think of a mathematical simplified model in which the relevant equations are:

$$A^{}=A^2,B^{}=22B^4,$$
(25)
in which the mixed term (physically due to stretched walks that still make an extensive number of contacts) is suppressed. This corresponds to balancing the free energy of a stretched coil with that of a compact globule. In this simplified treatment the re-entrance is present, as can be expected from naife estimates of the ground state entropies of the stretched and compact state, and the transition is first order.

Finally, we argue that the exponent $`a`$ in Eq.17 is given by $`a=\frac{\nu _\theta }{\varphi }`$ as anticipated above. Near $`T_\theta `$ at zero force the free energy behaves as $`F(T)(T_\theta T)^{2\alpha =\frac{1}{\varphi }}`$. For $`f\stackrel{>}{}0`$, on the other hand, at $`T=T_\theta `$, the free energy behaves as $`F(f)f^{\frac{1}{\nu _\theta }}`$. Consequently, we get $`a=\frac{\nu _\theta }{\varphi }`$ in Eq. 17. This is in agreement with the result found in $`d=3`$ in Ref. where the exponents both take their mean field values ($`1/2`$), and also with the mean field analysis in Ref. in which $`\varphi =1`$ and consequently $`a=\nu `$ ($`1/2`$ in the ideal case treated in ).

In conclusion, we presented an exact calculation of the phase diagram of a SAW in presence of a compacting contact energy and a stretching force. We deem it is interesting because it can be analyzed exactly. A mean field like treatment gives a reentrant boundary and a first order transition. Both these predictions are not confirmed by the exact treatment, which gives a critical line with a zero temperature negative slope and a second order transition. This has been explained in a simple way by analyzing a simplified version of the recursion relations, analogous to the real space RG equations. Whether or not the critical line in the hypercubic lattice shows reentrance is therefore not yet clear and appears to be an intriguing question. The 3DSG has $`D_f=2`$ and so the most natural comparison is with the Monte-Carlo simulations performed in Ref., which indeed gives a second order transition, even though we cannot be sure that the order of the transition is the same in the two-dimensional real and fractal lattices. The behaviour of the critical force near the theta point has been found and a general argument, valid also for hypercubic lattices, has been given in agreement with the result found here and also with the result found in $`d=3`$ in.

This work was supported by cofin2001.


text:
# A Laboratory Plasma Experiment for Studying Magnetic Dynamics of Accretion Discs and Jets

## 1 Introduction

The accretion disc occupies a leading role in astrophysics, figuring prominently in young stellar objects (YSO), binary star systems, and active galactic nuclei (AGN). An unsolved mystery is the origin of highly collimated bipolar jets and episodic flares associated with accretion discs. It was proposed some time ago that magnetic field dynamics can supply the necessary jet formation and collimation mechanisms (Blandford 1976; Lovelace 1976). Magnetohydrodynamic (MHD) simulations have shown that jets are a natural consequence of a rotating disc in the presence of a magnetic field (e.g. Shibata & Uchida 1985). Jet structure, such as “knots,” and also episodic behavior are observed in the simulations (Ouyed & Pudritz 1997; Goodson, Böhm & Winglee 1999; Nakamura, Uchida & Hirose 2001). The details of these models are unlikely to be tested by observations anytime in the near future. Therefore, data from laboratory experiments could be very useful for this purpose. It should be noted that astrophysical jets have been compared theoretically with plasma guns (Contopoulos 1995), and that plasma experimentalists have interpreted coaxial gun plasma flows in the context of astrophysical jet morphology (Caress 1996).

This work describes a new plasma gun based laboratory experiment and initial results which should give insight into the magnetic dynamics of accretion discs and jets. This experiment is the first to use a plasma gun which explicitly simulates the geometry and topology of a magnetically-linked star-disc system by using a co-planar disc-annulus electrode setup. The experiment reveals (1) the formation and helical structure of a magnetically-driven collimated plasma, similar to proposed models of astrophysical jets, and (2) plasma detachment which may be relevant for disc winds and flares and field-line opening in disc coronae. The resulting magnetic topology depends on the details of magnetic helicity injection, namely the force-free state eigenvalue $`\alpha `$ imposed at the boundary. The plasmas in this work satisfy the MHD criteria globally ($`S1`$, $`\rho _\mathrm{i}<<L`$, and $`V_\mathrm{A}<<c`$, where $`S`$ is the Lundquist number, $`\rho _\mathrm{i}`$ the ion gyro-radius, $`L`$ the plasma scale length, and $`V_\mathrm{A}`$ the Alfvén speed). It should be noted that $`S10^3`$ in the experiment, similar to that in MHD numerical simulations.

This paper is organized as follows. Section II gives a qualitative description of how the magnetic dynamics of an accretion disc system are modeled in the laboratory. Section III describes the experimental setup. Section IV presents initial experimental results, and Sec. V provides a discussion. Finally, a summary is given in Sec. VI.

## 2 Modeling disc magnetic dynamics in the laboratory

For non-relativistic accretion disc dynamics, such as in YSO’s and some binary star systems, global magnetic field evolution can be described using MHD. The ideal MHD Ohm’s law,

$$𝐄+𝐔\times 𝐁=0$$
(1)
implies that an electric field $`𝐄`$ will be induced due to the sheared toroidal rotation of the disc in the presence of a magnetic field $`𝐁`$. The $`E_\mathrm{R}`$ component creates an electric potential drop $`V`$ between the central object and the accreting disc, and therefore magnetic helicity $`K`$ will be injected into the disc magnetosphere at the rate $`\mathrm{d}K/\mathrm{d}t=2\psi V`$ (e.g. Bellan 2000), where $`\psi `$ is the poloidal magnetic flux emanating from the central object.

In the laboratory, a planar coaxial disc-annulus electrode setup simulates the topology and approximate geometry of the astrophysical star-disc system. However, rather than rotating the annulus to induce an $`E_\mathrm{R}`$ as in a real accretion disc, an electric potential $`V_{\mathrm{gun}}`$ is applied across the electrodes instead. The background poloidal $`𝐁`$, which corresponds to the field of the central object, is generated with an external coil. Again through Eq. (1), the background $`𝐁`$ and the applied voltage induce an $`𝐄\times 𝐁`$ toroidal rotation in the plasma, achieving a magnetic helicity injection rate $`\mathrm{d}K/\mathrm{d}t=2\psi _{\mathrm{gun}}V_{\mathrm{gun}}`$, where $`\psi _{\mathrm{gun}}`$ is the poloidal flux intercepting the inner disc electrode. This setup is commonly used to create spheromak plasmas (Rosenbluth & Bussac 1979; Jarboe 1994; Bellan 2000).

Therefore, the most basic magnetic dynamics near a star-disc system, namely the winding up of poloidal field lines generated by the central object by disc rotation, can be simulated in the laboratory by applying a voltage across coaxial electrodes in the presence of a background poloidal field. Of course there are limitations to the laboratory setup. For example, Keplerian sheared rotation and the existence and placement of a star-disc co-rotation radius cannot be exactly replicated. An experiment cannot hope to capture all the nuances and range of behaviors in all accretion disc systems. However, it can provide experimental data on the most fundamental plasma configurations which can result from a simple magnetically driven system.

## 3 Experimental setup

Figure 1 shows a side-view schematic of the experimental setup. A planar coaxial gun is installed on one end of a large vacuum chamber evacuated to $`1.5\times 10^7`$ Torr. The chamber is large compared to the plasma so as not to affect the plasma evolution. The gun setup, enlarged in Fig. 2, includes (1) an inner electrode consisting of a 20.3 cm diameter copper disc (attached to end of blue re-entrant port), (2) an outer electrode consisting of a 50.8 cm outer diameter copper annulus (in green), (3) an external solenoid (in red) to produce a poloidal bias magnetic flux $`\psi `$ linking the inner and outer electrodes, and (4) gas lines to deliver fast puffs of neutral hydrogen gas to the desired position of breakdown (adjacent to the electrodes). The gap between disc and annulus is 0.635 cm; plasma breakdown does not occur in the gap because the pressure-distance product there is too small to satisfy the Paschen breakdown criteria (Bellan 2000). A cylindrical coordinate system ($`R`$, $`\theta `$, $`Z`$) is utilized, with the origin at the center of the inner electrode and $`+Z`$ defined to be away from the electrode (toward the left in Fig. 2). There are eight gas injection holes on each electrode, distributed uniformly along the $`\theta `$ direction of the inner and outer electrodes. The bias field is characterized by the parameter $`\psi _{\mathrm{gun}}`$, which is the total initial bias flux intercepted by the inner electrode; sample contours of $`\psi `$ are shown in Fig. 2.

The plasma formation sequence is as follows. After the bias field and gas puffs are introduced, negative polarity high-voltage is applied to the inner electrode by discharging a 120$`\mu `$F capacitor bank through an ignitron; the outer electrode is maintained at vacuum chamber ground. In the presence of the high-voltage, the neutral hydrogen breaks down. The optimum path for breakdown is along vacuum field lines linking the inner and outer electrodes as shown in Fig. 2. The gun voltage $`V_{\mathrm{gun}}`$ drives a current $`I_{\mathrm{gun}}`$ between the electrodes, twisting up the purely poloidal vacuum field and producing a $`B_\theta `$. In an accretion disc, it is the disc rotation which achieves the same effect. The initally discrete flux tubes expand and eventually merge into an axisymmetric configuration. Throughout this process, magnetic helicity is being injected from the gun into the plasma.

Experimental results to date have centered around time-resolved global imaging of the plasma evolution. Images shown here were taken using a Cooke Corp. HSFC-PRO multiple-frame charge-coupled device (CCD) camera, which takes up to eight images per plasma discharge. The camera view is through the 20.3 cm window (labeled in Fig. 1) such that the gun electrodes appear on the right hand side of each frame. A false-color table is applied to the images for ease of viewing. Typically, the plasma is observed in unfiltered visible light. However, by using filters, it has been verified that most of the light emission is from neutral hydrogen line transitions. The exposure time of each frame is 20 ns and the interframe time is set typically to 1.5 $`\mu `$s, which is on the order of an Alfvén transit time. Additionally, a triple Langmuir probe is utilized to measure localized values of $`n_\mathrm{e}`$, $`T_\mathrm{e}`$, and floating potential. A Rogowski coil encircling the ceramic break (see Fig. 2) measures $`I_{\mathrm{gun}}`$, and a high-voltage probe measures $`V_{\mathrm{gun}}`$ at the mouth of the re-entrant port (shown in blue in Fig. 2). Typical experimental parameters are as follows: $`V_{\mathrm{gun}}4`$–6 kV, $`I_{\mathrm{gun}}70`$–130 kA, $`\psi _{\mathrm{gun}}0.5`$–2 mWb, $`B0.2`$–1 kG, $`T_\mathrm{e}T_\mathrm{i}5`$–20 eV, $`n_\mathrm{e}10^{14}`$ cm-3, and global $`\beta 2\mu _0nk(T_\mathrm{e}+T_\mathrm{i})/B^20.02`$$`0.1`$.

## 4 Initial Experimental results

Because the plasma is relatively low-$`\beta `$, it is reasonable as a first approximation to consider it as nearly force-free with plasma currents nearly parallel to the magnetic field, i.e.

$$\times 𝐁\alpha 𝐁.$$
(2)
Equation (2) implies that $`\alpha `$ is constant along field lines but not necessarily across them. Integrating Eq. (2) over the gun surface, it can be shown that the force-free state eigenvalue imposed by the gun is $`\alpha _{\mathrm{gun}}=\mu _0I_{\mathrm{gun}}/\psi _{\mathrm{gun}}.`$ Experimentally, both $`I_{\mathrm{gun}}`$ and $`\psi _{\mathrm{gun}}`$ can be adjusted to achieve a range of $`\alpha _{\mathrm{gun}}`$ values. Depending on the peak value of $`\alpha _{\mathrm{gun}}`$, several distinct plasma configurations are identified.

### 4.1 Collimated plasma: an analog for disc jets

For low values of $`\alpha _{\mathrm{gun}}`$, the formation of a collimated plasma is observed, as shown in Fig. 3. In this plasma, $`\alpha _{\mathrm{gun}}`$ peaks at approximately 66 m-1. In Fig. 3(a), gas breakdown has just occurred along eight discrete paths, each of which follows vacuum magnetic field lines and terminates at gas injection holes. In frame (b), the discrete arcs have expanded and begun coalescing. Magnetic reconnection is expected to occur along the $`Z`$-axis as the discrete flux tubes coalesce; the particularly intense light emission there may be due to higher plasma density resulting from compression. In frames (c)–(d), a central column forms and expands in the $`Z`$-direction. In frames (e)–(g), the central column extends in length and persists for several Alfvén transit times. In frame (h), the central column begins to break up.

The bright structures in the images are expected to correlate with magnetic topology and current flow; this was verified previously in a similar experiment using direct measurements of magnetic field (Yee & Bellan 2000). This correlation is consistent with (1) the light emission being predominantly from neutral hydrogen atoms excited by current-carrying electrons and (2) the fact that the plasma is low-$`\beta `$ and expected to be in a nearly force-free state in which $`𝐉`$ is nearly parallel to $`𝐁`$. It should be noted also that H–H+ charge-exchange time is estimated to be very fast ($`1`$ $`\mu `$s) in this experiment, and therefore light emission may be representative of plasma ion dynamics. Filamentary structures can be seen inside the column, suggesting complex field topologies within the globally collimated structure.

The cross section of the central column is clearly observable in Fig. 3. If there are no collimating forces, clearly the length of the column would not be as elongated as observed, nor would the cross section be so uniform along the $`Z`$-direction. This result verifies that the necessary magnetic structure and collimation for a disc jet can arise from magnetic forces associated with helicity injection.

### 4.2 Helical instability: an analog for jet structure

Higher $`\alpha _{\mathrm{gun}}`$ results in plasmas with unstable central columns, as shown in Fig. 4, in which $`\alpha _{\mathrm{gun}}`$ peaks at approximately 71 m-1. In Fig. 4(e), a highly nonlinear helical perturbation appears in the collimated plasma. By altering the CCD camera timing, it is possible to determine the characteristic growth time of the instability, which is about 1.7 $`\mu `$s, similar to the characteristic Alfvén time. The instability appears to be an ideal current-driven kink mode.

The Kruskal-Shafranov condition for ideal kink stability in a cylindrical column is (e.g. Bateman 1978)

$$q_{\mathrm{edge}}=2\pi aB_\mathrm{Z}/LB_\theta >1,$$
(3)
where $`a`$ is the column radius and $`L`$ the column length. Assuming that both the initial $`\psi _{\mathrm{gun}}`$ and the instantaneous $`I_{\mathrm{gun}}`$ are fully contained inside the plasma column, and using the definition of $`\alpha _{\mathrm{gun}}`$ and $`\psi _{\mathrm{gun}}`$, it is straightforward to show that the condition for stability is

$$\alpha _{\mathrm{gun}}L<4\pi .$$
(4)
Figure 5 plots $`\alpha _{\mathrm{gun}}`$ versus $`L`$ for a collection of central columns taken from both different discharges and different times from within the same discharge. Stable (triangles), marginally unstable (squares), and kinked (diamonds) columns are plotted in $`\alpha _{\mathrm{gun}}`$$`L`$ space, along with the Kruskal-Shafranov stability threshold (dashed line), showing good agreement between experiment and theory. This result indicates that plasma instabilities, in this case an ideal kink, can give rise to macroscopic structure within jets.

### 4.3 Detached plasma: an analog for disc flares

Still higher values of $`\alpha _{\mathrm{gun}}`$ result in detached plasmas, as shown in Fig. 6, in which $`\alpha _{\mathrm{gun}}`$ peaks at approximately 129 m-1. In this discharge, the plasma appears to detach from the electrodes in frame (d) and then propagates along the $`Z`$ direction at a speed of approximately $`6\times 10^4`$ m/s, a fraction of the estimated $`V_\mathrm{A}`$. It is likely that a spheromak configuration is formed here; this was verified previously with direct measurement of $`𝐁`$ on an experiment with a non-planar source but similar helicity injection (Yee & Bellan 2000). This result supports the idea of field line reconnection above accretion discs, which can lead to disc winds and also episodic high energy flaring.

## 5 Discussion

The results above show three distinct plasma configurations having accretion disc characteristics. All three configurations result from the same plasma formation process, the only difference being the peak value of the parameter $`\alpha _{\mathrm{gun}}=\mu _0I_{\mathrm{gun}}/\psi _{\mathrm{gun}}`$. Figure 7 illustrates this dependence by placing different plasmas into $`I_{\mathrm{gun}}`$$`\psi _{\mathrm{gun}}`$ parameter space, with detachment at larger $`\alpha _{\mathrm{gun}}`$, attached columns at lower $`\alpha _{\mathrm{gun}}`$, and kinked columns near $`\alpha _{\mathrm{crit}}60`$–70 m-1.

This $`\alpha _{\mathrm{gun}}`$ dependence has many implications. Most importantly, it suggests that the plasma configurations associated with accretion discs and jets are related to Taylor relaxation theory (Taylor 1986), which is a description of how a plasma evolves as magnetic energy is minimized subject to the constraint of constant magnetic helicity. This process can be cast as a variational problem, and it can be shown that the resulting magnetic field configuration satisfies Eq. (2) with uniform $`\alpha `$. In the case of a driven plasma like an astrophysical jet, $`\alpha `$ will be non-zero, and magnetic helicity will flow from regions of high to low $`\alpha `$, which tends to “relax” the plasma toward uniform $`\alpha `$ (Bellan 2000). The evolution of many laboratory plasmas, including spheromaks, can be understood in terms of Taylor relaxation. For example, it has been shown that a threshold $`\alpha _{\mathrm{crit}}`$ at the source must be exceeded in order to form the closed-field configuration of a spheromak (e.g. Yee & Bellan 2000). This property is not surprising since analytic solutions of Eq. (2) show that $`𝐁`$ transitions from sinh-like to sine functions above a critical $`\alpha `$ (Bellan 2000).

Well-known models of astrophysical jets have considered separately the roles played by poloidal (Blandford & Payne 1982) and toroidal fields (Contopoulos 1995) in jet formation. However, the present work shows that jet structure should result when the correct value of $`\alpha `$ occurs at the disc boundary, resulting in a non-arbitrary combination of poloidal and toroidal fields in the disc corona. Likewise, deformation of the jet structure will occur under different choice of $`\alpha `$. In the experiment, the observed instability is consistent with an ideal current-driven kink. This instability has been observed in numerical simulations of AGN radio jets (Nakamura, Uchida & Hirose 2001) and proposed as the responsible mechanism for wiggled structure in the jets.

## 6 Summary

A plasma gun based laboratory experiment approximates the boundary conditions and topology of a star-disc system, and it is shown experimentally that magnetic helicity injection with these boundary conditions leads naturally to both collimated plasmas and detached plasmas, suggestive of disc jets and flares, respectively. The onset of a helical instability in the plasma column is shown to be consistent with the Kruskal-Shafranov condition for an ideal current-driven kink. The magnetic topology depends on the force-free state eigenvalue $`\alpha _{\mathrm{gun}}`$ imposed at the plasma gun. It is argued that Taylor relaxation provides a useful description of the magnetic structures of accretion discs and jets. These results demonstrate experimentally that the concept of magnetically-driven jets is a viable one. More quantitative characterization of the observed plasmas is planned. This will include direct measurements of the magnetic field via insertable probes, ion flow velocity along the collimated plasma via Doppler spectroscopy, and plasma pressure profiles via Langmuir probes.

## acknowledgments

The authors are grateful to J. Hansen, S. Pracko, C. Romero-Talamás, and F. Cosso for technical assistance. This work was supported by a U.S. Dept. of Energy (US-DoE) Fusion Energy Postdoctoral Fellowship and US-DoE contract no. DE-FG-03-98ER54461.