text:
# Path integral Monte Carlo simulations of silicates

## I Introduction

The study and analysis of materials properties of crystalline silicates are important since these systems are used in many industrial processes and they occur also in many natural rocks. Many interesting effects have been found by experimental techniques at temperatures well below the Debye temperature. In this temperature range quantum effects like zero point motions and corresponding delocalizations of atoms are important which have to be taken into account in serious theoretical studies. Usually this is done by lattice dynamics theories, in the framework of the harmonic or quasiharmonic approximation. However, both near second-order structural phase transitions and quite generally at higher temperatures the accuracy of this approach is sometimes uncertain , and methods that work at all conditions would be desirable.

PIMC simulations  enable us to analyse the crystal low temperature thermal properties. In principle, this method yields exact quantum-statistical averages (apart from statistical errors) and reduces to classical statistical averages at high temperatures. In general the agreement with experimental data is better compared to classical computations. For $`\beta `$–cristobalite it turns out that even at temperatures as high as 600 K only with PIMC a good agreement with experimental findings is obtained. The negative thermal expansion in $`\beta `$-quartz however may be explained by a classical modelling.

The outline of this paper now is as follows: in Sec. II we define the models that are studied, in Sec. III we briefly review the PIMC technique, while Sec. IV describes our results in detail. We conclude in Sec. V with a brief summary.

## II The model

Starting from the diamond lattice of silicon solids and placing between each neighbour pair of Si atoms an O atom, we arrive at a silicate which exists in nature in cubic symmetry: $`\beta `$\- cristobalite.

For the computation of SiO2–structures a variety of potentials are studied in the literature, for a review see ref. . Concerning the computational time it is of great advantage that apparently two-body potentials like the TTAM–potential  and the BKS–potential  describe the system properties at least as good as three body potentials. Thus most of the investigations presented here have been done with the TTAM– and the BKS– potentials. Comparative studies have been done as well with other potentials like the two–body contribution of a complicated potential obtained via ab-initio computations  (“core-shell”-potential).

In case of the TTAM and the BKS potentials the interactions between the Si and O– ions are modelled by pair potentials of the form:

$$\mathrm{\Phi }(r_{ij})=\frac{Q_iQ_je^2}{r_{ij}}+A_{ij}\mathrm{exp}(B_{ij}r_{ij})\frac{C_{ij}}{r_{ij}^6}$$
(1)
Here $`e^2=\frac{1602.19}{48.8542\pi }eV\AA `$. In Eq.(1) each ion carries a partial charge: A Si- ion has the charge $`Q=2.4`$ and an O- ion has the charge $`Q=1.2`$.

The TTAM potential has been developed by Tsuneyuki et.al. . The parameters in Eq.(1) are listed in table I. The parameters in Eq.(1) for the BKS- potential  are tabulated in table II.

In general complex materials with several components interacting with long ranged Coulomb interactions have to be handled by appropriate numerical procedures such as the Ewald summation technique. In the numerical simulation the computation of the Coulomb contribution to the energy requires a summation not only over the particle indices but also over all periodic boxes:

$$E_{coul}=\frac{1}{2}\underset{\stackrel{}{n}}{}\underset{\genfrac{}{}{0pt}{}{i,j}{(ij\mathrm{for}\stackrel{}{n}=\stackrel{}{0})}}{}\frac{Q_iQ_je^2}{|\stackrel{}{r}_{ij}+\stackrel{}{n}|}$$
(2)
The sum in Eq.(2) is over all periodic images shifted by $`\stackrel{}{n}`$ and all particle indices except $`i=j`$ for $`\stackrel{}{n}=\stackrel{}{0}`$. For numerical details of the computation of the Coulomb energy and how to use a relatively simple and efficient algorithm we refer to Ref.

Both potentials have in common that at small interparticle distances the atomic repulsion is neglected and rather a divergency in the potentials at small distances appears. In order to avoid a “Coulomb-collapse” of the system, at small distances we have to include an additional repulsive term. This term can be modelled as a harmonic contribution  suitable for molecular dynamics simulations or as a purely hard sphere repulsion. In this study for the sake of simplicity we follow the latter prescription, the hard sphere distance is chosen as the distance of the first potential maximum in Eq. (1).

## III Path integral Monte Carlo simulations

Canonical averages $`<A>`$ of an observable $`A`$ in a system defined by the Hamiltonian $`=E_{kin}+V_{pot}`$ of $`N`$ particles in a volume $`V`$ are given by:

$$=Z^1Sp[Aexp(\beta )].$$
(3)
Here $`Z=Sp[exp(\beta )]`$ is the partition function and $`\beta =1/k_BT`$ is the inverse temperature. Utilising the Trotter–product formula,

$$exp(\beta )=\underset{P\mathrm{}}{lim}(exp(\beta E_{kin}/P)exp(\beta V_{pot}/P))^P,$$
(4)
we obtain the path integral expression  for the partition function:

$$Z(N,V,T)=\underset{P\mathrm{}}{lim}(\frac{mP}{2\pi \beta \mathrm{}^2})^{3NP/2}\underset{s=1}{\overset{P}{}}d\{𝐫^{(s)}\}\mathrm{exp}[\frac{\beta }{P}(\underset{k=1}{\overset{N}{}}\frac{mP^2}{2\mathrm{}^2\beta ^2}(𝐫_k^{(s)}𝐫_k^{(s+1)})^2+V_{pot}(\{𝐫^{(s)}\}))]$$
(5)
Here, $`m`$ is the particle mass, integer $`P`$ is the Trotter number and $`𝐫_k^{(s)}`$ denotes the coordinate of particle $`k`$ at Trotter-index $`s`$, and periodic boundary conditions apply, the particle with Trotter-index $`P+1`$ is the same as the particle with Trotter-index $`1`$. This formulation of the partition function allows us to perform Monte Carlo simulations for increasing values of $`P`$ approaching the true quantum limit for $`P\mathrm{}`$. Note that Eq. (5) does not take into account any quantum- mechanical exchange between particles (for atoms as heavy as Si and O this is an excellent approximation, though it would not work for He crystals).

Thermal averages in the ensemble with constant pressure $`p`$ are given via the corresponding partition function

$$\mathrm{\Delta }(N,p,T)=_0^{\mathrm{}}𝑑V\mathrm{exp}[\beta pV]Z(N,V,T)$$
(6)
. In order to make our results comparable with those of experiments, our PIMC simulations have been performed in the constant pressure ensemble. Most of our simulations  have been done with $`N=684`$ particles ($`\beta `$–cristobalite), $`N=576`$ ($`\beta `$–quartz), $`p=0`$, and $`P`$– values up to 100. A typical PIMC data point in figure 2 required about 1200 CPU hours on a CRAY–T3E (single processor).

## IV Results

### A Finite size effects and pair potential sensitivity

In order to analyse effects due to the finite system size simulations for the system in the $`\beta `$– cristobalite structure have been done using the three potentials described above and $`N=648`$ particles as well as with the BKS– potential and $`N=192`$ particles. In Figs. 1 and 2 the potential energy and the systems box volume is shown.

The energy obtained by the core-shell potential (in the two-body reduced form) is significantly smaller than the energy obtained by the other two potentials. The volume obtained by the core-shell potential is much larger than the experimental data  and the thermal expansion coefficient is negative over the entire temperature range studied- in contrast to the experimental results for $`\beta `$-cristobalite.

For the BKS– potentials the successful description of SiO2-melt properties  using a cut-off distance of 5.5$`\AA `$ for the short ranged interactions suggested the choice of this cutoff- distance in the present study as well. From figure 1 however it is apparent that the volume decreases strongly with increasing temperature- in contrast to the experimental data. For the ordered crystal it turned out that a larger cutoff value (12 $`\AA `$) is more adequate for the description of the system properties.

The energies obtained by the TTAM- potential and by the BKS- potential (for both system sizes) agree well. In case of the volume however the results for the BKS- potential agree for both system sizes within numerical scatter and show a qualitatively similar temperature dependence as the volume obtained by the TTAM- potential. The quantitative values of the volume for the TTAM- potential however are much larger than the experimental data, whereas the volume obtained by the BKS- potential agrees well with the experimental data for temperatures above 700 K.

In summary the best agreement with the experimental data can be obtained with the BKS- potential - at least in a temperature range above 700 K. For smaller temperatures apparently none of the classical computations can successfully describe the experimental data.

### B Quantum simulations for $`\beta `$\- cristobalite

#### 1 Potential energy and volume

A different picture emerges if results of quantum simulations are considered in addition (BKS- potential, $`N=648`$), see figure 2. At temperatures below 700 K the quantum results for the potential energies deviate from the classical results to higher values, but for higher temperatures they agree with the classical values. The volume of the simulation box obtained by the PIMC simulations (Fig. 2) agree within numerical scatter with the experimental data for all temperatures studied- in contrast to the classical simulations. Despite the apparently high temperatures of 700 K quantum effects are very important for the thermal properties of $`\beta `$-cristobalite and thus quantum effects should not be neglected in simulations of real materials. This supports the assumption that quantum effects may still play an important role above the Debye temperature which is at about 500 K for SiO2.

#### 2 Structural details

According to a model of Wyckoff for $`\beta `$\- cristobalite the O- atoms are located in the middle of a bond connecting two Si- atoms, the Si–O–Si bond angle is 180o and the Si-O bond length 1.54 $`\AA `$. In order to correct these values to obtain the actually observed angle (between 140o and 150o) and bond lengths (Si-O: about 1.61 $`\AA `$) several models have been suggested in which the O- positions deviate from the ideal middle places . The question which of these models describes the properties of $`\beta `$\- cristobalite best is not yet decided by experimental methods.

In figure 3 we show the bond length distribution function for the different atom pairs for three temperatures using the BKS potential and in figure 4 the average bond lengths as a function of temperature. Note that due to the choice of a finite size cubic simulation box with periodic boundary conditions, $`\beta `$-cristobalite is stable (or metastable, respectively) in the simulation over the entire temperature range, unlike experiment.

In the PIMC simulations (with P=30) the bond lengths have larger values at low temperatures compared to the classical case (P=1) (Fig. 4). The quantum results for the Si-O bond length are in much closer agreement with the literature value of 1.61 $`\AA `$. The Si-O as well as the O-O bond lengths scale approximately linearly with the temperature, the Si-Si bond length decrease is stronger in parallel to the temperature dependence of the volume. This may be due to the fact that the positions of the Si-atoms fix the lattice constant of $`\beta `$\- cristobalite and the O- atoms have a much larger mobility in the crystal lattice. In figure 4 we find again the different behaviour of the Si-Si bonds compared to the other bonds. The temperature dependence of the Si-O and the O-O bonds have the following features:

* The width of the distributions increases with temperature, the height decreases.
* The average values increase with the temperature.
* At the same temperature the distributions obtained by PIMC (P=30) have larger average values than the classical distributions.
* The difference in the distributions and the average values between the classical and the quantum results decrease with the temperature.

All these effects are present in the Si-Si bonds as well, but there is a significant difference: In Figs. 3 and 4 we find that for the other two bond types the differences between the classical and the quantum values at the temperature $`T`$ are always smaller than the difference between the classical values at $`T`$ and the next higher temperature studied. For the Si-Si bond the quantum distributions at 400 K are closer to the distributions at 1600 K (which do not differ very much any more) than to the classical distributions at 400 K. This shows that in particular the Si-Si bond properties are significantly determined by quantum effects.

In order to obtain more information on the long-range molecular structure, the radial distribution functions $`g(r)`$ were computed up to a distance of 10 $`\AA `$ for Si-Si, Si-O, and O-O bonds. The knowledge of $`g(r)`$ allows the calculation of the so-called total correlation function $`T(r)`$, which is a superposition of the various $`g(r)`$’s.  $`T(r)`$ has been measured recently by neutron diffraction.  The measurements allowed to accurately estimate the real Si-O bond lengths and it was concluded that the typical Si-O bond length is considerably longer than the Si-O bond length in the so-called ideal $`\beta `$-cristobalite geometry. In Fig. (5a), a comparison between the experimentally observed and the calculated function $`T(r)`$ is shown. The agreement between the two curves is at least semi quantitative and the conclusions drawn in Ref.  could be confirmed with the exception of the width of the first Si-O peak and the first O-O peak. These peaks are broadened in experiment with respect to simulations, which include quantum effects. The broadenings might be due to lattice defects or to internal surfaces in the experimental samples. Both kind of defects were absent in our simulations. Similar results on $`T(r)`$ were reported previously by purely classical simulations  on the basis of the TTAM potential.

Fig. (5b) allows to relate the features in $`T(r)`$ to the various $`g(r)`$’s. It is interesting to note that the peak in $`T(r)`$ at about $`5\AA `$ is due to a simultaneous maximum in $`g_{\mathrm{O}\mathrm{O}}(r)`$ and $`g_{\mathrm{Si}\mathrm{O}}(r)`$, while the peak in $`T(r)`$ at $`6.25`$ $`\AA `$ does not have a corresponding peak in any $`g(r)`$. The large local maximum in $`T(r)`$ at about $`9\AA `$ is then located at a position where all $`g(r)`$’s have a local maximum as well. The shape of $`T(r)`$ is particularly sensitive to details of the potentials at large distances, e.g., cutting off the short-range part of the potential at $`4.5\AA `$ alters $`T(r)`$ significantly for $`r8\AA `$.

### C Quantum simulations for $`\beta `$-quartz

The $`\beta `$\- quartz structure differs from the $`\beta `$\- cristobalite structure by the positions of the atoms in the unit cell. $`\beta `$\- cristobalite has a cubic structure, $`\beta `$\- quartz has a hexagonal structure and at a temperatures below 846 K a trigonal $`\alpha `$\- quartz is stable . For structural details of the $`\beta `$\- quartz see appendix A. The piezoelectricity (and pyroelectricity) of $`\beta `$-quartz is at least reduced for $`\beta `$\- cristobalite due to symmetry reasons.

As for $`\beta `$\- cristobalite the simulations have been done for Trotter values of P=1 and P=30 at temperatures between 600 K and 1600 K. Generalizing Eqs. (3),(5) and (6), in our constant pressure simulations a hydrostatic pressure $`p`$ has been applied and fluctuation trials of the box lengths in $`x`$-, $`y`$\- and $`z`$\- directions have been attempted independently keeping the angles between the box axes constant. Lattice constants ($`a`$ and $`c`$) have been obtained from the resulting average box geometry. The experimentally observed transition at 846 K is suppressed in the simulations since we keep our box in rectangular shape. A trigonal $`\alpha `$\- quartz structure thus cannot be stabilised, rather a metastable, undercooled $`\beta `$-quartz is present in the simulations. Above the phase transition temperature, the structure in the simulation agrees with the structure deduced from the experimental data.

In figure 6 the potential energy and lattice constants of $`\beta `$-quartz are compared to experimental data . The energy has a linear temperature dependence in the classical case and slightly higher values in the quantum simulations, the difference to the classical values increases the lower the temperature. Apparently in this temperature region quantum effects become important- as can be expected due to a Debye- temperature of 500 K for SiO2. The temperature dependence of the volume only qualitatively agrees with the experimental data. In figure 6 we note the qualitatively different behaviour between the experiments and the simulation data at temperatures close to the $`\alpha \beta `$ transition at 846 K. As explained above this different behaviour is caused by the box shape constraint in the simulation.

As in the experiment the lattice constant in a- direction is approximately temperature independent at high temperatures whereas in c- direction a negative thermal expansion coefficient is found. In the case of the lattice constant $`c`$, the PIMC results are in somewhat better agreement with the experimental data than the classical results. Since for $`\beta `$-cristobalite it turned out that by using the literature value (5.5 $`\AA `$) for the cutoff the agreement with the experimental data is not as good as without the cutoff (see above), we chose the maximum cutoff value which is possible in our box geometry. The deviation of 1% (c-direction) and 0.3% (a- direction) from the experimental values may be due to the choice of parameters in the interaction potential, in particular to the cutoff of the short ranged interactions in the BKS- potential. Probably a restriction to smaller cutoff values would result in a better agreement with experimental data.

A typical simulation data point in figure 6 also required a computational effort of about 1200 CPU hours on a CRAY-T3E computer (in units of single-processor time).

Surprisingly the c- lattice constant obtained by a classical simulation decreases slightly with increasing temperature. This behaviour- which is present in the experimental data as well- may thus be explained by a purely classical argument. In the literature the ‘rigid unit modes’  are suggested as being responsible for the negative thermal expansion coefficient. For the sake of a more complete discussion of our results, we recall this argument here: The idea is that in a covalently bonded crystal structure one can identify units which are essentially rigid, because the length of the covalent bonds are essentially rigid, but that at certain atomic positions where the rigid units are linked together, they may rotate relative to each other. These modes may be visualised schematically as in figure 7 for two dimensional structures. The quadratic bases in this example may correspond to SiO4 tetrahedron in $`\beta `$\- quartz. If the angle stabilising forces in such units are significantly stronger than the forces at the common O- atoms connecting such tetrahedron, thermal excitation leads to a rotation of the complete bases. In contrast to the positive thermal expansion coefficient caused by the anharmonicities of the interatomic interactions such collective modes may cause a negative thermal expansion coefficient. The possible mechanism is briefly sketched: By tilting of the unit cell by the angle $`\mathrm{\Theta }`$ the system volume $`A=A_0`$ (see Fig. 7) changes to $`A=A_0\mathrm{cos}^2(\mathrm{\Theta })`$. For small angles and harmonic rotational oscillations with frequency $`\omega `$ of the units (i.e. tetrahedra) with moment of inertia $`I`$ the resulting thermal average of $`A`$ is :

$$A(T)=_T=A_0(1<\mathrm{\Theta }^2>_T)=A_0(1\frac{k_BT}{I\omega ^2})$$
(7)
The thermal average of the tilt angle $`<\mathrm{\Theta }^2>_T`$ increases with temperature and thus leads to a decrease of system volume. At low temperatures this effect may more than compensate the ‘normal’ thermal expansion caused by the motion of single atoms such that the thermal expansion coefficient is negative.

#### 1 Structural details $`\beta `$-quartz

In order to get more information on the molecular structure, the radial distribution functions $`g(r)`$ were computed for $`\beta `$-quartz, again up to a distance of $`10\AA `$. The correlation function $`T(r)`$, which follows from the various $`g(r)`$’s, is shown in Fig. (8a) together with $`g_{\mathrm{Si}\mathrm{Si}}`$, $`g_{\mathrm{Si}\mathrm{O}}`$, and $`g_{\mathrm{O}\mathrm{O}}`$ in Fig. (8b).

It is instructive to compare Fig. (8) for $`\beta `$-quartz with the corresponding figure for $`\beta `$-cristobalite, Fig. (5). The location of the peaks in $`g(r)`$ are nearly identical for nearest neighbor Si-Si, Si-O, and O-O peaks, only the peaks are broadened in the case of $`\beta `$-quartz, which can entirely be attributed to the larger temperature of the $`\beta `$-quartz simulation. It can be concluded that the tetrahedra formed by SiO4 units are basically identical for both lattice structures. While the radial distribution functions differ markedly between the two lattices for distances $`r>3.5\AA `$, there is a much less pronounced, yet noticeable, effect in $`T(r)`$. E.g., $`T(r)`$ shows a shoulder for $`\beta `$-quartz at $`r3.7\AA `$ which is absent in the simulation of $`\beta `$-cristobalite.

### D “Size” of the quantum particles

By determining the average radius of the quantum-particles across the quantum chain their “size” can be mapped as a function of temperature. In figure 9 we compare the simulation results for these radii in the various lattice structures with the thermal de-Broglie wave lengths for free O- and Si- atoms. Apparently this quantity does not depend much on the particular lattice structure. This may be due to the fact that the quantum chain radius is determined by the close surrounding of the particles and the potentials. In both crystalline structures studied here the basic unit is a SiO4 tetrahedron and the BKS potential is used. The de-Broglie wave lengths for the free particles are much larger than the radii of the particles of the crystals. This may be due to the stronger binding and localisation of the latter.

## V Summary

In this paper we present the results of our PIMC simulations for crystalline silicates. With PIMC simulations we are able to analyse the thermal properties of crystals at temperatures at which quantum effects have to be considered as well as anharmonic effects. In general the agreement with experimental data is better compared to classical computations. For $`\beta `$–cristobalite it turns out that even at temperatures as high as 600 K only with PIMC a good agreement with experimental findings is obtained. The negative thermal expansion in $`\beta `$-quartz however may be explained by classical effects.

We conclude that for a realistic materials modelling the simulation by path integral Monte Carlo techniques is a useful method which should be used in future studies of thermal properties of other crystals as well.

###### Acknowledgements.

C.R thanks the SCHOTT– Glaswerke for support and P.N. thanks the DFG for support (Heisenberg foundation and SFB 513). The computations were carried out on the $`CRAYT3E`$ of the $`HLRZ`$ at Jülich.
## VI Appendix A

In most simulations a rectangular simulation box is used. It is thus of advantage to map the hexagonal system onto a larger rectangular box such that the periodic order results in the hexagonal structure.

The lattice structure for $`\beta `$\- quartz is sketched in Fig. 10 and may be described as follows. In the space group P6222 the unit vectors are:

$`\stackrel{}{e}_1`$ $`=`$ $`(4.9977,0,0)`$ (8)
$`\stackrel{}{e}_2`$ $`=`$ $`(2.49885,4.328135,0)`$ (9)
$`\stackrel{}{e}_3`$ $`=`$ $`(0,0,5.4601)`$ (10)
Thus the lattice constants are:

$`a`$ $`=`$ $`4.9977\AA `$ (12)
$`b`$ $`=`$ $`\sqrt{3}a`$ (13)
$`c`$ $`=`$ $`5.4601\AA `$ (14)
The coordinates of the particles in a unit cell are given in table III, the Cartesian coordinates in table IV.

Due to the differences in box lengths one has to reduce the effect of different surface areas on the simulation results. This is done by matching sufficiently many unit cells such that the simulation box is approximately cubic.


format: text

text:
# The Relative Star Formation Histories of Spiral Bulges and Elliptical Galaxies in the Hubble Deep Fields

## 1 Introduction

The origin of galactic bulges is an important and unsolved issue in our understanding of galaxy evolution. The hitherto established view is that galactic bulges form at high redshifts through early dissipationless collapse (Eggen, Lynden-Bell, & Sandage 1962, hereafter ELS). This is based principally on the evidence for old stellar populations concentrated in the bulge of our own Galaxy and contrasts with more recently-developed hierarchical galaxy formation models (Kauffmann et al 1993, Baugh et al 1998) where elliptical galaxies form from the merger of early disk systems which can, in turn, continue to accrete gas to form a two component spiral galaxy.

A robust prediction of all hierarchical models is that spiral bulges should, on average, contain older stars than their associated disks which form by subsequent accretion. Moreover, statistically at a given redshift, bulges should be older and redder than field ellipticals which predominantly form from the merger of previously created spirals. Importantly, these conclusions should remain valid regardless of the particular cosmological model or initial power spectrum which governs the rate of assembly of massive galaxies. As such, a comparison of the relative colors of ellipticals and spiral bulges offers a remarkably simple, but powerful, test of hierarchical assembly models.

A third alternative for the origin of stellar bulges proposes their manufacture via various instabilities of pre-existing disks (see the recent review by Combes 1999). Of greatest interest here is the suggestion that bulges may form through secular processes (those whose timescale are longer than the dynamical time), mostly through interactions and the evolution of galactic bars. Galactic bars can be transformed into bulges through a vertical heating of the inner disk via resonant scattering of stellar orbits by the bar potential. Numerical work also suggests that bulges may form via bar destruction through the growth and subsequent collapse of bar instabilities in cold rotating disks. Evidence that such secular activity may play an important role in the formation of bars is growing (Kormendy 1993; Kuijken & Merrifield 1995; Norman, Sellwood, & Hasan 1996) and there is rather good dynamical evidence that most bulges in very late-type galaxies are not true bulges, but rather vertically heated central disks (Kormendy 1993). This suggests that the relative importance of the three bulge formation processes may be a strong function of Hubble type.

Most investigators of bulge formation have sought to infer their origin by studying the stellar populations or dynamical properties of local galaxies. Of particular importance in this regard are extinction-free near-infrared imaging observations of the Milky Way bulge (Dwek et al 1995) and resolved studies of associated stellar populations (eg. Rich 1997). These studies and those possible for nearby galaxies have emphasized a much greater diversity in color and stellar history than had hitherto been believed (Wyse et al 1997) although the effects of dust and other factors remain unresolved (Peletier et al 1999).

Through the Hubble Deep Field (HDF) campaigns (Williams et al 1996, 1999), high resolution imaging data is now available with sufficient signal-to-noise to enable the direct study of bulge growth and evolution in situ to redshifts $`z1`$. Given the above discussion, it is clear that a detailed observational comparison of the resolved stellar populations of field spirals and ellipticals as a function of redshift will offer important insight into some of the most fundamental aspects of galaxy formation. Indeed, only via intermediate redshift data can the epoch corresponding to the formation of the bulge component be determined independently of the age of its stars.

In an earlier paper (Abraham et al. 1999a) we investigated the importance of secular processes in forming bulges by determining the proportion of barred spiral galaxies in the Hubble Deep Fields as a function of redshift, finding a tentative decline in the fraction of barred galaxies beyond redshift $`z0.5`$. The physical mechanisms responsible for such a deficit of barred spirals at high redshift remain unclear. Possibilities include dynamically hotter (or increasingly dark-matter dominated) disks or an enhanced efficiency in bar destruction at high redshift. Regardless, the absence of bars at high redshift does seem to pose an observational challenge for models where bulges form via secular activity initiated by bars. Such activities may be episodic (Sellwood 1999), or confined to low-luminosity late-type spirals (Kormendy 1993) which are too faint at high redshifts to have been included in our earlier investigation.

In a precursor of this study (Abraham et al. 1999b), the dispersion in the internal colors of a sample of high-redshift HDF spirals with redshifts was used to investigate the relative colors and star-formation histories of bulges and disks at high redshifts. This work concluded that bulges are redder and older than their surrounding disks in essentially all luminous high-redshift spirals. However, while bulges appear redder than their surrounding disks out to high redshifts, no investigation of the relative colors of bulges and ellipticals as a function of redshift was undertaken at that time.

Such a comparison forms the basis of this paper. Here we continue our HDF investigations concerning the origin of galactic bulges by measuring the relative colors of spiral bulges and ellipticals as a function of photometric and spectroscopic redshift. Using automated morphological classifications, our study includes a total sample of 95 spirals and 60 early-type galaxies with $`I_{814}`$$`<`$24 mag in the Northern and Southern Hubble Deep Fields. An additional aspect of this work is the augmentation of our optical data by near-infrared NICMOS observations (Dickinson 2000) of a subset of our data, extending coverage into the $`J_{F110W}`$ and $`H_{F160W}`$ photometric bands in the Northern Hubble Deep Field. These infrared observations are neither as deep nor as well-sampled as the optical data in the HDF-N, but provide a number of interesting constraints on the conclusions derived from optical photometry alone.

A plan of this paper follows. In §2 we describe our overall sample selection and the techniques used to assign redshifts to these galaxies. In §3 we describe our technique for estimating the relative colors of ellipticals and bulges via the measurement of a central aperture $`V_{606}I_{814}`$ color. We quantify the uncertainty in these estimates for the bulges in our sample by examining individual cases and applying a similar technique to the local spiral galaxy sample of de Jong (1996). We discuss the bulge-elliptical comparison in the context of hierarchical models in §4 and summarize our main conclusions in §5.

## 2 The HDF Sample

As in Abraham et al (1999a,b), we take advantage of the unique nature of both Hubble Deep Fields (HDFs, Williams et al 1996, 1999). The long exposures in these observations guarantees substantially improved signal to noise in several passbands over earlier multi-color HST imaging surveys such as the Medium Deep Survey (Ratnatunga et al 1999). This is particularly important in the reliable separation of ellipticals from spiral bulges with low surface brightness disks. The 0.04 arcsec/pixel sampling of the drizzled optical data makes isolation of the bulge comparatively easy, as described below (§3). By combining both HDFs, a larger sample is produced which is essential for statistical comparisons such as the color evolution of bulges and ellipticals at a given redshift.

For the Northern HDF field, we supplement our optical data with near-infrared NICMOS data. (The characteristics of these data are described in Dickinson 2000, and are only summarized here.) Although the HDF-N was observed in the near-infrared from the ground in several programs, the depth and angular resolution of these data ($``$1 arcsec FWHM) are a poor match to that of the optical WFPC2 images. The present NICMOS HST data was obtained by mosaicing the complete HDF-N with a mean exposure of 12600s per filter in F110W (1.1$`\mu `$m) and F160W (1.6$`\mu `$m). The sensitivity varies over the field of view but the mean depth is $`AB`$26.1 at 10$`\sigma `$ within a 0.7 arcsec diameter aperture. The drizzled NIC3 PSF has a FWHM = 0.22 arcsec and is primarily limited by the NIC3 pixel scale.

The $`I_{814}`$ =24 mag limit adopted for the present morphological study is a conservative one. For example, Ratnatunga et al (1999) claim reliable morphological classifications are possible to $`I_{814}`$ =26 mag in the northern HDF. An important consideration here is the reliable isolation of the bulge component at high redshifts. We defer discussion of this point to §3. However, we note at this point that bulges becomes difficult to reliably detect beyond a redshift $`z`$1. As a significant fraction of galaxies beyond $`I_{814}`$$`>`$24 have redshifts beyond unity, this is an additional justification for our adopted magnitude limit.

The baseline photometric catalogs for the present study are the publicly available SExtractor optical catalogs produced by the Space Telescope Science Institute (STScI). Throughout the present paper, we will adopt IAU-format object designations and coordinates keyed to these catalogs. To $`I_{814}=24`$ mag, both HDF images contain 242 galaxies of all types. Morphological classification of this sample was undertaken using the automated classification technique (based on measures of central concentration and rotational asymmetry) described in Abraham et al. (1996a,b), supplemented by visual classifications determined by two of the authors (RSE and RGA). This is the same strategy adopted in our analysis of bar structure in the HDF fields (Abraham et al 1999a). The galaxies in our sample are sufficiently bright that there is no ambiguity in matching galaxies to their counterparts in the corresponding NICMOS image of the northern HDF.

Table 1 summarizes the properties of the full sample. The table records integrated and central optical and near-IR colors as defined later in the paper, quantitative structural measures (central concentration, rotational asymmetry, axial ratio, and automated classification), and visual classifications according to two of the authors (RSE and RGA). Visual classifications are on the MDS numerical system: -2=star, -1=compact, 0=E, 1=E/S0, 2=S0, 3=Sab, 4=S, 5=Scdm, 6=Ir, 7=peculiar, 8=merger, 9=defect. A measure of the uncertainties for the purpose in hand (selecting spirals and ellipticals) can be gained from comparisons of both visual classification schemes with that determined automatically.

Automated classification results in a sample of 95 spirals and 60 early-type systems, the remainder being irregulars, mergers or peculiars. Visual classification by RSE (RGA) results in a total of 96 (123) spirals and 74 (64) early-types. In the present paper visually classified intermediate type systems, ie. S0/a galaxies, were categorized as early-type. (Nine galaxies in total were classed as S0/a systems, and we will demonstrate in §4 that our results do not depend sensitively on how they are partitioned).

The RSE classes are in excellent agreement with those determined automatically for spirals but a 23% excess of ellipticals is found. Conversely, RGA finds good agreement with the automated scheme for ellipticals but has a 29% excess of spirals. This level of discrepancy is consistent with earlier comparisons between independent visual classifications of faint galaxies (Abraham et al. 1996b). Concerning the apparent discrepancies for spirals, it should be noted that the automated classifier is tuned to categorize weakly distorted systems and very late-type spirals as peculiar. The excess of visually-classified ellipticals is thought to arise because some visual classifiers weight global symmetry more than central concentration (Marleau & Simard 1999).

The most important consideration here is the selection of reliable samples of ellipticals and well-defined spirals (the latter from which the bulge properties can be determined). In the case of the ellipticals, we follow precisely the selection criteria of Abraham et al (1996,1999a) for the HDF-N and its equivalent procedure for the HDF-S (Menanteau et al 2000). Discussion of possible contamination of these elliptical samples by compact galaxies can be found in Menanteau et al (2000). In the present paper no attempt is made to separate S0s from ellipticals.

As described below, we will restrict our spiral sample to those non-peculiar systems with prominent bulges. For this subset both visual and automated classifications give virtually identical results so far as locating spirals is concerned. Regular spirals and those with visible bulges were ascertained by visual inspection by both observers. However, one of us (RGA) examined those spirals he classed in more detail for blending effects and these more detailed visual morphological classifications (summarized in Table 1) form the basis of the remainder of the analysis. The principal goal here is to enable selection of various (presumed “clean”) subsets of spiral bulges for later analysis. We will show in §4 that the conclusions of this paper do not depend sensitively on the exclusive use of any one of these subsets.

For the RGA-classed spirals, the column denoted “bulge?” indicates the presence of a bulge component. The column denoted “blend?” indicates whether nearby companions or optical superpositions overlap with a galaxy image at a 2$`\sigma `$ isophotal threshold on the $`I_{814}`$ -band image. Blends are categorized as either “severe” or “minor”. Severe blends are cases where the overlapping images cannot be reliably disentangled. The distinction between a “severe” blend and a merger is obviously rather subjective. In some cases tidal tails make the latter classification obvious, while conversely “severe blending” is appropriate when contamination originates from, e.g. a diffraction spike presenting photometric difficulties. The “minor” blend category encompasses systems where the overlapping images can be cleanly separated by raising the $`2\sigma `$ isophotal threshold, used to define the galaxy above the sky background, by a small amount (typically an additional 1–3$`\sigma `$ above the sky noise level). The tabulated photometric parameters for galaxies flagged as minor blends in Table 1 correspond to the properties determined at this higher threshold. Throughout the remainder of this paper these minor blends will be treated in the same manner as isolated galaxies. The flag, “pec?” indicates the visual presence of an optical peculiarity that is insufficient to warrant classifying the galaxy as a peculiar, eg. a spiral that could reasonably be classed as Spec in local catalogs.

Figure 1 (Plate 1) presents a montage, ordered by $`I_{814}`$ magnitude, of the 68 isolated non-peculiar systems classified as spirals by RGA whose bulge components are clearly visible, ie. systems selected from Table 1 using the following criteria: (RGA = Sab, Sbc, or Disk)bulge? = yespec? = no, and blend? $``$ severe. We note that the prominent bulges and fairly tight spiral structure visible on many of the spirals shown Figure 1 indicates that the majority of systems are earlier than type Scd. This is of course consistent with our initial pre-selection designed to isolate spirals with prominent bulges.

Galaxy redshifts are listed on each panel of the montage. Spectroscopic redshifts become progressively fewer fainter than $`I_{814}`$ =22.5 mag in the HDF-N and are still rather scarce in the HDF-S. We have therefore augmented the spectroscopic redshift estimates with photometric redshifts based on the publicly available catalog produced by Gwyn (http://astrowww.phys.uvic.ca/grads/gwyn/pz/index.html). These photometric redshifts are based on optical photometry only and do not include broad-band infrared colors. However, the summary of results in Hogg et al. (1999) suggests that infrared photometry does not significantly improve the accuracy of the photometric redshifts in the HDF for the $`z<1`$ range of importance in this work.

## 3 Photometric Decomposition

Many studies of high-redshift bulges have attempted to decompose galaxian light into components by fitting analytical models to the surface brightness distributions of galaxy images (Schade et al 1995, Lilly et al 1998, Marleau & Simard 1998, Ratnatunga et al 1999). In this procedure, all pixels in the image are used to determine simultaneously the global structural parameters for all components of the galaxy, such as the scale lengths, ellipticities, and characteristic surface brightnesses for the bulge and disk. By integrating over the results from such fits, the colors of galactic components can be inferred.

Obtaining robust model fits using this technique can be challenging and, at a more fundamental level, the approach assumes the validity of canonical fitting laws and a close correspondence between stellar populations and global galactic structure (a relationship that is likely to vary reasonably strongly with Hubble type). Both of these may be poor assumptions, particularly for later Hubble types (Kormendy 1993). By utilizing all pixels in the image to determine model parameters this technique clearly optimizes signal-to-noise, but of course this advantage can also prove a liability, as a small number of pixels affected by features not incorporated into the models (such as dust) can have a global impact on the fits. Furthermore, the morphological irregularity of high redshift galaxies suggests that fitting smooth light distributions to distant galaxies may introduce spurious results if care is not taken to first isolate underlying smooth features from superposed asymmetrical structures. This idea is very difficult to implement in a scale-free manner, although some attempts based on enforcing rotational symmetry seem quite encouraging (Schade et al. 1995; Lilly et al. 1998), and alternative fitting methods which attempt to minimize the effects of irregular features by adopting particularly robust minimization strategies (eg. simulated annealing) also show promise (Marleau & Simard 1998).

In the present paper, we choose to adopt a simpler approach which foregoes any attempt to determine the global structural parameters of the galaxies in our sample, in order to focus specifically on bulge stellar populations in a robust non-parametric manner. We estimate bulge colors by measuring the colors of the innermost regions of galaxy images where the bulge component dominates. This approach is conceptually similar to that adopted by Peletier et al. (1999) in determining bulge colors for local spirals using HST NICMOS observations, although the details (described below) differ.

Ideally an aperture selected to locate the bulge should be defined so as to be linked physically to the intrinsic properties of the galaxy. In this way correlations with, e.g. redshift, disk surface brightness and luminosity would be minimized. However, such an approach suffers from many of the difficulties associated with a full decomposition (see above) and, in practice, when observing with an instrument (WFPC-2) with a fixed pixel size, turns out to be very difficult to arrange. A metric aperture cannot easily be chosen to simultaneously avoid disk contamination at low z whilst containing an adequate number of HST pixels (from signal/noise considerations) at large redshift.

Instead we adopted an aperture size which scales with the isophotal size of the galaxy. This has the merit of great simplicity, the argument being that any systematic effects can be explored through careful simulations based on nearby galaxies. A potential drawback of this approach is that the nuclear colors so defined may be contaminated by disk light in a manner that is a function of redshift (assuming disks are generally bluer than bulges at high redshifts, cf. Abraham et al. 1999b). We show below from detailed simulations that, for the nuclear regions of low-inclination spirals with morphological $`T`$-types $`T<6`$ (ie. Sc and earlier), bulge light dominates over disk light out to quite high redshifts in both $`V_{606}`$ and $`I_{814}`$ filters, so that the central $`V_{606}I_{814}`$ color traces bulge colors with an RMS accuracy of $`\sigma _{VI}<0.15`$ mag out to $`z<0.7`$.

Further tests we introduce to examine possible systematic effects are the analysis of subsets of our sample of spirals with prominent bulges defined in the previous section, restricted according to disk inclination and comparisons between optical and optical-infrared colors since disk contamination should be minimal in the latter case.

### 3.1 Methodology

In this section we describe the procedure used to analyze our primary WF/PC2 optical dataset. The analysis of our supplementary NICMOS observations of the HDF-N field proceeded along essentially identical lines. However, it is important to make clear at this stage that the under-sampling of the NIC3 detector relative to WF/PC2 observations, and the rather uncertain photometric calibration of the NICMOS instrument, make it substantially harder to apply the techniques described in this section to the NICMOS observations. For these reasons (and because the NIC3 data are restricted to HDF-N), in the present paper we will treat the NICMOS data as an interesting adjunct to our optical data that is mostly useful for testing and elaborating the conclusions derived at from data obtained at optical wavelengths, and defer further discussion of the NICMOS data until §5.4.

Our analysis is based on elliptical apertures defined using second order moments obtained from the $`I_{814}`$-band image of each galaxy in Table 1. These were used to define an ellipse whose area corresponds to that of the galaxy above a 2$`\sigma `$ isophotal threshold. Large aperture (‘total’) colors (see Table 1) were defined using this procedure applied to the $`I_{814}`$ image. For each galaxy a circular aperture was then defined with a radius of 5% of the semi-major axis length of the galaxy’s ellipse, and this inner aperture was then used to determine the $`V_{606}I_{814}`$ colors of the central portion of the galaxy.

The 5% isophotal radius aperture size was chosen according to the following considerations. Clearly the aperture should be as small as possible in order to isolate bulge light with a minimum of contamination from disk light. However the aperture must be large enough to encompass several pixels from signal/noise considerations. For a minimum of 10 pixels to be averaged, the aperture radius must be at least 2 pixels corresponding to an angular size of $``$0.08″using the dithered 0.04″pixels of the Version 2 HDF images released by STScI. Assuming $`H_0`$=70 km sec<sup>-1</sup> Mpc, $`\mathrm{\Omega }_M`$=0.3 and $`\mathrm{\Omega }_\mathrm{\Lambda }`$=0.7 (the cosmology adopted throughout the paper), 0.08″corresponds to 0.64 kpc at $`z`$=1.

The location of each aperture was automatically adjusted by the measurement algorithm so as to be positioned on the local centroid of the galaxy image. The central and total $`V_{606}I_{814}`$ colors for all galaxies in our sample (including ellipticals) are given in Table 1. The photometric zero points are given in Williams et al. (1996).

Initially we considered defining the bulge interactively for each galaxy by “growing” the bulge region out from the nucleus on a pixel-by-pixel basis, isolating the inner pixels of the galaxy whose resolved pixel-by-pixel color distribution was homogeneous at some pre-determined level. This methodology was found to work well for bright galaxies. However, as was found in Abraham et al. (1999b), the limiting magnitude for this approach in the HDF data is around $`I_{814}`$ =23.2 mag, motivating us to attempt to improve the signal-to-noise of the measurement by binning together central pixels using central apertures in order to reach a practical magnitude limit of $`I_{814}`$ =24 mag. Thus in the present paper the pure pixel-by-pixel approach is set aside in favor of a central aperture. It is important, nonetheless, to understand whether a 5% aperture radius allows us to chromatically isolate a typical bulge.

We can test this by examining the pixel-by-pixel optical colors in a number of bright galaxies as an independent check. Figure 2 shows how the bulge pixel colors are distributed with respect to those of the outer disk in 3 bright spirals with fractional aperture radii corresponding to 2%, 5% and 7% of the semi-major axis length. As can be seen, in all cases the bulge defined by a 5% radius is remarkably homogeneous and red compared to the pixels outside. In fact the relative colors of the spiral bulge and the elliptical in the third and fourth rows (seen at virtually identical spectroscopically-determined redshifts of $`z=1.013`$ and $`z=1.016`$) foreshadows our major result: the bulge in the spiral system is significantly bluer than the typical elliptical.

### 3.2 Possible Contamination of Bulge Colors by Disk Light

The suitability of deriving bulge colors using a constant fractional aperture can be further assessed by considering the spirals in the local sample of de Jong (1996; hereafter DJ96). In this section we utilize analytic profile fits to the DJ96 sample photometry and estimate, on their basis, the likely disk contamination expected when the DJ96 galaxies are artificially placed at various redshifts and analysed with our central aperture technique.

DJ96 presents fundamental disk and bulge parameters in $`BVRIHK`$ bands for 86 nearby low—intermediate inclination spiral galaxies as a function of Hubble T-type. The suitability of any currently-available local reference sample for the purposes of comparing with deep HST data is debatable, but the resolution and rest-frame depth of the data presented in DJ96 corresponds reasonably closely to that for deep images obtained with HST at high redshifts. The reader is referred to Lilly et al. 1998 (who also adopt DJ96 as the local calibration sample for a study of high redshift spirals) for a further discussion of this point.

It is also questionable whether an $`r^{1/4}`$-law profile or an exponential profile provides the most suitable description of local bulges. Following Frankston & Schild (1976) and Andredakis & Sanders (1994), DJ96 convincingly argue that local bulges are best fit by $`r^{1/4}`$ profiles in early-type spirals, and by exponential profiles in late-type systems, and presents fits based on both canonical profiles as a comparison. In order to undertake a more rigorous test of our methodology, and in order to estimate our measurement errors conservatively, in the present paper we will adopt the exponential-law profile fits given in DJ96. The sharp nuclear spike of the $`r^{1/4}`$ profile makes determining bulge colors on the basis of the technique described in the present paper more accurate.

Figure 3 illustrates the difference between the $`UV`$ and $`BR`$ bulge colors determined by 5%-isophotal radius aperture photometry and those determined by integrating over the DJ96 model fits. As de Jong (1996) does not tabulate model fits in $`U`$-band, we have assumed a constant rest-frame $`UB=0.12`$ mag for the disk component, and constant rest-frame $`UB=0.50`$ mag for the bulge component, corresponding to the mean colors for late and early-type galaxies in Table 1 of van den Bergh 1998. Note that $`U`$ and $`B`$ bands correspond roughly to rest-frame $`V_{606}`$ at redshifts of $`z0.4`$ and $`z0.7`$, respectively, so the left and right panels of the figure illustrate the expected errors introduced into our $`V_{606}I_{814}`$ bulge color measurements at these redshifts. On the basis of this figure, at $`z0.4`$ we can expect the RMS uncertainties on our bulge color estimates to be $`\sigma _{VI}<0.1`$ for spirals earlier than Sc, and $`\sigma _{VI}<0.15`$ for essentially all spirals earlier than type Sdm. At higher redshifts, blue rest-frame disk light may contaminate the bulge light to a larger degree, but even at $`z0.7`$ the RMS contamination is only $`\sigma _{VI}0.2`$ for spirals earlier than Sc, which we will show below is substantially smaller than the typical color difference between ellipticals and bulges at the redshifts of interest. Using NICMOS $`J_{110}H_{160}`$ colors, disk contamination is expected to be negligible at $`z<1`$.

### 3.3 Aperture Effects for the Elliptical Galaxies

In order to facilitate the simplest comparison between the colors of spiral bulges and ellipticals, we will utilize integrated (i.e. large aperture) colors for the field ellipticals. As described earlier, large aperture colors were derived within ellipses defined by the second order moment of each galaxy in the $`I_{814}`$ band. For completeness, Table 1 also lists 5% aperture photometry for these sources as this contains important information on the extent to which ellipticals represent a homogeneous population.

Abraham et al (1999b) found that as many as a third of the HDF ellipticals to $`I_{814}`$ = 23.2 have cores bluer than those of the galaxies’ periphery suggesting bursts of star formation involving at least a few percent of the galactic mass in the previous few Gyr. Because of this, the use of 5% aperture colors for the ellipticals would complicate any comparison with our sample of spiral bulges.

A complete discussion of the internal optical colors of the HDF elliptical population is presented in a separate paper (Menanteau et al 2000). Nonetheless it is interesting to examine the color inhomogeneities from Table 1 briefly in order to set a scale for the color variations we later discuss for spiral bulges. For the HDF ellipticals in the present study, the mean offset in color in the sense (aperture)-(total) is $``$0.08 mag with $``$15% showing no difference at all. Only 10% of the sample has offsets greater than 0.2 mag.

## 4 Relative Optical Colors of Bulges and Ellipticals

In this section we will investigate the optical colours of our complete sample of bulges and ellipticals, and defer consideration of the infrared subset of these data until §5.3.

The $`V_{606}I_{814}`$ colors of ellipticals and spiral bulges as a function of photometric redshift are compared in Figure 4. This figure also shows the predicted colors of a simple passively-evolving stellar population based on the spectral synthesis models of Charlot & Bruzual (1999). This baseline model prediction corresponds to a single 1 Gyr burst at high redshift ($`z_F`$=3) and assumes solar metallicity with a Scalo Initial Mass Function. In Figure 4 the symbol sizes are keyed to central concentration and it is seen that most spirals in the sample are indeed prominently nucleated early-types, as expected from our initial pre-selection of spirals with prominent bulges.

The prominent red locus defined by some (but not all) of the ellipticals in Figure 4 shows a redshift-dependent trend consistent with earlier work (Schade et al 1996, 1999, Brinchmann et al 1998) and, with some exceptions, is in excellent agreement with our simple passively-evolving model. As discussed by Abraham et al (1999b), a significant fraction of ellipticals lie blueward of the well-defined red-locus. Many of these early-type systems have bluer nuclei and are thought to be remants of recent merging providing good evidence for continued formation of spheroid galaxies through mergers of late-type galaxies (Menanteau et al 2000, Brinchmann & Ellis 2000).

Contrary to expectations, the bulges in our sample of regular, isolated spirals show a systematic blueward scatter relative to the red early-type locus. Whereas both bulges and ellipticals show a dispersion in color at a given redshift, the distributions are characteristically different. A much smaller fraction of the elliptical population is blue although there is strong hint of a bimodality in the distribution, consistent perhaps with bursts or star formation associated with recent merging (Menanteau et al 2000). By contrast, few of the bulges occupy the ‘red envelope’ region of the color - redshift relation as expected if they are intrinsically older than their elliptical counterparts. Moreover their distribution of colors is less bimodal and occupies the intermediate region, consistent on the one hand perhaps with less extreme evolutionary trends but on the other hand ones that are spread throughout the population.

Before accepting the surprising result that bulges are, as a population, systematically bluer than ellipticals in the mean, it is important to investigate possible sources of systematic error. The main source of systematic error is likely to be contamination from underlying disk light. As described in §3.2, the offset between typical bulge colors and the colors of red ellipticals is most significant for $`z0.5`$. At higher redshifts the bulk of the effect is unlikely to be explained fully by disk contamination unless the DJ96 is very unrepresentative of the high redshift spiral population. Other systematic effects can be best studied by introducing a measure of the color scatter of the bulges with respect to our baseline passively evolving high-redshift collapse model. We define a simple parameter characterizing the offset between a galaxy’s observed color and the predicted color of our baseline high-redshift burst model:

$$\delta (V_{606}I_{814})=(V_{606}I_{814})_{bulge}(V_{606}I_{814})_{passive}$$

In Table 2 we present $`\delta (V_{606}I_{814})`$, the median value of $`\delta (V_{606}I_{814})`$, and the corresponding RMS dispersion, $`\sigma _{\delta (V_{606}I_{814})}`$, for various subsamples of spiral bulge populations.

The distribution for the main sample (i.e. that plotted in Figure 4) is compared with that for the ellipticals in Figure 5 clearly illustrating the points made earlier. For this sample we obtain $`\delta (V_{606}I_{814})`$=-0.525 and $`\sigma _{\delta (V_{606}I_{814})}`$=0.266. The median colors are substantially bluer than their elliptical counterparts, and it is apparent they fill in the regime between the passive elliptical population and the bluest elliptical outliers. Particularly striking is the almost complete absence of systems as red as the passive-evolution track which provides such a good description for the majority of the early-type population.

Two important conclusions can be drawn from Figure 4 and Table 2. Firstly, as described in the §3.2, we estimate typical contamination of bulge colors from disk light to be $`V_{606}I_{814}`$$`0.15`$ mag for most early—intermediate spirals at redshifts $`z<0.5`$. The median offset for the main sample given in Table 2 is therefore substantially larger ($`\times `$3) than the expected RMS disk contamination for the galaxies in our spiral sample in this redshift range. At higher redshifts the significance of the blue bulge colors is hard to determine from optical data, but Figure 3 suggests that around a third of the spirals earlier than type Sbc at $`z<0.7`$ should show little disk contamination ($`V_{606}I_{814}`$ offsets less than 0.1 mag), although the optical colors of some systems will be significantly contaminated.

Secondly, in Table 2, the $`\delta (V_{606}I_{814})`$ color offsets and dispersions remain similar to the fiducial values for our main sample even if broadly different selection criteria are used to define the galaxy subsamples. Similar results are obtained when we: (a) restrict consideration only to large systems where the methodology is expected to be most effective; (b) restrict consideration only to low-inclination ($`\frac{b}{a}>0.5`$) systems; (c) restrict consideration only to subsets of galaxies with confirmed spectroscopic redshifts. The conclusions that most bulges are systematically bluer than the reddest ellipticals with a wide dispersion of colors is true regardless of the detailed prescription used to define the spiral subsamples.

## 5 Interpretation

Our principal result is that bulges are, statistically, optically bluer than the reddest ellipticals and show a large dispersion in their rest-frame colors. This result is insensitive to disk contamination, certainly for $`z<0.5`$ and probably to higher redshifts since only a very small proportion ($`<`$5%) of intermediate redshift bulges have colors as red as those ellipticals with occupy a typical passively-evolving track defined according to a single burst of star formation at high redshift.

In considering the interpretation of these results in this section, we must take great care to distinguish between three reasonable working definitions for the age of a bulge component: (i) optical luminosity-weighted age, (ii) structural age (corresponding to the epoch at which the bulge component is formed morphologically), and finally (iii) the age of the oldest stellar population in the bulge. In the context of our observations these are not merely semantic distinctions — in fact understanding which of these ages is being probed by various aspects of our data is central to linking our observations with the predictions of theoretical models.

Our results clearly indicate that the optical luminosity-weighted ages of bulges to at least $`z=0.5`$ are younger than those of the reddest ellipticals. Restriction of our sample to more robust (e.g. face-on) sub-samples does not change the strength of this conclusion. The pixel-by-pixel color distributions of representative spirals presented in Figure 2 gives a clear illustration of the overall result.

### 5.1 Constraints on Secular Models

We now consider the implications of younger luminosity-weighted ages for bulges from the viewpoint of models where bulges grow through secular evolution (Combes 1999). The paucity of bulges seen with optical colors close to the evolutionary track for a single burst of star formation seriously constrains a dominant population of old passively-evolving bulges. Similarly, the optical data is hard to reconcile with the notion that bulge formation proceeds secularly within the disk without associated star-formation in early–intermediate-type spirals. Further constraints on the secular picture depend upon whether the star-formation we infer from blue colors is associated with the formation epoch of the bulge itself, or with “pollution” by young stars upon a pre-existing, older, and morphologically established bulge component.

In the extreme case, if we assume that the blue colors are associated with the initial formation of the bulge, then conceivably some bulges formed at low redshift from the secular evolution of stellar disks with bars as an intermediate stage. Ellipticals could still form continuously from the merger of young disks at high redshift in which case their predominantly red colors in the redshift range 0$`<z<`$1 would be compatible with growth in a low density or $`\mathrm{\Lambda }`$-dominated cosmology (Kauffmann & Charlot 1998, Menanteau et al 1999, 2000). Given the small proportion of bulges as red as ellipticals in our redshift range, most bulges would have to have formed surprisingly recently and one would then expect a significant population of barred spiral precursors at higher redshift in contrast to the observations (Abraham et al 1999a). Interestingly, there is some evidence (c.f. Figure 2) that when bars are present, they do have redder colors than those of their bulges they contain as expected in the secular picture. However, numerically it seems unlikely that most present day bulges were formed in this way unless the number of $`z>`$0.5 barred systems has been seriously underestimated (c.f. Bunker 1999), or bulge formation is episodic, as suggested by Sellwood (1999).

### 5.2 Constraints on Hierarchical Models

Next we proceed to discuss the implications of our results in the context of hierarchical galaxy formation models. Figure 6 illustrates the difference in observed $`V_{606}I_{814}`$ colors expected for ellipticals, spiral bulges, and S0 galaxies in the $`\mathrm{\Lambda }`$CDM semi-analytical prescription of Baugh et al. (in preparation), as a function of both redshift and apparent magnitude. Contrary to what is observed, the hierarchical picture predicts a much tighter color distribution for both bulges and ellipticals with the former being $``$0.1-0.2 mag redder than ellipticals at high redshifts, and roughly similar colors at low redshift $`z0.4`$. In contrast to these predictions, we find that at $`z<0.5`$ (ie. where possible systematic effects introduced by disk contamination can be neglected) bulges are significantly bluer than ellipticals.

Whilst suggestive of a major shortcoming in the predictions of the hierarchical model, the severity of this failure depends once again on whether the luminosity-weighted ages probed can be associated with the true formation epoch. An interpretation of our data broadly consistent with the hierarchical framework would be that most bulges are indeed produced from minor mergers of earlier disk systems (and hence that the underlying stellar populations in bulges are older than ellipticals) but that other physical effects that are missing from the models also play an important role in stimulating central star-formation after the bulk of the bulge mass is already in place (eg. secular processes that stimulate a modest amount of ongoing star formation in the central regions).

Although our observations do not necessarily require a major revision of the merger-driven astrophysics underlying the hierarchical models, a significant puzzle is the greater scatter in the colors of bulges compared to those in ellipticals (for which similar merger-driven activities are presumably occurring). The observed differences might be understood in terms of the timescales involved (see discussion below), as well as by the preferential morphological selection of ellipticals as relaxed remnants (Menanteau et al 1999).

### 5.3 Models Incorporating Secondary Star Formation

The discussion in §5.1 and §5.2 suggests that the major issue for understanding the ramifications of our optical observations for both secular and hierarchical models is the relationship between the luminosity-weighted ages of bulges and their structural ages. The key to establishing this would seem to be the measurement of the third “age measure” referred to earlier: namely establishing the epoch of formation and mass associated with the initial bulge-formation event.

In Figure 4 we show optical color-redshift trajectories for passively-evolving stellar populations upon which we have superposed 15% of the stellar mass, either in the form of a short-lived starburst or as an extended star-formation episode (modeled as exponential star-formation with a 5 Gyr e-folding timescale). The latter star-formation activity may be associated with a quasi-continuous infall of gas. We illustrate models corresponding to the onset of star-formation activity at redshifts $`z=0.2,0.4,0.6,0.8,`$ and $`1.0`$.

A notable feature of the optical data is the bimodal nature of the deviations from the quiescent evolution track exhibited by the early-type galaxy population. Starburst models do not remain blue for much longer than the time associated with the active star-formation phase. At high redshifts, our observations probe sufficiently far into the ultraviolet that the observed $`V`$-band flux is dominated by current star-formation even for a moderate (10–15% by mass) burst. Unsurprisingly a “drip-feed” scenario results in a weaker but longer-lived blue phase suggesting that our present data can be used to place crude constraints on the duty-cycle of star-formation activity for both the bulges and early-type galaxies in our sample. The locus of very blue early-type systems combined with a paucity of intermediate-color early-type systems suggests that star-formation in high-redshift ellipticals is associated with rapid bursts whereas the broader color distribution observed for the bulges may be suggestive of a low-level, more extended star-formation.

As a quantitative example of the last point, consider two simple models: a passively evolving 5 Gyr old stellar population with roughly constant star-formation which then suffers a relatively minor nuclear starburst, and a pure starburst model with no underlying old stellar population. To be definitive, assume the starburst has an exponential star-formation history with an e-folding timescale of 100 Myr, and forms 10% of the mass of the underlying old population in the first model. Immediately after the burst phase, the composite (old + young) stellar population model at $`z=1`$ is only 0.12 mag redder in $`VI`$ than the pure burst model — at high redshifts, the young stellar population effectively masks the presence of any underlying old population in optical bands. However, the corresponding figure for $`IH`$ is 0.51 mag. Extending the present methodology to the infrared allows one to determine, in an individual case, whether a blue bulge is associated with initial bulge formation or incremental bulge growth.

### 5.4 Constraints from NICMOS Observations

To make further progress we now turn to the resolved near-infrared (IR) data made possible through NICMOS observations of the HDF-North. As stated earlier, analysis of these data proceeded in the same manner as the analysis of drizzled WF/PC2 optical data. However, because the sampling of the NICMOS data is relatively poor, care is needed to ensure precise co-registration of the optical and IR data. Prior to our analysis all optical data was first convolved with the $`H_{160}`$-band PSF and re-binned to match the characteristics of the NIC3 dataset. The entire analysis process described earlier was then repeated for the convolved and rebinned optical data as well as for the NICMOS data (using, once again, the $`I_{814}`$-band images to define all central apertures). In order to test whether the coarser resolution of the NIC3 data was a major source of error, we constructed a new version of Figure 4 using the convolved and resampled $`V_{606}`$ and $`I_{814}`$ images, and found that the trends remained similar.

Even though only half our sample has been imaged (and with coarser resolution) with HST in the infrared, the benefits NICMOS offers are considerable. Firstly, disk contamination is expected to be less important at longer wavelengths, so if we see the same trends as in Figure 4, then this provides very strong confirmation of differences in the evolutionary behaviour of bulges and ellipticals as discussed above. Secondly, the addition of the IR photometry can break some of the degeneracy in the modelling based on optical colors discussed above since stars of different main sequence timescales are involved. As emphasized earlier, these advantages are offset by the poorer angular resolution and sensitivity of the NICMOS data (as well as by possible uncertainties in the modelling of stellar population changes at near-IR wavelengths, as described by Charlot, Worthey, & Bressan 1996).

Figure 7 shows the $`J_{110}H_{160}`$ color-redshift relation for the HDF-N targets plotted in Figure 4 with deep IR imaging. The model tracks discussed earlier in the context of Figure 4 are also plotted for the IR photometric bands. As mentioned, the IR models are mainly used to indicate general trends and we have displaced the $`J_{110}H_{160}`$ BC96 model plotted in Figure 7 by +0.1 mag in order to ensure the passive evolution track lies closer to the locus defined by the redder E/S0 systems. We note that a $`J_{110}H_{160}`$$`=0.1`$ mag color offset is quite compatible with the systematic effects seen by Charlot, Worthey, & Bressan 1996, who investigated the agreement in the predictions of independently-derived infrared models. Despite these caveats, the qualitative benefit of IR models is clear; timescales for synchronization of colours following a burst are much shorter in the IR and, unlike the case with optical colours, there is a very significant difference in the size of the color excursions arising from secondary star formation depending upon whether the bursts are intense and short-lived or of a gradual “drip-feed” nature.

With only half the HDF galaxies available, small number statistics become a significant complication, and consequently we remain cautious regarding the interpretation of the IR sample. However the general impression from a comparison of Figures 4 and 7 is that the optical/IR trends are very similar to intermediate redshifts ($`z`$0.6) but become somewhat different at high redshift ($`z0.6`$). At intermediate redshifts, in both Figures 4 and 7 a dominant red E/S0 locus can be identified, although a minority of the E/S0 population lies well blueward of the overall relation. Clearly the photometric error is larger in the IR data. Furthermore, for $`z`$0.6, most of the spiral bulges have optical and IR colors intermediate between these extreme cases. However, at high redshifts, a significant fraction of infrared bulges now lie alongside the red ellipticals, and these are largely those in high central concentration spirals (large open symbols). Indeed many $`z0.6`$ bulges appear to be slightly redder than their elliptical counterparts.

Assuming (as justified earlier) that disk contamination does not significantly complicate the interpretation of Figure 4 ($`\mathrm{\S }`$3.2) to redshifts $`z`$0.5 and that this is completely negligible at all redshifts in Figure 7, then the combined data suggest that bulge building activity at higher redshifts is likely to be predominantly operating in a “short-burst” mode, rather than as part of a more extended “drip feed” process. Comparison with the model tracks suggests that the more extended recovery timescale for optical colors following a low-mass, short burst could well result in the majority of $`z0.6`$ bulges lying near the red locus in the IR while at while being significantly offset in optical color.

A more interesting interpretation, though tentative given the small samples involved, would arise if the transition in the bulge IR color distribution seen in Figure 7 at or around $`z`$0.6 is significant. This might imply that optically blue bulges are more common at $`z0.6`$ than at $`z0.6`$ which could be understood if this redshift corresponded to one where a fundamental change in the bulge-building history of spirals occurred. This epoch has been claimed to be that at which barred spiral systems become fairly common in deep imaging data (Abraham et al. 1999). An increase in the fraction of blue bulges at this redshift could provide further evidence that the formation of Hubble’s “tuning fork” occurs at or around this epoch, assuming the excess of $`z<0.6`$ blue bulges corresponds to the onset of secular bulge-building processes. This would strengthen the link between bar formation and bulge-building (either through bar-driven gas flows or bar dissolution). It is perhaps interesting that a greater proportion of blue bulges at $`z0.6`$ occur in low central concentration galaxies than is the case at $`z0.6`$, since the importance of secular activity is probably a strong function of Hubble type, as described earlier.

### 5.5 Future Work

The statistical relationship between bulge color and bulge-to-disk ratio as a function of redshift may also yield information on the epoch of initial bulge formation. If the structural parameters of the bulge can be reliably disentangled from the underlying disk, then high-redshift data that probes in situ formation from a structural standpoint can be used to isolate the epoch of bulge formation independently from the age of its stars. Such a study would also allow us to extend the present work in a number of interesting ways. In the present paper we have effectively decoupled measures of galactic structure from measures of stellar population, and hence we are unable to address quantitatively the important question of whether the enhanced blue colors of spiral bulges are associated with position on the Hubble sequence. Since the best evidence from local data is that secular activity is most closely associated with late-type spirals, studies linking the relative colors of bulges and ellipticals with measures of growth in bulge-to-disk ratio are required. Secular growth can also be examined directly in high-redshift galaxy data by kinematical observations that will soon become feasible with intermediate resolution spectrographs on 10m-class telescopes. Such dynamical data can be used to determine whether the bulk motions of high-redshift bulges have a closer kinematical connection to disk rotation or to thermal motions.

Coupling the present analysis to structural bulge-to-disk decompositions will also allow investigations of the extent to which the relatively blue bulges seen in our data may be the manifestation of a colour-magnitude relation (CMR) in field ellipticals, similar to that seen in well-studied cluster samples (Sandage & Visvanathan 1978). No CMR relation was found in high redshift HDF field ellipticals by Kodama, Bower, & Bell (1999), although its presence may have been masked by uncertainties arising from poorly-understood photometric redshifts.

Assuming that the slope of the field CMR at high redshifts resembles that seen in local cluster data, then the CMR has a completely negligible impact on our infrared bulge colors, and on the optical colors at $`z<0.5`$. For example, the slope of the $`BR`$ CMR in Coma is $`<0.03`$ mag<sup>-1</sup> (Milvang-Jensen & Jorgensen 2000), and presumably the $`VR`$ slope is even smaller. At higher redshifts, where the $`I_{814}`$ filter enters the rest-frame $`V`$-band, the CMR may begin to make optical bulge colors measurably bluer, although the effect probably represents a relatively small contribution to the systematic error budget. For example, in order for a Coma-like CMR to explain the whole of the present trend for blue bulges would require the bulges in our sample to be 6–8 mag fainter than typical ellipticals (assuming a Coma-like CMR of $`\mathrm{\Delta }(BV)`$=0.05 mag<sup>-1</sup>), which is clearly incompatible with our visual pre-selection isolating relatively early-type spirals with large bulge-to-disk ratios. Furthermore a CMR is already incorporated into the hierarchical predictions given in Figure 6. Nevertheless, it will be interesting to compare the absolute magnitudes of the bulges in the present sample with the corresponding absolute magnitudes for ellipticals at similar redshifts, in order to establish the extent to which cluster CMR data can be extended to spiral bulges in the field.

It is notable that the trends shown in Figures 4 and 7 are also seen when 5% apertures (instead of total apertures) are used to define the colors of the E/S0 population. This is not surprising, since Appendix 1 of Menanteau et al. (1999) indicates that the mean metallicity of an elliptical interior to its effective radius is reasonably close to the metallicity at the effective radius. This calculation is also illuminating when assessing the possible importance of metallicity gradients in biasing of the colors of spiral bulges. If the metallicity gradients in bulges are similar in strength to those in the luminous ellipticals studied by Arimoto et al. (1997), then the very slow drop-off in the mean metallicity of light interior to an aperture larger than the bulge effective radius implies that our 5% apertures should constitute a fair sample of the total bulge light.

Finally, we seek reassurance that the models outlined in Figures 4 and 7 are consistent with the colors and scatter observed locally for bulges and ellipticals. Surprisingly, however, there is no precise measure of the local homogeneity of bulges (Wyse et al 1997). Peletier et al (1999), in a particularly detailed optical-infrared analysis of 20 nearby systems with HST, draw attention to the redder colors of the inner 100-200 pc some component of which they attribute to dust obscuration. (We note that the presence of nuclear dust at this level, if confirmed, would further strengthen the results of the present paper.) On scales corresponding to bulge effective radii, where dust effects are expected to be negligible, Peletier et al (1999) find mean bulge colors comparable to those of Coma ellipticals. However, imprecise photometric conversions were necessary to make this comparison, rendering the precision of this statement no better than $`\pm `$0.2 mag in $`BI`$. The color dispersion in the bulges themselves is typically also $``$0.15-0.2 mag in $`BI`$. Clearly, improved local data are needed to establish the relative colors of bulges and ellipticals in the local Universe.

## 6 Conclusions

A robust prediction of hierarchical scenarios that associates morphological transformations with mergers and disk growth is that bulges must be statistically redder and older than ellipticals. In order to test this prediction, we have collated a magnitude limited sample of spirals and ellipticals from both Hubble Deep Fields and compared the aperture-based $`V_{606}I_{814}`$ color-redshift relation for the central bulges of spirals selected in various ways with those of field ellipticals.

In detail we find the following:

* At $`z0.6`$, where we demonstrate that contamination from underlying disk light is minimal, bulges are significantly bluer than elliptical galaxies. Very few bulges follow traditionally-accepted passive evolutionary tracks. At optical wavelengths most bulges are 0.3-1.0 mag bluer than ellipticals at comparable redshifts, with a significant color dispersion at a given redshift, in marked contrast to the expectations of recent semi-analytic predictions. Statistically speaking, bulges remain optically bluer than ellipticals at all redshifts $`z0.6`$ to the limits of our data, and contamination by bluer disk light and/or aperture mismatches is unlikely to explain the bulk of this effect.
* Although at first sight our data suggest that bulges are recently formed, we show that secondary star formation involving bursts of 15% by mass superimposed upon a pre-existing old stellar population can readily reproduce the color dispersion observed. We discuss two general scenarios: short-lived bursts and extended events possibly associated with continuous infall of gas.
* Near-infrared NICMOS observations of a sub-sample of our dataset agree with the optical trends at $`z0.6`$, but at higher redshifts most infrared bulges appear to be as red as most ellipticals. Model comparison with the optical trends suggests that the duty cycle for the optical blueing of high-redshift bulges is quite short. Most importantly, however, our data firmly rule out the traditional picture whereby bulges form at high redshift and evolve passively as miniature ellipticals. Although somewhat speculative at this stage, the divergence observed in the IR properties of bulges and ellipticals at $`z0.6`$ may suggest that other physical processes become important in forming galactic bulges, for example secular evolution in stellar disks linked to the formation of galactic bars.
* Much work needs to be done at both high and low redshifts in order to establish the star-formation history of bulges in various galaxy types and their environments. It is clear that such studies are central to testing hierarchical models for the morphological evolution of galaxies. In this paper we have shown significant differences between the distribution of colors in bulges and ellipticals which either seriously challenge our theoretical views on bulge formation or, more likely, reveal the necessity for including additional physics into these models, corresponding to periods of extended star formation in the inner regions of spiral galaxies.

## 7 Acknowledgements

We thank Carlton Baugh, Shaun Cole and Carlos Frenk for their enthusiastic support of this project and Carlton Baugh for generously giving us semi-analytic predictions in the context of our observations. We thank Stephen Gwyn for kindly making his photometric catalog available to us, and for matching the object identifications in his catalog to those given in the STScI HDF catalog. Reynier Peletier is thanked for making the results from his NICMOS bulge progam known to us in advance of publication. We acknowledge useful discussions with Jarle Brinchmann, Ray Carlberg, Roger Davies, Gerry Gilmore, Piero Madau, Ron Marzke, Felipe Menanteau, Mike Merrifield, Rachel Somerville and Chuck Steidel.


text:
# Percolation and Critical Behaviour in SU(2) Gauge Theory The work has been supported by the TMR network ERBFMRX-CT-970122 and the DFG under grant Ka 1198/4-1.

## 1 INTRODUCTION

The search for a geometrical description of phase transitions has stimulated a lot of interest over the last decades. Such descriptions are simple and elegant and allow to justify the relationship between critical phenomena and geometry which seems to play a key role in important aspects of critical behaviour, like the universality of the critical indices.

In particular, percolation theory  represents an ideal framework for the required formulation. Percolation takes place when geometrical clusters, formed by elementary objects of some system, give rise to an infinite network (percolation cluster) which spans the whole system. The percolation phenomenon turns out to have astonishing analogies with ordinary second order thermal phase transitions: $``$ the behaviour of the percolation variables at criticality is also described by simple power laws, with corresponding exponents; $``$ the values of the exponents, related to each other by simple scaling relations, are fixed only by the number of space dimensions of the system at study, regardless of its structure and of the special type of percolation process one considers. The spontaneous symmetry breaking of the Ising model can be indeed described by means of percolation of site-bond clusters . Since the critical behaviour of $`SU(2)`$ pure gauge theory is in the same universality class as the Ising model , an equivalent percolation picture for the confinement-deconfinement transition of $`SU(2)`$ has been recently proposed . Unfortunately, the approach suggested in  is valid only in the strong coupling limit and the problem of finding a general description remains open.

In this work, we will show that it is possible to define the clusters we need by constructing suitable effective theories for $`SU(2)`$ which admit an equivalent percolation picture. The effective models are Ising-like spin models with only spin-spin interactions, and we will show that they provide a working percolation description of the critical behaviour of $`(3+1)`$-$`d`$ $`SU(2)`$ both in the strong coupling limit ($`N_\tau =2`$) and for $`N_\tau =4`$, a case which approaches the weak coupling region.

## 2 PERCOLATION AND CRITICAL BEHAVIOUR IN THE ISING MODEL

The first percolation studies on the Ising model concerned the ordinary magnetic domains, i.e. clusters formed by nearest-neighbour spins of the same sign. In two dimensions such clusters happen indeed to percolate at the thermal critical point $`T_c`$ . Nevertheless, the values of the critical exponents differ from the corresponding Ising values . In three dimensions, the magnetic domains of the spins oriented in the direction of the magnetization percolate at any temperature; the domains formed by the spins opposite to the magnetization percolate for $`TT_p`$, with $`T_pT_c`$ .

It was then clear that the magnetic domains were too big to reproduce the thermal behaviour of the Ising model. The fact that two nearest neighbour spins of the same sign are always connected to each other implies that the clusters fail to reproduce the correct spin correlation because of geometrical effects. For this reason, Coniglio and Klein  suggested that two like-signed nearest neighbouring spins belong to the same cluster with a certain bond probability

$`p=1\mathrm{exp}(2J/kT)`$ (1)
($`J`$ is the Ising coupling, $`T`$ the temperature). The clusters introduced by Coniglio and Klein are thus site-bond clusters: their percolation transition is indeed equivalent to the magnetization transition of the Ising model.

As a matter of fact, this result is valid for a wide class of models. For example, if a theory is characterized just by ferromagnetic spin-spin interactions, the percolation picture of Coniglio and Klein can be trivially extended by introducing a bond between each pair of interacting spins, and a relative bond probability

$`p_i=1\mathrm{exp}(2J_i/kT),`$ (2)
where $`J_i`$ is the coupling associated to the $`i`$-interaction of the theory.

This generalization is at the basis of our work.

## 3 POLYAKOV LOOP PERCOLATION IN SU(2) GAUGE THEORY

Let us now consider the case of $`SU(2)`$ gauge theory. We have again a $`Z(2)`$ symmetric variable, the Polyakov loop $`L`$, whose thermal lattice average is the order parameter of the confinement-deconfinement transition. Therefore, the $`SU(2)`$ configurations are similar to the Ising ones, with the Polyakov loop playing the role of the spin variable. Also here we will have regions in which $`L`$ has the same sign, analogous to the magnetic domains of the Ising model. However, to identify the physical clusters, as we have seen above, we need to introduce some bond probability. This is a difficult problem, since the $`SU(2)`$ lattice action cannot be expressed in terms of the Polyakov loop $`L`$. In order to derive the correct bond weights, in fact, it is necessary to know how the ”gauge spins”, i.e. the Polyakov loops, interact with each other.

One way to face the problem is to approximate $`SU(2)`$ by means of effective theories which are easy to handle. In  one derived an effective theory from a strong coupling expansion of $`SU(2)`$. That method thus cannot be applied to the weak coupling case.

To find a general procedure, we shall construct the effective theory starting not from the $`SU(2)`$ Lagrangian, but from the Polyakov loop configurations. If it were possible to find an effective model which reproduces well the Polyakov loop configurations and has an equivalent percolation formulation, the problem would be solved. The Polyakov loop, unlike the Ising spins, is a continuous variable, but we shall consider the Ising-projected configurations, i.e. the distributions of the signs of the Polyakov loops. This is done assuming that the $`Z(2)`$ symmetry is the only relevant feature at the basis of the critical behaviour. Our ansatz for the effective theory will be an Ising-like model with just spin-spin interactions. Its Hamiltonian $`(s)`$ is

$`(s)=J_1{\displaystyle \underset{NN}{}}s_is_jJ_2{\displaystyle \underset{NNN}{}}s_ks_letc.,`$ (3)
where the distance between coupled spins increases progressively starting from the simple nearest-neighbour ($`NN`$) case ($`NNN=`$next-to-nearest, and so on). We have seen in the previous section that, if $`J_i>0i`$, such a model has a simple percolation picture a la Coniglio-Klein.

The relationship between $`SU(2)`$ and the effective model is established through the equation

$`e^{(s)/kT}={\displaystyle [dU]\underset{𝐧}{}\delta [s_𝐧,sgn(L_𝐧)]e^{S_{SU2}}},`$ (4)
where $`L_𝐧`$ is the value of the Polyakov loop at the spatial point n and $`S_{SU2}`$ the $`SU(2)`$ lattice action.

The couplings of the effective theory are calculated by solving a set of Schwinger-Dyson equations, a method which gives good results in Monte Carlo renormalization group studies of field theories . If the couplings turn out to be all positive, we can use them to calculate the corresponding bond weights we need to build the clusters in the original Polyakov loop configurations. Because of the expression (2) for the bond weights, with this approach only the configurations of the signs of the Polyakov loops are relevant for the cluster building.

## 4 RESULTS FOR $`(3+1)`$-$`d`$ SU(2)

We performed simulations of $`(3+1)`$-$`d`$ $`SU(2)`$, focusing on the cases $`N_\tau =2`$ and $`N_\tau =4`$ ($`N_\tau `$ is the number of lattice spacings in the temperature direction). We found that the two cases can be indeed approximated by models like (3), with 15 and 19 ferromagnetic couplings, respectively. The distances between the spins in the longest-range interactions, in lattice units, are $`\sqrt{17}`$ for $`N_\tau =2`$ and $`\sqrt{27}`$ for $`N_\tau =4`$. The relative strength of such interactions compared to the nearest-neighbour couplings is about $`1/1000`$ in both cases, which shows how rapidly the couplings decrease with the distance between the spins.

Both for $`N_\tau =2`$ and for $`N_\tau =4`$ we used four different lattice sizes, in order to be able to extract the results by applying finite size scaling techniques.

To locate the critical point of the percolation transition, we made use of the percolation cumulant. Such variable is determined simply by counting, for a fixed lattice size and a value of the $`SU(2)`$ coupling $`\beta =4/g^2`$, how many configurations contain a spanning cluster. This number, divided by the total number of configurations, is directly a scaling function, analogous to the Binder cumulant in continuous thermal phase transitions .

Figs. 1 and 2 show the behaviour of the percolation cumulant curves as a function of $`\beta `$, for $`N_\tau =2`$ and $`4`$, respectively. The curves cross remarkably at the same point, within errors, in good agreement with the thermal thresholds, indicated within one standard deviation by the dashed lines.

The percolation cumulant gives the critical point $`\beta _c`$ but also the critical exponent $`\nu `$. In fact, if we rewrite the percolation cumulant as a function of $`(\beta \beta _c)L^{1/\nu }`$ ($`L`$ is the linear spatial dimension of the lattice), we should get the same curve for each lattice size. Fig. 3 and 4 show the rescaled percolation cumulant using for $`\beta _c`$ the position of the crossing points and for the exponent $`\nu `$ the 3D Ising value. In both cases, the curves coincide.

To complete the investigation of the geometrical transition, we determined the ratio of critical exponents $`\gamma /\nu `$ from a standard finite size scaling analysis. We were not able to evaluate the ratio $`\beta /\nu `$, because of large fluctuations of the percolation order parameter, the percolation strength $`P`$, around the critical coupling. The final results are listed in Tables 1 and 2, where we give, for comparison, also the critical indices of the thermal transition. The agreement is good in both cases.

## 5 CONCLUSIONS

We have shown that the deconfinement transition of SU(2) gauge theory in 3+1 dimensions can be described in terms of Polyakov Loop percolation. The procedure we used is, necessarily, approximate. However, in comparison with the method illustrated in , it has the advantage to be applicable also to the weak coupling limit of $`SU(2)`$.