text:
# THE H𝛼 LUMINOSITY FUNCTION OF THE GALAXY CLUSTER ABELL 521 AT 𝑧 = 0.25Based on data collected at the Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.

## 1 INTRODUCTION

Recent observations of clusters of galaxies at intermediate-redshift ($`0.1z0.5`$) demonstrated that environmental processes play an important role on both the star-formation activity and the morphology of galaxies (e.g., Dressler et al. 1997; Balogh et al. 1998; Poggianti et al. 1999; Moss & Whittle 2000). The galaxy evolution in such clusters may be explained as the result of the increased activity of the field galaxies at more distant redshift (Lilly et al. 1996), modulated by changing the infall rate onto the clusters (Bower 1991). However, little information on the global star formation history in such intermediate-redshift clusters has been obtained even to date.

In many previous spectroscopic studies such as the Canadian Network for Observational Cosmology (CNOC) and MORPHS surveys (e.g., Abraham et al. 1996; Poggianti et al. 1999), the \[O ii\]$`\lambda 3727`$ emission line has often been used to investigate the star formation activity in cluster galaxies. However, the H$`\alpha `$ emission line provides a more reliable indicator of star formation because it is less sensitive both to dust extinction and to metallicty effects (Kennicutt 1998). Recently, an H$`\alpha `$ survey has been made for nearby clusters of galaxies Abell 1367 and Coma (Iglesias-Páramo et al. 2002). These data are utilized to derive the H$`\alpha `$ luminosity function (LF) of the cluster member galaxies. Optical spectroscopic surveys are also useful in studying H$`\alpha `$ emission-line properties of cluster member galaxies (Balogh et al. 2002; Couch et al. 2001). However, such spectroscopic sample are selected by using optical broad-band color properties, being independent from their H$`\alpha `$ emission-line properties. Therefore, these data sets may miss a part of H$`\alpha `$-emitting galaxies that are too faint to be selected in a continuum-magnitude-limited sample. In order to determine the shape of the H$`\alpha `$ LF unambiguously, it is required to carry out a deeper H$`\alpha `$ imaging survey.

Motivated by this, we performed our deep H$`\alpha `$ imaging survey of a cluster of galaxies, Abell 521 at $`z=0.25`$ which is a relatively rich (richness class 1) cluster; note that its Bautz-Morgan Type is III (Abell, Corwin, & Olowin 1989). Recent studies of this cluster showed that the projected galaxy density distribution in the central region has a very anisotropic morphology, as it exhibits two high-density filaments crossing in an X-shaped structure at the barycenter of the cluster (Arnaud et al. 2000; Maurogordato et al. 2000; Ferrari et al. 2003). The observed high velocity dispersion of the cluster ($``$ 1325 km s-1) may be attributed to the presence of those substructures. Since this velocity dispersion is also high compared to that expected from the temperature of the X-ray gas ($`kT=6.3`$ keV, Arnaud et al. 2000), it is suggested that this cluster is undergoing strong dynamical evolution driven by mergers among substructures. In this paper, we present observational properties of H$`\alpha `$-emitting galaxies in Abell 521 and then derive the H$`\alpha `$ LF for the first time. Throughout this paper, magnitudes are given in the AB system (Fukugita, Shimasaku, & Ichikawa 1995). We adopt $`H_0=50`$ km s-1 Mpc-1, $`q_0=0.5`$, and $`\mathrm{\Lambda }=0`$ cosmology, for comparison with previous results in the literature.

## 2 OBSERVATIONS AND DATA REDUCTION

The imaging data of Abell 521 were obtained with Suprime-Cam (Miyazaki et al. 2002) on the 8.2 m Subaru Telescope atop Mauna Kea during two runs (2001 October and 2002 February). We used the narrow-band filter, $`NB816`$ ($`\lambda _0=8150`$Å and $`\mathrm{\Delta }\lambda =120`$Å); the central wavelength corresponds to a redshift of $`0.24`$ for H$`\alpha `$ emission. We also used broad-band filters, $`B`$, $`V`$, $`R_\mathrm{C}`$, $`i^{}`$, and $`z^{}`$. A journal of our observations is summarized in Table 1.

We used IRAF and the mosaic-CCD data reduction software, NEKOSOFT (Yagi et al. 2002), to reduce and combine the individual CCD data. Photometric and spectrophotometric standard stars used for the flux calibration are as follows; SA98 (Landolt 1992) for the $`B`$ data, SA95 and SA98 (Landolt 1992) for the $`V`$ data, SA95, SA113, and PG0231 (Landolt 1992) for the $`R_\mathrm{C}`$ data, GD50 (Oke 1990) for the $`i^{}`$ data, and GD50, GD108, and PG1034+001 (Stone 1996) for the $`NB816`$ data. Since the zero point of Landolt (1992) is based on Vega, we convert the $`B`$, $`V`$, and $`R_\mathrm{C}`$ magnitudes to the AB system assuming monotonic colors of $`B(\mathrm{Vega})B(\mathrm{AB})=0.140`$, $`V(\mathrm{Vega})V(\mathrm{AB})=0.019`$, and $`R_\mathrm{C}(\mathrm{Vega})R_\mathrm{C}(\mathrm{AB})=0.169`$ (Fukugita et al. 1995). For $`i^{}`$ and $`NB816`$, we multiplied the atmospheric, instrumental, and filter transmittance, mirror reflectivity and CCD quantum efficiency with the SED of standard stars to calculated absolute flux. The $`z^{}`$ data were calibrated by using the magnitude of the quasar SDSSp J104433.04–012502.2 (Fan et al. 2000), which were obtained in the same observing night (see Fujita et al. 2003). Since any quasar is a potentially variable object, the photometric calibration of the $`z^{}`$ data may be more unreliable than those of the other band data. In order to estimate the uncertainty, we compared the colors of stars in our area with Galactic stars given by Gunn & Stryker111Bruzual-Persson-Gunn-Stryker (BPGS) Library (ftp://ftp.stsci.edu/cdbs/cdbs2/grid/bpgs) and applied the correction to our $`z^{}`$ data. The combined images for the individual bands were aligned and smoothed with Gaussian kernels to match the PSF ($`=1\mathrm{}.1`$). The final images cover an area of 32 $`\times `$ 20 arcmin2, which corresponds to 9.4 $`\times `$ 5.8 Mpc2 at $`z=0.25`$.

Source detection and photometry were performed by using SExtractor version 2.2.2 (Bertin & Arnouts 1996). The $`NB816`$ image was used as a reference image for detection, and photometry was performed for each band image. We used a criterion that a source has at least a 13-pixel connection at $`2\sigma `$ level. For each object, we measured the magnitude within a fixed aperture of a $`3\mathrm{}`$ ($``$ 15 kpc at $`z=0.25`$) diameter for colors and the magnitudes in all 6 band images. However, we used an adaptive aperture of diameter $`2.5r_\mathrm{K}`$ for measurement of total magnitudes, where $`r_\mathrm{K}`$ is the Kron radius (Kron 1980). Our $`3\sigma _{\mathrm{rms}}`$ within a $`3\mathrm{}\varphi `$ aperture in each band is given in Table 1.

In order to correct for the Galactic extinction, we calculated the extinction at each band (Cardelli et al. 1989), assuming that $`A_B=0.324`$ (Schlegel et al. 1998, NED222NASA/IPAC Extragalactic Database http://nedwww.ipac.caltech.edu/). The extinction at each band, $`A_\lambda `$, is given in Table 1.

## 3 SELECTION OF H$`\alpha `$ EMITTERS

We defined the off-band continuum flux of objects as $`i^{}z^{}0.6i^{}+0.4z^{}`$, which was computed from a linear interpolation between the effective frequencies of the $`i^{}`$ and $`z^{}`$ bands. Taking the scatter in the $`i^{}z^{}NB816`$ color into account, emission-line objects were selected with the following criteria,

$`i^{}z^{}NB816`$ $`>`$ $`\mathrm{max}(0.1,3\sigma \mathrm{error}\mathrm{of}i^{}z^{}NB816),`$ (1)
$`i^{}z^{}`$ $`<`$ $`24.9.`$ (2)
The color of $`i^{}z^{}NB816=0.1`$ corresponds to an equivalent width of $`12`$Å in the rest frame at $`z=0.247`$. The magnitude limit corresponds to the $`3\sigma _{\mathrm{rms}}`$ of the $`i^{}z^{}`$ band. These criteria are shown by the solid, dashed, and dotted lines in Figure 1. Then we found 324 objects that satisfy the above criteria; note that they are brighter than the limiting magnitude at each band.

A narrow band imaging can potentially detect galaxies with different emission lines at different redshifts. Emission lines detected in our narrow-band imaging survey are, for example, H$`\alpha `$ ($`z0.24`$), \[O iii\]$`\lambda \lambda 4959,5007`$ ($`z0.63`$), H$`\beta `$ ($`z0.68`$), and \[O ii\]$`\lambda 3727`$ ($`z1.18`$). In order to distinguish H$`\alpha `$ emitters at $`z0.24`$ from other emission-line objects, we investigate their broad-band color properties. In Figure 2, we show the diagram between $`BV`$ and $`Vi^{}`$ for the 324 emission-line objects. We also show the model galaxy colors that estimated by using the population synthesis model GISSEL96 (Bruzual & Charlot 1993). The star formation rate of model galaxies is proportional to $`\mathrm{exp}(t/\tau )`$, where we use $`\tau =1`$ Gyr. Changing the age of model galaxies, we generate the SEDs of five types of galaxies \[$`t=`$ 1, 2, 3, 4, and 8 Gyr for SB (young starburst), Irr, Scd, Sbc, and E, respectively\]. Then we find that H$`\alpha `$ emitters at $`z0.24`$ can be selected by adopting the following color criterion,

$$BV>0.75(Vi^{})+0.25.$$
(3)
As shown in Figure 2, this criterion appears to isolate H$`\alpha `$ emitters from other emitters at different redshifts. We then obtain a sample of 165 H$`\alpha `$ emitters.

In order to examine whether or not this color criterion works well more quantitatively, we investigate color properties of the spectroscopically confirmed member galaxies in Abell 521. The previous spectroscopic surveys identified 125 member galaxies in this cluster (Ferrari et al. 2003; Maurogordato et al. 2000). These members are identified in our catalog. We found that 123 objects satisfy the color criterion based on our photometry. The remaining two objects were rejected in our selection procedure because they have neighboring galaxies, causing large photometric errors. As a further check, we examined the spatial distribution of the galaxies which are selected by the above color criterion. In Figure 3, the excess of surface number density by the cluster members is clearly seen in the area within the radial distance of 2 Mpc from the center of Abell 521; note that the center is taken at the position of the main X-ray cluster (Arnaud et al. 2000).

In order to investigate the detection completeness of H$`\alpha `$ emitters, we performed a detection test on simulated artificial images. We prepared model galaxies with broad-band colors of late type galaxies (SB, Irr, Scd, and Sbc) at $`z=0.24`$ generated by GISSEL96 model. Since the equivalent width of H$`\alpha `$ emission is generally different form galaxy to galaxy, we put the excess of $`NB816`$ to $`i^{}z^{}`$ continuum as a parameter. We studied a color range of $`0.1<i^{}z^{}NB816<1.1`$ with five bins of $`\mathrm{\Delta }(i^{}z^{}NB816)=0.2`$. Since the total magnitude of $`NB816`$ is another parameter, we studied a magnitude range of $`18.0<NB816<24.0`$ with twelve bins of $`\mathrm{\Delta }(NB816)=0.5`$. We prepared 500 model galaxies in each set of $`NB816`$ and $`i^{}z^{}NB816`$. The light distribution of model galaxies is generated as the exponential law by IRAF ARTDATA. Their sky positions, half-light radius (1 to 7 kpc), and ellipticities (0.3 to 1.0) were randomly set and put into the CCD data together with Poisson noises. After smoothing model-galaxy images to match to the seeing size, we tried to detect them using SExtractor and to select as H$`\alpha `$ emitters with the same procedure as that used for the observed data. The detectability of the model emitters is shown in Figure 4 as a function of $`NB816`$ magnitude. We found that the detectability is higher than 80% for objects with $`NB816<22.0`$, but lower than 30% for objects with $`23.0<NB816<24.0`$.

Next, we examined the contamination by emission-line objects at other redshifts (e.g., \[O iii\]$`\lambda \lambda 4959,5007`$ at $`z0.63`$) in the same manner as that for the detection completeness estimate. We found the contamination to be only less than 5% even for the faint objects with $`23.0<NB816<24.0`$. In conclusion, our color criterion allows us to select reliable H$`\alpha `$ emitters at $`z0.24`$ for sources with $`NB816<23.0`$.

## 4 H$`\alpha `$ LUMINOSITY FUNCTION

### 4.1 Emission-Line Equivalent Width

At first we derive the emission-line equivalent width of the individual H$`\alpha `$ emitters. Since the pass band of our narrow-band filter is too wide to separate \[N ii\]$`\lambda \lambda 6548,6584`$ from H$`\alpha `$, the equivalent with derived here is that contributed from the three emission lines. We adopt the same method as that used by Pascual et al. (2001) to calculate H$`\alpha `$ \+ \[N ii\]$`\lambda \lambda 6548,6584`$ (hereafter H$`\alpha `$ \+ \[N ii\]) equivalent width. The flux density in each filter can be expressed as the sum of the line flux and the continuum density (the line is covered by both filters),

$`f_{NB}`$ $`=`$ $`f_\mathrm{C}+(F_\mathrm{L}/\mathrm{\Delta }NB),`$ (4)
$`f_{i^{}z^{}}`$ $`=`$ $`f_\mathrm{C}+0.6(F_\mathrm{L}/\mathrm{\Delta }i^{}),`$ (5)
where $`f_\mathrm{C}`$ is the continuum flux density, $`F_\mathrm{L}`$ the line flux, $`\mathrm{\Delta }NB`$ and $`\mathrm{\Delta }i^{}`$ the $`NB816`$ and $`i^{}`$ band filter effective widths, respectively, and $`f_{NB}`$ and $`f_{i^{}z^{}}`$ the flux density in each filter. Then the equivalent width can be expressed as follows;

$$EW_{\mathrm{obs}}(\mathrm{H}\alpha +[\text{N II}])=\frac{F_\mathrm{L}}{f_\mathrm{C}}=\mathrm{\Delta }NB\left[\frac{f_{NB}f_{i^{}z^{}}}{f_{i^{}z^{}}0.6f_{NB}(\mathrm{\Delta }NB/\mathrm{\Delta }i^{})}\right],$$
(6)
where $`\mathrm{\Delta }i^{}=1535`$Å for the $`i^{}`$ band and $`\mathrm{\Delta }NB=120`$Å for the $`NB816`$ band.

### 4.2 H$`\alpha `$ Luminosity

In order to obtain the H$`\alpha `$ luminosity, we have to correct for the flux contribution from the \[N ii\] doublet. We also have to apply correction for internal extinction. For the two corrections, we adopt the flux ratio of H$`\alpha `$/\[N ii\] = 2.3 (Kennicutt 1992) and $`A_{\mathrm{H}\alpha }=1`$ mag as a typical internal extinction at $`z0.2`$ (Tresse & Maddox 1998). Applying these corrections, we estimate the H$`\alpha `$ luminosity, $`L(\mathrm{H}\alpha )`$, from the total magnitude for each objects. We then obtain the star formation rate, $`SFR`$, using the following relation (Kennicutt 1998),

$$SFR(M_{}\mathrm{yr}^1)=7.9\times 10^{42}L(\mathrm{H}\alpha )(\mathrm{erg}\mathrm{s}^1).$$
(7)
The net H$`\alpha `$ flux depends on the source redshift because of the fixed transmittance of our $`NB816`$ filter. In order to correct for this effect, we adopt the following assumptions. (1) The velocity (redshift) distribution of the cluster galaxies is a Gaussian probability density function $`P(z)`$ around a mean velocity of 74019 km s-1 with a dispersion of 1325 km s-1 (Ferrari et al. 2003). (2) The equivalent width distribution of the H$`\alpha `$ emitters in the rest frame is expressed as the following exponential functional form; $`Q(EW)=A\mathrm{exp}(EW/B)`$, where $`A`$ and $`B`$ are free parameters. Then, the rest frame equivalent width probability of the H$`\alpha `$ emitters at redshift $`z`$ can be expressed as $`P(z)\times Q(EW)`$. The best-fit parameters are obtained so as to fit the distribution of $`i^{}z^{}NB816`$.

We apply the correction for the incomplete velocity coverage of the $`NB816`$ filter. Since the number of the H$`\alpha `$ emitters excluded from the filter can be calculated from the value of the equivalent width, we correct the equivalent width distribution function that derived above. Therefore, using the equivalent width distribution function, we obtain the corrected H$`\alpha `$ luminosity distribution.

### 4.3 Statistical Weights

Since the excess of surface galaxy density by the cluster member is clearly seen in the area within 2 Mpc radius (see Figure 3), it seems safe to use the H$`\alpha `$ emitters located in the region; note that there are 53 H$`\alpha `$ emitters in this region. Although the $`NB816`$ filter corresponds to the redshifted H$`\alpha `$ emission line at $`z=0.242\pm 0.009`$, the mean redshift of the cluster is $`z=0.247`$. Therefore, our H$`\alpha `$ emitter sample may suffer contamination from field H$`\alpha `$ emitters, in particular from field galaxies in the foreground. In order to correct for this effect, we estimate the field contamination by the following method.

In order to obtain a sample of field galaxies, we use an area of 3 Mpc outside from the cluster center, taking account that the virial radius is estimated to be 2.74 Mpc (Girardi & Mezzetti 2001). However, since a subclump is seen to the North-East at about 4.5 Mpc from the cluster center, we exclude this area from our analysis. The remaining sky area is used to estimate the surface number density of field galaxies as a function of $`NB816`$ magnitude. Since Abell 521 shows elongated galaxy distribution (Ferrari et al. 2003), we use the local surface number density, which is calculated from 10 nearest neighbors, rather than the radial surface number density distribution. In this way, all H$`\alpha `$ emitters in the cluster region are weighted by the luminosity-dependent factor, $`w_\mathrm{l}`$ and the density-dependent factor, $`w_\mathrm{d}`$, after renormalizing $`w_\mathrm{d}`$ so that the sum of the total galaxy weights, $`w_\mathrm{l}\times w_\mathrm{d}`$ is equal to the sum of $`w_\mathrm{l}`$.

### 4.4 H$`\alpha `$ Luminosity Function

We derive the H$`\alpha `$ luminosity function (LF) for the cluster region. In order to normalize the H$`\alpha `$ LF to a proper volume, we use the sphere volume within 2 Mpc radius. The H$`\alpha `$ LF for our sample is well fitted by the Schechter (1976) function;

$$\varphi (L)\mathrm{d}L=\varphi ^{}(L/L^{})^\alpha \mathrm{exp}(L/L^{})\mathrm{d}(L/L^{}).$$
(8)
The size of the bins is taken to be $`\mathrm{\Delta }\mathrm{log}L(\mathrm{H}\alpha )=0.5`$. The number counts in the individual bins are given in Table 2, and Figure 5 shows the H$`\alpha `$ LF for Abell 521 and its Schechter fit. After correcting for the detection completeness (see section 3), we obtain the best-fit parameters for the cluster membership within 2 Mpc radius are

$`\alpha `$ $`=`$ $`0.75\pm 0.23,`$
$`\varphi ^{}`$ $`=`$ $`10^{0.25\pm 0.20}\mathrm{Mpc}^3,`$
$`L^{}`$ $`=`$ $`10^{42.03\pm 0.17}\mathrm{erg}\mathrm{s}^1.`$
The inset panel in Figure 5 shows the $`\chi ^2`$ error contours corresponding to the 1 and 2$`\sigma `$ levels. It can be seen that $`\alpha `$ and $`L^{}`$ values are highly correlated.

## 5 DISCUSSION

In Figure 5, we also show the H$`\alpha `$ LFs of field galaxies in the local Universe (Gallego et al. 1995), at $`z0.2`$ (Tresse & Maddox 1998), at $`z0.24`$ (Fujita et al. 2003), and those of the nearby clusters of galaxies (Iglesias-Páramo et al 2002). As for the nearby clusters of galaxies, Iglesias-Páramo et al. (2002) derived the H$`\alpha `$ LFs for Abell 1367 and Coma (richness class 2). The physical areas analysed by them are nearly the same as that of Abell 521. Although they applied a slightly different extinction correction; $`A_{\mathrm{H}\alpha }=1.1`$ mag for type Scd or earlier and $`A_{\mathrm{H}\alpha }=0.6`$ for type Sd or later (Boselli et al. 2001), they obtained $`\alpha 0.7`$ and $`L^{}10^{41.25}`$ erg s-1 for these clusters. The inset panel in this figure shows the best-fit parameters of Abell 521 and the others. There is the significant difference in H$`\alpha `$ LFs between the clusters and the fields. It is found that the faint end slope, $`\alpha `$, of the H$`\alpha `$ LF of Abell 521 is flatter than that of the field H$`\alpha `$ LFs at the same redshift. Therefore, although the sample analysed here is not so large, it appears that the flatter faint end slope of the H$`\alpha `$ LF is commonly seen in the clusters of galaxies between $`z0`$ and $`z0.2`$.

The most important finding in this study is that the characteristic luminosity, $`L^{}`$ of Abell 521 is $``$ 6 times (or $``$ 2 mag) brighter than those of the local cluster H$`\alpha `$ LFs, and nearly the same as those of the field H$`\alpha `$ LFs. The evolution of the field H$`\alpha `$ LF from $`z0.2`$ to $`z0`$ appears closely related to the $`\varphi ^{}`$ parameter rather than to the $`L^{}`$ and $`\alpha `$ parameters. However, as for the cluster H$`\alpha `$ LFs, the $`L^{}`$ parameter appears to evolve strongly from $`z0.2`$ to $`z0`$. Furthermore, we investigate the following galaxy clusters for which spectroscopic H$`\alpha `$ emitter surveys were already made; Abell 1689 at $`z=0.18`$ (Balogh et al. 2002) and AC 114 at $`z=0.31`$ (Couch et al. 2001). Since their observations are based on spectroscopy, it is difficult to compare with our H$`\alpha `$ LF straightforwardly. Therefore, we compare the H$`\alpha `$ LF of Abell 521 with the other clusters, only using the H$`\alpha `$ luminous galaxies ($`L^{}10^{41.5}`$ erg s-1). It should be noted that the same dust extinction correction (i.e., $`A_{\mathrm{H}\alpha }=1`$ mag) is applied for these clusters, and the H$`\alpha `$ luminosity from spectroscopy may be underestimated due to aperture bias.

Figure 6 shows the surface density of H$`\alpha `$ luminosity distribution of each cluster; note that the surface density is used instead of the volume density which is used in Figure 5, to compare simply with previous results. Although there are some galaxies whose $`SFR`$s exceed $`SFR4M_{}`$ yr-1 in Abell 521 and Abell 1689, no such galaxies are found in Abell 1367 and Coma. This $`SFR`$ is about the same level as that observed in the Milky Way (Rana 1991). This result implies that Abell 521 and Abell 1689 contain more active star-forming galaxies than the nearby clusters. However, AC 114 has nearly the same property in the H$`\alpha `$ LF as those of the nearby clusters in spite of its higher redshift.

This evolution may simply reflect the Butcher-Oemler (BO) effect, which is the increase in the fraction of blue galaxies, presumably star-forming galaxies, in clusters with increasing redshift (Butcher & Oemler 1984). In fact, the fraction of blue galaxies in our H$`\alpha `$ emitters by the definition of Butcher & Oemler (1984) is $``$ 85%. The blue fraction, $`f_\mathrm{B}`$, in Abell 1689 is comparable to that of Abell 521 ($`f_\mathrm{B}0.17`$, see for details Appendix) and much higher than those in the nearby rich clusters such as Coma ($`f_\mathrm{B}=0.03`$). However, AC 114 and Abell 1367 also have a large fraction of blue galaxies ($`f_\mathrm{B}0.20`$). Therefore, the star formation activity in clusters of galaxies may not be understood solely by the BO effect.

These results suggest that some physical processes work in their star formation activity, independent of the blue fraction of each cluster. We investigated what parameter of cluster correlates with the existence of high $`SFR`$ galaxies. One possible parameter is the dynamical status, as pointed out by Balogh et al. (2002). Although Abell 1689 appears to show a round shape in X-ray, the observed radial velocity distributions of the cluster shows substructures. On the other hand, AC 114 appears to be relaxed dynamically, although the cluster member galaxies shows an elongated morphology (Balogh et al. 2002). Therefore, it seems important to investigate the radial velocity distributions for the H$`\alpha `$ emitters of each cluster (Figure 7). In this figures, we use the relative radial velociy, $`\mathrm{\Delta }v_\mathrm{r}`$, to the systemic velocity, $`v_{\mathrm{r},0}`$, of each cluster, and $`v_{\mathrm{r},0}`$ is shown in the left corners of each figure. For Abell 521, Abell 1367, and Coma, their redshift are unknown from their imaging survey, since the H$`\alpha `$ emitters are detected using the narrow-band filters. Therefore, their distributions are shown using available spectroscopic data in the literature (Ferrari et al. 2003; Iglesias-Páramo et al. 2002; Cortese et al. 2003). The vertical dashed lines in this figures show the velocity coverage corresponding to the $`FWHM`$ of the narrow-band filters. For Abell 521, in particular, since some H$`\alpha `$ emitters with $`\mathrm{\Delta }v_\mathrm{r}>0`$ km s-1 may be excluded from the filter transmittance, we also show the velocity distribution of galaxies identified with other emission lines (\[O II\] or H$`\beta `$) taken from Ferrari et al. (2003). It is found that this radial velocity distribution is consistent with that of the H$`\alpha `$ emitters with $`\mathrm{\Delta }v_\mathrm{r}<0`$ km s-1. Therefore, we expect that it also traces that of the H$`\alpha `$ emitters with $`\mathrm{\Delta }v_\mathrm{r}>0`$ km s-1. In this figure, the radial velocity distribution for the H$`\alpha `$ emitters in Abell 521 is inconsistent with a Gaussian distribution which we assumed in section 4.1. However, since the spectroscopic sample is rather small, it is dangerous to apply this distribution to our sample. Even if this distribution is valid, it does not change the result that $`L^{}`$ of Abell 521 is much brighter than those of the nearby clusters.

In Figure 7, the radial velocity distributions in Abell 1367 and Coma tend to be concentrated around $`\mathrm{\Delta }v_\mathrm{r}0`$ km s-1. Therefore, it is suggested that these galaxies are associated with a galaxy population in the cluster core, probably a virialized system. For AC 114, the radial velocity distribution suggests that the most H$`\alpha `$ emitters are associated with the structure centered at $`\mathrm{\Delta }v_\mathrm{r}+2000`$ km s-1. Even if they comprise a substructure infalling on to the cluster center, the star-forming activity seems to have already been ceasing by the environment effect (Kodama et al. 2001; Gómez et al. 2003). Abell 521 and Abell 1689 show that the radial velocity distributions show several peaks in a wide velocity range, and there are few H$`\alpha `$ emitters at the mean cluster redshifts. And for Abell 521, they are not thought to be field contamination but be associated with the cluster, since their number density is high at the center of the cluster (Figure 3). The H$`\alpha `$ emitters in the two clusters are likely to be in the process of accretion of the field population. If this is the case, we may explain why their star-forming activities are higher than those of the other clusters. It is thus suggested that the star formation activity in clusters seems to be strongly related to the dynamical status of the clusters.

We would like to thank the Suprime-Cam team and the Subaru Telescope staff for their invaluable help, and T. Hayashino for his technical help. We also an anonymous referee for his/her useful comments and suggestions. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. The data reduction is partly done using computer system at the Astronomical Data Analysis Center of the National Astronomical Observatory of Japan. This work was financially supported in part by the Ministry of Education, Culture, Sports, Science, and Technology (Nos. 1004405, 10304013, and 15340059). TN and MA are JSPS fellows.

## Appendix A COMMENTS ON THE BLUE FRACTION

In order to derive the fraction of blue galaxies, $`f_\mathrm{B}`$, in Abell 521, we followed the procedure of Butcher & Oemler (1984). We calculated the fraction of galaxies that are 0.2 mag bluer than the color-magnitude sequence of early-type galaxies in the rest-frame $`BV`$ color within a magnitude limit, $`M_V=20`$, and that are enclosed within $`R_{30}`$, the radius which contains 30% of the cluster population.

The blue fraction, $`f_\mathrm{B}`$, was calculated by subtracting the number count expected for the field from that obtained for the cluster sample. In order to obtain the value of $`f_\mathrm{B}`$, we used the $`VR_\mathrm{C}`$ color index, which is very close to $`BV`$ rest frame at the cluster redshift. The magnitude limit, $`M_V=20`$, in the rest frame corresponds to $`R_\mathrm{C}20.8`$ at cluster redshift, and a difference of $`\mathrm{\Delta }(BV)=0.2`$ mag is $`\mathrm{\Delta }(VR_\mathrm{C})0.17`$ mag. The value of $`R_{30}`$ is determined from the surface number density profile (Figure 3). We found $`R_{30}=3.1`$ arcmin ($`0.9`$ Mpc), being very similar to the typical radius of $`R_{30}1`$ Mpc for a sample of rich clusters at $`0.2z0.4`$ (Kodama & Bower 2001). As a result, we obtain $`f_\mathrm{B}=0.17\pm 0.03`$.


format: text

text:
# Density functional theory and demixing of binary hard rod-polymer mixtures

## Abstract

A density functional theory for a mixture of hard rods and polymers modeled as chains built of hard tangent spheres is proposed by combining the functional due to Yu and Wu for the polymer mixtures \[J. Chem. Phys. 117, 2368 (2002)\] with the Schmidt’s functional \[Phys. Rev. E 63, 50201 (2001)\] for rod-sphere mixtures. As a simple application of the functional, the demixing transition into polymer-rich and rod-rich phases is examined. When the chain length increases, the phase boundary broadens and the critical packing fraction decreases. The shift of the critical point of a demixing transition is most noticeable for short chains.

The addition of the nonadsorbing polymer to a mono-disperse suspension of colloidal particles can lead to a phase separation due to depletion interactions Asakura54  arising from a tendency of the system to reduce the volume excluded to the polymer coils. One of the simplest theoretical models taking into account this phenomenon is the so-called Asakura-Oosawa (AO) model of colloid-polymer mixtures Vrij76  where the ideal polymer coils (modeled as spheres) can freely interpenetrate each other but the polymer-colloid and colloid-colloid interactions are of the hard sphere type. Initial studies of such systems focused on the bulk phase behavior Gast83 ; Lekkerkerker92 , however recently developed density functional theory (DFT) for the AO model Schmidt00  initiated investigations of inhomogeneous colloid-polymer mixtures. When brought close to a hard wall, such mixtures may develop a sequence of layering transitions in the partial wetting regime prior to a transition to complete wetting Brader02a ; Dijkstra02 .

Similar mechanism of fluid-fluid phase separation can be found if other mesoscopic particles such as hard rods are used as depletant agents. Bolhuis and Frenkel (BF) Bolhuis94  used computer simulations and free volume theory Lekkerkerker92  to study bulk phase behavior of mixtures of colloidal hard spheres and vanishingly thin hard rods. They found a surprisingly good agreement (c.f. Fig. 3 in Bolhuis94 ) between the free volume theory and Gibbs ensemble Monte Carlo simulations. Although the vanishing thickness and volume constitutes a significant simplification (e.g. it rules out the isotropic-nematic transition), the BF model is relevant for some experimental systems e.g. mixtures of silica spheres and silica coated bohemite rods Koenderink99 ; Vliegenthart99 . In order to describe inhomogeneous hard rod-hard sphere systems Schmidt Schmidt01  proposed a density functional theory which in its structure closely resembles the AO functional Schmidt00 . The functional for the BF model incorporates the exact low-density limit and yields the same equation of state as in Ref. Bolhuis94 . Moreover, entropic surface phase transitions found previously in model colloid-polymer mixtures close to a hard wall were also recently encountered in hard rod-hard sphere mixtures Roth03 . This further demonstrates deep similarities between the two models.

The aim of the present work is to construct a functional for a mixture of vanishingly thin hard rods and polymers. Such systems can be regarded as simple microscopic models of the liquid crystal-polymer mixtures. The functional is constructed by combining the Schmidt’s functional for the BF model with the Yu and Wu Yu02  functional for mixtures of polymeric fluids. To make this conjecture we take the advantage of the fact that both functionals underlie the fundamental measure theory (FMT) of Rosenfeld Rosenfeld89 . As a simple application we investigate bulk phase diagrams resulting from the proposed theory.

Consider a mixture of hard, vanishingly thin needles (species $`N`$) of length $`L`$ and polymers (species $`P`$) modeled as chains composed from $`M`$ tangentially bonded hard-sphere segments of diameter $`\sigma `$. The hard-sphere monomers building up the chains are freely jointed i.e. they can adopt any configuration as long as it is free of the intermolecular and intramolecular overlap. The interaction potential between needles $`V_{NN}=0`$ for all separations, while the pair potential between a polymer segment and a hard rod, $`V_{PN}`$, and between two polymer segments, $`V_{PP}`$, is of a hard-core type i.e. is infinite if a pair of objects overlap and zero otherwise. The grand potential of such system as a functional of local densities of polymers $`\rho _P(𝐑)`$ and needles $`\rho _N(𝐫,𝝎)`$ can be written as

$`\mathrm{\Omega }[\rho _P(𝐑),\rho _N(𝐫,𝝎)]=F_{id}[\rho _P(𝐑),\rho _N(𝐫,𝝎)]+`$ (1)
$`F_{ex}[\rho _P(𝐑),\rho _N(𝐫,𝝎)]+{\displaystyle 𝑑𝐑\rho _P(𝐑)(V_{ext}^{(P)}(𝐑)\mu _P)}`$
$`+`$ $`{\displaystyle 𝑑𝐫\frac{d𝝎}{4\pi }\rho _N(𝐫,𝝎)(V_{ext}^{(N)}(𝐫,𝝎)\mu _N)},`$
where $`𝐑(𝐫_1,𝐫_2,\mathrm{},𝐫_M)`$ denotes a set of coordinates describing the segment positions, $`𝝎`$ describes the orientation of the rod, $`V_{ext}^{(P)}(𝐑)`$, $`\mu _P`$, $`V_{ext}^{(N)}(𝐫,𝝎)`$ and $`\mu _N`$ are the external and the chemical potentials for polymers and rods, respectively. The ideal part of the free energy is known exactly

$`\beta F_{id}[\rho _P(𝐑),\rho _N(𝐫,𝝎)]=\beta {\displaystyle 𝑑𝐑\rho _P(𝐑)V_b(𝐑)}+`$ (2)
$`{\displaystyle 𝑑𝐑\rho _P(𝐑)[\mathrm{ln}(\mathrm{\Lambda }_P^3\rho _P(𝐑))1]}+`$
$`{\displaystyle 𝑑𝐫\frac{d𝝎}{4\pi }\rho _N(𝐫,𝝎)[\mathrm{ln}(\mathrm{\Lambda }_N^3\rho _N(𝐫,𝝎))1]},`$
where $`\mathrm{\Lambda }_P`$ and $`\mathrm{\Lambda }_N`$ are the thermal wavelengths for polymers and needles, respectively. The total bonding potential $`V_b(𝐑)`$ is a sum of bonding potentials $`v_b`$ between the segments $`V_b(𝐑)=_{i=1}^{M1}v_b(|𝐫_{i+1}𝐫_i|)`$ and for the tangential hard spheres this contribution satisfies $`\mathrm{exp}[\beta V_b(𝐑)]=_{i=1}^{M1}\delta (|𝐫_{i+1}𝐫_i|\sigma )/4\pi \sigma ^2`$.

Within the framework of the fundamental measure theory the excess free energy density $`\mathrm{\Phi }`$ is expressed as a simple function of the weighted densities $`n_\alpha ^{(i)}`$. We assume that $`\mathrm{\Phi }`$ can be represented as a sum of the orientation-independent polymer contribution, $`\mathrm{\Phi }_{POL}`$, and the orientation-dependent polymer-needle contribution, $`\mathrm{\Phi }_{PN}`$. For the orientation-independent contribution we use an extension of the Rosenfeld’s FMT to the polymeric fluids Yu02 , where the polymer excess free energy density $`\mathrm{\Phi }_{POL}`$ is assumed to be a functional of only segment densities $`\rho _{PS}(𝐫)`$ defined as

$$\rho _{PS}(𝐫)=\underset{i=1}{\overset{M}{}}\rho _{PS,i}(𝐫)=\underset{i=1}{\overset{M}{}}𝑑𝐑\delta (𝐫𝐫_i)\rho _P(𝐑),$$
(3)
where $`\rho _{PS,i}(𝐫)`$ is the local density of the segment $`i`$ of the polymer. Following Yu02  we assume that $`\mathrm{\Phi }_{POL}`$ can be split into the hard-sphere $`\mathrm{\Phi }_{HS}`$ contribution resulting from the hard-sphere repulsion of polymer segments and $`\mathrm{\Phi }_P`$ contribution arising due to the chain connectivity. Taking into account the assumptions mentioned above we have $`\mathrm{\Phi }=\mathrm{\Phi }_{PN}+\mathrm{\Phi }_{HS}+\mathrm{\Phi }_P`$. The total excess free energy of the inhomogeneous system is obtained via the integration of the excess free energy density

$$\beta F_{ex}=𝑑𝐫\frac{d𝝎}{4\pi }\mathrm{\Phi }_{PN}(\{n_\alpha ^{(i)}\})+\mathrm{\Phi }_{HS}(\{n_\alpha ^{(P)}\})+\mathrm{\Phi }_P(\{n_\alpha ^{(P)}\}).$$
(4)
For the hard-sphere part $`\mathrm{\Phi }_{HS}`$ we choose the original Rosenfeld expression Rosenfeld89

$`\mathrm{\Phi }_{HS}(\{n_\alpha ^{(P)}\})`$ $`=`$ $`n_0^{(P)}\mathrm{ln}(1n_3^{(P)})+`$ (5)
$`{\displaystyle \frac{n_1^{(P)}n_2^{(P)}𝒏_{V1}^{(P)}𝒏_{V2}^{(P)}}{1n_3^{(P)}}}`$ $`+`$ $`{\displaystyle \frac{(n_2^{(P)})^33n_2^{(P)}𝒏_{V2}^{(P)}𝒏_{V2}^{(P)}}{24\pi (1n_3^{(P)})^2}}.`$
The polymer weighted densities $`n_\alpha ^{(P)}(𝐫)`$ are evaluated via spatial convolutions

$$n_\alpha ^{(P)}(𝐫)=𝑑𝐫^{}\rho _{PS}(𝐫^{})w_\alpha ^{(P)}(𝐫𝐫^{}),$$
(6)
where the weight functions $`w_\alpha ^{(P)}(𝐫)`$ are given by

$`w_3^{(P)}(𝐫)=\mathrm{\Theta }({\displaystyle \frac{\sigma }{2}}|𝐫|)`$ , $`w_2^{(P)}(𝐫)=\delta ({\displaystyle \frac{\sigma }{2}}|𝐫|),`$ (7)
$`𝒘_{V2}^{(P)}(𝐫)={\displaystyle \frac{𝐫}{|𝐫|}}\delta ({\displaystyle \frac{\sigma }{2}}|𝐫|)`$ , $`w_1^{(P)}(𝐫)={\displaystyle \frac{w_2^{(P)}(𝐫)}{2\pi \sigma }},`$ (8)
$`w_0^{(P)}(𝐫)={\displaystyle \frac{w_2^{(P)}(𝐫)}{\pi \sigma ^2}}`$ , $`𝒘_{V1}^{(P)}(𝐫)={\displaystyle \frac{𝒘_{V2}^{(P)}(𝐫)}{2\pi \sigma }}.`$ (9)
The contribution due to the chain connectivity is evaluated using the Wertheim’s first-order perturbation theory for a bulk fluid Wertheim87  and extended (using FMT-style weighted densities) by Yu and Wu to inhomogeneous systems Yu02

$$\mathrm{\Phi }_P(\{n_\alpha ^{(P)}\})=\frac{1M}{M}n_0^P\zeta \mathrm{ln}[y_{hs}(\sigma ,\{n_\alpha ^{(P)}\})],$$
(10)
where $`\zeta =1𝒏_{V2}^{(P)}𝒏_{V2}^{(P)}/(n_2^{(P)})^2`$, while $`y_{hs}`$ is connected with the Carnahan-Starling expression for the contact value of the radial distribution function of a hard-sphere mixture

$$y_{hs}(\sigma ,\{n_\alpha ^{(P)}\})=\frac{1}{1n_3^{(P)}}+\frac{n_2^{(P)}\sigma \zeta }{4(1n_3^{(P)})^2}+\frac{(n_2^{(P)})^2\sigma ^2\zeta }{72(1n_3^{(P)})^3}$$
(11)
We note here that both $`\mathrm{\Phi }_{HS}`$ and $`\mathrm{\Phi }_P`$ are independent on the local density of rods.

To specify the polymer-needle contribution we make use of Schmidt’s DFT for hard rod-sphere mixtures Schmidt01  and write the excess free energy density due to vanishingly thin needles as

$$\mathrm{\Phi }_{PN}(\{n_\alpha ^{(i)}\})=n_0^{(N)}\mathrm{ln}(1n_3^{(P)})+\frac{n_1^{(N)}n_2^{(PN)}}{1n_3^{(P)}}.$$
(12)
The equation above implies that the polymer-needle contribution to the excess free energy stems solely from the hard-core repulsion between the needle and the hard sphere forming the polymer segment.

The needle weighted densities, $`n_\alpha ^{(N)}`$, are obtained through spatial convolutions of the needle local density and the corresponding orientation-dependent weight functions

$$n_\alpha ^{(N)}(𝐫,𝝎)=𝑑𝐫^{}\rho _N(𝐫^{},𝝎)w_\alpha ^{(N)}(𝐫𝐫^{},𝝎),\alpha =0,1,$$
(13)
while the “mixed” polymer segment-needle weighted density, $`n_2^{(PN)}`$, is obtained via spatial convolution of the polymer segment density and an orientation-dependent weight function

$$n_2^{(PN)}(𝐫,𝝎)=𝑑𝐫^{}\rho _{PS}(𝐫^{})w_2^{(PN)}(𝐫𝐫^{},𝝎).$$
(14)
Corresponding weight functions are given as Schmidt01

$`w_1^{(N)}(𝐫,𝝎)`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle _{L/2}^{L/2}}𝑑l\delta (𝐫+𝝎l),`$
$`w_0^{(N)}(𝐫,𝝎)`$ $`=`$ $`{\displaystyle \frac{1}{2}}[\delta (𝐫+𝝎l)+\delta (𝐫𝝎l)],`$
$`w_2^{(PN)}(𝐫,𝝎)`$ $`=`$ $`2|𝒘_{V2}^{(P)}(𝐫)𝝎|.`$ (15)
This completes the prescription for the functional. The present theory reduces to the Schmidt’s functional Schmidt01  if $`M=1`$ and to Yu and Wu functional Yu02  if the density of rods $`\rho _N=0`$. Note also that, similarly to Schmidt01 , the functional is linear in the local density of rods.

Although the proposed DFT is intrinsically designed to study inhomogeneous systems, here we restrict ourselves to the isotropic bulk phases where the density of both species are constant. In this case the scalar weighted densities become proportional to the bulk densities of needles and polymers, while the vector weighted densities vanish. Consequently from Eq. 6 we obtain $`n_\alpha ^{(P)}=\xi _\alpha ^{(P)}M\rho _P=\xi _\alpha ^{(P)}\rho _{PS}`$, with $`\xi _3^{(P)}=\pi /6\sigma ^3`$, $`\xi _2^{(P)}=\pi \sigma ^2`$, $`\xi _1^{(P)}=\sigma /2`$ and $`\xi _0^{(P)}=1`$. Likewise Eq. 13 yields $`n_\alpha ^{(N)}=\xi _\alpha ^{(N)}\rho _N`$ where $`\xi _1^{(N)}=L/4`$ and $`\xi _0^{(N)}=1`$. Finally Eq. 14 leads to $`n_2^{(PN)}=\xi _2^{(P)}\rho _{PS}`$. After inserting the above expressions for the weighted densities into the Eq. 5 and Eqs. 10-12 the total free energy per unit volume, $`\mathrm{\Phi }_v`$, is evaluated as $`\mathrm{\Phi }_v=\mathrm{\Phi }_{HS}+\mathrm{\Phi }_P+\mathrm{\Phi }_{PN}+\beta ^1\rho _P[\mathrm{ln}(\rho _P\mathrm{\Lambda }_P^3)1]+\beta ^1\rho _N[\mathrm{ln}(\rho _N\mathrm{\Lambda }_N^3)1]`$. From $`\mathrm{\Phi }_v`$ the pressure, $`p`$, and the chemical potentials of both species are easily calculated

$$\beta p=\mathrm{\Phi }_v+\underset{j=P,N}{}\rho _j\frac{\mathrm{\Phi }_v}{\rho _j},\beta \mu _j=\frac{\mathrm{\Phi }_v}{\rho _j}.$$
(16)
Under appropriate conditions a mixture of polymers and hard rods undergoes entropically driven demixing transition to polymer-rich (rod-poor) and polymer-poor (rod-rich) phases. The coexsisting equilibrium densities (binodals) were obtained by solving simultaneously equations for the equality of pressures and chemical potentials in two phases. The spinodals delimiting the regions stable against fluctuations of density and composition were evaluated from the condition $`det[^2\mathrm{\Phi }_v/\rho _i\rho _j]=0`$, $`i,j=P,N`$. The critical points were evaluated from

$`s^3{\displaystyle \frac{^3\mathrm{\Phi }_v}{\rho _P^3}}+3s^2{\displaystyle \frac{^3\mathrm{\Phi }_v}{\rho _P^2\rho _N}}`$ $`+`$ $`3s{\displaystyle \frac{^3\mathrm{\Phi }_v}{\rho _P\rho _N^2}}+{\displaystyle \frac{^3\mathrm{\Phi }_v}{\rho _N^3}}=0,`$
$`s{\displaystyle \frac{^2\mathrm{\Phi }_v}{\rho _P\rho _N}}`$ $`/`$ $`{\displaystyle \frac{^2\mathrm{\Phi }_v}{\rho _P^2}}.`$ (17)
The above equation arises from the fact that at the critical point the tie-line connecting the coexisting densities becomes tangent to the spinodal line Rowlinson59 .

In Fig. 1 we show examples of binodals (solid lines), spinodals (dashed lines) and critical points (black circles) resulting from the present theory. The phase diagrams plotted in the polymer packing fraction, $`\eta _P=\frac{\pi }{6}\sigma ^3M\rho _P`$ versus scaled needle reservoir density, $`\rho _N^{(r)}L^2\sigma `$ representation (see the inset) were evaluated for systems with $`M=1,10,100`$ for constant $`q=L/\sigma =1`$. For the special case $`M=1`$ (the uppermost diagram) the system reduces to the BF model note1 .

As the chain length is increased, the phase boundary becomes more asymmetric and the critical point moves towards lower polymer packing fractions and towards lower reservoir needle densities. Another interesting feature observed in Schmidt01  is that when the polymer packing fraction-actual needle density representation instead of the polymer packing fraction-reservoir needle density representation is chosen, the critical points for different size ratios $`q`$ are located on the spinodal of the system with $`q=1`$. We find the same behavior also in the present model i.e. for $`M>1`$. In the main plot of Fig. 1 we show phase diagram for a mixture of 10-mers and hard rods. The critical points for mixtures with different size ratios (open diamonds, from top to bottom for $`q=20`$, 15, 8, 5 and 0.1) lie on the spinodal for the system with $`q=1`$ (dashed line).

Let us now consider the limit of very long polymer chains. Intuitively one could argue, that in this regime the phase behavior of a needle-polymer mixture should barely depend on the rod length to polymer-segment diameter ratio $`q=L/\sigma `$ because the physical dimensions of the polymer e.g. the gyration radius, $`R_g`$ become much bigger than the rod elongation. This scenario is captured within the present theory. In Fig. 2 we show the phase diagram for a mixture of hard-sphere monomers ($`M=1`$) and hard rods for the size ratio $`q=1`$. The critical points for different size ratios $`q`$ (black diamonds) are located on the spinodal. Open diamonds denote the critical points for the systems with $`M=2`$, 5, 20, 100, and 10000, respectively. As the chain length increases the critical points for the systems with different size ratios but the same $`M`$ come closer together and for $`M=10000`$ they virtually merge (on the figure scale) into one point. In the limit $`M\mathrm{}`$ the critical scaled needle density tends to the value 2.5526 for all $`q`$ while the critical polymer packing fraction tends to zero. Similar limiting behavior was found in mixtures of spherical colloids and polymers with excluded volume interactions Bolhuis02 ; Sear02 ; Paricaud03 .

In conclusion, in this work we propose a density functional theory for a mixture of vanishingly thin hard rods and polymers modeled as chains built of hard tangent spheres. The functional is constructed by combining the functional due to Yu and Wu for polymer mixtures Yu02  with Schmidt’s functional for hard rod-sphere mixtures Schmidt01 . The proposed theory predicts a demixing transition similar in its nature to that observed for sphere-rod systems. The present functional is well suited to study inhomogeneous systems. It would be of interest to consider a fluid-fluid interface or to investigate surface phase transitions such as entropic wetting or layering that have been discovered in colloid-polymer and colloid-rod mixtures. It is also straightforward to incorporate the Onsager limit Brader02b  of the needle contribution to the functional generating thus a more sophisticated model of inhomogeneous liquid crystal-polymer mixtures. Some of these topics are already under study in our laboratory.

This work has been supported by KBN of Poland under the Grant ”Wpływ samoorganizacji na równowagi fazowe w płynach złożonych”.


text:
# Pentaquark spectrum in string dynamics

## Abstract

The masses of $`uudd\overline{s}`$ and $`uudd\overline{d}`$ pentaquarks are evaluated in a framework of the Effective Hamiltonian approach to QCD using the Jaffe-Wilczek $`[ud]^2\overline{q}`$ approximation. The mass of the $`[ud]^2\overline{s}`$ state is found to be $`400`$ MeV higher than the observed $`\mathrm{\Theta }^+(1540)`$ mass.

Recently several collaborations reported the observation of a very narrow peak in the $`K^+N`$ invariant mass distribution -. The reported masses have been clustered around $`1540`$ MeV, with widths less than $`25`$ MeV. These results are in remarkable agreement with a chiral soliton model prediction  of a $`\mathrm{\Theta }^+`$ state at $`1530`$ MeV and a width less than $`15`$ MeV 111The problem of the $`\mathrm{\Theta }^+`$ width is crucial since the widths beyond the few-MeV level seem to be excluded by the $`K^+N`$ phase-shift analysis  and total cross sections ..

In the sense of the quark model such a state is manifestly exotic – it is not a simple three-quark state but rather the pentaquark $`uudd\overline{s}`$ state. If so, the discovery of $`\mathrm{\Theta }^+`$ provides an opportunity to refine our understanding of nonperturbative quark dynamics at low energy. The $`\mathrm{\Theta }`$-hyperon has hypercharge $`Y=2`$ and third component of isospin $`I_3=0`$. The apparent absence of the $`I_3=+1`$, $`\mathrm{\Theta }^{++}`$ in $`K^+p`$ argues against $`I=1`$, therefore it is usually assumed the $`\mathrm{\Theta }`$ to be an isosinglet, although other suggestions have been made in the literature . The other quantum numbers are not established yet. The uncorrelated quark models, in which all quarks are in the same spatial wave function, naturally predict the ground state energy of a $`J^P=\frac{1}{2}^{}`$ pentaquark to be lower than that of a $`J^P=\frac{1}{2}^+`$ one. Several suggestions were made to reverse this order ,  but no quantitative microscopic evaluations have been performed yet. The QCD sum rules predict a negative parity $`\mathrm{\Theta }^+`$ of mass $`1.5`$ GeV, while no positive parity state was found . The lattice QCD study also predicts that the parity of the $`\mathrm{\Theta }(1540)`$ is most likely negative ,.

Jaffe and Wilczek proposed  that for $`\mathrm{\Theta }(1540)`$ with $`I^P=0^+`$ and other $`q^4\overline{q}`$ baryons the four quarks are bound into two scalar, singlet isospin diquarks combined with the antiquark into a color singlet. The problem of diquarks has a long history and sound physical motivation . Both gluon and pion exchanges in the singlet spin and isospin $`ud`$ state result in the strong short-range attraction which may keep $`u`$ and $`d`$ quarks close together thus forming almost point-like diquark.

The purpose of this letter is to test this interpretation quantitatively using the Effective Hamiltonian (EH) approach in QCD . An attractive feature of this formalism is that it contains the minimal number of input parameters: current (or pole) quark masses, the string tension $`\sigma `$ and the strong coupling constant $`\alpha _s`$, and does not contain fitting parameters as e.g. the total subtraction constant in the Hamiltonian. Therefore the EH method is the best way to compute the masses of yet unknown quark systems like pentaquarks.

In order to illustrate the accuracy of the method, we quote in Table 1 a few EH results for baryons compared with the lattice ones and experiment. One can observe that the accuracy of the EH method for the three-quark systems is $`100`$ MeV or better. We expect the same accuracy for the diquark-diquark-(anti)quark system considered in .

The EH for the three constituents has the form

$$H=\underset{i=1}{\overset{3}{}}\left(\frac{m_i^2}{2\mu _i}+\frac{\mu _i}{2}\right)+H_0+V,$$
(1)
where $`H_0`$ is the kinetic energy operator, $`V`$ is the sum of the perturbative one-gluon exchange potentials and the string potential which is proportional to the total length of the string, i.e to the sum of the distances of (anti)quark or diquarks from the string junction point. The dynamical masses $`\mu _i`$ are expressed in terms of the current masses $`m_i`$ from the condition of the minimum of the hadron mass $`M_H^{(0)}`$ as function of $`\mu _i`$ 333Technically, this is done using the auxiliary field approach to get rid of the square root term in the Lagrangian , . Applied to the QCD Lagrangian, this technique yields the EH for mesons or baryons depending on auxiliary fields $`\mu _i`$. In practice, these fields are finally treated as $`c`$-numbers determined from (2).:

$$\frac{M_H^{(0)}(m_i,\mu _i)}{\mu _i}=0,M_H^{(0)}=\underset{i=1}{\overset{3}{}}\left(\frac{m_i^2}{2\mu _i}+\frac{\mu _i}{2}\right)+E_0(\mu _i),$$
(2)
$`E_0(\mu _i)`$ being eigenvalue of the operator $`H_0+V`$. The physical mass $`M_H`$ of a hadron is

$$M_H=M_H^{(0)}+\underset{i}{}C_i,$$
(3)
where the constants $`C_i`$ have the meaning of the constituents self energies and are expressed in terms of string tension $`\sigma `$ . The self-energy is due to constituent spin interaction with the vacuum background fields and equals zero for any scalar diquark.

In the quark model five quarks are connected by seven strings (Fig. 1a). In the diquark approximation the short legs on this figure shrink to points and the five-quark system effectively reduces to the three-body one (Fig. 1b), studied within the EH approach in ,. Consider a pentaquark consisting of two identical diquarks with current mass $`m_{[ud]}`$ and antiquark with current mass $`m_{\overline{q}}`$ ($`q=d,s`$). The three-body problem is conveniently solved in the hyperspherical formalism . The wave function $`\psi (𝝆,𝝀)`$ expressed in terms of the Jacobi coordinates $`𝝆`$ and $`𝝀`$ can be written in a symbolical shorthand as

$$\psi (𝝆,𝝀)=\underset{K}{}\psi _K(R)Y_{[K]}(\mathrm{\Omega }).$$
(4)
where $`Y_{[K]}`$ are eigen functions (the hyperspherical harmonics)  of the angular momentum operator $`\widehat{K}(\mathrm{\Omega })`$ on the 6-dimensional sphere: $`\widehat{K}^2(\mathrm{\Omega })Y_{[K]}=K(K+4)Y_{[K]}`$, with $`K`$ being the grand orbital momentum. For *identical* diquarks, like $`[ud]^2`$, the lightest state must have a wave function antisymmetric under diquark space exchange. There are two possible pentaquark wave functions antisymmetric under diquark exchange, one (with lower energy) corresponding to the total orbital momentum $`L=1`$, and the second one (with higher energy) corresponding to $`L=0`$. For a state with $`L=1,l_\rho =1,l_\lambda =0`$ the wave function in the lowest hyperspherical approximation $`K=1`$ reads

$$\psi =R^{5/2}\chi _1(R)u_1(\mathrm{\Omega }),u_1(\mathrm{\Omega })=\sqrt{\frac{8}{\pi }}\mathrm{sin}\theta (𝒏_𝝆)_z,$$
(5)
where $`R^2=𝝆^2+𝝀^2`$. For a state with $`L=0,l_\rho =1,l_\lambda =1`$ the wave function in the lowest hyperspherical approximation $`K=2`$ is

$$\psi =R^{5/2}\chi _2(R)u_2(\mathrm{\Omega }),u_2(\mathrm{\Omega })=\frac{4}{\sqrt{\pi ^3}}\mathrm{sin}\theta \mathrm{cos}\theta (𝒏_𝝆𝒏_𝝀).$$
(6)
The Schrödinger equation written in terms of the variable $`x=\sqrt{\mu }R`$, where $`\mu `$ is an arbitrary scale of mass dimension which drops off in the final expressions, reads:

$$\frac{d^2\chi _K(x)}{dx^2}+2\left[E_0+\frac{a_K}{x}b_Kx\frac{(K+\frac{3}{2})(K+\frac{5}{2})}{2x^2}\right]\chi _K(x)=0.,$$
(7)
with the boundary condition $`\chi _K(x)𝒪(x^{5/2+K})`$ as $`x0`$ and the asymptotic behavior $`\chi _K(x)\mathrm{Ai}((2b_K)^{1/3}x)`$ as $`x\mathrm{}`$. In Eq. (7)

$`a_K`$ $`=R\sqrt{\mu }{\displaystyle V_{\text{Coulomb}}(𝒓_1,𝒓_2,𝒓_3)u_K^2𝑑\mathrm{\Omega }},`$ (8)
$`b_K`$ $`={\displaystyle \frac{1}{R\sqrt{\mu }}}{\displaystyle V_{\text{string}}(𝒓_1,𝒓_2,𝒓_3)u_K^2𝑑\mathrm{\Omega }},`$
where

$$V_{\text{Coulomb}}(𝒓_1,𝒓_2,𝒓_3)=\frac{2}{3}\alpha _s\underset{i<j}{}\frac{1}{r_{ij}},$$
(9)
and explicit expression of $`V_{\text{string}}(𝒓_1,𝒓_2,𝒓_3)`$ in terms of Jacobi variables is given in . Note, that in the Y-shape in Fig. $`1b`$, the strings meet at $`120^{}`$ in order to insure the minimum energy. This shape moves continuously to a two–legs configuration where the legs meet at an angle larger than $`120^{}`$.

The mass of the $`\mathrm{\Theta }^+`$ obviously depends on $`m_{[ud]}`$ and $`m_s`$. The current masses of the light quarks are relatively well-known: $`m_{u,d}0`$, $`m_s170`$ MeV. The effective mass of the diquark $`m_{[ud]}`$ is basically unknown. In principle, this mass could be computed dynamically. Instead, one can estimate $`m_{[ud]}`$ from the nucleon mass calculations in the quark-diquark approximation. In this case, neglecting for simplicity Coulomb-like interaction, one obtains for the mass of the nucleon:

$$M_N=M_N^{(0)}\frac{2\sigma }{\pi \mu _d},$$
(10)
where $`M_N^{(0)}`$ and $`\mu _d`$ are defined from the minimum condition,

$$\frac{M_N^{(0)}}{\mu _{[ud]}}=\frac{M_N^{(0)}}{\mu _d}=0.$$
(11)
In Eq. (11)

$$M_N^{(0)}=\frac{m_{[ud]}^2}{2\mu _{[ud]}}+\frac{\mu _{[ud]}+\mu _d}{2}+E_0,E_0=\gamma \left(\frac{\sigma ^2}{2}\frac{\mu _{[ud]}+\mu _d}{\mu _{[ud]}\mu _d}\right)^{1/3}$$
(12)
with $`\gamma =2.338`$ being the first zero of the Airy function: $`\mathrm{Ai}(\gamma )=0`$. Equating $`M_N`$ to the experimental value of the $`N\mathrm{\Delta }`$ center of gravity ($`M_N=1.085`$ GeV), we obtain that $`m_{[ud]}`$ varies in the interval $`0m_{[ud]}300`$ MeV, when $`\sigma `$ varies in the interval $`0.15\text{ GeV}^2\sigma 0.17\text{ GeV}^2`$. The last value of $`\sigma `$ is preferred by meson spectroscopy, while the former one follows from the lattice calculations of Ref. . The lowest $`uudd\overline{s}`$ pentaquark corresponds to the case when one fixes the effective mass of a diquark $`m_{[ud]}=0`$. In what follows, we use $`m_{[ud]}=0`$, $`\sigma =0.15\text{ GeV}^2`$, but explicitly include the Coulomb–like interaction between quark and diquarks.

The calculated masses of $`[ud]^2\overline{s}`$ pentaquarks are given in Table 2. For illustration we show also the masses of $`[ud]^2\overline{d}`$ pentaquarks which have been identified in  with the otherwise perplexing Roper resonance. The states with $`K=2`$ are always higher than those with $`K=1`$ by $``$ 300 MeV, the difference being mainly due to the kinetic energy term in (7). Increasing $`\alpha _s`$ up to $`0.6`$ (the value used in the Capstick-Isgur model ) decreases the $`[ud]^2\overline{s}`$ mass by $``$ $`100`$ MeV. The hyperfine interaction due to $`\eta `$ exchange between diquarks and strange antiquark can lower the $`\mathrm{\Theta }^+`$ energy by $`100`$ MeV or less . Recall that the accuracy of our approach for pentaquarks is expected to be $`100`$ MeV. Since the calculated value is the lower bound on the pentaquark mass, we conclude that the lightest $`uudd\overline{s}`$ state is still $`400`$ MeV higher than the observed $`\mathrm{\Theta }^+(1540)`$ state. Note that $`[ud]^2\overline{d}`$ pentaquarks lie $`100`$ MeV lower than $`[ud]^2\overline{s}`$ pentaquarks. This is the consequence of the different current masses $`m_{\overline{d}}`$ and $`m_{\overline{s}}`$. Also note, that for the $`\frac{1}{2}^{}`$ pentaquark the $`[ud][d\overline{s}]u`$ configuration with $`K_{\text{min}}=0`$ is more preferable giving the mass $`1.8`$ GeV.

In conclusions, we have presented the first dynamical calculations of the pentaquark masses within the Jaffe–Wilczek approximation. Our predictions are universal in the sense that they use only two parameters inherent to QCD: the string tension $`\sigma `$ and the strong coupling constant $`\alpha _s`$. These results do not require any additive mass shift used in constituent quark models. Our predictions can be viewed as lower bounds for the masses of pentaquark states. We showed that the $`\mathrm{\Theta }^+`$ mass is much higher than the observed one. This implies that $`\mathrm{\Theta }^+(1540)`$ can not be explained in terms of the standard string interaction between quarks in the diquark–diquark–antiquark picture. Note that the chiral interactions (e.g. pion exchanges) are automatically taken into account inside diquarks, while are not possible between diquarks and strange antiquark. Therefore a drastic modification of the present results requires a completely new dynamics, either the chiral soliton type, as in , or other quark clusters, like $`q^4\overline{q}`$  or $`q^2\overline{s}q^2`$ .

The authors are grateful to K.G.Boreskov for careful reading of the manuscript and valuable suggestions. This work was supported by RFBR grant $`\mathrm{\#}`$ 03-02-17345 and PRF grant for leading scientific schools $`\mathrm{\#}`$ 1774.2003.2.


text:
# The fluid–fluid interface in a model colloid–polymer mixture: Application of grand canonical Monte Carlo to asymmetric binary mixtures

(October 17, 2003)

## Abstract

We present a Monte Carlo method to simulate asymmetric binary mixtures in the grand canonical ensemble. The method is used to study the colloid–polymer model of Asakura and Oosawa. We determine the phase diagram of the fluid–fluid unmixing transition and the interfacial tension, both at high polymer density and close to the critical point. We also present density profiles in the two–phase region. The results are compared to predictions of a recent density functional theory.

In colloid experiments, hard–sphere–like systems can be realized in which the corresponding phase behavior is of purely entropic origin. An example is the fluid–fluid unmixing transition that is observed in solutions of colloids and non–adsorbing polymers asakura1954a ; vrij1976a . This transition is due to a depletion effect asakura1954a ; vrij1976a ; vrij1997a ; roth2000a ; aarts2002a . Each colloidal particle is surrounded by a depletion zone from which polymers are excluded. When two colloids are close together, their depletion zones may overlap, thereby increasing the free volume available to the polymers. This results in an anisotropic pressure exerted by the polymers onto the colloids which gives rise to an effective attraction between the colloids. Such depletion forces also occur in mixtures of large and small hard spheres, but in this case it is still debated whether phase separation occurs bhs .

A simple model for colloid–polymer mixtures was first introduced by Asakura and Oosawa asakura1954a  and later independently by Vrij vrij1976a . In this model (the so–called AO model) colloids and polymers are treated as spheres with respective radii $`R_\mathrm{c}`$ and $`R_\mathrm{p}`$. Hard sphere interactions are assumed between colloid–colloid (cc) and colloid–polymer (cp) pairs, while polymer–polymer (pp) pairs can interpenetrate freely. This yields the following pair potentials for the AO model:

$`u_{\mathrm{cc}}(r)`$ $`=`$ $`\{\begin{array}{cc}\mathrm{}\hfill & \text{for }r<2R_\mathrm{c}\hfill \\ 0\hfill & \text{otherwise,}\hfill \end{array}`$ (3)
$`u_{\mathrm{cp}}(r)`$ $`=`$ $`\{\begin{array}{cc}\mathrm{}\hfill & \text{for }r<R_\mathrm{c}+R_\mathrm{p}\hfill \\ 0\hfill & \text{otherwise,}\hfill \end{array}`$ (6)
$`u_{\mathrm{pp}}(r)`$ $`=`$ $`0,`$
where $`r`$ is the distance between two particles. The polymers thus represent ideal polymer coils with a radius of gyration $`R_\mathrm{p}`$ which can be realized experimentally in a $`\theta `$ solvent.

The AO model has been the subject of many studies in the framework of density functional theories (DFT) brader2000 ; schmidt2002a ; brader2002a . In particular, these studies yielded the phase diagram for a wide range of colloid to polymer size ratios. Moreover, in the case of the fluid–fluid unmixing transition, they predicted interfacial tensions that are roughly one thousand times lower than those for simple liquids, in agreement with experiments hoog1999a ; chen2000 . A drawback of the latter DFTs is that they are mean field theories and thus cannot recover the 3D Ising critical behavior observed, for instance, in a recent experiment on a real colloid–polymer mixture chen2000 .

Monte Carlo simulations are also well suited to study the phase behavior of colloid–polymer mixtures. Recent simulations in the Gibbs ensemble were performed to determine phase diagrams of the AO model dijkstra1997a  and also of a model that considers non–zero interactions between the polymers bolhuis2002a . In the latter study even quantitative agreement with experiments was obtained. However, interfaces are absent in the Gibbs ensemble panagiotopoulos1987a , so these simulations do not enable investigations close to the critical point, nor investigations of the interface in the two–phase region.

The general problem in the simulation of asymmetric binary mixtures (such as the AO model) is that, by displacing or inserting a large particle, overlap with a number of small particles will likely result. Such overlaps are generally unfavorable and will be rejected in the majority of cases. If no special steps are taken, one ends up displacing mainly small particles while the large particles remain essentially frozen. For asymmetric binary mixtures in the canonical ensemble a number of specialized algorithms have been developed that circumvent this problem biben1996a ; buhot1998a . Unfortunately, it is difficult to obtain the surface tension accurately in the canonical ensemble because it must be derived from the rather small anisotropy of the pressure tensor in that case: long wavelength interfacial fluctuations are hard to equilibrate in the canonical ensemble varnik2000a .

In this letter we present a grand canonical Monte Carlo method that enables direct simulations of asymmetric binary mixtures. We use the method to study the fluid–fluid unmixing transition in the AO model. Our method consists of collective Monte Carlo moves in conjunction with an umbrella sampling technique that was recently developed by Virnau and Müller virnau2003a . This way, we are able to calculate the phase diagram of the AO model close to the critical point. At the same time we can calculate the interfacial tension $`\gamma `$ in the two–phase region using the method of Binder binder1982a . We present a comparison of our data to recent DFT results brader2002a ; schmidt2002a ; schmidt2003a  and we show that the AO model displays 3D Ising critical behavior.

The grand canonical ensemble requires that the temperature $`T`$, the volume $`V`$ and the respective chemical potentials $`\{\mu _\mathrm{p},\mu _\mathrm{c}\}`$ of the small (polymer) and large (colloid) particles are fixed. In a standard grand canonical move one tries to insert or remove a particle using Metropolis sampling landau2000a . This approach is not efficient for asymmetric binary mixtures because the insertion of large particles (colloids) is severely hampered by the presence of small particles (polymers). The method that we present is aimed to circumvent this problem. The main idea is to not transfer particles one–at–a–time, but to swap clusters of small particles instead. This stimulates the formation of voids, in which a large particle can be inserted without producing overlap.

Figs. 1A and 1B demonstrate the insertion of an additional large particle into an asymmetric binary mixture currently containing $`N_\mathrm{c}`$ large particles. First, a point is selected randomly in the mixture and a sphere with radius $`\delta `$ and volume $`V_\delta =4\pi \delta ^3/3`$ is drawn around it (see Fig. 1A). Let $`n_\mathrm{p}`$ denote the number of small particles inside the sphere: a particle is inside the sphere when the coordinates of its center are inside the sphere. Next, we choose a uniform random integer $`n_\mathrm{r}`$ from the interval $`0n_\mathrm{r}<m`$, with $`m`$ an integer that will be specified later. If $`n_\mathrm{r}>n_\mathrm{p}`$ the move is rejected but if $`n_\mathrm{r}n_\mathrm{p}`$, $`n_\mathrm{r}`$ small particles are randomly selected from the sphere out of the $`n_\mathrm{p}`$ present (these particles are shaded grey in Fig. 1A). The $`n_\mathrm{r}`$ selected particles are then removed from the mixture and a large particle is inserted at the center of the sphere (see Fig. 1B). The new configuration is accepted with probability:

$$A_+=\mathrm{min}[1,\frac{z_\mathrm{c}V}{N_\mathrm{c}+1}\frac{(n_\mathrm{p})!}{(n_\mathrm{p}n_\mathrm{r})!}\frac{\mathrm{e}^{\beta \mathrm{\Delta }E}}{(z_\mathrm{p}V_\delta )^{n_\mathrm{r}}}],$$
(7)
with $`\mathrm{\Delta }E`$ the potential energy difference between the initial and the final configuration and $`\{z_\mathrm{c},z_\mathrm{p}\}`$ the fugacity of large and small particles, respectively. The fugacity $`z`$ is related to the chemical potential via $`z=\mathrm{exp}(\beta \mu )`$ with $`\beta =1/(k_BT)`$ and $`k_B`$ the Boltzmann constant.

The reverse move is illustrated in Figs. 1C and 1D. First, a large particle is selected at random and a sphere with radius $`\delta `$ is drawn around the center of this particle. Next, a uniform random integer $`n_r`$ from the interval $`0n_r<m`$ is chosen followed by the selection of $`n_r`$ random locations inside the sphere. These random locations are marked as crosses in Fig. 1C. The selected large particle is now removed from the mixture and $`n_\mathrm{r}`$ small particles are placed on the locations selected before (the newly inserted particles are shaded grey in Fig. 1D). The new configuration is accepted with probability:

$$A_{}=\mathrm{min}[1,\frac{N_\mathrm{c}}{z_\mathrm{c}V}\frac{(n_\mathrm{p})!(z_\mathrm{p}V_\delta )^{n_\mathrm{r}}}{(n_\mathrm{p}+n_\mathrm{r})!}\mathrm{e}^{\beta \mathrm{\Delta }E}],$$
(8)
the notation being the same as in Eq. (7).

It is straightforward to show that the acceptance probabilities $`A_+`$ and $`A_{}`$ enforce detailed balance, which ensures that the algorithm is not statistically biased newman1999a . The algorithm is also ergodic because single large particles have a finite probability of being inserted anywhere in the system with one move. Similarly, a small particle can be inserted anywhere via a combination of moves: for example, by the insertion of a large particle followed by the removal of the same large particle.

In order to apply the method to the AO model, the parameters $`\delta `$ and $`m`$ still need to be specified. We use $`\delta =R_\mathrm{c}+R_\mathrm{p}`$ which is just big enough to contain one colloid in a sea of polymers. The integer $`m`$ must be chosen large enough to allow for the formation of voids. In the pure polymer phase, the polymer density equals $`z_\mathrm{p}`$ because the polymers behave like an ideal gas. The insertion sphere will then contain $`z_\mathrm{p}V_\delta `$ polymers on average and thus we choose $`m`$ slightly above this value. With this choice of $`\delta `$, the insertion of a colloid can only succeed if all polymers are removed from the insertion sphere $`V_\delta `$. This will in general happen a fraction $`1/m`$ of the time. To boost the acceptance rate, we choose to remove all polymers from $`V_\delta `$ when we attempt to insert a colloid, provided their number does not exceed $`m`$: moves that attempt to remove more than $`m`$ polymers are rejected. To maintain detailed balance, the acceptance probabilities $`A_+`$ and $`A_{}`$ must be multiplied by $`1/m`$ and $`m`$, respectively.

The phase diagram of the AO model can be expressed in terms of the reduced polymer packing fraction $`\eta _\mathrm{p}^\mathrm{r}(4\pi /3)z_\mathrm{p}R_\mathrm{p}^3`$ as a function of the colloid packing fraction $`\eta _\mathrm{c}(4\pi /3)R_\mathrm{c}^3N_\mathrm{c}/V`$. This is analogous to the temperature–density phase diagram for simple fluids. In case of the AO model, $`\eta _\mathrm{p}^\mathrm{r}`$ plays the role of inverse temperature and $`\eta _\mathrm{c}`$ that of order parameter. In order to determine the coexistence curve of the fluid–fluid unmixing transition, we calculate for a given value of $`\eta _\mathrm{p}^\mathrm{r}`$ the probability distribution $`P(\eta _\mathrm{c})`$. This is the probability of observing a mixture with colloid packing fraction $`\eta _\mathrm{c}`$. If phase separation occurs $`P(\eta _\mathrm{c})`$ is bimodal: the peak at low $`\eta _\mathrm{c}`$ corresponds to the colloid vapor phase, the peak at high $`\eta _\mathrm{c}`$ to the colloid liquid phase, and the region in between is the phase–separated regime. To ensure phase coexistence, the colloid fugacity $`z_\mathrm{c}`$ is tuned such that the area under both peaks is equal landau2000a .

A crucial point in our simulation is the use of a new biased sampling technique called successive umbrella sampling. This technique was recently developed by Virnau and Müller virnau2003a . Combination of our grand canonical Monte Carlo move and successive umbrella sampling enables us to sample $`P(\eta _\mathrm{c})`$ also in regions where, due to the free energy barrier separating both phases, $`P(\eta _\mathrm{c})`$ may be very low.

In the following we consider an AO mixture with a polymer to colloid size ratio $`qR_\mathrm{p}/R_\mathrm{c}=0.8`$. Note that all lengths are given in units of $`R_\mathrm{c}1`$. From previous studies we expect fluid–fluid phase separation for $`q=0.8`$ in a wide range of $`\eta _\mathrm{p}^\mathrm{r}`$. The simulations were performed in an elongated box with aspect ratio $`1/2`$ and periodic boundary conditions. We have performed simulations for two different volumes of the simulation box: $`V_1=4652.4`$ and $`V_2=9314.9`$. When the colloid packing fraction reaches $`0.45`$, the smaller volume contains $`500`$ colloids and the larger volume $`1000`$ colloids.

Fig. 2 shows the phase diagram for the two system sizes $`V_1`$ and $`V_2`$. Comparison of the two data sets shows that finite–size effects are relatively small, even close to the critical point. Also included in Fig. 2 is the phase diagram obtained from recent DFT schmidt2002a ; brader2002a . We observe that DFT underestimates the location of the critical point by about 30%. More importantly, DFT yields the typical mean–field parabolic shape of the binodal (critical exponent $`\beta =1/2`$) while the simulation yields the expected flatter binodal ($`\beta 0.325`$; Ising model universality class luiten2001a ). We will demonstrate below that the simulation indeed displays Ising critical behavior.

The probability distribution $`P(\eta _\mathrm{c})`$ can also be used to extract the interfacial tension $`\gamma `$ between the coexisting colloid vapor and colloid liquid phases. To this end one uses a formula that was first derived by Binder binder1982a :

$$\gamma \underset{L\mathrm{}}{lim}\frac{F_\mathrm{L}}{2L^2}=\underset{L\mathrm{}}{lim}\frac{1}{2L^2}\mathrm{ln}\left[\frac{P^{\mathrm{max}}(\eta _\mathrm{c})}{P^{\mathrm{min}}(\eta _\mathrm{c})}\right],$$
(9)
with $`P^{\mathrm{max}}(\eta _\mathrm{c})`$ and $`P^{\mathrm{min}}(\eta _\mathrm{c})`$ the value of $`P(\eta _\mathrm{c})`$ at its maxima and its minimum, respectively, and $`L`$ the length of the simulation box parallel to the interface. An example distribution is shown in the inset of Fig. 3 for $`\eta _\mathrm{p}^\mathrm{r}=1.225`$. Note that the presence of a flat region between the two peaks is important for an accurate estimate of $`\gamma `$. This is enforced in the simulation by using an elongated box.

Fig. 3 shows the reduced interfacial tension $`\gamma ^{}4R_\mathrm{c}^2\gamma `$ as a function of the difference between the colloid packing fractions in the coexisting phases, as obtained from our simulation together with a recent DFT result. Again, the DFT result deviates from the simulation by about 30%. This demonstrates that the rather perfect agreement of DFT with experimental data as claimed in Ref. brader2002a  is coincidental. According to our simulation data, the values of $`\gamma ^{}`$ for the AO model underestimate experimental data significantly which shows that more sophisticated models are required to describe colloid–polymer mixtures on a quantitative level.

We have carried out additional simulations to calculate colloid density profiles in the two–phase region. These simulations were performed using $`1156`$ and $`2889`$ colloidal particles at a colloid packing fraction of $`0.13`$. Periodic boundary conditions were again used and the aspect ratio of both boxes was $`1/3`$. This corresponds to a box length of $`L_\mathrm{z}=69.46`$ and $`L_\mathrm{z}=94.27`$ for the smaller and larger system, respectively (here $`L_\mathrm{z}`$ is the box length perpendicular to the interface). Three density profiles for the larger system are shown in the inset of Fig. 4. We have estimated the (10%–90%)-interfacial width $`W`$ from these profiles by fitting the profiles to a hyperbolic tangent. Since the interfacial width is expected to diverge with the same exponent as the bulk correlation length near the critical point, one expects $`W\left(\eta _\mathrm{p}^\mathrm{r}\eta _{\mathrm{p},\mathrm{crit}}^\mathrm{r}\right)^\nu `$, with $`\nu =0.63`$ corresponding to 3D Ising critical behavior luiten2001a . Fig. 4 shows that the data for both system sizes is consistent with this power law. We can also infer from Fig. 4 that $`W`$ depends strongly on $`L_z`$, even far away from the critical point. This is most likely due to capillary waves binder2000  and care has to be taken when comparing $`W`$ from this simulation to interfacial widths obtained in experiments or analytical theories.

In summary, we have presented a grand canonical Monte Carlo method which is well suited to simulate asymmetric binary mixtures. It is particularly powerful when combined with a re–weighting scheme: both the phase diagram and the surface tension can be obtained accurately in that case. We have used the method to determine the coexistence line of the fluid–fluid transition in the AO model with high accuracy. We have also presented new analysis of the interface between coexisting phases in the AO model, in particular estimates of the interfacial tension.

We are grateful to the Deutsche Forschungsgemeinschaft for support (TR6/A5) and to K. Binder, M. Müller, M. Schmidt and P. Virnau for stimulating discussion.


text:
# Noise-induced synchronization for phase turbulence

## Abstract

Phase turbulence is suppressed by applying common noise additively to the Kuramoto-Sivashinsky type equation, and the noise-induced phase synchronization is realized. The noise strength necessary for the suppression of phase turbulence is evaluated theoretically.
PACS: 05.45.+b, 05.40.+j, 02.50-r

Recently, various noise effects to nonlinear systems have been studied. The response of a bistable system or an excitable system to a periodic force is enhanced by the noise effect. The stochastic resonance improves signal detection by the superposed noise . Noise-enhanced entrainment among coupled oscillators is found experimentally in Belousov-Zhabotinsky reactions . Frequency locking of noise-sustained oscillations is found in coupled excitable systems . A small amount of noise may change a chaotic trajectory into a rather regular trajectory, and it is called noise-induced order . Common noise may induce complete synchronization for uncoupled chaotic oscillators . It is called noise-induced synchronization. In this paper, we apply common noise for a modified equation of the Kuramoto-Sivashinsky equation. Without the common noise, the model equation exhibits phase turbulence. Noise-induced synchronization occurs and a spatially uniform state is observed owing to the common noise.

The model equation is written as

$$\varphi _t=\omega r\mathrm{sin}\varphi \mu \varphi _{xx}\varphi _{xxxx}+\varphi _x^2+\xi (t),$$
(1)
where $`\varphi (x,t)`$ is a phase variable, $`\omega >r`$ is assumed, and $`\xi (t)`$ represents spatially uniform Gaussian white noise satisfying

$$\xi (t)\xi (t^{})=2D\delta (tt^{}).$$

The model equation is equal to the original Kuramoto-Sivashinsky equation, when $`r=D=0`$. The Kuramoto-Sivashinsky equation exhibits phase turbulence for $`\mu >0`$, owing to the long-wavelength instability of the uniform state. For nonzero $`r`$, the phase motion of each oscillator is not uniform. This type of phase oscillator model has been used to study cooperative phenomema in coupled limit cycle oscillators . Equation (1) represents a model of spatially coupled phase oscillators. Without the noise term, the spatially uniform rotation is unstable and a phase turbulence is observed similarly to the original Kuramoto-Sivashinsky equation.

Firstly, we study the stability of the spatially uniform state under the common noise. The spatially uniform state $`\varphi (x,t)=\varphi _0(t)`$ satisfies the Langevin equation

$$d\varphi _0/dt=\omega r\mathrm{sin}\varphi _0+\xi (t).$$
(2)
The corresponding Fokker-Planck equation is

$$\frac{P}{t}=\frac{}{\varphi _0}\{(\omega r\mathrm{sin}\varphi _0)P\}+D\frac{^2P}{\varphi _0^2}.$$
(3)
For $`D=0`$, the stationary probability distribution is

$$P_s(\varphi _0)\frac{1}{\omega r\mathrm{sin}\varphi _0}.$$
(4)
For nonzero D, the stationary probability distribution is expressed as

$$P_s(\varphi _0)=\mathrm{exp}\{(\omega \varphi _0r+r\mathrm{cos}\varphi _0)/D\}P_s(0)\left\{1+\frac{(e^{2\pi \omega /D}1)_0^{\varphi _0}e^{(\omega \psi r\mathrm{cos}\psi )/D}𝑑\psi }{_0^{2\pi }e^{(\omega \psi r\mathrm{cos}\psi )/D}𝑑\psi }\right\},$$
(5)
where $`P_s(0)`$ is determined from the normalization condition $`_0^{2\pi }P_s(\varphi _0)𝑑\varphi _0=1`$. This stationary probabiliy distribution approaches a uniform distribution $`P_s=1/2\pi `$ for $`D\mathrm{}`$. The stationary distribution for $`\omega =0.8,r=0.5`$ and $`D=0.7`$ is displayed in Fig. 1(a). The peak of the probability distribution is located at $`\varphi \pi /4<\pi /2`$.

The linearized equation of the Fourier amplitude $`\varphi _k`$ with wave number $`k`$ around the uniform state for Eq. (1) satisfies an equation

$$d\varphi _k/dt=(r\mathrm{cos}\varphi _0(t)+\mu k^2k^4)\varphi _k,$$
(6)
where $`\varphi _k=1/\sqrt{L}_0^L\varphi (x)e^{ikx}𝑑x`$, and $`\varphi _0(t)`$ obeys Eq. (2). The long-time average of $`r\mathrm{cos}\varphi _0(t)+\mu k^2k^4`$ determines the stability for the perturbation with wave number $`k`$. For $`D=0`$, the average value of $`\mathrm{cos}\varphi _0`$ is zero, since the distribution (4) is symmetric around $`\varphi =\pi /2`$. For $`D\mathrm{}`$, the average value of $`\mathrm{cos}\varphi _0(t)`$ is also zero, since the distribution is uniform. For general parameter values of $`D`$, the average value of $`r\mathrm{cos}\varphi _0(t)`$ is negative, since the probability distribution has a peak below $`\pi /2`$ as shown in Fig. 1(a), and $`\mathrm{cos}\varphi _0`$ is positive for $`0<\varphi _0<\pi /2`$. We have calculated the average value of $`r\mathrm{cos}\varphi _0(t)`$ using the stationary distribution (5) for various $`D`$’s for $`\omega =0.8`$ and $`r=0.5`$, and plotted the results in Fig. 1(b). The curve of $`r\mathrm{cos}\varphi _0`$ has a characteristic form of the resonant one, although it is not easily expected from the distribution (5). The average value of $`r\mathrm{cos}\varphi _0`$ approaches 0 for $`D0`$ and $`D\mathrm{}`$ and it takes a minimum value -0.08 at $`D0.7`$. The average value of $`r\mathrm{cos}\varphi _0(t)+\mu k^2k^4`$ takes a maximum $`r\mathrm{cos}\varphi _0+\mu ^2/4`$ at $`k=\sqrt{\mu /2}`$, when $`k`$ is changed. The uniform state is linearly stable if $`r\mathrm{cos}\varphi _0+\mu ^2/4<0`$. At $`\omega =0.8`$, $`r=0.5`$ and $`D=0.7`$, the uniform state is stable if $`\mu `$ is smaller than 0.4. For $`\mu =0.5`$, $`\mu k^2k^4`$ takes a maximum value 0.0625 at $`k=0.5`$, therefore, the uniform state is stable for $`0.3<D<1.58`$. It implies that the common noise of intemediate strength can suppress long-wavelength fluctuations and stabilize the uniform state. Common noise of intermediate strength is most effective for the noise-induced synchronization, which is analogous to the stochastic resonance.

We have performed numerical simulation to see the noise-induced synchronization for the modified Kuramoto-Sivashinsky equation. We have used the finite difference method with time step 0.0001 and space step 160/512 for a numerical simulation, and the periodic boundary conditions with system size 160 is used. The parameter values are $`\omega =0.8`$, $`r=0.5`$ and $`\mu =0.5`$. The noise intensity $`D`$ was changed as a control parameter. Figure 2 displays snapshot profiles of $`\varphi (x)`$ at $`D=0,0.2,0.7`$ and 2.5. At $`D=0`$, the spatio-temporal chaos is observed. At $`D=0.7`$, the spatially uniform state appears. It corresponds to the noise-induced synchronization. At $`D=0.2`$ and $`D=2.5`$, the spatial profile is not uniform, but fairly regular. The wavelength of the most dominant Fourier mode is 160/12, and the corresponding wavenumber is 0.471, which is close to $`\sqrt{\mu /2}=0.5`$. The amplitude of spatial fluctuations changes intermittently in time. The intermittent time sequence of the Fourier amplitude $`|\varphi _k(t)|`$ with $`k=0.471`$ for $`D=0.2`$ is shown in Fig. 4(a). The amplitudes of the spatial fluctuations seem to differ very much for $`D=0.2`$ and 2.5 in Figs. 2(b) and (d), but it is a problem of timing to take snapshots. The average amplitudes of the spatial fluctuations do not differ so much for the two parameter values of $`D`$, as shown in Figs. 3(b) and (d).

Figure 3 displays the average of the Fourier amplitude: $`P_w(k)=|\varphi _k|`$ for the four parameters of the noise intensity. For $`D=0`$, the spectrum has a peak at $`k\sqrt{\mu /2}`$ and increases near $`k0`$. This is a characteristic spectrum of phase turbulence. For $`D=0.7`$, the uniform state is stable and the spectrum is almost zero. For $`D=0.2`$ and $`D=2.5`$, the spectrum has strong peaks at $`k\sqrt{\mu /2}`$ and its higher harmonics, but the long-wavelength fluctuations are suppressed in contrast to the case of $`D=0`$. The suppression of the long-wavelength fluctuations makes the spatial profiles of $`\varphi (x)`$ regular, as shown in Fig. 2.

The time evolution is intermittent, and $`\varphi (x,t)`$ is close to the uniform state $`\varphi _0(t)`$, when the amplitude of spatial fluctuations becomes sufficiently small. Then, the Fourier amplitudes obey Eqs. (2) and (6) approximately. The variable $`\varphi _0(t)`$ is a random varibale, since Eq. (2) is a Langevin equation. Equation (6) represents a random multiplicative process, since the growth rate $`r\mathrm{cos}\varphi _0(t)+\mu k^2k^4`$ for the Fourier amplitude with wave number $`k`$ is randomly fluctuating. If the average value of the growth rate is positive (negative), the Fourier amplitude generally tends to grow from zero (decay to zero). Even if the average value of the growth rate is positive, there are chances that the growth rate keeps negative in a certain time interval. The Fourier amplitude takes very small values for the time interval. Inversely, even if the average value of the growth rate is negative, there are chances that the growth rate keeps positive in a certain time interval. The Fourier amplitude takes non-small values for the time interval. This is a mechanism that a random multiplicative process induces an intermittent time sequence . When the amplitude of spatial fluctuations becomes large, the nonlinear term acts and the growth of the spatial fluctuations is suppressed. Figures 4(a) and (c) display time sequences of Fourier amplitudes $`|\varphi _k|`$ with (a) $`k=0.471`$ and (c) $`k=0.942`$ at $`D=0.2`$. The average growth rate is positive for wave number $`k=0.471`$ and the average growth rate is negative for wave number $`k=0.942`$. For both wave numbers, intermittent time sequences appear. Figures 4(b) and (d) display histograms of the probability distributions $`P(|\varphi _k|)`$ of $`|\varphi _k|`$ for (a) $`k=0.471`$ and (b) $`k=0.942`$. The probability distributions obey approximately a power law for a rather large range of $`k`$, which is characteristic of multiplicative stochastic processes . The exponent of the power-law distribution is nearly -0.74 for $`k=0.471`$ and -0.93 for $`k=0.942`$. The exponents change with wave numbers. Similar type intermittent behaviors are also observed for $`D=2.5`$.

To summarize, we have found noise-induced synchronization for the modified Kuramoto-Sivashinsky equation. The phase turbulence is suppressed by the common noise. The mechanism for the noise-induced synchronization is clear in this model equation. That is, the stationary phase distribution of the uniform state is deformed by the noise and the uniform state is stabilized against the long-wavelength fluctuations. For $`r=0`$ as in the original Kuramoto-Sivashinsky equation, this type of noise-induced synchronization does not occur. The growth rate of the Fourier amplitude with the long wavenumber is decreased especially for an intermediate range of the noise sterngth. The average growth rate is positive for a finite interval of wave numbers, and the spatial profile becomes rather regular owing to the common noise. The fluctuation of the growth rate induces an intermittent time sequence and a power-law distribution for the Fourier amplitude has been observed.


text:
# Depairing currents in superconductor / ferromagnet Nb/CuNi trilayers close to Tc.

## 1 Introduction

The issue of the proximity effect between a superconductor and a ferromagnet is the focus of much current research. In principle, the superconductivity in the S layer is suppressed due to the pair breaking of the Cooper pair by the exchange field $`h_{ex}`$ experienced in the ferromagnet. In the F-layer, a superconducting order parameter can still exist, but it becomes spatially modulated because the electrons which create a Cooper pair belong to different spin subbands. This leads to several effects observed experimentally, including oscillatory behavior of the superconducting transition temperature Tc as function of ferromagnetic layer thickness in F/S bi- or multilayers , and $`\pi `$-junctions in S/F/S Josephson junctions . A different mechanism to influence the superconducting order parameter in the S-layer by the adjacent F-layers is by varying the relative directions of the magnetization (and therefore $`h_{ex}`$) in the two F-layers. It was predicted that, when the thickness of superconductor dS is of the order of the superconducting coherence length $`\xi _S`$, the superconducting transition temperature is higher when the two magnetization directions are antiparallel than when they are parallel to each other . The difference can even be enhanced by tuning the thickness of the F-layers $`d_F`$ such that $`d_F`$ $`\xi _F/2`$ (with $`\xi _F`$ the coherence length in the ferromagnet) . A first observation of these so-called spin-switch effects was recently reported in ref. . Basically, spin switching comes about when the Cooper pair samples different directions of $`h_{ex}`$ simultaneously, which raises the question whether similar effects can be observed if the F-layers contain a domain structure rather than a homogeneous magnetization . Two recent experiments reported an influence of the domain state on the transport properties of the superconductor. In Nb/Co bilayers it was found that the superconducting transition temperature Tc and critical current density Jc (defined using a voltage criterion) are higher at the coercive field of the Co layer ; a similar effect on $`J_c`$ was observed in Co/Nb/Co trilayers, again using a voltage criterion . Here we present results of a similar study, which differs in two ways from the ones cited above. We use F/S/F trilayers with S = Nb and F = Cu0.46Ni0.54, which is a very weak ferromagnet with a saturation magnetization $`M_s`$ of about 0.1 $`\mu _B`$ ($`\mu _B`$ is the Bohr magneton)  to be compared to about 1 $`\mu _B`$ for Co. The experiment we perform is measuring the depairing current by a pulsed-current method. This quantity reflects the superconducting order parameter directly, and we showed before that it can be used in S/F systems . We find that at least close to Tc the critical current of the trilayer as function of an in-plane magnetic field peaks at the coercive field Hco of the F-layer, which means that superconductivity is enhanced in the domain state.

## 2 Sample preparation and magnetic properties

Single F-films and F/S/F trilayers were DC-magnetron sputtered in an ultra high vacuum system with base pressure of about 10-9 mbar and sputtering argon pressure about 6\*10-3 mbar. A Cu1-xNix target was used with $`x`$ = 0.50, which yielded a slightly different Ni concentration in the samples of $`x`$ =0.54. We will call this alloy CuNi. The thickness of the F-layers was chosen to be 9 nm, which is about 0.5 $`\xi _F`$ . The S-layer thickness was chosen at 23 nm in order to have the superconducting transition temperature around 4 K .

Electron beam lithography was used to pattern the samples for 4-point measurements, with the bridge between the voltage contacts 2 $`\mu `$m wide and 20 $`\mu `$m long. The patterned samples were Ar-ion etched. Unstructured samples were measured by SQUID-magnetometry in order to determine Ms, Hco and the Curie temperature TCurie. Fig. 1 shows the typical dependence of the magnetization $`M`$ on temperature $`T`$ for a single F-layer, using a standard procedure: the sample was magnetized at 5 K to its saturation in a field of 0.7 T, then the field was set to a small value (typically of about 0.1 mT), and $`M`$(T) was measured up to 300 K and back down to 5 K. TCurie was defined at the temperature where irreversibility sets in when cooling down. The dependence of $`M`$ on applied field $`H_a`$ at 5 K of a single F-layer is given in Fig. 2a. It shows a hysteresis loop typical for a ferromagnet with Hco about 25 mT. To determine possible differences in Hco between structured and unstructured samples caused by the sample shape anisotropy, we measured the magnetoresistance R(Ha) of a structured samples at 5 K. With Ha in the plane of the film, but perpendicular to the bridge and therefore the current, R(Ha) behaves as shown in Fig. 2b. Two peaks are found at Hco, due to the anisotropic magnetoresistance effect. The slight difference in aspect ratio $`A`$ for the structured ($`A`$ = 10) and the unstructured ($`A`$ = 5) samples does not lead to a significant difference in Hco.

## 3 Superconducting properties and discussion

For the determination of the critical current Ic at different temperatures a pulsed-current method was used. Current pulses of about 3 ms with growing amplitude were sent through the sample below Tc = 4.0 K (shown in the onset of Fig. 3). Each pulse was followed by a long pause of about 7 s. The voltage response of the system was observed on a oscilloscope triggered for the time of a single pulse. To improve the signal resolution a differential amplifier combined with low-noise band filters was used. Two typical current (I)- voltage (V) characteristics are shown in Fig. 3, taken at temperatures 3.9 K (reduced temperature t = T / Tc = 0.98) and 3.7 K (t = 0.93). One can see a clear jump from the superconducting into the normal state at Ic, which we designate the depairing current . For all samples, down to t = 0.98 a small onset voltage was observed below Ic, probably because of vortex motion. In order to be sure that this has no influence on the depairing current, the sample temperature was probed during every current pulse. A significant increase was found only when Ic was reached.

The dependence of Ic as function of Ha at T = 3.8 K is shown in Fig. 4. Upon lowering Ha from the positive high field side, Ic rises and reaches its maximum value at a negative value of about -20 mT, quite close to Hco. It then decreases again. When increasing Ha from a large negative field, a similar maximum occurs at +20 mT, with clearly hysteretic behavior. Clearly, the suppression of superconductivity in the S-layer by the proximity effect is less for the case when a domain structure is present in the outer F-layers, which is similar to the findings in refs. . Note that the result does not imply the existence of a coupling between the F-layers : the variation of exchange field directions in the domain walls might be sufficient to reduce the suppression. The scatter in the data points in Fig. 4 may be caused by the changing domain structure; either a varying suppression of the order parameter or an induced vortex state  could result in a slightly different current distribution at each point. The main result of a peak in Ic around Hco is also found for the case of thick ($``$ 50 nm) ferromagnetic layers. This seems to exclude the possibility that effects of inducing an inhomogeneous order parameter in the F-banks play a role, at least in this temperature regime.

## 4 Acknowledgements

This work is part of the research program of the ’Stichting voor Fundamenteel Onderzoek der Materie (FOM)’, which is financially supported by NWO.