text:
# A universal angular momentum profile for galactic halos

## 1 Introduction

00footnotetext: Hubble Fellow

The origin of the distribution of mass and angular momentum in disk galaxies is still an open issue, despite its long history.

The archetypical model of Eggen, Lynden-Bell & Sandage (1962) provides a useful framework, in which a disk galaxy (the Milky Way) forms by centrifugal support of a collapsing gas cloud. Mestel (1963) added the assumption that the specific angular momentum of each mass element, $`j`$, is conserved during the collapse, and demonstrated that the final mass profile of the disk could then be related to the initial mass and angular momentum distribution of the gas. The key is that the cumulative mass with specific angular momentum less than $`j`$, $`M(<j)`$, is preserved during the collapse. Crampin & Hoyle (1964) (also Innanen 1966; Freedman 1970) then realized that the observed exponential density profiles of disk galaxies in circular motions are consistent with the angular momentum distribution of a hypothetical uniform sphere in solid body rotation,

$$M(<j)=M_{\mathrm{tot}}\left[1(1j/j_{\mathrm{max}})^{3/2}\right]$$
(1)
where $`j_{\mathrm{max}}=\omega R^2`$ with $`\omega `$ the fixed angular velocity and $`R`$ the radius of the sphere, and $`M_{\mathrm{tot}}`$ the total mass of the sphere.

Fall & Efstathiou (1980) re-examined this question in the context of the more modern view according to which disk galaxies form from contracting gas within extended dark-matter halos (White & Rees 1978). They added the assumption that the gas and dark matter were initially well-mixed such that the distribution of $`j`$ in the disk is equal to that of the halo. Rather than assuming an angular momentum profile, they implicitly deduced it from the constraint that the final surface density profile of the disk is exponential. The implied disk scale radius, $`r_\mathrm{d}`$, is related to the total pre-collapse angular momentum of the halo, roughly $`r_d\lambda `$. Here, $`\lambda `$ is the dimensionless spin parameter (e.g. Peebles 1969)

$$\lambda \frac{J\left|E\right|^{1/2}}{GM^{5/2}},$$
(2)
where $`J`$,$`E`$ and $`M`$ are the total angular momentum, energy and mass of the system, and $`G`$ is Newton’s constant.

Angular momentum is presumably acquired by the dark matter (and gas) through tidal interactions with neighboring objects (Peebles 1969). Barnes & Efstathiou (1987) (following Efstathiou & Jones 1979; Efstathiou & Barnes 1983; Zel’dovich & Novikov 1983) utilized cosmological N-body simulations with a CDM power spectrum to make detailed predictions for halo angular momentum structure and statistics. In particular they found that the distribution of spin parameters is roughly log normal, with a median value of $`\lambda 0.05`$. They also confirmed that the angular momentum of pre-collapse halos grows roughly linearly with time, as predicted by linear tidal-torque theory (Doroshkevich 1970; White 1984). Several numerical and analtyical investigations have followed, which confirmed the nature of the $`\lambda `$ distribution, addressed the alignment of the spin vectors between neighboring galaxies, the (lack of) dependence of halo spin on environment, and the effect of major mergers on halo spin (e.g., Frenk et al. 1988; Heavens & Peacock 1988; Zurek, Quinn & Salmon 1988; Warren et al. 1992; Catelan & Theuns 1996; Cole & Lacey 1996; Lemson & Kauffmann 1997; Gardner 2000).

Armed with this global information about the initial state of the dark matter and gas mixture, several authors have used the above disk formation scenario to predict observable quantities. Blumenthal et al. (1986) addressed the rotation curves of disk galaxies in the context of the CDM hierarchical formation scenario, assuming $`10\%`$ of the mass in the gas component (Blumenthal et al. 1984). They pointed out that the non-dissipative dark halo should react to the dissipative gas infall by considerable contraction, and showed that this process can be approximated using an adiabatic invariant. This coupling between the halo and the disk “conspires” to produce a continuity between the disk-dominated and halo-dominated regions of the rotation curves, as observed (Burstein & Rubin 1985).

Flores et al. (1993) explored the dependence of the final rotation curves on a range of assumed parameters, including the measured spread in $`\lambda `$. Dalcanton, Spergel, & Summers (1997) included a more realistic initial halo density profile (Hernquist 1990) and explored the properties of both high surface brightness and low surface brightness galaxies. Mo, Mao, & White (1998a; 1998b, 1999) performed a similar examination within the context of various cosmologies, assuming that the dark halos followed the density profile advocated by Navarro, Frenk, & White (1996, 1997; NFW) and compared their predicted galaxy properties to both local and high-redshift data. van den Bosch (1998; 2000) and van den Bosch & Dalcanton (2000) extended the approach to include a bulge component and investigated the origin of the Hubble sequence of galaxy types and other aspects of galaxy formation.

The main commonality in all of these investigations is that the variation in observable galaxy properties depends almost exclusively on the initial total spin parameter of the halo hosting the galaxy. However, all of the predictions rest on the implicit assumption that, if specific angular momentum is conserved, the initial angular momentum profiles of the halos produce disks that are exponential in form. 111An exception is Dalcanton et al. (1997), who assumed that the initial $`M(<j)`$ profile was that of a uniform solid body rotator, which produces disk density profiles that are similar to exponential (Ryden & Gunn 1987) but only over about two disk scale lengths (Olivier, Primack, & Blumenthal, 1991).

Within this simple framework for disk formation, the key missing ingredient is the actual distribution of angular momentum in halos. Partially motivated by this, we set out to measure the $`j`$ profiles directly using a statistical sample of halos in a high resolution simulation, and to estimate the intrinsic scatter in these distributions. The long-term goal is to learn what the detailed form of halo angular momentum profiles can teach us about the origin of galaxy types, disk surface brightness profiles, and bulge properties.

An additional strong motivation for this study comes from the results of hydrodynamic simulations of disk galaxy formation in a cosmological context. Unlike the simplified, monolithic-collapse framework discussed above, the hydrodynamic treatments model gas cooling within the hierarchical growth of dark halos, and typically find that the resulting disks are significantly smaller than real galactic disks (Navarro & Steinmetz 2000; Weil, Eke, & Efstathiou 1998; Navarro & Steinmetz 1997; Navarro, Frenk, & White 1995). The problem apparently arises because most of the mass is accreted through mergers with objects whose gas component has already cooled and clumped, causing a large fraction of their angular momentum to be transferred to the dark halo or transported outwards. Since angular-momentum transport is likely to always be from the inside out, the simple case in which $`j`$ is conserved during the collapse, which we can study in detail with high accuracy, can serve as a useful limiting case. This will provide a better basis for addressing the more complex physical processes required for solving the angular-momentum problem seen in the hydro simulations.

In §2 we describe our methods, including the simulation and halo finding, the measurement of the spin parameter, the associated errors, and the measurement of the angular-momentum profile. In §3 we discuss the measured (universal) mass distribution of angular momentum in halos. In §4 we address the spatial distribution of angular momentum, its alignment, cylindrical symmetry, and profile in spherical shells. In §5 we use both tidal-torque theory and an estimate of angular momentum transfer due to minor mergers to explore the origin of the angular-momentum profile. In §6 we briefly describe the implications for the disk surface density profile and the “angular momentum catastrophe”. We discuss our results and conclude in §7. An Appendix is devoted to testing our results with a simulation having $`8`$ times the mass resolution, and investigating $`M(<j)`$ profiles of halos at high redshift.

## 2 Method

### 2.1 Halos in the ART simulation

Large cosmological N-body simulations have reached the stage where detailed structural properties of many dark-matter halos can be resolved simultaneously. A method which allows the required force and mass resolution is the Adaptive Refinement Tree (ART) method (Kravtsov, Klypin, & Khokhlov 1997) which implements successive refinements of the spatial grid and time steps in high density environments. The simulations based on the ART code provide a statistical sample of halos, covering a wide mass range, with the resolution appropriate for studying angular-momentum profiles.

We have used the ART code to simulate the evolution of collisionless dark matter within the “standard” low-density flat $`\mathrm{\Lambda }`$CDM model ($`\mathrm{\Omega }_\mathrm{m}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`h=0.7`$, and $`\sigma _8=1.0`$ at $`z=0`$). The simulation followed the trajectories of $`256^3`$ particles within a periodic box of comoving size $`60h^1\mathrm{Mpc}`$ from redshift $`z=40`$ to the present. A basic $`512^3`$ uniform grid was used, and six refinement levels were introduced in the regions of highest density, implying a dynamic range of $`32,000`$. The formal resolution of the simulation is thus $`f_{\mathrm{res}}=1.8h^1\mathrm{kpc}`$, and the mass per particle is $`m_p=1.1\times 10^9h^1M_{}`$.

Our halo finding algorithm (Bullock et al. 2000) is based on the Bound Density Maxima technique (Klypin & Holtzman 1997) and has been specifically designed to identify halos (and subhalos) in such a high resolution simulation. In the current study of angular momentum, we limit ourselves to ‘distinct’ halos, which do not reside within a larger host halo, in the mass range $`10^{12}10^{14}h^1M_{}`$. In the Appendix, we use a higher resolution simulation to extend our mass range down to $`10^{11}h^1M_{}`$.

A brief description of the halo finder is as follows: After finding all the maxima in the smoothed density field of the simulation, we unify overlapping maxima, define a minimum number of particles per halo (say 50), and iteratively find the center of mass of a sphere about each of the remaining maxima. We compute the spherical density profile about each center and identify the halo virial radius $`R_\mathrm{v}`$ inside which the mean overdensity has dropped to a value $`\mathrm{\Delta }_\mathrm{v}`$, based on the top-hat spherical infall model. For the family of flat cosmologies ($`\mathrm{\Omega }_\mathrm{m}+\mathrm{\Omega }_\mathrm{\Lambda }=1`$), the value of $`\mathrm{\Delta }_\mathrm{v}`$ can be approximated by (Bryan & Norman 1998) $`\mathrm{\Delta }_\mathrm{v}(18\pi ^2+82x39x^2)/(1+x)`$, where $`x\mathrm{\Omega }_m(z)1`$. In the $`\mathrm{\Lambda }`$CDM model used in the current paper, $`\mathrm{\Delta }_\mathrm{v}`$ varies from about 180 at $`z1`$ to $`\mathrm{\Delta }_\mathrm{v}340`$ at $`z=0`$.

We then fit each halo density profile with a universal functional form. We adopt the NFW profile,

$$\rho _{\mathrm{NFW}}(r)=\frac{\rho _s}{(r/r_s)\left(1+r/r_s\right)^2},$$
(3)
with the two free parameters $`r_s`$ and $`\rho _s`$. An equivalent pair of parameters is, for example, the virial mass $`M_\mathrm{v}`$ and the concentration parameter $`c_\mathrm{v}R_\mathrm{v}/r_s`$. Using this fit, we iteratively remove unbound particles from each modeled halo, and unify every two halos that overlap in their $`r_s`$ and are gravitationally bound.222Removal of unbound particles is a new feature allowed by the high resolution of our simulation; it could not have been properly implemented in earlier studies. The modeling of the halos with a given functional form allows us to assign to them characteristics such as a virial mass and radius, and to estimate sensible errors for these quantities.

The halo finding is 100% complete for halos of more than $`150`$ particles, $`M1.5\times 10^{11}h^1M_{}`$ (see Sigad et al. 2000). For our purpose here, we limit the sample to halos of more than $`1000`$ particles.

### 2.2 Angular-momentum errors

A first step towards measuring the distribution of specific angular momentum within each of the halos is measuring the global spin parameter of each halo. The angular momentum of a halo of $`N`$ particles is defined by

$$\text{J}=m_i\underset{i=1}{\overset{N}{}}𝐫_𝐢\times 𝐯_𝐢,$$
(4)
where $`𝐫_𝐢`$ and $`𝐯_𝐢`$ are the position and velocity of the $`i`$th particle with respect to the halo center of mass. In principle, given the spin parameter $`\lambda `$, the value of the global specific angular momentum, $`J/M`$, can be determined by using an assumed energy content for the halo in Eq. (2). In practice, however, this is not a straightforward procedure. For example, the energy of a halo in a crowded region is somewhat ambiguous because it depends on the environment. A related difficulty arises when only a sub-volume of the virial sphere is concerned, e.g., within the cooling radius. We therefore define an alternative and more practical spin parameter by

$$\lambda ^{}\frac{J}{\sqrt{2}MVR},$$
(5)
given the angular momentum $`J`$ inside a sphere of radius $`R`$ containing mass $`M`$, and where $`V`$ is the halo circular velocity at radius $`R`$, $`V^2=GM/R`$. This spin parameter reduces to the standard $`\lambda `$ when measured at the virial radius of a truncated singular isothermal halo. The value of $`\lambda ^{}`$ turns out to be robust to the choice of an outer radius; we use below the virial radius $`R_\mathrm{v}`$, but we find a very similar distribution of $`\lambda ^{}`$ values when using $`R_\mathrm{v}/2`$ instead. This allows modelers to freely scale our results using any desired outer radius appropriate for the problem at hand.

The spin parameters $`\lambda ^{}`$ and $`\lambda `$ are in fact very similar for typical NFW halos. They are related by $`\lambda ^{}=\lambda |E_{SIS}/E_{NFW}|^{1/2}\lambda f(c_\mathrm{v})^{1/2}`$, where $`E_{SIS}`$ and $`E_{NFW}`$ are the energies of virialized halos with isothermal and NFW density profiles respectively, and $`f(c_\mathrm{v})[2/3+(c_\mathrm{v}/21.5)^{0.7}]`$ (Mo et al. 1998) 333The exact expression is $`f(c_\mathrm{v})=0.5c_\mathrm{v}[(1+c_\mathrm{v})^212(1+c_\mathrm{v})\mathrm{ln}(1+c_\mathrm{v})]/[c_\mathrm{v}(1+c_\mathrm{v})\mathrm{ln}(1+c_\mathrm{v})]^2`$.. The function $`f(c_\mathrm{v})`$ is about unity for typical concentrations, $`c_\mathrm{v}10`$.

The distribution of $`\lambda ^{}`$ over the $`500`$ halos in our sample is shown in Fig. 2.1. It is well fit by a log-normal distribution,

$$P(\lambda ^{})=\frac{1}{\lambda ^{}\sqrt{2\pi }\sigma }\mathrm{exp}\left(\frac{\mathrm{ln}^2(\lambda ^{}/\lambda _0^{})}{2\sigma ^2}\right),$$
(6)
with best fit values $`\lambda _0^{}=0.035\pm 0.005`$ and $`\sigma =0.5\pm 0.3`$. Not surprisingly, the distribution of $`\lambda ^{}`$ is very similar to the known distribution of $`\lambda `$ (Barnes & Efstathiou 1987). The distribution of $`\lambda `$ values for our halos has best fit values $`\lambda _0=0.042\pm 0.006`$ and $`\sigma =0.5\pm 0.35`$.

To highlight a small difference between $`\lambda ^{}`$ and $`\lambda `$, Fig. 2.2 displays the values of $`\lambda ^{}`$ versus $`c_\mathrm{v}`$ for our sample of simulated halos, and shows a weak correlation. Given that $`\lambda `$ is known to be uncorrelated with $`c_\mathrm{v}`$ (NFW), this correlation arises from $`f(c_\mathrm{v})`$ in the above relation between the two spin parameters (as demonstrated in the figure) 444In agreement with the NFW result, we also find no correlation between $`c_\mathrm{v}`$ and the traditional spin parameter $`\lambda `$..

The error on the measurement of $`𝐉`$ in the simulations is of primary importance for the interpretation of the results and, here, for the specific task of fitting a functional form to the angular-momentum profile. Obtaining a precise measurement of $`𝐉`$ in halos is difficult because the signal-to-noise is relatively small due to the fact that the local coherent velocity of the halo contributing to the angular momentum is small compared the velocity dispersion. The variable $`\lambda ^{}`$ itself roughly characterizes the signal to noise (per particle) at any radius $`R`$, Eq. (5), which is typically at the level of a few percent.

In order to obtain a crude estimate of the error, which is probably an upper limit, we singled out several large halos with $`N_i`$ simulated particles ($`>2\times 10^4`$), and measured their angular momentum using a series of sparsely sampled realizations of $`N`$ particles down to $`1\%`$ of $`N_i`$, with 100 random realizations for each value of $`N`$. The average and standard deviation of the measured value $`J`$ relative to the original value obtained with all the particles, $`J_i`$, as a function of $`N`$, are shown in Fig. 2.2 for three halos with different values of $`\lambda ^{}`$. We find that for a subsample of $`N`$ particles the standard deviation $`\sigma _J`$ scales like

$$\frac{\sigma _J}{J}=\sqrt{1/N+1/(25\lambda ^2N)}\frac{0.2}{\lambda ^{}\sqrt{N}}.$$
(7)
This scaling relation is shown as error bars in Fig. 2.2, and it fits rather well the measured scatter. The first term in Eq. (7) corresponds to the general Poisson sampling error and the second term reflects the noise due to velocity dispersion. The second term typically dominates because $`\lambda ^{}`$ is small, and thus the spin parameter itself directly influences our ability to measure $`J`$ accurately, with the number of particles needed for a fixed accuracy scaling like $`\lambda ^2`$. For $`\lambda ^{}=0.03`$, even $`50\%`$ accuracy requires $`N150`$.

The above estimate of the error probably tends to overestimate the true error because it includes small-scale noise that is not present when a smaller number of particles is used in the actual dynamical simulation of the halo. On the other hand, our procedure above does not mimic the error made in following the true dynamics of the smooth halo when a smaller number of particles is used — an effect which tends to cause an underestimate of the error. An accurate evaluation of the error in $`J`$ requires a direct comparison of the same halo as simulated with different resolutions, and is clearly sensitive to the $`N`$-body code used. In our analysis here we limit ourself to the crude estimate, Eq. (7). In the Appendix we use a simulation with $`8`$ times the mass resolution in order to check for systematic differences in the derived $`M(<j)`$ profiles, and find none. This analysis gives us confidence that mass resolution does not affect our results in any severe way.

### 2.3 Measuring the angular-momentum profile

Our goal is to determine how the angular momentum is distributed in the halo. We wish to compute the mass distribution of specific angular momentum, $`M(<j)`$, and also the spatial distribution of $`j`$.

We first compute the total $`𝐉`$ for each halo, and let it define the $`z`$ direction. Then we subdivide the spherical volume of each halo into many spatial cells, each containing many particles, and compute the specific angular momentum in each cell, projected along the $`z`$ direction. In general we use the symbol $`j`$ to refer to this $`z`$-projected angular momentum. In some instances we also investigate the non-projected, absolute magnitude of angular momentum in each cell, which we refer to as $`|j|`$. When comparisons between the two measurements are made, we explicitly add the subscript $`z`$ to the projected value, $`j_z`$.

The cell geometry and spatial distribution were designed to maximize the number of particles in each cell while at the same time sampling as much of the $`j`$ distribution as possible. The vector character of angular momentum naturally introduces a preferred axis, convolved with the general spherical symmetry of the mass distribution. The cells are defined using the usual spherical coordinates ($`r`$,$`\theta `$,$`\varphi `$) about the halo center. Each of our cells spans the full $`2\pi `$ range in $`\varphi `$, and they span the range of ($`r/R_\mathrm{v}`$, $`\mathrm{sin}\theta `$) from ($`0,0`$) to ($`1,1`$). The radial shells are spaced such that each contains the same number of particles. The shells are then divided into three azimuthal cells of equal volume between $`\mathrm{sin}\theta =0`$ and $`1`$. Positions with the same $`r`$ $`\mathrm{sin}\theta `$ above and below the equatorial plane belong to the same cell. For all the $`M(<j)`$ profiles studied in detail and presented in the figures below, we adjusted the number of radial bins such that each cell contains roughly $`N_p`$ particles, where $`N_p=500`$ if $`M_\mathrm{v}>5.5\times 10^{12}h^1M_{}`$ and $`N_p=0.1M_\mathrm{v}/m_p`$ for $`M_\mathrm{v}<5.5\times 10^{12}h^1M_{}`$.

We construct $`M(<j)`$ profiles for each halo by ranking the cells by their $`j`$ values, and counting the cumulative mass in cells with angular momentum less than $`j`$. Because $`j`$ corresponds to a projected component, it is possible for a cell to have a negative value of $`j`$. The anti-alignment needed for a cell to have a negative $`j`$ value is rare, but when one occurs, we remove the cell completely from the ranked list used to construct the $`M(<j)`$ profile. About $`5\%`$ of halos have a significant amount of their total mass ($`10\%`$) contained in negative $`j`$ cells. We do not consider these anti-aligned halos in our $`M(<j)`$ analysis.555We include all halos when we present the distribution of half-mass alignment cosines in §4.1

We assign an error to the $`j`$ in each cell using Eq. (7). Every halo analyzed has $`N>1000`$ particles within its virial radius, corresponding to $`M_\mathrm{v}10^{12}h^1M_{}`$. We have $`600`$ halos that meet our requirements. In practice, we determined the shape of the functional form which fits $`M(<j)`$ by examining the group of ($`200`$) halos with $`N>6000`$ particles, for which the profiles are determined with greater accuracy. We extended the mass range down to $`1000`$ particles in order to determine whether more typical galactic-sized halos have a similar profile. Although these smaller halos have larger errors on their fit parameters, they seem to obey the same $`M(<j)`$ profile.

The main limitation of this analysis is mass resolution. Because we need $`100`$ particles per cell in order to obtain a reasonable measure of $`j`$, our determination of $`M(<j)`$ is only reliable down to $`10\%`$ for $`M_\mathrm{v}<5.5\times 10^{12}h^1M_{}`$ and $`(5\times 10^4m_p/M_\mathrm{v})\%`$ for $`M_\mathrm{v}>5.5\times 10^{12}h^1M_{}`$.

## 3 Mass distribution of angular momentum

### 3.1 A universal angular-momentum profile

We find that the specific angular-momentum profiles of halos in the simulation are well-fit by the following two parameter function:

$$M(<j)=M_\mathrm{v}\frac{\mu j}{j_0+j},\mu >1.$$
(8)
The profile has an implicit maximum specific angular momentum $`j_{\mathrm{max}}=j_0/(\mu 1)`$. It is roughly a power-law for $`jj_0`$, and flattens out for $`jj_0`$. As can be seen by its relation to $`j_{\mathrm{max}}`$, the quantity $`\mu `$ ($`>1`$) acts as a shape parameter: for $`\mu 1`$, $`M(<j)`$ is a pure power law, while $`\mu 1`$ means that only half the mass falls within the power-law regime and the bend is pronounced.

The pair of parameters $`\mu `$ and $`j_0`$ fully defines the angular momentum distribution of the halo. The global spin parameter is related to $`\mu `$ and $`j_0`$ via

$$j_0b(\mu )=\sqrt{2}V_\mathrm{v}R_\mathrm{v}\lambda ^{},$$
(9)
where

$$b(\mu )_0^1\frac{m}{\mu m}𝑑m=\mu \mathrm{ln}(1\mu ^1)1.$$
(10)
Once the value of $`\lambda ^{}`$ is known, Eq. 8 is a one-parameter fitting function. In fitting the $`M(<j)`$ profiles for each halo, we first measured $`\lambda ^{}`$ and then obtained the best fit $`M(<j)`$ with the constraint that $`\mu `$ and $`j_0`$ return the measured $`\lambda ^{}`$. The pair $`(\lambda ^{},\mu )`$ is perhaps the most useful pair of parameters for characterizing the halo angular momentum because $`\lambda ^{}`$ is the conventional global measure and $`j_0`$ can be explicitly determined from $`\lambda ^{}`$ and $`\mu `$ via Eq. (9), while $`\mu `$ cannot be determined explicitly from the other two parameters without iteratively solving Eq. (10). Alternatively, we provide the following approximation for the iterative solution of Eq. (10):

$$\mu (b)1+\frac{0.5}{[\mathrm{exp}(1.09b)1]},$$
(11)
which is accurate to within $`8\%`$ in $`\mu 1`$ over the range $`\mu 1=0.0110`$.

Figure 1 shows $`M(<j)`$ examples for several of our halos. In the left panel we have normalized the profiles by $`M_\mathrm{v}`$ and $`j_{\mathrm{max}}`$ in order to illustrate how different values of $`\mu `$ affect the distribution. As $`\mu `$ approaches its minimum value of $`1`$, a larger fraction of the halo mass is spinning slowly compared to $`j_{\mathrm{max}}`$. Larger values of $`\mu `$ imply a more uniform $`j`$ distribution. The right panel shows the same four halos, now normalized to their best-fit $`j_0`$ and $`M_0`$ values, where $`M_0=M(<j_0)=\mu M_\mathrm{v}/2`$. Notice how remarkably the halos follow the characteristic functional shape.

The spread in profile shapes can be described by the distribution of $`\mu 1`$ ($`=j_0/j_{\mathrm{max}}`$). Figure 4.1 shows the distribution of $`\mathrm{log}_{10}(\mu 1)`$ for all of our halos (histogram) along with a Gaussian distribution with the same mean and standard deviation: $`\mathrm{log}_{10}(\mu 1)=0.6`$, $`\sigma =0.4`$. The implied $`90\%`$ range is $`\mu 10.061.0`$.

### 3.2 Correlations between parameters

Although the spin parameter $`\lambda ^{}`$ and the shape parameter $`\mu `$ clearly measure two different aspects of the angular-momentum distribution, they are not necessarily uncorrelated. Figure 4.1 shows the joint distribution of these two parameters for our simulated halos. We observe a weak but significant linear correlation between $`\mathrm{log}(\mu 1)`$ and $`\mathrm{log}\lambda ^{}`$, with a Pearson correlation coefficient of $`r0.23`$ (corresponding to a probability $`p10^6`$ for no correlation). Thus, high-spin systems tend to have more evenly distributed (power-law) $`M(<j)`$ profiles than low-spin systems, but the scatter about this trend is large.

Does the angular-momentum distribution correlate with the mass-density distribution? We know that the global spin characterized by $`\lambda ^{}`$ does not correlate strongly with the halo mass (Barnes & Efstathiou 1987) and correlates only weakly with the concentration parameter (Fig. 2.2). Figure 4.1 shows $`\mu 1`$ versus these mass parameters. We detect a marginal anti-correlation with mass, with $`r0.142`$ ($`p4\times 10^3`$ for no correlation), and an even less significant anti-correlation with $`c_\mathrm{v}`$, of $`r0.055`$ ($`p0.26`$ for no correlation).

## 4 Spatial distribution of angular momentum

In the previous section we explored the distribution by mass of the specific angular momentum within each halo. This does not tell us much about how well the angular momentum is aligned throughout the halo or how angular momentum is distributed spatially. This section is explicitly devoted to these issues. In §4.1 we address the question of alignment, in §4.2 we investigate the cylindrical symmetry of the angular-momentum distribution, and in §4.3 we explore the angular-momentum profile in spherical shells.

### 4.1 Alignment

Here we address the question of how well the smoothed angular momentum is aligned throughout the halo volume.

For our first alignment statistic, we divide the particles in each halo into an inner half-mass sphere and an outer half-mass shell and measure the cosine of the angle between inner and outer angular-momentum vectors, $`\mathrm{cos}\theta _{1/2}`$. The distribution of measured alignment values is shown in the top panel of Figure 4.1.

Before we attempt to interpret these results, it is essential to keep in mind the large errors involved in determining angular momentum direction vectors. Even for a halo with an intrinsically well-aligned angular momentum distribution ($`\mathrm{cos}\theta _{1/2}1`$), measurement errors in each of the six directional components involved in the alignment cosine will tend to drive the measured value towards zero, and, if large enough, will lead to a false measure of misalignment. The distribution in the top panel can therefore serve as a lower bound for the true alignment.

In order to obtain a conservative estimate of the distribution of true alignment cosines, we use Eq. (7) to assign errors to each component of $`𝐉`$ in both the inner and outer half-mass regions and perform a standard propagation of errors to obtain an estimate for the measurement error on $`\mathrm{cos}\theta _{1/2}`$. We add the estimated error in each case to the measured alignment value in order to obtain a “corrected” distribution for $`\mathrm{cos}\theta _{1/2}`$. This is shown in the lower panel of Figure 4.1. As argued above, we believe Equation 7 is an overestimate for the measurement error, so our corrected alignment cosines represent a plausible upper limit for the alignment. The true intrinsic distribution probably lies somewhere between the corrected and uncorrected distributions. If we make this assumption, then the two distributions allow us to conclude that between $`70\%`$ and $`90\%`$ of the halos are aligned to a greater degree than $`\mathrm{cos}\theta _{1/2}=0.7`$, although there is also a tail of significantly misaligned halos.

If the angular momentum within a halo is poorly aligned, one might expect a reflection of this to appear in its $`M(<j_z)`$ profile. There is a hint for such a trend in Fig. 4.1, where the distribution of $`\mathrm{cos}\theta _{1/2}`$ for all the halos is compared to the distribution for low $`\mu `$ halos of $`\mu <1.1`$. The latter clearly show a higher fraction of misaligned halos. Figure 4.1 shows the corrected $`\mathrm{cos}\theta _{1/2}`$ values versus the shape parameter $`\mu `$ for each halo, and indeed there is a significant trend, with $`r0.36`$ ($`p10^{14}`$). 666We also find a trend with the measured value of $`\mathrm{cos}\theta _{1/2}`$ and $`\lambda ^{}`$; however, this seems to be due to the correlation between $`J`$ errors and the spin. The correlation is not apparent when we use our corrected estimate for $`\mathrm{cos}\theta _{1/2}`$.  The misaligned halos tend to be associated with low-$`\mu `$ profiles — that is, profiles that deviate significantly from power laws, with a relatively larger amount of mass in the tails of their $`j_z`$ distributions. This result suggests that there are common aspects to the origin of misalignment and that of the non-powerlaw nature of the $`M(<j)`$ profiles.

Since halos with misaligned angular momentum distributions may be less likely to host large disk galaxies, this result should be kept in mind when modeling galaxy formation within halos with low-$`\mu `$ angular momentum profiles.

Fig. 4.1 addresses the angular-momentum alignment in halos in another way, independent of any specific spatial symmetry. Here we show the $`M(<j)`$ profiles for four different halos, comparing the distributions of both $`j_z`$ and $`|j|`$ in cells throughout the halo volume. We see that in most cases, the two kinds of profiles are similar; for at least half the mass they differ by less than a factor of two. This indicates that the angular momentum is reasonably well aligned throughout the halo.

The slight misalignment is expressed in a trend that the $`|j|`$ distributions systematically resemble a power-law more than the $`j_z`$ distributions. This is demonstrated in Fig. 4.1, which shows the values of $`\mu 1`$ obtained using $`|j|`$ versus those obtained using $`j_z`$. The correlation is obvious, with $`r0.75`$. The value of $`\mu 1`$ is typically about a factor of 3 higher for the $`|j|`$ profiles. The fact that this factor is somewhat higher for low values of $`\mu `$, reflects the correlation seen in Fig. 4.1 between alignment and $`\mu `$.

### 4.2 Cylindrical symmetry

Cylindrical symmetry is the natural symmetry in the spatial distribution of angular momentum as a result of tidal-torques or a sequence of mergers. In order to explore the degree of cylindrical symmetry, we have again divided the halos into cells, as described in §2.3, except that we no longer explicitly impose symmetry about the equatorial plane and rotational symmetry in $`\varphi `$. Specifically, we define the cells as outlined in §2.3, but then, for each value of $`r\mathrm{sin}\theta `$, we further divide the cell particles into four regions: above and below the equatorial plane, and split by rotation angle $`\varphi =0\pi `$ and $`\varphi =\pi 2\pi `$. Because each cell now contains roughly one fourth of the particles, we increased our minimum halo mass threshold by a factor of four to $`4\times 10^{12}h^1M_{}`$, only for this section of our analysis.

Figure 2 shows a map of the $`j`$ distribution in cells in four of our halos. Except for the halo depicted in the lower-left panel, the symmetry of the angular-momentum distribution tends to be more cylindrical than spherical, with a weak dependence on $`\theta `$ at a constant distance from the global spin axis, $`r\mathrm{sin}\theta `$. On the other hand, the angular momentum of the halo in the lower-left panel tends to be concentrated in the $`z=0`$ plane.

Another way of visualizing the degree of cylindrical symmetry of halos is shown in Figure 4.2. It shows the values of $`j`$ in the same cells and for the same four halos displayed in Figure 2, now plotted as a function of the distance from the angular momentum axis, $`r/R_\mathrm{v}\mathrm{sin}\theta `$, and distinguished by different symbols according to their average distance from the equatorial plane, $`|z/R_\mathrm{v}|`$. The different symbol types should be well-mixed if a halo is cylindrically symmetric. This is the case for three out of the four halos shown, except in the lower-left panel, as expected from Figure 2.

Although the angular momentum distributions in most of the halos tend to qualitatively show a certain degree of cylindrical symmetry, there is some indication that, for fixed $`r\mathrm{sin}\theta `$, mass near the equatorial plane tends to have more specific angular momentum than corresponding mass in the poles. In order to quantify the extent to which this is true, we compare the quantity $`j_{z=0}/j_{\mathrm{pole}}`$, defined as the ratio of specific angular momentum in cells near the equatorial plane ($`|z/R_v|=00.35`$) to that in cells near the poles ($`|z/R_v|=0.651`$). For each halo, we determine $`j_{z=0}/j_{\mathrm{pole}}`$ by averaging over all values of $`r\mathrm{sin}\theta `$ for which the two $`z`$ ranges overlap. For example, in Figure 4.2, $`j_{z=0}/j_{\mathrm{pole}}`$ can be estimated for each halo by comparing the average value of $`j`$ for the solid squares ($`|z/R_v|0.17`$) to that for the open circles ($`|z/R_v|0.83`$) at fixed values of $`r\mathrm{sin}\theta `$. Figure 4.2 shows the histogram of this ratio for all of our halos. The mean is $`\mathrm{log}_{10}(j_{z=0}/j_{\mathrm{pole}})0.13`$, and the standard deviation is $`0.18`$.

In order to test whether our specific choice of cell geometry would bias this measure, we also generated 500 realizations of halos with perfect cylindrical symmetry ($`j_{z=0}/j_{\mathrm{pole}}=1`$). They each had NFW density profiles and $`M(<j)`$ profiles in the form of Equation 8. The resulting $`j_{z=0}/j_{\mathrm{pole}}`$ histogram for the artificial catalog is shown by the shaded historgram in Figure 4.2. The sharp peak near a ratio of $`1`$ indicates that the bias due to our cell geometry is much smaller than the detected deviation from cylindrical symmetry. Recall however that the average deviation from cylindrical symmetry is small, about 35%.

### 4.3 Profile in spherical shells

Although, as demonstrated in §4.2, the specific angular momentum distribution in halos is closer to being cylindrically symmetric than spherically symmetric, an exploration of how $`j`$, averaged over shells, behaves as a function of radius is of particular interest because this is the profile that is most directly addressed by theory (see below).

Barnes & Efstathiou (1987) found in their early work an indication for a rough $`j(r)r`$ behavior. For our sample of well-resolved halos, we indeed find that a power-law: $`j(r)r^\alpha `$, with $`\alpha `$ roughly distributed over the halos like a Gaussian: $`\alpha =1.1\pm 0.3`$ is a good approximation.

In order to connect more closely to the mass distribution profiles discussed in previous sections, and also to allow a more direct connection with theory, we focus on characterizing the spherical profile as a function of the mass contained within the radius $`r`$, $`M=M(<r)`$. We find that the power-law approximation describes the profiles reasonably well,

$$j(M)M^s,$$
(12)
with $`s`$ roughly distributed over the halos like a Gaussian: $`s=1.3\pm 0.3`$.

Figure 4.2 shows three examples of the spherical $`j(M)`$ profiles along with their best-fit power-laws. Notice that although the power-laws serve as a good general characterization of the $`j`$ behavior with $`M`$, the profiles do not always follow a power-law form in detail. In many cases, the profiles show an upwards bend in the outer shells, reminiscent of the characteristic bend in the low-$`\mu `$ mass distribution, $`M(<j)`$. In some cases, $`j`$ is not always monotonic with $`M`$. So, although the $`s`$ values provide a useful characteristic of the spherical distribution of $`j`$, some information is lost under this approximation.

Shown in Fig. 4.2 is the distribution of $`s`$ values, versus halo mass. Despite the large scatter, there is an anti-correlation between slope and mass, with $`r0.49`$, ($`p10^{27}`$). As we will discuss in §5, a simple calculation based on linear tidal-torque theory provides a hint at understanding the typical slope value $`s1`$, and predicts a mild anti-correlation between $`s`$ and the halo mass, in agreement with the simulations.

The spherical profile $`j(M)`$ and the mass distribution $`M(<j)`$ are expected to be related. If the mass distribution is spherically symmetric and the angular-momentum distribution is cylindrically symmetric, this relation can be spelled out explicitly. In this case the inverse of $`M(<j)`$ is the spatial profile $`j_{\mathrm{cyl}}(M)`$ in growing cylinders about the total spin axis. A power law with a bend for $`M(<j)`$ translates to a corresponding power law with a weaker bend in the opposite direction for $`j(M)`$. We thus expect a correlation between the shape parameter $`\mu `$ of $`M(<j)`$ and the slope $`s`$ of the spherical profile $`j(M)`$. A high $`s`$ and a low $`\mu `$ both correspond to halos with a high fraction of mass spinning slowly relative to $`j_{\mathrm{max}}`$. Figure 4.2 shows $`s`$ and $`\mu `$ for our halos. An anti-correlation is evident, as expected, with $`r0.36`$ ($`p10^{13}`$). The scatter reflects deviations from the global spatial symmetries. It is likely, for example, that partial misalignment and pockets of low-angular momentum material at large radii contribute to the scatter. By comparing Figures 4.2 and 4.1 it is evident that the scatter in the $`s`$-$`\mu `$ relation is large enough to wash away any mass dependence of $`\mu `$ that would be implied by the weak mass dependence of $`s`$.

## 5 On the origin of the profile

In this section we explore the possible origin of the angular momentum profile based on two different pictures for how specific angular momentum is acquired in halos. In the first picture, angular momentum is built up shell by shell, and is modeled using linear theory with no further angular momentum transfer within the halo. In the second picture, the angular momentum profile arises solely from angular momentum transfer from one or more satellite halo merger events. Since halo mass accretion histories typically reflect some combination of relatively quiescent mass accretion as well as more pronounced mergers (where $`j`$ transfer is likely), the $`j`$ profiles of individual halos may reflect a complicated combination of the two processes explored here. Indeed the two processes are intrinsically linked, since, in the hierarchical framework, all mass accretion can be treated as the accretion of smaller objects. As shown below, the profiles calculated from linear theory as well as those calculated using $`j`$ transfer from mergers produce a range of $`\mu `$ values, in qualitative agreement with what is observed in the simulations.

### 5.1 TTT and EPS

Linear tidal-torque theory (TTT) provides a hint for the origin of the detected power-law-like $`j`$ profile, as follows. TTT (Doroshkevich 1970; White 1984) implies that the angular momentum gained by a halo at time $`t`$ before its turn-around is

$$J_i(t)=a(t)^2\dot{D}(t)ϵ_{ijk}T_{jl}I_{lk},$$
(13)
where the time growth is from some fiducial initial time, $`a(t)`$ is the expansion factor at $`t`$, $`D(t)`$ is the linear growth factor, $`I_{lk}`$ is the inertia tensor of the proto-halo at the initial time and $`T_{jl}`$ is the tidal tensor at the halo center at the initial time, smoothed on the halo scale. This is based on assuming the Zel’dovich approximation for the velocities inside the proto-halo, and a 2nd-order Taylor expansion of the potential. We show elsewhere (Porciani, Dekel & Hoffman 2000) that the standard scaling relation of TTT should be slightly modified; it should read

$$jD(t_\mathrm{c})^{3/2}\sigma (M)M^{2/3},$$
(14)
where $`j`$ is the specific angular momentum of the halo, $`t_\mathrm{c}`$ is the turn-around time of the halo, and $`\sigma (M)`$ is the rms density fluctuation on scale $`M`$ at the initial time. If the perturbation turned around while the cosmology was still Einstein-de Sitter, then $`jt_\mathrm{c}`$777 The $`j`$ dependence on the time of collapse is $`t_\mathrm{c}`$ rather than $`t_\mathrm{c}^{1/3}`$, because the tidal part of the deformation tensor, the source of angular momentum, is actually independent of the initial density of the perturbation, which determines the collapse time. The shear depends instead on the rms fluctuation $`\sigma `$, which involves an implicit mass dependence.

To understand the general power-law behavior of the angular-momentum profile we apply Eq. (14). We can assume that mass is accreted in shells and that the turn-around time of each shell is determined by $`\overline{\delta }(M)D(t_\mathrm{c})1`$, where $`\overline{\delta }(M)`$ is the mean density inside $`M`$. For a Gaussian field, the typical density fluctuation profile about a random point scales like (Dekel 1981) $`\delta (r)\xi (r)`$, where $`\xi (r)`$ is the linear two-point correlation function. This is accurate to a few percent also around a high peak (Bardeen et al. 1996, Fig. 8). Thus, for a power-law power spectrum $`P_kk^n`$ and a flat universe, we obtain the following spin profile within each halo:

$$j(M)M^{2/3+(3+n)/3}.$$
(15)
This implies $`1<s<4/3`$ for $`2<n<1`$, the range appropriate for $`10^{12}10^{14}h^1M_{}`$ halos in a CDM spectrum, in pleasant agreement with our finding for the simulated halos (Fig. 4.2).

The tidal torque theory can be combined with the Extended Press-Schechter formalism (EPS, Bond et al. 1991; Lacey & Cole 1993) to obtain a more detailed model of $`j`$-distribution, as well as quantitative predictions for values of $`s`$ and $`\mu `$ and their mass dependence.

Let us assume for simplicity that halo mass is acquired via accretion of material with specific angular momentum given by Equation 14 and that the direction of angular momentum of each mass shell is perfectly aligned.

We can readily estimate the slope $`s`$ of the $`j(M)M^s`$ distribution if we assume that the turn-around time for each mass shell is approximated by the time at which the mass was first accreted onto the halo888The time of accretion is expected to be proportional to the turnaround time; see e.g. Sugerman et al. (2000). by comparing the expected $`j`$ value at any two shells encompassing masses $`M_1`$ and $`M_2`$,

$$s\frac{\mathrm{ln}[j(M_2)/j(M_1)]}{\mathrm{ln}(M_2/M_1)}.$$
(16)
Equation 14 provides the required value of $`j`$ associated with each mass shell given its “accretion time”, $`t_c(M)`$. Considering, for example, $`M_\mathrm{v}`$ and $`M_\mathrm{v}/2`$, we obtain

$$s(M_\mathrm{v})1.44\mathrm{ln}\left(2^{2/3}D[t_{1/2(M_\mathrm{v})}]^{3/2}\frac{\sigma (M_\mathrm{v})}{\sigma (M_\mathrm{v}/2)}\right),$$
(17)
where $`t_{1/2}`$ is the time at which the halo of mass $`M_\mathrm{v}`$ first accreted half of its mass, $`t_{1/2}(M_\mathrm{v})t_c(M_\mathrm{v}/2)`$. We adopted for all halos, by the requirement of virialization, $`t_c(M_\mathrm{v})=t_0`$ (the current age of the universe), and set $`D(t_0)=1`$. The most probable value of $`t_{1/2}(M_\mathrm{v})`$ can be estimated via EPS (Lacey & Cole 1993, Eq. 2.26). The predicted $`s(M_\mathrm{v})`$ relation is added to Fig. 4.2 as the thick solid line. The agreement with the mean trend shown by the simulated halos is remarkable, especially since our model has no free parameters.

We can further elaborate the model by following the more detailed mass accretion history of halos using EPS merger trees. Specifically, we have used the method of Somerville & Kolatt (1999) to model mass growth and halo accretion histories for a random ensemble of dark matter halos formed in the $`\mathrm{\Lambda }`$CDM model and with masses in the mass range studied in the simulation. At each epoch, $`t_i`$, and time interval, $`\mathrm{\Delta }t_i`$, the model draws a mass $`\mathrm{\Delta }M_i`$ accreted by a halo of mass $`M_i(t_i)`$ using the EPS probability distribution. If we assume that this mass is accreted with specific angular momentum given by Eq. (14) and that there is no angular momentum loss, we can integrate the specific halo mass accretion history to get the $`M(<j)`$ distribution for the final halo at $`z=0`$. Figure 5.1 shows $`M(<j)`$ profiles calculated using this model along with the fits of the form $`M(<j)=M_\mathrm{v}\mu j(j+j_0)^1`$, Eq. (8), for three representative mass accretion histories from the ensemble of realizations. One can see that the model reproduces the general shape of the $`M(<j)`$ profile with values of $`\mu `$ that are similar to those of simulated halos.

For a given halo mass, the differences in merger histories lead to scatter in the resulting values of $`\mu `$, with the range of $`\mu 10.12.5`$. The absence of small $`\mu `$ ($`1.1`$) may result from our simplified model assumption of perfect alignment. Indeed, the range of these modeled $`\mu `$ values is similar to the range obtained for our simulated halos when $`|j|`$ profiles are considered instead of $`j_z`$ (see Fig. 4.1).

We can also reproduce both the range of values and the trend with mass of the slope $`s`$ if we force a power-law fit to the model $`j(M)`$ profiles. The halos with low (high) $`\mu `$ in this model are those that accrete large fractions of their mass at early (late) epochs.

Interestingly, the value of $`\mu `$ for a given halo in this model strongly depends on the halo’s formation redshift, $`z_f`$, defined as the redshift where the mass of the most massive progenitor was half of the final mass for the first time \[$`\mu 1(1+z_f)^2`$\]. In light of this expectation, it is somewhat surprising that we find no noticeable trend with redshift for $`\mu `$ parameters in halos of fixed mass (see Appendix) but we are currently investigating whether a trend is evident using individual halo merger histories (Wechsler et al. 2001). If there is a trend with formation time, it may have interesting implications for galaxy formation.

### 5.2 Profile due to minor mergers

A simple toy calculation for how angular momentum is being deposited in a large halo during a minor merger seems to provide another hint for the origin of the characteristic $`j`$ profile. This is studied in more detail using toy models and full $`N`$-body simulations by Dekel & Burkert (2000).

Consider a fixed halo and an incoming satellite halo of mass profiles $`M(r)`$ and $`m(\mathrm{})`$ respectively. Assume that as it moves in, the satellite is losing mass outside a tidal radius $`\mathrm{}_\mathrm{t}`$, which we approximate at $`r`$ to be determined by

$$\frac{m(\mathrm{}_\mathrm{t})}{\mathrm{}_\mathrm{t}^2}=\frac{\mathrm{}_\mathrm{t}\kappa (r)}{r^3},\kappa (r)=\left[2M(r)r\frac{\mathrm{dM}}{\mathrm{dr}}\right].$$
(18)
For host halos and satellite halos with power-law mass profiles ($`Mr^\alpha ,\alpha <2`$), the above approximation implies that the total mass of a satellite halo located at a radius r of its host is $`m(r)m[\mathrm{}_t(r)]M(r)`$.

Suppose that the satellite is spiraling in due to dynamical friction roughly along circular orbits, and that the mass and $`j`$ are deposited locally. Then the spherical $`j(r)`$ profile is obtained by averaging over shells,

$$4\pi r^2\rho (r)j(r)=m(r)\frac{\mathrm{d}[\mathrm{rV}_\mathrm{c}(\mathrm{r})]}{\mathrm{dr}}+\frac{\mathrm{dm}(\mathrm{r})}{\mathrm{dr}}rV_\mathrm{c}(r).$$
(19)
The two terms on the right hand side reflect angular-momentum transfer due to the slowdown of the satellite by dynamical friction and the direct tidal stripping of mass from the satellite halo, respectively. For an isothermal host and satellite, we obtain

$$j(M)M.$$
(20)
Fig. 5.2 shows the expected $`j(M)`$ profile based on this toy model for an NFW halo swallowing a satellite, using three alternative recipes for mass loss by the satellite: (a) $`m(r)M(r)`$, (b) a satellite with an isothermal profile, and (c) a satellite of fixed mass, for comparison. The general shape of $`j(M(<r))`$ in the outer decade of mass, when realistic mass loss (a or b) is considered, is reminiscent of the profiles detected in the cosmological simulation (see Fig. 4.2).

Full $`N`$-body simulations of mergers, spanning a range of halo and collision parameters, confirm the robust production of such characteristic profiles (Dekel & Burkert 2000). The resulting $`j(M)`$ from a sequence of minor mergers is expected to be a sum of similar contributions, each projected onto the direction of the total net angular momentum (perhaps determined by the most major merger involving a large fraction of the final mass).

## 6 Disk Structure

We now explore certain possible implications of our results on the angular-momentum distribution in the halos for the formation of galactic disks. We follow the limiting case assumption that the specific angular momentum is conserved during gas infall into a centrifugally supported disk, and compare the resultant gas density profile to an exponential disk. Also, since the angular momentum structure of halos is governed by at least two parameters — say $`\mu `$, in addition to $`\lambda ^{}`$ — we specifically explore how variations in the shape parameter may affect the disk characteristics.

Suppose that each dark-matter halo has a specific angular momentum profile $`M(<j)`$ given by Eq. (8). Given $`\mu `$, $`M(<j)`$ is determined by the spin parameter via $`j_0=\sqrt{2}\lambda ^{}V_\mathrm{v}R_\mathrm{v}b(\mu )^1`$. Assume that initially the specific angular momentum distribution of the gas mirrors that of the halo, $`M_{\mathrm{gas}}(<j)=fM_\mathrm{v}(<j)`$, where $`f`$ is the mass fraction of the halo that ends up as cool baryons in the disk, $`M_{\mathrm{disk}}=fM_\mathrm{v}`$.

After infall into a disk while preserving $`j`$, the specific angular momentum of a gas element that ends up in a circular orbit of radius $`r`$ is $`j(r)\sqrt{GM(r)r}`$, where $`M(r)`$ is the total mass of DM and baryons within $`r`$. All the gas with specific angular momentum less than $`j(r)`$ will wind up interior to $`r`$. Using Eq. (8), we obtain for the disk mass profile

$$m_\mathrm{d}(r)f\mu M_\mathrm{v}\frac{j(r)}{j_0+j(r)},j(r)<j_{\mathrm{max}}.$$
(21)
In order to estimate the implied disk structure, and to gain a qualitative understanding of how the value of $`\mu `$ may affect the disk size and density profile, we assume first that the dark halo does not react to the infall of gas and approximate the total mass distribution of the system to be that of an isothermal sphere, $`M(r)r`$. In this case $`j(r)rV(r)=rV_\mathrm{v}`$. The final mass distribution of the disk is then

$$m_\mathrm{d}(r)f\mu M_\mathrm{v}\frac{r}{r_\mathrm{d}+r},r<r_{\mathrm{max}}.$$
(22)
Here, $`r_\mathrm{d}\sqrt{2}\lambda ^{}R_\mathrm{v}b(\mu )^1`$ and $`r_{\mathrm{max}}=r_\mathrm{d}/(\mu 1)`$. The implied surface density profile is

$$\mathrm{\Sigma }_\mathrm{d}(r)=\frac{f\mu M_\mathrm{v}}{2\pi }\frac{r_\mathrm{d}}{r(r_\mathrm{d}+r)^2},r<r_{\mathrm{max}}.$$
(23)
In our simplified model, the surface density vanishes beyond $`r_{\mathrm{max}}`$.

Several implications are apparent from the derived disk surface-density profile of Eq. (23). First, the distribution is qualitatively similar to an exponential disk — it falls rapidly at large $`r`$, and flattens towards the origin. However, in detail, there are significant differences. At large radii, $`rr_d`$, the fall off $`\mathrm{\Sigma }_\mathrm{d}(r)r^3`$ is slower than exponential. Near the origin, $`rr_d`$, the profile is steeper than exponential, $`\mathrm{\Sigma }_\mathrm{d}(r)r^1`$ versus $`\mathrm{const}`$. As expected, the scale radius $`r_d`$ scales linearly with $`\lambda ^{}`$. However, the scale radius is also a strongly decreasing function of $`\mu `$. For a fixed spin parameter $`\lambda ^{}`$, as $`\mu `$ varies over its $`90\%`$ range about the mean ($`\mu 1=0.061.0`$), $`r_\mathrm{d}`$ varies by about a factor of $`5`$. This is similar to the variation of $`r_\mathrm{d}`$ as a function of $`\lambda ^{}`$ for fixed $`\mu `$ (a factor of $`5`$ over the 90% range of $`\lambda ^{}`$).

Figure 6 shows the disk surface-density profile derived using a more accurate numerical calculation. Here we assume that the dark halo initially followed an NFW density profile with $`c_\mathrm{v}=14`$. The $`M(<j)`$ profile is characterized by $`\lambda ^{}=0.04`$ and $`\mu =1.25`$, and the disk mass fraction is assumed to be $`f=0.03`$. The disk profile is obtained by self-consistently solving for the disk circular velocity at $`r`$ and taking into account the adiabatic contraction of the halo in response to the baryonic infall (Blumenthal et al. 1986). Shown in comparison is the surface density profile obtained under the same assumptions regarding the halo structure and $`j`$ conservation except that the given $`M(<j)`$ distribution is replaced with the assumption that the final disk profile is exponential. A similar result is obtained by assuming that the initial $`M(<j)`$ is that of a uniform-density sphere in solid body rotation.

Figure 6 shows how the derived disk profile varies as a function of $`\mu `$. As expected, the central concentration of the disk is a decreasing function of $`\mu `$.

As anticipated, compared to the exponential disk (and solid body rotator) the derived disk profiles have core and tail excesses. Only halos with the highest $`\mu `$ values ($`\mu 2`$, corresponding to $`5\%`$ of the halos) produce surface density profiles that begin to resemble an exponential disk.

In comparison with the observed light profiles of disks, the derived gas surface density profiles deviate in their general shape and in particular in their central concentration. The general deviation of shape is less worrisome because the final stellar profile may not necessarily mirror the gas profile. As Lin & Pringle (1987) and Olivier et al. (1991) argued, viscous transport processes can act to produce exponential stellar disks under rather general assumptions about the nature of the initial gas disk. However, since typical viscous processes would tend to amplify the central density even further, the indicated central density excess is of some concern.

This result is reminiscent of the more severe problem of excessively small disk sizes seen in hydrodynamic cosmological simulations (e.g., Navarro & Steinmetz 1997, 2000). That we find a disk central-density excess even under the assumption that angular momentum is conserved during the collapse indicates that the solution to the problem may not be trivial, because any simple angular-momentum transport mechanism would drive the angular-momentum from the gas to the halo and from the inside out and could only worsen the discrepancy. The solution may require that the baryons somehow obtain an angular momentum distribution that is biased relative to that of the dark matter towards high-$`j`$.

Alternatively, if the resultant profiles are interpreted to represent both a bulge and a disk component, with the slowly rotating core associated with the bulge, then the central concentration problem may be resolved. Testing this hypothesis will require more detailed modeling. In particular, it will be necessary to match the observed distribution of galaxy bulge-to-disk ratios as well as disk scale lengths.

Another way out may be associated with the fact that some of the approximations adopted here may break down when the $`j`$ distribution is significantly misaligned throughout the halo, which indeed preferentially happens for low-$`\mu `$ halos (§4.1).

## 7 conclusions

We have studied the angular momentum structure within dark-matter halos in the mass range $`10^{12}10^{14}h^1M_{}`$, using a sample of $`600`$ halos in an N-body simulation of the $`\mathrm{\Lambda }`$CDM cosmology. We examined in detail the more massive of these halos, for which the angular momentum is measured more accurately, and then verified that the more typical galaxy-sized halos have a similar angular momentum distribution.

Our primary result is that the mass distribution of specific angular momentum in halos obeys a roughly universal form, which is well fit by the two-parameter function $`M(<j)=\mu M_\mathrm{v}j/(j_0+j)`$. For a fixed value of the spin parameter $`\lambda ^{}`$, the parameter $`\mu `$ determines the shape of the distribution, with high-$`\mu `$ corresponding to a pure power-law and a smaller contrast between the tails of the distribution. The distribution of $`\mu 1`$ is roughly log normal, with the $`90\%`$ range spanning $`\mu 1=0.061.0`$. Compared to an idealized uniform sphere in solid-body rotation, the simulated halos tend to have more of their mass in the tails of the distribution, especially at small $`j`$. The shape parameter is only weakly correlated with the spin parameter, $`\mu 1(\lambda ^{})^\delta `$ with $`\delta =2\pm 1.5`$ and a correlation coefficient $`r0.23`$.

Most halos have well-aligned angular momentum vectors throughout their volume. Between 70% and 90% have alignment cosines greater than $`0.7`$ between the inner and outer half masses. However, at least 10% of the halos have significant misalignments, which may be interesting for several reasons. First, they tend also to be halos with low $`\mu `$, indicating that the direction of J in accreted material, or the halo merger history, plays a role in determining the shape of the $`M(<j)`$ profile. In addition, it is unlikely that the standard simple picture of disk formation could be valid within a severely misaligned halo, so perhaps these objects typically host spheroidal stellar components. Marginally misaligned halos may play a role in the formation of galactic warps (e.g., Dekel & Schlosman 1983).

The spatial distribution of halo angular momentum tends to be more cylindrically symmetric than spherically symmetric. At a fixed distance from the rotation axis, mass near the equatorial plane typically has only about $`35\%`$ more specific angular momentum than corresponding mass near the poles.

The mean $`j`$ is spherical shells encompassing mass $`M`$ is well-fit by a power law, $`jM^s`$. The power $`s`$ is distributed like a Gaussian with a mean of $`s=1.3`$ and standard deviation $`\sigma =0.3`$.

We pointed out two possible explanations for the origin of the universal angular-momentum profile. The first is based on applying a corrected version of the linear tidal-torque theory (Porciani & Dekel 2000) to extended Press-Schechter mass accretion histories. The second is based on the nonlinear process of angular momentum transfer from satellite orbits in halo mergers (see also Dekel & Burkert 2000). Since halo mass growth arises as a combination of merger events and relatively quiescent mass accretion, the origin of the profile in individual halos may reflect some combination of the two processes. Each of these processes seems to produce $`j`$ profiles that are similar in form to those observed in simulated halos, and thus may provide a starting point for a deeper understanding of the origin of halo angular-momentum structure.

Finally, we have started to explore the implications of our universal $`M(<j)`$ distribution in the context of the archetypical model of galactic disk formation, namely, adiabatic baryonic infall to a rotationally supported disk while conserving angular momentum in every mass element. We find that the implied surface-density profiles, which vary as a function of both $`\mu `$ and $`\lambda ^{}`$, deviate significantly from an exponential disk for all but the largest values of $`\mu `$. The resultant surface profiles are more extended than exponential at large $`r`$, and are overly concentrated at small $`r`$. Since the observed light profiles of disks are closer to exponential, our result indicates that the disk formation process cannot be fully understood within this simplified picture.

We mentioned that the general deviation from an exponential surface density at moderate radii might be overcome in this picture if a viscous transport mechanism is included (Lin & Pringle 1987). As shown by Olivier et al. (1991), such a mechanism tends to produce exponential stellar distributions regardless of the initial gas profile. However, most of the expected $`j`$ transport typically acts to increase the central mass density, so the deviation at small radii may be a significant problem. Even though we have assumed that the baryons experience no angular-momentum loss, we ended up with an “angular momentum problem” similar to that detected in hydrodynamic simulations (e.g., Navarro & Steinmetz 1997, 2000, where it is associated with $`j`$ transport from the gas to the dark matter and from inside out). This problem is most severe for halos with small $`\mu `$, so the tendency for these halos to have misaligned angular-momentum distributions, which makes them less-likely hosts of large disk galaxies, may partly ease the problem in the limiting case where angular-momentum is conserved. However, even if we focus on halos that tend to be well aligned, with $`\mu 1.1`$, the central densities remain higher than those in exponential disks.

A solution to this problem may be to associate the derived central mass concentrations with bulges. It is indeed possible that the central regions of disks with profiles as in Fig. 6 are unstable to self gravity. More detailed modeling is needed for testing the viability of this solution. In particular, one may be worried about the model over-predicting bulge to disk ratios and in particular under-predicting the number of bulge-less galaxies. On the other hand, such an interpretation may shed new light on the very origin of the Hubble sequence, since both $`\mu `$ and $`\lambda ^{}`$ may play a role in determining disk scale lengths, surface brightnesses, and bulge fractions.

## Acknowledgments

This work has been supported by grants from the US-Israel BSF, the Israel SF, and NASA and NSF at UCSC and NMSU. JSB was supported by NASA LTSA grant NAG5-3525 and NSF grant AST-9802568. Support for AVK was provided by NASA through Hubble Fellowship grant HF-01121.01-99A from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555. CP was supported by a Golda Meir fellowship. We thank Andi Burkert, Stephan Corteau, Sandra Faber, Ricardo Flores, Savvas Koushiappas, Ari Maller, Yair Sigad, Rachel Somerville, Frank van den Bosch, Maya Vitvitska, David Weinberg, and Risa Wechsler for useful discussions.

Since our error estimates are only approximate, it is important to carry an independent test of the validity of our $`M(<j)`$ profiles in view of the finite mass resolution of the simulation. The test utilizes the output from another $`\mathrm{\Lambda }`$CDM simulation with identical cosmological parameters and the same number of particles as used in our main simulation, but now in a box of side $`30h^1\mathrm{Mpc}`$, namely half the original size. The mass resolution is thus $`8`$ times higher, so a statistical comparison with halos analyzed from our lower resolution simulation at fixed mass will provide a useful test for the effects of mass resolution. The smaller box simulation also allows us to extend our analysis down to much smaller galaxy-mass halos ($`M_\mathrm{v}10^{11}h^1M_{}`$). Unfortunately, the higher resolution simulation was stopped at $`z=1.7`$, so we cannot compare directly to the $`z=0`$ results presented in the main paper. However, a comparison performed at high $`z`$ can perfectly serve our purpose. This also allows us the opportunity to check for any evolution in the $`M(<j)`$ profiles with redshift.

Figure Acknowledgments shows the best-fit $`\mu `$ parameters as a function of halo mass for halos from our main simulation at $`z=3`$. The distribution of $`\mu `$ values is very similar to that seen at $`z=0`$ (Fig. 4.1), so there is no indication of evolution. Figure Acknowledgments also shows $`\mu `$ versus $`M_\mathrm{v}`$ for halos from the higher resolution simulation at the same time. Although the overlapping region in mass is rather small, there is no indication of any offset between the two simulations. The halos have a similar distribution of shape parameters over about three orders of magnitude in mass, $`10^{11}10^{14}h^1M_{}`$.

We conclude that mass resolution does not seem to limit our analysis, and also that the $`M(<j)`$ profiles of halos show very little variation as a function of redshift and mass.


format: text

text:
# Probability-Changing Cluster Algorithm for Potts Models

## Abstract

We propose a new effective cluster algorithm of tuning the critical point automatically, which is an extended version of Swendsen-Wang algorithm. We change the probability of connecting spins of the same type, $`p=1e^{J/k_BT}`$, in the process of the Monte Carlo spin update. Since we approach the canonical ensemble asymptotically, we can use the finite-size scaling analysis for physical quantities near the critical point. Simulating the two-dimensional Potts models to demonstrate the validity of the algorithm, we have obtained the critical temperatures and critical exponents which are consistent with the exact values; the comparison has been made with the invaded cluster algorithm.

PACS numbers: 05.50.+q, 64.60.Ak, 75.10.Hk

Cluster algorithms  have been successfully used to overcome slow dynamics in the Monte Carlo simulation. Swendsen and Wang (SW)  applied the Kasteleyn-Fortuin (KF)  representation to identify clusters of spins. The problem of the thermal phase transition is mapped onto the geometric percolation problem in the cluster formalism . Quite recently, based on the cluster formalism, the multiple-percolating clusters of the Ising system with large aspect ratio have been studied .

Machta et al.  proposed another type of cluster algorithm, which is called the invaded cluster (IC) algorithm; this algorithm samples the critical point of a spin system without a priori knowledge of the critical temperature. It is in contrast with the usual procedure that one makes simulations for various parameters to determine the critical point. The IC algorithm has been shown to be efficient in studying various physical quantities in the critical region, but the ensemble is not necessarily clear. Moreover, it has a problem of “bottlenecks”, which causes the broad tail in the distribution of the fraction of the accepted satisfied bonds .

In this Letter, extending the SW algorithm , we propose a new algorithm of tuning the critical point automatically. The basic idea of our algorithm is that we change the probability of connecting spins of the same type, $`p=1e^{J/k_BT}`$, in the process of the Monte Carlo spin update, where $`J`$ is the exchange coupling. We decrease or increase $`p`$ depending on the observation whether the KF clusters are percolating or not percolating. This simple negative feed-back mechanism together with the finite-size scaling (FSS) property of the existence probability (also called the crossing probability) $`E_p`$, the probability that the system percolates, leads to the determination of the critical point. Since our ensemble is asymptotically canonical as $`\mathrm{\Delta }p`$, the amount of the change of $`p`$, becomes 0, the distribution functions of physical quantities obey the FSS; as a result, we can determine critical exponents using the FSS analysis.

Let us explain the procedure for our probability-changing cluster (PCC) algorithm in detail. As an example, we consider the ferromagnetic $`q`$-state Potts model  whose Hamiltonian is given by

$$=J\underset{<i,j>}{}(\delta _{\sigma _i,\sigma _j}1),\sigma _i=1,2,\mathrm{},q,$$
(1)
and for $`q`$ = 2 this corresponds to an Ising model. The procedure of Monte Carlo spin update is as follows:

1. Start from some spin configuration and some value of $`p`$.
2. Construct the KF clusters using the probability $`p`$, and check whether the system is percolating or not. Update spins following the same rule as the SW algorithm, that is, flip all the spins on any KF cluster to one of $`q`$ states.
3. If the system is percolating (not percolating) in the previous test, decrease (increase) $`p`$ by $`\mathrm{\Delta }p(>0)`$.
4. Back to the process 2.

After repeating the above processes, the distribution of $`p`$ for our Monte Carlo samples approaches the Gaussian distribution of which mean value is $`p_c(L)`$; $`p_c(L)`$ is the probability of connecting spins, such that the existence probability $`E_p`$ becomes 1/2. The width of the distribution depends on the choice of $`\mathrm{\Delta }p`$ in the process 2, and becomes 0 in the limit of $`\mathrm{\Delta }p0`$. We should note that $`p_c(L)`$ depends on the system size $`L`$, and $`E_p`$ follows the FSS near the critical point,

$$E_p(p,L)X(tL^{1/\nu }),t=(p_cp)/p_c,$$
(2)
where $`p_c`$ is the critical value of $`p`$ for the infinite system ($`L\mathrm{}`$), and $`\nu `$ is the correlation-length critical exponent. (As for the FSS of $`E_p`$, see Ref. , for example.) We can estimate $`p_c`$ from the size dependence of $`p_c(L)`$ using Eq. (2) and, in turn, estimate $`T_c`$ through the relation $`p_c=1e^{J/k_BT_c}`$. We have chosen the value of $`E_p`$ which gives $`p_c(L)`$ as 1/2 because it is the simplest. We may modify the update process such that this value is different from 1/2.

A comment should be made here on the choice of criterion to determine percolating. Machta et al.  used both the extension rule and the topological rule for their stopping condition. The former rule is that some cluster has maximum extent $`L`$ in at least one of the $`d`$ directions in $`d`$-dimensional systems. The latter rule is that some cluster winds around the system in at least one of the $`d`$ directions. We can use any rule to determine percolating, but FSS functions for physical quantities, therefore $`p_c(L)`$, depend on the rule.

There is one free parameter in our algorithm; we may choose the difference $`\mathrm{\Delta }p`$ in the process 3. In the limit of small $`\mathrm{\Delta }p`$ we obtain the canonical ensemble, but it takes a long time to equilibrate for small $`\mathrm{\Delta }p`$. Practically, we may start with rather large $`\mathrm{\Delta }p`$, and switch to smaller $`\mathrm{\Delta }p`$ with monitoring the trail of the values of $`p`$. Small steps of preparation are enough for equilibration.

In order to show the validity of the present method, we have made simulations for the $`2d`$ ferromagnetic 2-state Potts model (Ising model) and 3-state Potts model. We have treated the systems with linear sizes $`L`$ = 64, 128, 256, and 512. We start with $`\mathrm{\Delta }p=0.01`$, and gradually decrease $`\mathrm{\Delta }p`$ to the final value. We have chosen this final value of $`\mathrm{\Delta }p`$ as $`1/(20\times L^2)`$; the steps for preparation are 10,000 for the largest size ($`L=512`$). After reaching the final small value of $`\mathrm{\Delta }p`$, we have taken 100,000 (200,000) Monte Carlo samples in the case of $`q=2`$ ($`q=3`$) with keeping $`\mathrm{\Delta }p`$ as constant. From each bond configuration, we have made 5 (10) spin configurations in order to get better statistics for magnetization in the case of $`q=2`$ ($`q=3`$). We have performed several runs for each size, and have checked the statistical errors.

Let us show the results of the $`2d`$ Ising model ($`q=2`$). In Fig. 1, we plot the size dependence of $`T_c(L)`$ for both rules in units of $`J/k_B`$. We have determined $`p_c(L)`$ from the average of $`p`$, and calculated $`T_c(L)`$ through the relation $`p_c(L)=1e^{J/k_BT_c(L)}`$. In this plot, as an illustration, we have used the known value of $`\nu `$ for the $`2d`$ Ising model ($`\nu =1`$). Using the least square method, we estimate $`T_c`$ as 1.1344(2) (1.1346(2)) for the extension (topological) rule, which is consistent with the exact value, $`[\mathrm{ln}(1+\sqrt{2})]^1`$ = 1.1346. Here, the number in the parenthesis denotes the uncertainty in the last digit. We have used the known value of $`\nu `$ but we may treat $`\nu `$ as an unknown parameter to be determined. Assuming the FSS relation, $`T_c(L)=T_c+aL^{1/\nu }`$, which is derived from Eq. (2), we may follow the three-parameter ($`T_c,1/\nu ,a`$) non-linear fitting procedure. Then, we have obtained ($`T_c,1/\nu `$)=(1.1345(2), 1.00(4)) for the extension rule, and (1.1344(2), 1.04(4)) for the topological rule. Both estimates of $`T_c`$ and $`\nu `$ are consistent with the exact values.

We plot the distribution of $`p`$, $`f(p)`$, for $`L=64`$ with the topological rule, as an example, in the inset of Fig. 1, which shows that the distribution of $`p`$ is sharply peaked at 0.58196 with the standard deviation of 0.0004 for our choice of the final $`\mathrm{\Delta }p`$. For comparison, we also show $`f(p)`$ for the IC algorithm of the same size with the same rule; different scales are used for the vertical axis in the inset to express two quite different data. We notice that the distribution of $`p`$ for the IC algorithm is far broader. A simple linear analysis yields that the width of the distribution for the PCC algorithm is proportional to $`\sqrt{\mathrm{\Delta }p/a}`$, where $`a`$ is the value of $`dE_p/dp`$ at $`p_c(L)`$; using the FSS we expect $`aL^{1/\nu }`$. It should be noted that we have obtained the expected $`\mathrm{\Delta }p`$\- and $`L`$-dependence for the width of the distribution of $`p`$. We use the average value of $`p`$ for the estimate of $`p_c(L)`$. Performing 10 runs, we have estimated $`p_c(64)`$ as $`0.581954\pm 0.000013`$; the statistical errors are very small.

The resulting energy histogram, $`f(E)`$, of our PCC algorithm for $`L=64`$ with the topological rule is given in Fig. 2. The energy histogram obtained by the constant-temperature calculation using the SW algorithm is shown by a solid curve in Fig. 2; the temperature has been chosen as $`T_c(L)`$ determined by the PCC algorithm. The energy histogram of the PCC algorithm is indistinguishable from that by the constant-temperature calculation because of the sharp distribution of $`p`$. Thus, we may say that our ensemble is actually canonical for small enough $`\mathrm{\Delta }p`$. The energy histogram of the IC algorithm, which is also given in Fig. 2, has broad tails for both high-energy and low-energy sides. Although our ensemble is almost canonical, there are deviations in physical quantities, in principle; the variance of energy, $`E^2E^2`$, becomes larger with non-zero $`\mathrm{\Delta }p`$, for example. We have checked the $`\mathrm{\Delta }p`$-dependence of the systematic deviation for large $`\mathrm{\Delta }p`$. However, the deviation of the variance of energy is smaller than the statistical error, 2 % for $`L=64`$, with our choice of $`\mathrm{\Delta }p`$.

In order to estimate another critical exponent $`\beta `$, the magnetization exponent, we plot the average of the squared magnetization $`m^2`$ as a function of $`L`$ in logarithmic scale in Fig. 3. Since our Monte Carlo samples are sharply peaked at $`p=p_c(L)`$, in other words, at $`T=T_c(L)`$, we can use the FSS relation,

$$m^2_{T=T_c(L)}L^{2\beta /\nu },$$
(3)
for the estimate of $`\beta /\nu `$. From the slopes of the data for both rules, we have $`\beta /\nu `$ = 0.125(2) (0.126(2)) for the extension (topological) rule, which is consistent with the exact value, 1/8 (=0.125).

It is quite interesting to study the distribution function of physical quantities. We show the FSS plot of the distribution function $`p(m)`$ in Fig. 4, based on the FSS relation,

$$p(m)_{T=T_c(L)}L^{\beta /\nu }f(mL^{\beta /\nu }).$$
(4)
The scaling plot for the extension rule (a) and that for the topological rule (b) are given there. The data for various sizes are collapsed on a single curve. We have very good FSS behavior for both rules. Two rules give different $`tL^{1/\nu }`$ in Eq. (2) for $`E_p=1/2`$. It is easier to percolate for the extension rule compared with the topological rule. Therefore, $`T_c(L)`$ of the extension rule is higher than that of the topological rule, which results in the difference in the FSS functions for $`p(m)`$ between two rules.

Next turn to the case of the 3-state Potts model. The size dependence of $`T_c(L)`$ for both rules is shown in Fig. 5. We have used the known value of $`\nu `$ for this plot; the exponent $`\nu `$ for the $`2d`$ 3-state Potts model conjectured by the conformal field theory  is 5/6. We estimate the extrapolated value of $`T_c`$ as 0.99490(6) (0.99494(6)) for the extension (topological) rule from Fig. 5. This value is consistent with the exact value, $`[\mathrm{ln}(1+\sqrt{3})]^1`$ = 0.99497. The convergence is very good for the 3-state Potts model. It is in contrast with the situation of the IC algorithm , where the convergence is not good enough maybe due to the problem of “bottlenecks”.

We can estimate the critical exponent $`\beta /\nu `$ for the 3-state Potts model from the size dependence of $`m^2`$ as in the case of the Ising model. For the order parameter of the 3-state Potts model, we use the vector order parameter $`(m_x,m_y)`$. The $`x`$ and $`y`$ components of the vector order parameter are obtained from the three components, $`m_1,m_2`$ and $`m_3`$, as $`m_x=m_1(1/2)(m_2+m_3)`$ and $`m_y=(\sqrt{3}/2)(m_2m_3)`$. Using the FSS relation (3) for the 3-state Potts model, we have $`\beta /\nu `$ = 0.131(2) (0.134(2)) for the extension (topological) rule, which is again consistent with the conjectured value, 2/15 (=0.1333) . We have also found nice FSS behavior for the order-parameter distribution functions $`p(m)`$ of the 3-state Potts model as in the case of the Ising model. Here, $`m`$ stands for the absolute value of the vector order parameter. The details will be given in a separate paper.

Here we may remark on the computation time. We only need to modify the code of the SW algorithm slightly in order to change $`p`$. The typical computation time to get $`10^5`$ Monte Carlo samples for $`L=64`$ is 471 seconds using the Alpha 21164A (533 MHz) machine, which is about 30 % longer than that for the SW algorithm, 361 seconds. The most of time is spent on the usual procedure of the SW algorithm, that is, the assignment of the cluster and the cluster spin update. With a small cost of computation time, we can determine the critical point automatically without making simulations for various parameters. In contrast, it takes 958 seconds in the case of the IC algorithm to get $`10^5`$ samples for $`L=64`$ because one should check whether the system is percolating several times before getting one Monte Carlo sample.

To summarize, we have proposed a new cluster algorithm of tuning the critical point automatically. Our algorithm is the extension of the SW algorithm , but we change the probability of connecting spins $`p`$ in the process of Monte Carlo spin update. The resulting distribution of $`p`$ is sharply peaked at $`p_c(L)`$. We approach the canonical ensemble in the limit of small $`\mathrm{\Delta }p`$, which has been explicitly checked by the energy histogram. This is in contrast with the histogram of the IC algorithms , which has broad tails for both high-energy and low-energy sides. Thus, we can use the FSS analysis for the physical quantities. The estimated values of the critical temperatures and the critical exponents for the $`2d`$ Potts models are consistent with the exact values. In order to get more accurate estimate of the critical point and critical exponents, the FSS analysis employed by Ferrenberg and Landau  in a high-resolution Monte Carlo study is useful for the data obtained by the PCC algorithm; we extract $`\nu `$ first, and with $`\nu `$ determined quite accurately we can estimate $`T_c`$. Since the main purpose of the present Letter is to present a new and simple idea of the cluster algorithm, the refined data analysis including the corrections to FSS will be left to a subsequent study.

In the present study, we have shown the application of the PCC algorithm to the thermal phase transition of Potts models, but the idea is based only on the property of a percolation problem. Thus, of course, we can use this algorithm in the study of the geometric percolation problem, that is, the $`q=1`$ Potts model, to get the percolation threshold, $`p_c`$, and various critical properties.

Moreover, we can apply the PCC algorithm to any problem where the mapping to the cluster formalism exists. It is straightforward to apply this algorithm to the diluted Ising (Potts) model. The lack of self-averaging has been recently discussed for the three-dimensional diluted Ising model . The PCC algorithm is quite useful for the problem with the lack of self-averaging, where the distribution of $`T_c`$ due to the randomness is important, because we can tune the critical point of each random sample automatically. Another direction of application is the cluster quantum Monte Carlo simulation , and this problem is left to a future study.

We would like to thank N. Kawashima and M. Kikuchi for valuable discussions. This work was supported by a Grant-in-Aid for Scientific Research from the Ministry of Education, Science, Sports and Culture, Japan.