text:
# The Formation of the Hubble Sequence of Disk Galaxies: The Effects of Early Viscous Evolution

## 1 Introduction

The current picture of disk galaxy formation and evolution has as its basis the dissipative infall of baryons within a dominant dark halo potential well (White & Rees 1978). The collapse and spin-up of the baryons, with angular momentum conservation, can provide an explanation for many of the observed properties of disks, with the standard initial conditions of baryonic mass fraction $`F0.1`$ and dark halo angular momentum parameter $`\lambda 0.07`$ (Fall & Efstathiou 1980; Gunn 1982; Jones & Wyse 1983; Dalcanton, Spergel & Summers 1997; Hernandez & Gilmore 1998; Mo, Mao & White 1998; van den Bosch 1998). Galaxies such as the Milky Way which have an old stellar population in the disk, must, within the context of a hierarchical-clustering scenario, evolve through only quiescent merging/accretion, so as to avoid excessive heating and disruption of the disk (Ostriker 1990). Further, the merging processes with significant substructure cause angular momentum transport to the outer regions, which must somehow be suppressed to allow the formation of extended disks as observed (e.g. Zurek, Quinn & Salmon 1988; Silk & Wyse 1993; Navarro & Steinmetz 1997).

Thus here we adopt the simplified picture that disk galaxies form from smooth gaseous collapse to centrifugal equilibrium, within a steady dark halo potential. We discuss where appropriate below how this may be modified to take account of subsequent infall, or earlier star formation. Our model incorporates the adiabatic response of the dark halo to the disk infall, and we provide new, more general, analytic solutions for the density profile, given a wide range of initial density profiles and angular momentum distributions. We provide new insight into the ‘disk-halo’ conspiracy within the context of this model, demonstrating how an imperfect conspiracy is improved by the disk-halo interaction. We explicitly include subsequent viscous evolution of the gas disk to provide the exponential profile of the stellar disk, and develop analytic expressions that illustrate the process. The resulting radial inflow builds up the central regions of the disk and we investigate the properties of ‘bulges’ that may form as a consequence of instabilities of the central disk. We derive new constraints on the characteristic redshift of disk star formation. We obtain a simple relation connecting the initial conditions, such as spin parameter and baryonic mass fraction, to the efficiency of viscous evolution and star formation.

## 2 The Disk Galaxy Formation Model

In this section we shall derive the mass profiles of disk and halo after the collapse of the baryons. We shall follow earlier treatments of disk galaxy formation (e.g. Mo, Mao & White 1998) by assuming that the virialized dark halo, mixed with baryonic gas, is ‘formed’ – or at least assembled – at redshift $`z_f`$. This virialized halo has a limiting radius $`r_{200}`$ within which the mean density is $`200\rho _{crit}(z_f)`$, and contains a baryonic mass fraction $`F`$. Then

$$r_{200}=\frac{V_{200}}{10H(z_f)};M_{tot}=\frac{V_{200}^2r_{200}}{G}=\frac{V_{200}^3}{10GH(z_f)},$$
(1)
where $`H(z_f)`$ is the value of the Hubble parameter at redshift $`z_f`$, $`M_{tot}(z_f)`$ is the total mass within virialized radius $`r_{200}`$, and $`V_{200}`$ is the circular velocity at $`r_{200}`$.

The baryonic gas cools and settles into a disk, causing the dark halo to contract adiabatically (Blumenthal et al. 1986). The specific angular momentum distribution of the gas is assumed to be conserved during these stages. We shall include below the subsequent re-arrangement of the disk due to angular momentum transport. This we investigate by variation of the disk angular momentum distribution function, choosing an appropriate analytic functional form. Let $`m_d(r)`$ and $`m_h(r)`$ respectively denote the fraction of the total baryonic mass, and total dark mass, that is contained within radius $`r`$, and denote the baryonic mass angular momentum distribution function by

$$m_d\left(<j\right)=f\left(j/j_{max}\right),$$
(2)
where $`j_{max}`$ is the maximum specific angular momentum of the disk.

We will be requiring that the functional form, $`f(j/j_{max})`$, vary as the disk evolves, and it is convenient to introduce the notation $`\mathrm{}j/j_{max}`$ and define

$$c_f1_0^1f\left(\mathrm{}\right)𝑑\mathrm{},$$
(3)
which represents the area above the angular momentum distribution function curve $`f(\mathrm{})`$ for $`0\mathrm{}1`$. We will mimic the effects of viscous evolution by decreasing the value of $`c_f`$ in our evolving disk models in section 4 below.

In terms of this parameter the total disk angular momentum is:

$`J_d`$ $`=`$ $`FJ_{tot}=M_d{\displaystyle _0^{j_{max}}}j{\displaystyle \frac{dm_d}{dj}}𝑑j`$ (4)
$`=`$ $`M_dj_{max}\left(1{\displaystyle _0^1}f\left(\mathrm{}\right)𝑑\mathrm{}\right)=FM_{tot}j_{max}c_f,`$
with $`M_d=FM_{tot}`$. Thus $`j_{max}c_f`$ is the average specific angular momentum of the disk material.

The specific angular momentum of the disk material is assumed to follow that of the dark halo, but in general will not be a simple analytic function (e.g. Quinn & Binney 1992). For illustration, we adopt an analytic monotonic increasing function $`f(b,\mathrm{})`$ containing a free parameter $`b`$, with $`0b1`$. We require the initial angular momentum distribution to be scale free, representing the angular momentum distribution of the virialized halo, and will adopt $`f(b=0,\mathrm{})=\mathrm{}^n`$. We shall vary the value of the parameter $`b`$ to mimic the effects of viscous evolution on the angular momentum distribution.

The total energy is:

$$E_{tot}=\frac{ϵ_0GM_{tot}^2}{2r_{200}}=\frac{ϵ_0M_{tot}V_{200}^2}{2},$$
(5)
where $`ϵ_0`$ is a constant of order unity, depending on the dark halo density profile, and since it is constant for any specific halo model, we can take $`ϵ_0=1`$ without loss of generality. The spin parameter $`\lambda `$ is by definition

$$\lambda J_{tot}|E_{tot}|^{1/2}G^1M_{tot}^{5/2}.$$
(6)
Thus the mean specific angular momentum of the disk material may be expressed as

$$c_fj_{max}=\sqrt{2}\lambda V_{200}r_{200}.$$
(7)
Assuming spherical symmetry, the rotationally-supported disk has mass profile given by

$$m_d=f\left(j/j_{max}\right)=f\left(\frac{\sqrt{GM_{tot}\left(m_d\right)r\left(m_d\right)}}{j_{max}}\right)=f\left(\mathrm{}\right),$$
(8)
while the initial virialized halo mass profile is

$$g(R_{ini})=m_h(R_{ini})=M_{ini}(R_{ini})/M_{tot},$$
(9)
with $`Rr/r_{200}`$.

### 2.1 Constraints on the Final Dark Halo Profile and Mass Angular Momentum Distribution Function

The above equations describe the disk and halo just upon the settling of the gas disk to the mid-plane, prior to the subsequent adiabatic compression of the halo. A self-consistent calculation of the modified disk and halo density profiles may be made by consideration of the adiabatic invariance of the angular action, $`I_\theta v_\theta r𝑑\theta =\sqrt{GM_{tot}(r)r}`$, together with the assumption of no shell crossing (cf. Blumenthal et al. 1986).

Suppose a dark matter particle initially at $`r_{ini}`$ finally settles at $`r(m_d)`$, the radius within which the dark halo mass fraction is $`m_h`$. Then under adiabatic invariance the disk mass profile, $`m_d`$, and the halo mass profile, $`m_h`$, are related by:

$`GM_{tot}\left(m_d\right)r\left(m_d\right)`$ $`=`$ $`GM_{tot}\left(m_h\right)r\left(m_h\right)`$ (10)
$`=`$ $`GM_{ini}\left(m_h\right)r_{ini}\left(m_h\right)`$ (11)
$`=`$ $`GM_{tot}m_hr_{200}g^1\left(m_h\right),`$ (12)
where $`g^1(m_h)`$ is the inverse function of $`g(R)`$, the initial virialized halo mass profile.

Further manipulation of these relations is simplified by introduction of the parameter $`\xi `$, given by

$$\xi \frac{\sqrt{GM_{tot}r_{200}}}{j_{max}}=\frac{c_f}{\sqrt{2}\lambda }.$$
(13)
Generally $`\xi `$ is a quantity that is closely related to the overall disk collapse factor.

From equations (8) and (12), we have

$`m_d`$ $`=`$ $`f(\mathrm{}),`$ (14)
$`\mathrm{}`$ $`=`$ $`\xi \left(m_hg^1(m_h)\right)^{1/2},`$ (15)
where $`\mathrm{}`$ is the normalized specific angular momentum. Again, $`0\mathrm{}1`$, and $`\mathrm{}=1`$ corresponds to the maximum specific angular momentum of the disk, which occurs at the edge of the disk, equivalently at the disk cutoff radius. The fraction of the dark matter contained within the disk thus has a maximum value, $`m_{hc}`$, given by $`\mathrm{}=1`$ in the above equation, and for radii with $`m_hm_{hc}`$, $`m_d=1`$.

To illustrate the physical meaning of these parameters, consider the rigid singular isothermal halo, for which $`g(m_{hc})=m_{hc}=R_c`$. Then from equation (15), $`R_c=R(\mathrm{}=1)=1/\xi `$ which corresponds to the disk cutoff radius. Thus in this case, remembering that $`R`$ is the normalized radius, $`\xi =1/R_c`$ is the disk collapse factor.

Up to now we know the mass profile of disk and halo after collapse, in terms of the normalized specific angular momentum $`\mathrm{}`$, as given in equations (14) and (15), for given forms of the angular momentum distribution function, $`f`$, and initial virialized dark halo mass profile, $`g`$. Next we shall obtain the relation between $`\mathrm{}`$ and radius $`R`$, to complete the derivation of the mass profiles of disk and halo after collapse.

Returning to a general halo density profile, the total mass contained within the radius corresponding to $`m_h`$ is

$$M(m_h)=M_{tot}\left((1F)m_h+Fm_d\right).$$
(16)
From equations (9) - (12) and (16), we have

$`R`$ $`=`$ $`{\displaystyle \frac{GM_{ini}(m_h)r_{ini}(m_h)}{GM(m_h)r_{200}}}`$ (17)
$`=`$ $`{\displaystyle \frac{m_hg^1(m_h)}{(1F)m_h+Fm_d}}.`$ (18)
Introducing the radius variable $`x`$ and coefficient $`c_0`$ as

$`x`$ $``$ $`\xi (1F)R,`$ (19)
$`c_0`$ $``$ $`\xi F/(1F),`$ (20)
we may finally derive the functional dependences on $`\mathrm{}`$ of the disk mass $`m_d`$, of the halo mass $`m_h`$, and of the radius $`x`$:

$`m_d`$ $`=`$ $`f(\mathrm{}),`$ (21)
$`m_hg^1(m_h)`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{\xi ^2}},`$ (22)
$`x`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{\xi m_h(\mathrm{})+c_0f(\mathrm{})}}.`$ (23)
Thus $`\mathrm{}`$ can be thought of as a normalized radius. As we shall see later, $`c_0`$ is a measure of the compactness of the final collapsed disk due to the competition between the spin parameter $`\lambda `$ and the baryonic mass fraction $`F`$.

These equations (21) - (23) can be used to derive disk and halo properties for a free choice of virialized halo profile $`g(R)`$ and angular momentum distribution function $`f(\mathrm{})`$. Within the disk cutoff radius, with $`0\mathrm{}1`$, the disk surface density, circular velocity and the disk-to-dark mass ratio as function of radius $`\mathrm{}`$ or $`R`$ have the generic forms:

$`\mathrm{\Sigma }_d`$ $`=`$ $`{\displaystyle \frac{10H(z)FV_{200}}{2\pi G}}{\displaystyle \frac{1}{R}}{\displaystyle \frac{df}{d\mathrm{}}}{\displaystyle \frac{d\mathrm{}}{dR}},`$ (24)
$`V_c`$ $`=`$ $`{\displaystyle \frac{V_{200}\mathrm{}}{\xi R}},`$ (25)
$`{\displaystyle \frac{M_d(\mathrm{})}{M_h(\mathrm{})}}`$ $`=`$ $`{\displaystyle \frac{c_0f(\mathrm{})}{\xi m_h(\mathrm{})}},`$ (26)
where $`M_d(\mathrm{})=M_d(\mathrm{}=1)m_d(\mathrm{})=M_dm_d(\mathrm{})`$ and $`M_h(\mathrm{})`$ is defined similarly.

The circular velocity at radii beyond the disk cutoff, but within the halo, is given by:

$`V_c`$ $`=`$ $`{\displaystyle \frac{V_{200}\sqrt{m_hg^1(m_h)}}{R}},`$ (27)
$`R`$ $`=`$ $`{\displaystyle \frac{m_hg^1(m_h)}{(1F)m_h+F}}.`$ (28)
Armed with these relations, one may now look at various initial virialized halo density profiles and angular momentum distributions, and determine the allowed parameter space from observed properties of disk galaxies. It is convenient to adopt power-law approximations for the initial virialized halo mass profile and angular momentum distributions, such that $`m_h(\mathrm{})\mathrm{}^m`$ and $`f(\mathrm{})\mathrm{}^n`$ for small $`\mathrm{}`$. Figure 1 shows the location of various fiducial models in the plane of these power law indices $`m`$ and $`n`$; the value $`m=1`$ corresponds to the singular isothermal sphere, the value $`m=4/3`$ corresponds to the Hernquist (1990) and to the Navarro, Frenk & White (1997) profiles, while the value $`m=3/2`$ corresponds to the non-singular isothermal sphere with a constant-density core. These profiles span the range of dark-halo profiles suggested by theory, and plausibly consistent with observations. The shaded region is the allowed parameter space for these models, constrained by the surface-density profile, rotation curve and disk-to-halo central mass ratio. The line ABC is the maximum angular-momentum index consistent with a disk surface density profile in the central regions that declines with increasing radius, while the line DEF is the minimum angular-momentum index consistent with a finite value of the central circular velocity. The line BF denotes the maximum values of $`n`$ consistent with a non-zero central disk-to-halo mass ratio.

This power-law approximation has an initial virialized halo density profile at small radius as $`\rho _{h,ini}(R)R^{\frac{m}{2m}3}`$ (seen by solution of (22) for the form of $`g`$). Solving equations (21) - (23) above gives the corresponding profile after adiabatic infall. In the central region, where the disk dominates the gravitational potential (i.e. $`c_0f(\mathrm{})\xi m_h(\mathrm{})`$), the halo density profile is $`\rho _h(R)R^{\frac{m}{2n}3}`$. Note that in the region where the dark halo dominates the gravitational potential (i.e. $`c_0f(\mathrm{})\xi m_h(\mathrm{})`$), the halo density profile is essentially unaffected by the disk, as expected.

For the case $`n=m`$, the central halo density profile is unaffected since the disk mass density profile and the initial virialized halo density profile have the same dependence on $`\mathrm{}`$. The viable models within the shaded region have $`nm`$, so that the final halo density profile in the central, disk-dominated region should be steeper than its initial virialized profile in this region, not surprisingly. Thus the final halo profile for these models ranges from $`\rho _hR^{0.75}`$ to $`\rho _hR^2`$. It is interesting to note that $`\rho _hR^{0.75}`$, the outcome of an initial virialized halo with a constant density core responding to the settling of a disk with angular momentum index $`n=4/3`$, corresponds to the de-projected de Vaucouleurs central density profile. Thus provided the virialized halo does not have a declining density profile with decreasing radius, which is unphysical, the final dark halo cannot have a constant density core, at least in the very central region where the disk dominates, but should be cuspy.

### 2.2 The Singular Isothermal Sphere

The singular isothermal sphere provides a virialized halo density profile that is the most tractable analytically, and we can obtain some important scaling relations without having to specify the angular momentum distribution $`f(\mathrm{})`$; aspects of the analysis of this profile should hold in general, and provide insight.

The final disk and halo mass profiles are given by solution of:

$`x`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{1+c_0\frac{f(\mathrm{})}{\mathrm{}}}},`$ (29)
$`m_h`$ $`=`$ $`\mathrm{}/\xi ,`$ (30)
$`m_d`$ $`=`$ $`f(\mathrm{}),`$ (31)
where $`0\mathrm{}1`$ and $`c_0`$ and $`\xi `$, the parameters describing the compactness of the collapsed disk and its collapse factor, are defined above in equations (20) and (13).

The collapse factor, defined as the ratio of the pre-collapse radius $`r_{200}`$ to the cutoff disk radius $`r_c`$ at $`\mathrm{}=1`$ (to be distinguished from the final disk scale length) is

$$\frac{r_{200}}{r_c}=\xi (1F)(1+c_0)=\frac{c_f}{\sqrt{2}\lambda }\left(1F+\frac{c_f}{\sqrt{2}}\frac{F}{\lambda }\right),$$
(32)
where $`c_f`$, defined in equation (3), is a measure of the shape of the angular momentum distribution, and small values of $`c_f`$ mean steeply-rising angular momentum distributions.

Typical angular momentum distributions for disks are shown in Figure 2 (a,b) and discussed below. Here one should just note that a value $`c_f0.5`$ is reasonable. With this, and $`\lambda =0.06`$ and $`F=0.1`$, we obtain $`\xi =5.9`$ and $`c_0=0.66`$. The collapse factor as defined here is then about a factor of 9. For fixed angular momentum distribution $`c_f`$, the collapse factor depends not only on $`1/\lambda `$ but also on compactness $`c_0F/\lambda `$. This is consistent with the results of previous studies of the collapse factor in two extreme cases (Jones & Wyse 1983; Peebles 1993): if the final disk is so self-gravitating that the value of $`F/\lambda `$ corresponds to $`c_01`$, then the collapse factor is $`F/\lambda ^2`$; on the other hand if the final disk is sufficiently far from self-gravitating, with small $`F/\lambda `$ and $`c_01`$, the collapse factor is $`1/\lambda `$. The above collapse factor relationship is valid over the entire range of values for the parameters $`\lambda `$ and $`F`$ (and is the most general relation derived to date).

It should be noted that varying the normalized angular momentum distribution by varying $`c_f`$ changes the derived collapse factor; this is investigated further below.

Comparison between the sizes of observed disks and those predicted from such collapse calculations provides a constraint on the redshift at which the collapse happened (e.g. Mo, Mao & White 1998). Our disk cutoff radius for non-trivial $`f(\mathrm{})`$ may be expressed in terms of the initial conditions as :

$$r_c=\frac{\sqrt{2}\lambda V_{200}}{10H(z_f)c_f(1F+\frac{c_f}{\sqrt{2}}\frac{F}{\lambda })},$$
(33)
where $`z_f`$ is the ‘formation’ or assembly redshift, at which the halo is identified to have given mass and circular velocity, with no mass infall after this epoch. One can see from this relation that both $`\lambda `$ and $`F/\lambda `$ are equally important; previous determinations considered the baryonic mass fraction fixed (Dalcanton, Spergel & Summers 1997; Mo, Mao & White 1998).

The explicit inclusion of the parameter $`c_f`$ allows us to take account of disk evolution, as gas is transformed to stars. Adopting a disk cutoff radius at three disk scale-lengths (e.g. van der Kruit 1987), and choosing specific values of the present-day stellar disk scale-length $`r_d=3.5`$kpc, $`\lambda =0.06`$, $`V_{200}=200`$ kms<sup>-1</sup> and an Einstein-de-Sitter Universe with Hubble constant $`0.5<h<1`$, the above relation gives the formation redshift $`1.6<z_f<3.1`$ for $`F=0.05`$ and $`c_f=1/3`$; $`1.4<z_f<2.8`$ for $`F=0.1`$ and $`c_f=1/3`$; $`0.9<z_f<2.0`$ for $`F=0.05`$ and $`c_f=1/2`$; $`0.7<z_f<1.7`$ for $`F=0.1`$ and $`c_f=1/2`$.

Consistent with previous calculations, for fixed formation redshift smaller values of $`\lambda `$ can lead to smaller disk size; we have here explicitly demonstrated that larger $`F/\lambda `$ can also lead to this result. Lower values of $`c_f`$ also lead to higher redshift of formation. Thus for no viscous evolution i.e. there is no angular momentum redistribution during the evolution of the galactic disk and $`c_f`$ has a time-independent value, the formation redshift $`z_f`$ determined by assuming fixed initial $`\lambda `$ and $`F`$ and fixed disk size $`r_c`$ is smaller than would be determined if $`c_f`$ could be decreased (to mimic say viscous evolution). Lower values of $`c_f`$ for given total angular momentum content imply a larger disk scale-length; the effect of viscous evolution is to re-arrange the disk material so as to increase the disk scale-length. Typical values of viscosity parameters lead to a factor of 1.5 increase in disk scale-length in a Hubble time. As we demonstrate below, a smaller value for $`c_f`$ corresponds to larger bulge-to-disk ratio. This result is consistent with the results of Mo, Mao & White (1998): halo and disk formation redshift can be pushed to higher value when a bulge is included. These trends are general, and not tied to the specific model of the halo density profile.

One should bear in mind that the old stars in the local thin disk of the Milky Way have ages of at least 10 Gyr, and may be as old as the oldest stars in the Galaxy (Edvardsson et al. 1993); the age distribution of stars at other locations of the Galactic disk is very poorly-determined, but it is clear that a non-trivial component of the thin disk was in place at early times (at redshift $`z>2`$ for the cosmologies considered above). A common assumption in previous work is that the mass angular momentum distribution of the disk is that of a solid-body, rotating uniform density sphere, $`f(\mathrm{})=1(1\mathrm{})^{3/2}`$ (Mestel 1963). For this distribution, $`c_f=0.4`$, and one derives a low redshift of formation for a galaxy like the Milky Way (Mo, Mao & White 1998), which has difficulties with the observations.

#### 2.2.1 The Singular Isothermal Halo with Simple $`f(\mathrm{})`$

Analytic solutions to equations (29) - (31) can be obtained by assuming a simple monotonic increasing function $`f(b,\mathrm{})`$ containing one parameter $`b`$ with $`0b1`$. In order to avoid the situation where one obtains a trivial collapse factor due to a very small amount of disk material at very large radius, we restrict the shape of the mass angular momentum distribution function $`f(\mathrm{})`$ to avoid too shallow an asymptotic slope as $`f`$ approaching 1 when $`\mathrm{}`$ increases (see Figure 2). From Figure 1, the singular isothermal halo at point F requires $`f(\mathrm{})\mathrm{}`$ for $`\mathrm{}1`$. A simple form consistent with this is $`f(\mathrm{})=(1+b)\mathrm{}b\mathrm{}^2`$ with $`0\mathrm{}1`$ and $`0b1`$. The angular momentum parameter $`c_f=(3b)/6`$. Figure 2a shows this angular momentum distribution function with different values of the parameter $`b`$, compared with the mass angular momentum distribution of a solid-body, rotating uniform density sphere, $`f(\mathrm{})=1(1\mathrm{})^{3/2}`$ (Mestel 1963). The normalized total angular momentum $`c_f`$ corresponds to the area above each curve.

Within the disk, where $`0\mathrm{}1`$, the galactic disk surface density, circular velocity and the disk-to-dark mass ratio as a function of radius are:

$`\mathrm{\Sigma }_d(\mathrm{})={\displaystyle \frac{10V_{200}H(z)F(1F)^2\xi ^2}{\pi G}}`$
$`\times {\displaystyle \frac{(1+b2b\mathrm{})(1+c_0+c_0bc_0b\mathrm{})^3}{2\mathrm{}(1+c_0+c_0b)}},`$ (34)
$$V_c(\mathrm{})=V_{200}(1F)[1+c_0(1+bb\mathrm{})],$$
(35)
$$\frac{M_d(\mathrm{})}{M_h(\mathrm{})}=\frac{(1+bb\mathrm{})(3b)F}{6\sqrt{2}\lambda (1F)},$$
(36)
where

$$\mathrm{}=\frac{[1+c_0(1+b)]\xi (1F)R}{1+bc_0\xi (1F)R}.$$
(37)
From equation (27)-(28), the circular velocity at radii beyond the disk cutoff, but within the halo, is:

$$V_c(R)=\frac{V_{200}(1F)}{2}\left[1+\sqrt{1+\frac{4F}{(1F)^2R}}\right].$$
(38)
These results are plotted in Figure 3 (a,b,c) for $`\lambda =0.06`$, $`F=0.1`$, $`0\mathrm{}1`$. The different curves correspond to $`b=0,1/4,1/2,3/4,1`$. Larger values of $`b`$ yield larger disk cut-off radii; note that the circular velocities for points beyond the cut-off radius of a given model may be obtained by forming the envelope of the values for the cut-off radius for larger values of $`b`$. Varying the value of the parameter $`b`$ changes the angular momentum distribution function similar to the effects of viscous evolution; the ratio of disk mass to dark halo mass increases at small radius with increasing $`b`$, which can be interpreted as due to radial inflow of disk material.

#### 2.2.2  Halo Density Profile and Angular Momentum Combinations

For initial virialized halo profiles other than the singular isothermal sphere, the collapse factor has a more general form:

$$\frac{r_{200}}{r_c}=\xi (1F)(\xi m_{hc}+c_0),$$
(39)
where $`m_{hc}`$ is the solution of $`1=\xi ^2m_{hc}g^1(m_{hc})`$ and is the mass fraction of the dark halo that is contained within the cut-off radius of the disk. The trend of the dependence of the collapse factor on the values of $`\lambda `$ and $`F`$ remains the same as found above for the singular isothermal halo.

The range of viable models represented by the shaded region in Figure 1 can be investigated by the appropriate virialized halo profile $`g(R)`$ and angular momentum distribution function $`f(b,\mathrm{})`$ corresponding to the points E,B,D and C. The distribution of surface density, circular velocity and disk-to-halo mass ratio for these models, varying parameter $`b`$, are shown in Figures 4–7, by choosing fixed $`\lambda =0.06`$, $`F=0.1`$ (and halo core size $`c=4`$ if the halo has a core radius).

Model E: The results of a model corresponding to point E in Figure 1 are shown in Figure 4; this has a Hernquist halo profile $`g(R)=\frac{(1+c)^2R^2}{(1+cR)^2}`$ and $`f(b,\mathrm{})=(1+b)\mathrm{}b\mathrm{}^2`$. The different curves correspond to $`b=0,1/4,1/2,3/4,1`$.

Model B: The results of a model corresponding to point B in Figure 1 are shown in Figure 5; this has a Hernquist halo profile $`g(R)=\frac{(1+c)^2R^2}{(1+cR)^2}`$ and $`f(b,\mathrm{})=(1+10b)\mathrm{}^{4/3}10b\mathrm{}^{22/15}`$. This $`f(b,\mathrm{})`$ is shown in Figure 2b. For $`\mathrm{}1`$, $`f(b=0,\mathrm{})\mathrm{}^{4/3}`$; the $`22/15`$ index in the second term is determined by the requirement that $`f(b,\mathrm{})`$ should be a monotonic increasing function of $`\mathrm{}`$ for all values of $`0b1`$. The different curves correspond to $`b=0,0.1,0.3,0.6,1`$.

Model D: The results of a model corresponding to point D in Figure 1 are shown in Figure 6; this has a non-singular isothermal halo with a constant density core, $`g(R)=\frac{cRarctan(cR)}{carctan(c)}`$, and $`f(b,\mathrm{})=(1+b)\mathrm{}b\mathrm{}^2`$. Again the different curves correspond to $`b=0,1/4,1/2,3/4,1`$.

Model C: The results of a model corresponding to point C in Figure 1 are shown in Figure 7; this has a non-singular isothermal halo with constant density core, $`g(R)=\frac{cRarctan(cR)}{carctan(c)}`$, and $`f(b,\mathrm{})=(1+10b)\mathrm{}^{4/3}10b\mathrm{}^{22/15}`$. Again the different curves correspond to $`b=0,0.1,0.3,0.6,1`$.

As can be seen from the figures, the different halo profiles and angular momentum distributions produce disks with a variety of surface density profiles and rotation curves. As in Figure 3b, the circular velocities for points beyond the cut-off radius of a given model may be obtained by forming the envelope of the values for the cut-off radius for larger values of $`b`$. Thus if the disk is very compact, from equations (27)-(28) or equation (38), we find that the circular velocity beyond the edge of the disk tends to decreases with radius. The location of the edge of the disk depends on both $`\lambda `$ and $`F`$ (in addition to $`b`$). Thus we have shown that rotation curves should show an imperfect disk–halo ‘conspiracy’ if the disk is too compact or too massive. This is consistent with observations (Casertano & van Gorkom 1991).

Further, these figures demonstrate that with increasing $`b`$, the inner rotation curves become flat, and the transition between disk-dominated and halo-dominated sections of the rotation curve becomes more and more smooth with increasing $`b`$. As discussed above, and illustrated in Figure 2, a higher value of $`b`$ corresponds to a flatter specific angular momentum distribution function, and increasing $`b`$ mimics the effects of viscous evolution in transporting angular momentum. This indicates that viscous evolution can help the creation of an apparent disk-halo ‘conspiracy’.

As can be seen from equations (21) - (26), for given $`c_0`$ (the compactness parameter defined in equation (20)), or $`F/\lambda `$ ratio, within a given model of $`f(\mathrm{},b)`$ and virialized dark halo profile $`g(R)`$, the normalized properties of the disks formed for $`b=0`$ are very similar. In particular the disk surface density profile, rotation curve and disk-to-halo mass ratio profile scale similarly with $`\mathrm{}`$. For example, in the case of the singular isothermal sphere (model F), equations (34) - (37) show that for $`b=0`$ the normalized disks are identical for given $`F/\lambda `$ ratio. Thus $`F/\lambda `$ must be an important factor in distinguishing one disk from another. The overall normalization of the surface density is $`F/\lambda ^2`$, so that again $`F/\lambda `$ and $`\lambda `$ enter separately and are both important.

Note that here we are not insisting that the surface density profile of the gas disk so formed be exponential, unlike previous work (Mo, Mao & White 1998). We shall however appeal to viscous evolution tied to star formation to provide a stellar exponential disk. We now turn to this.

## 3 The Viscous Evolution and Star Formation

In this paper we aim to link viscous evolution within disks and the Hubble sequence of disk galaxies. One of the motivations for invoking viscous disks is that if the timescale of angular momentum transport via viscosity is similar to that of star formation, a stellar exponential disk is naturally produced independent of the initial gaseous disk surface density profile (Silk & Norman 1981; Lin & Pringle 1987; Saio & Yoshii 1990; Firmani, Hernandez & Gallagher 1996). Angular momentum transport and associated radial gas flows (both inwards and outwards) can also, as shown above, provide a tight ‘conspiracy’ between disk and halo rotation curves, and, as demonstrated below, provide a higher phase space density in bulges as compared to disks.

The star formation rate per unit area in a disk can be represented by a modified Schmidt law involving the dynamical time and the gas density (Wyse 1986; Wyse & Silk 1989). We shall use the form of the global star formation rate per unit area, $`\mathrm{\Sigma }_\psi `$, of Kennicutt (1998), based on his observations of the inner regions of nearby large disk galaxies:

$$\mathrm{\Sigma }_\psi =\alpha \mathrm{\Sigma }_{gas}\mathrm{\Omega }_{gas},$$
(40)
where $`\mathrm{\Sigma }_{gas}`$ is total gas surface density, $`\mathrm{\Omega }_{gas}`$ is the dynamical time at the edge of the gas disk, and the normalization constant, related to the efficiency of star formation, has the value $`\alpha =0.017`$ (Kennicutt 1998). Note that observations of the star-formation rates in the outer regions of disk galaxies suggest that it is actually volume density that should enter the Schmidt law, rather than surface density, and since many (if not all) gas disks flare in their outer regions, equation (40) will over estimate the star formation rate (Ferguson et al. 1998). This is beyond the scope of the present model, but should be borne in mind and will be incorporated in our future work.

For our models here, the edge of the initial gas disk is where $`\mathrm{}=1`$, $`R=R_c`$ and thus $`\mathrm{\Omega }_{gas}=\mathrm{\Omega }_c`$. For halo formation redshift $`z_f`$, we can obtain the relationship between global star formation timescale and the galaxy initial conditions in the general form using equation (1), (25) and (39):

$`t_{}^1`$ $`=`$ $`\alpha \mathrm{\Omega }_c=\alpha V_c(r_c)/r_c`$ (41)
$`=`$ $`10\alpha H(z_f)\xi (1F)^2(\xi m_{hc}+c_0)^2,`$
where again $`m_{hc}`$ is the solution of $`1=\xi ^2m_{hc}g^1(m_{hc})`$, and is the fraction of the dark halo mass that is contained within the cut-off radius of the disk. The parameters $`\xi `$ and $`c_0`$ are defined in equations (13) and (20); $`F`$ is the baryonic mass fraction in the initial density perturbation.

In the case of the singular isothermal halo, this relation has a simple form:

$$t_{}^1=10\alpha H(z_f)\xi (1F)^2(1+c_0)^2.$$
(42)
The gas consumption timescale will be longer than the characteristic star formation timescale due to the gas returned by stars during their evolution and death. For a standard stellar Initial Mass Function, $`t_g2.5t_{}`$ (e.g. Kennicutt et al. 1994).

This modified Schmidt law is based on observations of the inner regions of nearby large disk galaxies, and simple theoretical principles. Assuming it holds at all epochs allows one to estimate the properties of present-day disks from the initial conditions of earlier sections in this paper. Let us assume that the dark halo is fully virialized at redshift $`z_f`$. In keeping with the spirit of hierarchical clustering, let us allow for some star formation that could have taken place in the disk, from an earlier redshift $`z_i`$, and that the total mass of the system could increase until $`z_f`$ (although to maintain the thin disk, this accretion and merging must be only of low mass, low density systems). So for any time $`t`$ or redshift $`z`$ between $`z_i`$ and $`z_f`$,

$`{\displaystyle \frac{dM_g}{dt}}`$ $`=`$ $`F{\displaystyle \frac{dM_{tot}}{dt}}{\displaystyle \frac{dM_{}}{dt}},`$ (43)
$`{\displaystyle \frac{dM_{}}{dt}}`$ $`=`$ $`{\displaystyle \frac{M_g}{t_g(z)}},`$ (44)
where $`M_{tot}`$ is defined in equation (1) and $`M_{}`$ is mass locked up into stars. In an Einstein-de-Sitter Universe, $`H(z)=H_0(1+z)^{3/2}`$ and $`H(z)t=2/3`$. We have

$$\frac{dM_g}{dt}=A\frac{BM_g}{t},$$
(45)
where $`A=\frac{3FV_{200}^3}{20G}`$ and $`B=\frac{8}{3}\alpha \xi (1F)^2(1+c_0)^2`$, and $`\xi `$ and $`c_0`$ are expressed in terms of $`c_f`$, $`\lambda `$ and $`F`$ through equations (13) and (20) in section 2.1. Hence, from equation (1), $`M_{tot}\frac{A}{FH(z)}`$.

Identifying the dark halo to have fixed $`V_{200}`$, independent of redshift, leads to $`A`$ also being a constant, and thus the total mass grows as $`M_{tot}t`$. This differs from the standard solution of infall onto a point, $`M_{tot}t^{2/3}`$ (Gunn & Gott 1972).

The solution to the above equation is then

$$M_g=\frac{A}{1+B}t\left[1+B\left(\frac{t}{t_i}\right)^{(1+B)}\right],$$
(46)
where $`t_i`$ corresponds to the redshift $`z_i`$ of the onset of star formation. Thus at the halo formation redshift $`z_f`$, the disk gas fraction is

$$f_g(z_f)=\frac{1+B\left(\frac{1+z_f}{1+z_i}\right)^{3(1+B)/2}}{1+B},$$
(47)
Since $`B`$ depends on $`c_f`$, $`\lambda `$ and $`F`$, i.e. $`B\frac{c_f}{\lambda }(1F+\frac{c_f}{\sqrt{2}}\frac{F}{\lambda })^2`$, the value of the constant $`B`$ may be evaluated for $`\lambda =0.06`$ and various reasonable values of $`F`$ and $`c_f`$ as: for $`F=0.1`$, $`c_f=1/3`$, $`B=0.30`$; for $`F=0.05`$, $`c_f=1/3`$, $`B=0.23`$; for $`F=0.1`$, $`c_f=0.5`$, $`B=0.59`$; and for $`F=0.05`$, $`c_f=0.5`$, $`B=0.41`$. For fixed $`\lambda `$ and $`F`$, small values of $`B`$ correspond to small values of $`c_f`$, and hence flatter specific angular momentum distributions.

Thus for $`1+z_f\mathrm{}<2(1+z_i)`$, $`f_g(z_f)1/(1+B)`$, typically $`\mathrm{}>2/3`$.

We assume that after $`z_f`$, there is no further infall, and that the gas in the disk will be consumed with characteristic timescale $`t_g(z_f)`$. Thus the gas fraction of a typical disk galaxy at the present time is

$$f_g=f_g(z_f)\mathrm{exp}\left[\delta t(z_f)/t_g(z_f)\right],$$
(48)
where $`\delta t(z_f)`$ is the time interval between the halo formation redshift $`z_f`$ and present time $`z=0`$.

Assuming an Einstein-de-Sitter universe, we have $`\delta t(z_f)=t_0(1H_0/H(z_f))`$ and $`t_0H_0=2/3`$. Thus a typical value for the present gas fraction of disk galaxies is:

$$\mathrm{ln}\left(\frac{f_g(z_f)}{f_g}\right)=B\left[(1+z_f)^{3/2}1\right].$$
(49)
With the approximation $`f_g(z_f)=1/(1+B)`$, we have

$$(1+z_f)^{3/2}=1+\frac{\mathrm{ln}f_g\mathrm{ln}(1+B)}{B}.$$
(50)
For the Milky Way Galaxy, if we adopt $`5\times 10^9M_{}`$ for the atomic $`HI`$ gas, and $`1.3\times 10^9M_{}`$ for the molecular $`H_2`$ gas (Blitz 1996; Dame 1993), then with an estimate of $`46\times 10^{10}M_{}`$ for the total baryonic mass of the Milky Way depending on stellar exponential scale length (Dehnen & Binney 1998), we obtain gas fraction, the total gas mass includes $`24\%`$ helium mass, $`f_g15\%`$ or even higher depending on the mass model of the Galaxy. As mentioned earlier, for fixed $`\lambda `$ and $`F`$, small values of $`B`$ correspond to small values of $`c_f`$, and hence flatter specific angular momentum distributions. For small values of the parameter $`B`$, say $`B0.3`$, we obtain $`z_f2.4`$, while for large values of $`B`$, $`B0.6`$, we obtain $`z_f1.2`$. The larger value of $`z_f`$ is preferred, given what we know of the age distribution of stars in the local thin disk (e.g. Edvardsson et al.1993). The effect of viscous evolution is equivalent to choosing small $`c_f`$, i.e. small values of $`B`$. So the inclusion of viscous evolution can give relatively higher halo formation redshift, which is consistent with the results from the constraint on the redshift of formation that we obtained from considerations of the size of the disk.

It should be noted that for fixed halo ‘formation’ redshift $`z_f`$, the star formation timescale $`t_{}^1\frac{1}{\lambda }(1F+\frac{1}{2\sqrt{2}}\frac{F}{\lambda })^2`$. Again, both $`\lambda `$ and $`F/\lambda `$ are important. As we discussed in the previous section, the structure of the normalized disk depends strongly on $`F/\lambda `$ while the overall normalization depends strongly on $`\lambda `$ for fixed $`F/\lambda `$. As we show later, the bulge-to-disk ratio also depends on the importance of $`\lambda `$ and $`F/\lambda `$. Thus many aspects of the Hubble sequence of disk galaxies – star formation timescale, disk gas fraction, bulge-to-disk ratio – depend on both $`\lambda `$ and $`F/\lambda `$.

The star formation timescale derived above is independent of $`V_{200}`$, which at first sight is surprising, given that the Hubble sequence of disk galaxies has been interpreted as a sequence of star formation timescales, relative to collapse times (Sandage 1986), and observations show that the Hubble type of a disk galaxies is broadly correlated with the disk luminosity (Lake & Carlberg 1988; de Jong 1995). However, $`V_{200}`$ is not an easily-observed quantity.

The present-day luminosity of a disk in our model can be written as

$$L=\frac{FM_{tot}(1f_g)}{\gamma _{}}=\frac{F[1f_g(B,z_f)]V_{200}^3}{10GH(z_f)\gamma _{}},$$
(51)
with $`\gamma _{}`$ the current value of the mass-to-light ratio. Estimation of the predicted Tully-Fisher relation depends on what model parameter we use for the width of the HI line. Obviously if we identify this with $`V_{200}`$ we have a reasonable relationship provided that the coefficient $`\gamma _{tf}`$ is constant, i.e. $`F[1f_g(B,z_f)]/H(z_f)=\gamma _{TF}`$. This then requires that there should be a correlation between $`F`$ and $`z_f`$, in the sense that for large $`F`$, $`z_f`$ is large. Then we can see from equation (42) that this leads to the star formation timescale being small, implying more efficient viscous evolution, and leading to larger B/T ratio. Large $`z_f`$ may be correlated with small $`V_{200}`$ in the context of hierarchical-clustering cosmologies, and in that case, a short star formation time and large B/T ratio should be correlated with low luminosity, which is not consistent with observation. However, an interpretation that is compatible with observations is that high n-sigma fluctuations for fixed $`V_{200}`$ can form high luminosity disks with large B/T ratio.

The Tully-Fisher relationship is not a simple relation between luminosity and $`V_{200}`$, but depends on where the circular velocity $`V_c(R)`$ is measured (Courteau 1997). The relationship between $`V_{200}`$ and $`V_c(R)`$ obviously depends on the details of the halo density profile and the angular momentum distribution. From the rotation curves in fig. 3 to fig. 7, we can see that it is appropriate to choose for our estimate of $`V_c`$ the circular velocity at the cutoff radius of disk, adopted as 3 three scale lengths. From equation (35), we have $`V_c=V_{200}(1F)(1+c_0)`$, where $`c_0=c_fF/\sqrt{2}\lambda (1F)`$, as defined in equation (20), is the compactness of disk. Now the predicted Tully-Fisher relation is

$$L=\frac{F[1f_g(B,z_f)]V_c^3}{10GH(z_f)\gamma _{}(1+c_0)^3(1F)^3}.$$
(52)
Requiring that the coefficient $`\gamma _{TF}`$ in the Tully-Fisher relation be constant, i.e. $`\gamma _{TF}F[1f_g(B,z_f)]/H(z_f)(1+c_0)^3(1F)^3`$, then gives, from equation (42), $`t_{}^1c_0/(1+c_0)`$. So for small $`c_0`$, the star formation timescale is large, which can cause less efficient viscous evolution. So less efficient viscous evolution and small $`F/\lambda `$ will lead to small B/T ratio. Also small $`c_0`$ is correlated with large $`z_f`$ from the constancy of the coefficient in the Tully-Fisher relation. Similarly large $`z_f`$ may be correlated with small $`V_{200}`$ in the context of hierarchical clustering cosmology. So small $`c_0`$ and small $`V_{200}`$ will lead to small $`V_c`$ and lower disk luminosity. Thus this version of the predicted Tully-Fisher relation appears fully compatible with the observations. However, one should bear in mind that we have adopted a fixed constant of proportionality $`\alpha `$ in the star formation law, and this may well vary with global potential well depth (White & Frenk 1991) or local potential well depth (Silk & Wyse 1993).

A further test of the model is the relation between disk scale and circular velocity, and its variation with redshift; observations indicate that $`R_d/V_c`$ is smaller at high redshift, $`z1`$ (Vogt et al. 1996; Simard et al. 1999). In our model there is little change in total mass between redshifts of unity and the present, and so this evolution of disk size cannot be due to the halo mass growth, as had been proposed by Mao, Mo & White (1998). Instead in our model this is due to the different and changing scale lengths of gas and stars. In our model, due to early star formation, the gas fraction at $`z_f`$ is about $`f_g(z_f)1/(1+B)`$, typically $`2/3`$. It is natural that the gas component of the disk will have a larger scale length than the stellar component. The scale length of the gas component can increase with time due to viscous evolution while the scale length of stellar component can also increase with time due to the non-linear local star formation law (Saio & Yoshi 1990). So the stellar scale length at high redshift, when the gas fraction is large, should be much smaller than the stellar scale length at present time, when gas fraction is lower. The study of detailed evolution of gas and stellar component is beyond the scope of this paper; but the prediction of stellar size evolution of our model is qualitatively consistent with the observations.

The measured distribution of $`R_d/V_c`$ for the local disk galaxy sample is approximately peaked at $`R_d/V_c1.5`$ ($`R_d`$ in kpc, $`V_c`$ in kms<sup>-1</sup>) with a spread from -2 to -1 (Mao, Mo & White 1998; Courteau 1996). Our estimation of the disk size or scale length in Section 2 is valid for galaxies at the present time, when the gas fraction is small. Then assuming the disk cut off radius is three scale lengths, from equation (41), $`R_d/V_c=R_c/3V_c=\alpha t_{}/3`$, and the present $`R_d/V_c`$ is an indication of the galactic global star formation timescale; further, from the constancy of $`\gamma _{TF}`$, the coefficient of the Tully-Fisher relation in equation (52), we have

$$R_d/V_c=\alpha t_{}/3=\frac{G}{3}\gamma _{TF}\gamma _{}\frac{1+c_0}{c_0(1f_g)}.$$
(53)
Thus this ratio is also an indication of the disk compactness. Adopting the B-band mass-to-light ratio of our local disk $`\gamma _{}2.5M_{}/L_{}`$, and using the luminosity of our galaxy $`L_B3\times 10^{10}L_{}`$ and $`V_c220`$ kms<sup>-1</sup> to estimate $`\gamma _{TF}=L_B/V_c^3`$, we can obtain that $`\mathrm{log}(R_d/V_c)2.0+\mathrm{log}(\frac{1+c_0}{c_0(1f_g)})`$. Obviously $`\mathrm{log}(R_d/V_c)2`$ roughly corresponds to the predicted lower limit of the local sample, which is consistent with observations. The distribution of $`F`$ from 0.05 to 0.2 and the distribution of $`\lambda `$ from 0.03 to 0.12 will cause the value of the compactness parameter $`c_0`$ to spread from 0.1 to 2 approximately. Adopting typical values $`F=0.1`$, $`\lambda =0.06`$, $`c_f=1/3`$, the peak will be located at $`\mathrm{log}(R_d/V_c)=1.5`$ with $`c_0=0.44`$, which is again consistent with observations. The spread of $`R_d/V_c`$ is simply caused by the spread in the value of compactness parameter $`c_0`$.

## 4 The Formation of Bulges and the Hubble Sequence

Galaxies classification schemes based on morphology are the basic first step in understanding how galaxies form and evolve (van den Bergh 1998). The bulge-to-disk luminosity ratio is one of the three basic classification criterion for the Hubble sequence (Sandage 1961). However, the relation between bulge-to-disk ratio and Hubble type has a fair amount of scatter, some of which must be related to the difficulty of a decomposition of the light profile into bulge and disk, and the bulge-to-disk ratio is dependent on the band-pass used to define the luminosity (de Jong 1996). The current observational data show that bulges can be diverse and heterogeneous (Wyse, Gilmore & Franx 1997). Some share properties of disks and some are more similar to ellipticals. Models of bulge formation can be classified into several categories: the bulge is formed from early collapse of low angular momentum gas, with short cooling time and efficient star formation (Eggen et al. 1962; Larson 1976 – who invoked viscosity to transport angular momentum away from the proto-bulge; van den Bosch 1998); the bulge is formed from merging of disk galaxies (Toomre & Toomre 1972; Kauffmann, White & Guiderdoni 1993); the bulge is formed from the disk by secular evolution process after bar instability (Combes et al. 1990; Norman, Sellwood & Hasan 1996; Sellwood & Moore 1999; Avila-Reese & Firmani 1999); the bulge is formed from early dynamical evolution of massive clumps formed in disk (Noguchi 1999). However it is speculated that large bulges, which tend to have de Vaucouleur’s law surface brightness profiles, share formation mechanisms with ellipticals, while smaller bulges, which tend to be better fit by an exponential profile, are formed from their disks through bar dissolution. It should be noted that the significantly higher phase space density of bulges when compared to disks suggests that gaseous inflow should play a part in the instability (Wyse 1998), just as invoked in earlier sections of this paper, for other reasons.

Here we only consider this latter case, the formation of small bulges from the disk. Early studies of disk instabilities showed that Toomre’s local stability criterion $`Q>1`$ is also sufficient for global stability to axisymmetric modes (Hohl 1971; Binney & Tremaine 1987). It is known that the bar instability requires a similar condition (Hockney & Hohl 1969; Ostriker & Peebles 1973). Efstathiou, Lake & Negroponte (1982) used N-body techniques to study the global bar instability of a pure exponential disk embedded in a dark halo and proposed a simple instability criterion for the stellar disk, based on the disk-to-halo mass ratio. However, it has been argued recently that there is no such simple criterion for bar instability (Christodoulou, Shlosman & Tohline 1995; Sellwood & Moore 1999). Further, recent N-body simulations (Sellwood & Moore 1999; Norman, Sellwood & Hasan 1996) show that every massive disk form a bar during the early stages of evolution, but later the bar is destroyed by the formation of a dense central object, once the mass of that central concentration reaches several percent of the total disk mass. This can be understood in terms of the linear mode analysis work and nonlinear processes of swing amplifier and feed back loops (Goldreich & Lynden-Bell 1965; Julian & Toomre 1966). Toomre (1981) argued that the bar-instability can be inhibited in two ways: one is that a large dark halo mass fraction can reduce the gain of the swing amplifier, while the other is that feedback through the centre can be shut off by a inner Lindblad resonance (ILR). The dense central object can destroy the bar via the second of these mechanisms. However, most of these studies assume that the dark halo has a constant density core. On the contrary, we have shown that dark halo profile cannot have a constant density core after adiabatic infall, if one starts from physical initial conditions. It will be interesting to study bar formation under different disk-halo profiles in addition to the well-studied harmonic core.

We will for simplicity adopt the simple criterion of Efstathiou, Lake & Negroponte (1982), interpreted to determine the size of a bar-unstable region, with the radial extent of the bar, $`r_b`$ defined by

$$M_d(r_b)/M_h(r_b)\beta ,$$
(54)
with the value of the parameter $`\beta `$ chosen to fit observations. As we have argued in section 2, $`F/\lambda `$ is the important quantity determining the overall normalization of the disk surface density. From N-body simulations (Warren et al. 1992; Barnes & Efstathiou 1987; Cole & Lacey 1996; Steinmetz & Bartelmann 1995), the distribution of $`\lambda `$ can be fit by a log-normal distribution:

$$P(\lambda )d\lambda =\frac{1}{\sqrt{2\pi }\sigma _\lambda }exp\left[\frac{\mathrm{ln}^2(\lambda /\lambda _0)}{2\sigma _\lambda }\frac{d\lambda }{\lambda }\right],$$
(55)
where $`\lambda _0=0.06`$ and $`\sigma _\lambda =0.5`$ (this result is fairly independent of the slope of the power spectrum of density fluctuations).

What of the possible range of the baryonic fraction $`F`$? Some previous studies suggested that $`F\lambda `$ as an explanation of the disk-halo ‘conspiracy’ (Fall & Efstathiou 1980; Jones & Wyse 1983; Ryden & Gunn 1987; Hernandez & Gilmore 1998). Some interpretations of the observed Tully-Fisher relationship suggest that indeed $`F`$ is not invariant (McGaugh & de Blok 1998). Here we shall assume that the distribution of $`F`$, similar to $`\lambda `$, is log-normal, centreed at $`F_0=0.1`$ with $`\sigma _F=0.05`$. So $`F`$ is mainly within the range $`0.050.2`$.

We generated a Monte Carlo sample of disks, with fixed halo $`V_{200}`$ and halo formation redshift $`z_f`$, but with values of $`F`$ and $`\lambda `$ following the above distributions. Then for given virialized dark halo profiles $`g(R)`$ and angular momentum distribution $`f(b,\mathrm{})`$ in Figure 1, we can calculate the star formation timescale from equation (41). The parameter $`b`$ represents the efficiency of viscous evolution; on the assumption that the viscous timescale is equal to the star formation timescale, we can use a simple linear correlation between the value of $`b`$ and the star formation timescale. Then the value of parameter $`b`$ can be obtained. From equations (21), (22), (23), (26), (50) and (54), one can calculate the bulge-to-disk ratio and final disk gas fraction.

Thus we can plot B/T ratio versus disk gas fraction, for the whole Monte Carlo sample, confining the parameter space of $`\lambda `$ and $`F`$ to give small bulges with $`B/T<0.5`$. Larger values of B/T would correspond to such an unstable disk that the exercise is invalid. Disks that are too stable, having low $`F/\lambda `$ or large $`\lambda `$, will evolve little and probably end up as low surface brightness systems.

Figure 8a shows model F, representing the singular isothermal halo. Assuming the typical values $`\lambda =0.06`$ and $`F=0.1`$, then the choice of $`\beta 0.8`$ is required to allow the existence of bulges. This is just the value of $`\beta `$ used in the global bar instability criterion (Efstathiou, Lake & Negroponte 1982). For this model the disk-to-halo mass ratio varies little with radius (as shown in Figure 3a), so the B/T ratio is strongly dependent on $`\lambda /F`$. Only a small range of values of $`\lambda /F`$ is allowed, so as to not over-produce either bulge-less disks or completely unstable disks.

The relation between B/T ratio and gas fraction for this model is given in Figure 8b; there a general trend in the observed sense, with large scatter. The different curves in this plot correspond to different values of $`\lambda /F`$, which sets the overall trend.

Figure 9 (a,b) shows the equivalent plots of model E, representing a Hernquist profile halo model, which is probably a more realistic case. Here the the disk-to-halo mass ratio varies strongly with radius when approaching the centre, and the allowed parameter space for values of $`\lambda `$ and of $`F`$ that allowing the formation of bulges can be large. Overly-unstable disks are denoted by asterisk symbols in Figure 9a, and low surface-brightness disks are denoted by cross symbols. The relation between B/T ratio and gas fraction is similar to that for the isothermal halo.

### 4.1 Constraints from the Milky Way

The Milky Way bulge is reasonably well-fit by an exponential profile, with a scale-length approximately one-tenth that of the disk (Kent et al. 1991). The morphology of the bulge is consistent with some triaxiality (Blitz & Spergel 1991; Binney et al. 1997). Perhaps the Milky Way is a system in which the bulge has formed from the disk, through a bar instability? Observations show no evidence for a significant young or even intermediate-age stellar population in the field population of the Galactic Bulge (Feltzing & Gilmore 1999), despite their being ongoing star formation in the inner disk. This implies that if the bulge were formed from the disk through bar dissolution, only one such episode is allowed, and this should have happened at high redshift. In the context of the present model, the lower star formation rates and longer viscosity timescales of later times act to stabilize the system. However, it remains to be seen if the observed relative frequencies of bars, bulges and central mass concentration is consistent with the models of bar dissolution.

## 5 Summary

In the context of hierarchical clustering cosmology, the dark halo of a disk galaxy can be formed by quiescently merging small sub-halos into the primary dark halo, or by smoothly accreting matter into the dark halo. We derive the generic solution to the adiabatic infall model of disk galaxy formation pioneered by many authors (Mestel 1963; Fall & Efstathiou 1980; Gunn 1982; Faber 1982; Jones & Wyse 1983; Ryden & Gunn 1987; Dalcanton, Summer & Spergel 1997; Mo, Mao & White 1998; Hernandez & Gilmore 1998). Through exploring the allowed parameter space of dark halo profile and angular momentum distribution function, we show that the central halo density profile should be cuspy, with the power law index ranging from $`0.75`$ to $`2`$, in the central regions where the disk mass dominates.

Using a modified Schmidt law of global star formation rate, we derive a simple scaling relationship between the disk gas fraction and the assembly redshift. We explicitly allow a distribution in the values of the baryonic mass fraction, $`F`$, in addition of the distribution in values of the spin parameter $`\lambda `$. These two are found to play different role in determining the structural properties of the final disk, the star formation properties and bulge-to-disk ratio.

We mimic viscous evolution of disks by varying the specific angular momentum distribution of the disk, to redistribute angular momentum as a function of time. We derive a consistent picture of the formation of galaxies like the Milky Way, with old stars in the disk. Under the assumption that the viscous evolution timescale is equal to the star formation timescale, we can further combine the $`\lambda `$ and $`F`$ with the efficiency of angular momentum redistribution caused by viscosity. Assuming that small bulges are formed from their disks through bar dissolution, we can use the global bar instability condition to obtain bulge-to-total ratio, and explore the dependence on $`F`$, $`\lambda `$ and the viscous evolution efficiency.

The inclusion of viscous evolution has the merits of addressing several important issues: the conspiracy between disk and halo, the formation of the exponential profile of stellar disk, the high phase space density of bulges. We have presented an analytic treatment, to illustrate these points and identify areas of particular need for more work.

## Acknowledgments

We acknowledge support from NASA, ATP Grant NAG5-3928. BZ thanks Colin Norman, Jay Gallagher for helpful comments. RFGW thanks all at the Center for Particle Astrophysics, UC Berkeley, for their hospitality during the early stages of this work.


format: text

text:
# Anderson’s “Theorem” and Bogoliubov-de Gennes Equations for Surfaces and Impurities

## Abstract

In order to incorporate spatial inhomogeneity due to nonmagnetic impurities, Anderson  proposed a BCS-type theory in which single-particle states in such an inhomogeneous system are used. We examine Anderson’s proposal, in comparison with the Bogoliubov-de Gennes equations, for the attractive Hubbard model on a system with surfaces and impurities.

The procedure for examining surface and impurity effects on a microscopic level is by now well established. One uses a mean field-like decoupling, with potentials which are determined from self-consistency requirements. These potentials are then used in the effective Hamiltonian, which is numerically diagonalized. This process is continued until self-consistency is achieved. This is the essence of the Bogoliubov-de Gennes (BdG)  formalism.

An earlier proposal was suggested by Anderson , in which the single-particle problem is first diagonalized. Eigenvalues and eigenstates are obtained, with which one can formulate the BCS problem, but in a vector space associated with these eigenstates. In certain situations, the single-particle problem can be obtained analytically (open boundaries, for example ), or it can be obtained numerically with significantly less effort than required by the full BdG process. In these instances it would be advantageous to utilize the Anderson prescription. In this paper we report on some test cases to evaluate the Anderson prescription.

The BdG equations are well documented . In this work we utilize the attractive Hubbard model, with open boundaries, and with the possibility for single site impurity potentials. The resulting equations are :

$`E_nu_n(\mathrm{})={\displaystyle \underset{\mathrm{}^{}}{}}A_{\mathrm{}\mathrm{}^{}}u_n(\mathrm{}^{})+\mathrm{\Delta }_{\mathrm{}}v_n(\mathrm{})`$ (1)
$`E_nv_n(\mathrm{})={\displaystyle \underset{\mathrm{}^{}}{}}A_{\mathrm{}\mathrm{}^{}}v_n(\mathrm{}^{})+\mathrm{\Delta }_{\mathrm{}}^{}u_n(\mathrm{})`$ (2)
where

$$A_{\mathrm{}\mathrm{}^{}}=t\underset{\delta }{}\left(\delta _{\mathrm{}^{},\mathrm{}\delta }+\delta _{\mathrm{}^{},\mathrm{}+\delta }\right)+\delta _{\mathrm{}\mathrm{}^{}}\left(V_{\mathrm{}}\mu +ϵ_{\mathrm{}}\right).$$
(3)
The self-consistent potentials, $`V_{\mathrm{}}`$, and $`\mathrm{\Delta }_{\mathrm{}}`$, are given by

$`\mathrm{\Delta }_{\mathrm{}}=|U|{\displaystyle \underset{n}{}}u_n(\mathrm{})v_n^{}(\mathrm{})(12f_n)`$ (4)
$`V_{\mathrm{}}=|U|{\displaystyle \underset{n}{}}\left[|u_n(\mathrm{})|^2f_n+|v_n(\mathrm{})|^2(1f_n)\right],`$ (5)
where $`|U|`$ is the strength of the attractive interaction, the index $`n`$ labels the eigenvalues (there are $`2N`$ of them), the index $`\mathrm{}`$ labels the sites (1 through N), and the composite eigenvector is given by $`\left(\genfrac{}{}{0pt}{}{u_n}{v_n}\right)`$, of total length $`2N`$. The sums in Eqs. (4,5) are over positive eigenvalues only. The other physical parameters are the single-particle hopping, $`t`$, the single site impurity potentials, $`ϵ_{\mathrm{}}`$, and the chemical potential, $`\mu `$. The $`f_n`$ is the Fermi function, with argument $`\beta E_n`$, where $`\beta \frac{1}{k_BT}`$, with $`T`$ the temperature. The single site electron density, $`n_{\mathrm{}}`$, is given, through Eq. (5), by $`V_{\mathrm{}}=|U|\frac{n_{\mathrm{}}}{2}`$.

These equations are iterated to convergence, with results to be presented below.

The alternative Anderson formalism  first solves for the eigenvalues and eigenstates of the ‘non-interacting’ problem, i.e.,

$$E_n^0w_n(\mathrm{})=\underset{\mathrm{}^{}}{}A_{\mathrm{}\mathrm{}^{}}^0w_n(\mathrm{}^{}),$$
(6)
where

$$A_{\mathrm{}\mathrm{}^{}}^0=t\underset{\delta }{}\left(\delta _{\mathrm{}^{},\mathrm{}\delta }+\delta _{\mathrm{}^{},\mathrm{}+\delta }\right)\delta _{\mathrm{}\mathrm{}^{}}\left(\mu ϵ_{\mathrm{}}\right).$$
(7)
The $`N\times N`$ matrix equation (6) is solved for its eigenvalues $`E_n^0`$ and eigenvectors $`w_n`$. This amounts to determining the unitary matrix $`U_\mathrm{}n`$ that gives a basis for the electron operators

$$c_\mathrm{}\sigma ^{}=\underset{n}{}U_\mathrm{}n^{}\stackrel{~}{c}_{n\sigma }^{},$$
(8)
which diagonalizes the single-particle Hamiltonian. From this matrix we determine the transformed electron-electron interaction:

$$V_{nm,n^{}m^{}}=|U|\underset{\mathrm{}}{}U_\mathrm{}n^{}U_\mathrm{}m^{}U_\mathrm{}n^{}U_\mathrm{}m^{},$$
(9)
which now mediates the (generally off-diagonal) electron-electron interaction. The gap and number equations are derived in the usual way; they are in general complicated — the gap is a function of the quantum label $`n`$ and the chemical potential is shifted by an $`n`$-dependent quantity. Once these are obtained, we can transform back to real space, and examine the gap function or the electron density, for example, as a function of position.

Figure 1 illustrates the gap parameter obtained by the BdG formalism as a function of position, for all densities, in the case of open boundary conditions (OBC), in one dimension. Results in higher dimension will be very similar . Variations in the gap are strongest near the boundaries, as expected, and the Anderson prescription is reasonably accurate in reproducing the oscillations (not shown). In Fig. 2 we show the gap as a function of position for the case of a single impurity (at site 16) with a repulsive potential, with periodic boundary conditions (PBC). As expected, the gap is suppressed at this site, and once again, the Anderson prescription semi-quantitatively reproduces the BdG result.


text:
# A NEW MEASUREMENT OF THE MUON MAGNETIC ANOMALY

## 1 Physics Motivation

The magnetic anomaly of fermions $`a=\frac{1}{2}(g2)`$ describes the deviation of their magnetic g-factor from the value 2 predicted in the Dirac theory. This quantity has been measured for single electrons and positrons in Penning traps by Dehmelt and his coworkers to 10 ppb . Accurate calculations for $`a`$ of these two particles are possible to this level, which involve exclusively the ”pure” Quantum Electrodynamics (QED) of electron, positron and photon fields. The presently most accurate value for the fine structure constant $`\alpha `$  can be obtained from a comparison between experiment and theory, where it appears as an expansion coefficient. The high accuracy, to which QED calculations can be performed, is demonstrated by the compatibility of this value of $`\alpha `$ and the ones obtained in measurements based on the quantum Hall effect  or the number extracted from the precisely known Rydberg constant using an accurate determination of the neutron de Broglie wavelength and relevant mass ratios. Moreover, the agreement of $`\alpha `$ values determined from the electron magnetic anomaly and from the hyperfine splitting in the muonium atom  may be interpreted as the most precise reassurance of the internal consistency of QED, because the first case involves the theory of free particles whereas in the second case distinctively different bound state approaches need to be applied .

The anomalous magnetic moment of the muon $`a_\mu `$ is more sensitive by a factor of $`(m_\mu /m_e)^2410^4`$ to heavier particles, which appear virtually in loop graphs, and other than electromagnetic interactions. Such effects can be studied in a precise determination of $`a_\mu `$, because very high confidence in the validity of calculations of the dominating QED contribution arises from the excellent description of the electron magnetic anomaly and electromagnetic transitions in fundamental systems like, e.g. hydrogen and muonium atoms.

In a series of three experiments at CERN  $`a_\mu `$ could be measured to 7.2 ppm. This has verified the muons nature as a heavy leptonic particle and the proper description of its electromagnetic interactions to very high accuracy by QED. In the last of these measurements contributions arising from strong interactions, which amount to 57.8(7) ppm , could be verified. At BNL a new dedicated experiment has been designed to determine the muon magnetic anomaly $`a_\mu `$ with 0.35 ppm relative accuracy, meaning a 20 fold improvement over the previous approaches. At this level exists particular sensitivity to contributions arising from weak interaction through loop diagrams with W and Z bosons (1.3 ppm) . The experiment promises here a clean test of renormalization in weak interaction. The muon magnetic anomaly may also contain contributions from new physics . A variety of speculative theories can be tested which have been invented to extend the present standard model in order to explain some of the features which are described but not fundamentally understood yet. The spectrum of such theoretical models includes physics concepts like muon substructure, new gauge bosons, supersymmetry, an anomalous magnetic moment of the W boson, leptoquarks and violation of Lorentz and CPT invariance. Here a precise measurement of $`a_\mu `$ can be complementary to searches in high energy experiments and the sensitivity may even be higher.

## 2 The Brookhaven Muon g-2 Experiment

In the new experiment at the alternating gradient synchrotron (AGS) of BNL polarized muons are stored in a magnetic storage ring with highly homogeneous field $`B`$ and with weak electrostatic focussing. The difference $`\omega _a`$ of the spin precession and the cyclotron frequencies,

$$\omega _a=\omega _s\omega _c=a_\mu \frac{e}{m_\mu c}B,$$
(1)
is measured, with $`m_\mu `$ the muon mass and $`c`$ the speed of light. Positrons (electrons) from the weak decays $`\mu ^\pm e^\pm +2\nu `$ are observed. For relativistic muons the influence of a static electric field vanishes , if $`a_\mu =1/(\gamma _\mu ^21)`$ which corresponds to $`\gamma _\mu =29.3`$ and a muon momentum of $`p=3.09`$ GeV/c, where $`\gamma _\mu =1/\sqrt{1(v_\mu /c)^2}`$ and $`v_\mu `$ is the muon velocity. For sufficient accuracy of the electric field correction the average muon momentum $`p`$ needs to be within a few parts in 104 of magic momentum.

For a homogeneous field the magnet must have iron flux return and shielding. Because of the particular momentum requirement and in order to avoid strong magnetic saturation effects of the iron a device of 7 m radius was built. It has a C-shaped iron yoke cross section with the open side facing towards the center of the ring. It provides 1.4513 T field in a 18 cm gap. The magnet is energized by 4 superconducting coils carrying 5177 A current.

The magnetic field is determined by a newly developed narrow band magnetometer system which is based on pulsed nuclear magnetic resonance (NMR) of protons in water and vaseline. It has a capability for an absolute measurement to $`50`$ ppb. . The field and its homogeneity are continuously monitored by 380 NMR probes. They are distributed around the ring and they are embedded near the magnet poles in the walls of the Al vacuum tank. For mapping the field inside the storage volume a trolley carrying 17 NMR probes is run in regular intervals, typically twice a week. This device contains a fully computerized magnetometer built entirely from nonferromagnetic components. The field accuracy is derived from and related to a precision measurement of the proton gyromagnetic ratio in a spherical water sample . On average the field around the ring is homogeneous to 1 ppm (Fig.1). This has been achieved using mechanical shimming methods which include movable iron wedges in an air gap between the low carbon steel pole pieces and the magnet yoke as well as iron strips of adjusted width fixed to the neighbourhood of junctions between poles. A set of 60 electrical coils, which run on the surface of the pole pieces around the ring and which can be driven at individually different currents, allows the compensation of other than dipole components of the field. The absolute value of the field integral in the storage region is known at present to better than 0.5 ppm. There is a potential for a significant improvement in this figure. Field drifts are compensated using a set of 36 selected fixed NMR probes. Their average is kept within 0.1 ppm of the nominal value by regulating the main magnet power supply. To avoid large short term thermal effects, the magnet yoke has been dressed with passive thermal insulation material.

The weak muon focussing is provided by electrostatic quadrupole electrodes with 10 cm separation between opposite plates. They cover 43 % of the ring circumference. The electric field is applied by pulsing $`\pm 24.5`$ kV voltage for $`1.4`$ms duration to minimize electron trapping and avoid electrical breakdown. The storage volume diameter is defined by circular apertures to 9 cm.

Due to parity violation in the weak muon decay process the positrons (electrons) are emitted preferentially along (opposite) to the muon spin direction. This causes a time dependent variation of the spatial distribution of decay particles in the muon eigensystem which translates into a time dependent variation of the energy distribution in this experiment. Inside the ring the positrons(electrons) are observed in 24 shower detectors consisting of scintillating fibers embedded in lead. They have 13 readiation lengths thickness and an average resolution of $`\sigma `$/E = 6.8% at the nominal energy cut of E = 1.8 GeV which is applied in the analysis leading to the positron distribution shown in Fig. 2. All positron events are digitized individually in a custom waveform digitizer at 400 MHz rate and stored for analysis. The time standard of the detectors and the field measurement system is a single LORAN C receiver with better than $`10^{11}`$ long term stability.

The technical improvements over previous experiments at CERN include an azimuthally symmetric iron construction for the magnet with superconducting coils, a larger gap and higher homogeneity of the field, segmented positron (electron) detectors covering a larger solid angle and improved electronics. A major advantage is the two orders of magnitude higher primary proton intensity available at the AGS Booster at BNL. Further conceptually novel features are the NMR trolley, a superconducting static inflector magnet and direct muon injection onto storage orbits by means of a magnetic kicker. Previously pions had been introduced in to the ring some of which decayed into stored muons. The electrostatic quadrupoles at BNL have twice the field gradient compared to the CERN experiment and in addition the vacuum requirements are more relaxed due to a new design which minimizes electron trapping. The vacuum chamber is scalloped to avoid preshowering.

## 3 Present Status of results

By now data taking has been carried out for $`\mu ^+`$ in two extended periods. In the startup phase of the experiment in 1997  pion injection was used. The efficiency of this process was below the theoretical expectation, which is $`2510^6`$, resulting in $`10^3`$ stored muons per injection pulse (with $`510^{12}`$ protons from the AGS on target). This method is accompanied by a significant flash in the detectors caused by hadronic interactions of unused pions. The impact of this effect was minimized by gated photomultiplier operation. The data were useful only after 22-75 $`\mu `$s, depending on the detector position. The first result obtained in this way was $`a_{\mu ^+}=1165925(15)10^9`$(13 ppm).

Muon injection, which is employed regularly since 1998 with an efficiency of order 5%, gives about an order of magnitude more muons per injection pulse and largely reduces flash background. In addition, major improvements in the magnetic field homogeneity and stability were made and the detector efficiency was increased. A part of the new data ($``$ 4%) have already been completely analyzed and provides the preliminary value of $`a_{\mu ^+}=1165919(6)10^9`$ (5 ppm) where the uncertainty is dominated by statistics. Among the systematic errors the dominating contributions arise from positron pileup in the detectors, flashlets, i.e. the additional delivery of small bunches of protons after the AGS main pulse, and the field calibration. Combining all the measured values from CERN and BNL (Fig.3) yields $`a_\mu (expt)=1165921(5)10^9`$ (4 ppm). This agrees with the latest theoretical value  $`a_\mu (theor)=116591628(77)10^{11}`$ (0.66 ppm). The dominating error here arises from the knowledge of the hadronic part, which has been calculated using electron positron annihilation into hadrons and hadronic $`\tau `$-decays .

## 4 Perspectives

The data recorded up to now cover more than $`210^9`$ decay positrons. This leads to an expected statistical uncertainty at the 1 ppm level. Further data taking is in progress, now with typically $`4010^{12}`$ protons per AGS cycle which provides 10 pulses. The systematic errors are expected to sum up to a few 0.1 ppm. In order for the new muon g-2 experiment to reach its 0.35 ppm design accuracy, besides $`\omega _a`$ and the field also the muon mass respectively its magnetic moment needs to be known to 0.1 ppm or better (see eq.(1)). This has been achieved very recently by microwave spectroscopy of the Zeeman effect in the muonium atom ($`\mu ^+e^{}`$) ground state hyperfine structure, resulting in a measurement of the ratio of the muon magnetic moment to the proton magnetic moment $`\mu _\mu /\mu _p`$ to 120 ppb. (A comparison of the simultaneously obtained muonium ground hyperfine interval with QED theory may be interpreted in terms of an even more precise value of this quantity at 30 ppb.)

In minimal supersymmetric models, as a particular example of relevant speculative models, a contribution to $`a_\mu `$ of

$$\mathrm{\Delta }a_\mu (SUSY)/a_\mu 1.25\mathrm{ppm}\left(\frac{100GeV/c^2}{\stackrel{~}{m}}\right)^2\mathrm{tan}\beta ,$$
(2)
is expected, where $`\stackrel{~}{m}`$ is the mass of the lightest supersymmetric particle and $`\mathrm{tan}\beta `$ is the ratio of the vacuum expectation values for the two involved Higgs fields. At the projected accuracy for g-2, there is a sensitivity to large values of the latter parameter.

The experiment is planned for both $`\mu ^+`$ and $`\mu ^{}`$ as a test of CPT invariance. There is actual interest in view of the suggestion  to compare tests of CPT invariance in different systems on a common basis, i.e. by using the energies of the states involved. For fermion magnetic anomalies particles with spin down in an external field need to be compared to their antiparticles with spin up. The nature of g-2 experiments is such that they provide a figure of merit $`r=|a^{}a^+|\frac{\mathrm{}\omega _c}{mc^2}`$ for a CPT test, where $`a^{}`$ and $`a^+`$ are respective magnetic anomalies, and $`m`$ is the particle mass. For the past electron and positron measurements one has $`r_e1.210^{21}`$  which is a much tighter bound than from the neutral kaon system, were the mass differences between $`K^0`$ and $`\overline{K^0}`$ yield $`r_K110^{18}`$. An even more stringent CPT test arises therefore already from the past muon magnetic anomaly measurements were $`r_\mu 3.510^{24}`$. Hence, this may be viewed as the presently best known CPT test based on system energies. The BNL g-2 experiment allows to look forward to a 20 times more precise test of this fundamental symmetry.

According to the standard theory an elementary particle is not allowed to have a finite permanent electric dipole moment (edm) as this would violate CP and T symmetries, if CPT is assumed to be conserved. An edm of the muon would manifest itself in the g-2 experiment in a time dependent up down asymmetry of decay positrons which can be searched for along with the muon g-2 measurements. The BNL experiment is expected to provide one order of magnitude improvement over the present limit at $`1.0510^{18}`$ e cm. This is possible through proper segmentation of the detector packages. A further highly promising approach has been suggested as a dedicated follow on experiment. It is expected that based on the g-2 setup an experiment can be tailored to allow 5-6 orders of magnitude increase in sensitivity. It should be noted that a non standard model value of $`a_\mu `$ would call for a muon edm search, as both quantities are intimately linked in many theories, where their sizes are connected through a CP violating phase . Another possibility is using the magnet as spectrometer in which pion decays are observed for restricting the muon neutrino mass by a further factor of 20.

## 5 Acknowledgements

This work was supported in part by the U.S. Department of Energy, the U.S. National Science Foundation, the German Bundesminster für Bildung und Forschung, the Russian Ministry of Science and the US-Japan Agreement in High Energy Physics.


text:
# Quantum transport in ballistic conductors: evolution from conductance quantization to resonant tunneling

## Abstract

We study the transport properties of an atomic-scale contact in the ballistic regime. The results for the conductance and related transmission eigenvalues show how the properties of the ideal semi-infinite leads (i.e. measuring device) as well as the coupling between the leads and the conductor influence the transport in a two-probe geometry. We observe the evolution from conductance quantization to resonant tunneling conductance peaks upon changing the hopping parameter in the disorder-free tight-binding Hamiltonian which describes the leads and the coupling to the sample.

Mesoscopic physics  has changed our understanding of transport in condensed matter systems. The discovery of new effects, such as weak localization  or universal conductance fluctuations,  has been accompanied by rethinking of the established ideas in a new light. One of the most spectacular discoveries of mesoscopics is conductance quantization (CQ)  in a short and narrow constriction connecting two high-mobility (ballistic) two-dimensional electron gases. The conductance of these quantum point contacts as a function of the constriction width $`W\lambda _F`$ has steps of magnitude $`2e^2/h`$. New experimental techniques have allowed observation of similar phenomena  in metallic point contacts of atomic size. The Landauer formula  for the two-probe conductance

$$G=\frac{2e^2}{h}\text{Tr}(\mathrm{𝐭𝐭}^{})=G_Q\underset{n=1}{\overset{M}{}}T_n,$$
(1)
has provided an explanation of the stepwise conductance in terms of the number $`NM`$ of transverse propagating states (“channels”) at the Fermi energy $`E_F`$ which are populated in the constriction. Here $`𝐭`$ is the transmission matrix, $`T_n`$ transmission eigenvalues and $`G_Q=2e^2/h`$ is the conductance quantum. In the ballistic case $`(\mathrm{𝐭𝐭}^{})_{ij}`$ is $`\delta _{ij}`$, or equivalently $`T_n`$ is $`1`$. Further studies have explored CQ under a range of conditions.  They include geometry,  scattering on impurities,  temperature effects, and magnetic field.

In this paper we study the influence of the attached leads on ballistic transport ($`\mathrm{}>L`$, $`\mathrm{}`$ being elastic mean free path, $`L`$ being the system size) in a nanocrystal. We assume that in the two-probe theory an electron leaving the sample does not reenter the sample in a phase-coherent way. This means that at zero temperature phase coherence length $`L_\varphi `$ is equal to the length of the sample $`L`$. In the jargon of quantum measurement theory, the leads act as a “macroscopic measurement apparatus”. Our concern with the influence of the leads on conductance is therefore also a concern of quantum measurement theory. Recently, the effects of a lead-sample contact on quantum transport in molecular devices have received increased attention in the developing field of “nanoelectronics”.  Also, the simplest lattice model and related real-space Green function technique are chosen here in order to address some practical issues which appear in the frequent use of these methods  to study transport in disordered samples. We emphasize that the relevant formulas for transport coefficients contain three different energy scales (corresponding to the lead, the sample, and the lead-sample contact), as discussed below.

In order to isolate only these effects we pick the strip geometry in the two-probe measuring setup shown on Fig. 1.

The nanocrystal (“sample”) is placed between two ideal (disorder-free) semi-infinite “leads” which are connected to macroscopic reservoirs. The electrochemical potential difference $`eV=\mu _L\mu _R`$ is measured between the reservoirs. The leads have the same cross section as the sample. This eliminates scattering induced by the wide to narrow geometry  of the sample-lead interface. The whole system is described by a clean tight-binding Hamiltonian (TBH) with nearest neighbor hopping parameters $`t_{\mathrm{𝐦𝐧}}`$

$$\widehat{H}=\underset{𝐦,𝐧}{}t_{\mathrm{𝐦𝐧}}|𝐦𝐧|,$$
(2)
where $`|𝐦`$ is the orbital $`\psi (𝐫𝐦)`$ on the site $`𝐦`$. The “sample” is the central section with $`N_x\times N_y\times N_z`$ sites. The “sample” is perfectly ordered with $`t_{\mathrm{𝐦𝐧}}=t`$. The leads are the same except $`t_{\mathrm{𝐦𝐧}}=t_\text{L}`$. Finally, the hopping parameter (coupling) between the sample and the lead is $`t_{\mathrm{𝐦𝐧}}=t_\text{C}`$. We use hard wall boundary conditions in the $`\widehat{y}`$ and $`\widehat{z}`$ directions. The different hopping parameters introduced here have to be used to get the conductance at Fermi energies throughout the whole band extended by the disorder,  i.e. $`t_\text{L}>t`$. Thus, one has to be aware of the conductances we calculate in our analysis when engaging in such studies.

Our toy model shows exact conductance steps in multiples of $`G_Q`$ when $`t_\text{C}=t_\text{L}=t`$. This is a consequence of infinitely smooth (“ideally adiabatic” ) sample-lead geometry. Then we study the evolution of quantized conductance into resonant tunneling conductance while changing the parameter $`t_\text{L}`$ of the leads as well as the coupling between the leads and the conductor $`t_\text{C}`$. An example of this evolution is given on Fig. 2. The equivalent evolution of the transmission eigenvalues $`T_n`$ of channels is shown on Fig. 3. A similar evolution has been studied recently in one-atom point contacts.

The non-zero resistance is a purely geometrical effect

caused by reflection when the large number of channels in the macroscopic reservoirs matches the small number of channels in the lead. The sequence of steps ($`1,3,6,5,7,5,6,3,1`$ multiples of $`G_Q`$ as the Fermi energy $`E_F`$ is varied) is explained as follows. The eigenstates in the leads, which comprise the scattering basis, have the form $`\psi _𝐤\mathrm{sin}(k_ym_y)\mathrm{sin}(k_zm_z)e^{ik_xm_x}`$ at atom $`𝐦`$, with energy $`E=2t_\text{L}[\mathrm{cos}(k_xa)+\mathrm{cos}(k_ya)+\mathrm{cos}(k_za)]`$, where $`a`$ is the lattice constant. The discrete values $`k_y(i)=i\pi /(N_y+1)a`$ and $`k_z(j)=j\pi /(N_z+1)a`$ define subbands or “channels” labeled by $`(k_y,k_z)(i,j)`$, where $`i`$ runs from $`1`$ to $`N_y`$ and $`j`$ runs from $`1`$ to $`N_z`$. The channel $`(k_y,k_z)`$ is open if $`E_F`$ lies between the bottom of the subband, $`2t_\text{L}[1+\mathrm{cos}(k_ya)+\mathrm{cos}(k_za)]`$, and the top of the subband, $`2t_\text{L}[1+\mathrm{cos}(k_ya)+\mathrm{cos}(k_za)]`$. Because of the degeneracy of different transverse modes in 3D, several channels $`(k_y,k_z)`$ open or close at the same energy. Each channel contributes one conductance quantum $`G_Q`$. This is shown on Fig. 2 for a sample with $`3\times 3`$ cross section where the number of transverse propagating modes is $`M=9`$. In the adiabatic geometry, channels do not mix, i.e. the transmission matrix is diagonal in the basis of channels defined by the leads.

We compute the conductance using the expression obtained in the framework of Keldysh technique by treating the coupling between the central region and the lead as a perturbation.  This provides the following, Landauer-type, formula for the conductance in the non-interacting system

$`G`$ $`=`$ $`{\displaystyle \frac{2e^2}{h}}\text{Tr}\left(\widehat{\mathrm{\Gamma }}_L\widehat{G}_{1N_x}^r\widehat{\mathrm{\Gamma }}_R\widehat{G}_{N_x1}^a\right)={\displaystyle \frac{2e^2}{h}}\text{Tr}(\mathrm{𝐭𝐭}^{}),`$ (3)
$`𝐭`$ $`=`$ $`\sqrt{\widehat{\mathrm{\Gamma }}_L}\widehat{G}_{1N_x}^r\sqrt{\widehat{\mathrm{\Gamma }}_R}.`$ (4)
Here $`\widehat{G}_{1N_x}^r`$, $`\widehat{G}_{N_x1}^a`$ are matrices whose elements are the Green functions connecting the layer $`1`$ and $`N_x`$ of the sample. Thus only the block $`N_y\times N_z`$ of the whole matrix $`\widehat{G}(𝐧,𝐦)`$ is needed to compute the conductance. The positive operator $`\widehat{\mathrm{\Gamma }}_L=i(\widehat{\mathrm{\Sigma }}_L^r\widehat{\mathrm{\Sigma }}_L^a)=2\text{Im}\widehat{\mathrm{\Sigma }_L}>0`$ is the counterpart of the spectral function $`\widehat{A}=i(\widehat{G}^r\widehat{G}^a)`$ for the self-energy $`\widehat{\mathrm{\Sigma }}_L`$ introduced by the left lead. It “measures” the coupling of the open sample to the left lead ($`\widehat{\mathrm{\Gamma }}_R`$ is equivalent for the right lead).The Green operator is defined as the inverse of $`(E\widehat{H})`$ including the relevant boundary conditions. Instead of inverting the infinite matrix we invert only $`(E\widehat{H}_S)`$ defined on the Hilbert space spanned by orbitals $`|𝐦`$ inside the sample

$$\widehat{G}^r=(E\widehat{H}_S\widehat{\mathrm{\Sigma }}^r)^1,$$
(5)
where $`\widehat{H}_S`$ is TBH for the sample only. This is achieved by using the retarded self energy $`\widehat{\mathrm{\Sigma }}^r=\widehat{\mathrm{\Sigma }}_L^r+\widehat{\mathrm{\Sigma }}_R^r`$ introduced by the left $`(L)`$ and the right $`(R)`$ lead. In site representation Green operator $`\widehat{G}^{r,a}`$ is a Green function matrix $`\widehat{G}^{r,a}(𝐧,𝐦)=𝐧|\widehat{G}^{r,a}|𝐦`$. Equation (5) does not need the small imaginary part $`i0^+`$ necessary to specify the boundary conditions for the retarded or advanced Green operator $`\widehat{G}^{r,a}`$ because the lead self-energy $`(\widehat{\mathrm{\Sigma }}^a=[\widehat{\mathrm{\Sigma }^r}]^{})`$ adds a well defined imaginary part to $`E\widehat{H}_S`$. This imaginary part is related to the average time an electron spends inside the sample before escaping into the leads. The self-energy terms have non-zero matrix elements only on the edge layers of the sample adjacent to the leads. They are given  in terms of the Green function on the lead edge layer and the coupling parameter $`t_\text{C}`$

$`\widehat{\mathrm{\Sigma }}_{L,R}^r(𝐧,𝐦)`$ $`=`$ $`{\displaystyle \frac{2}{N_y+1}}{\displaystyle \frac{2}{N_z+1}}{\displaystyle \underset{k_y,k_z}{}}\mathrm{sin}(k_yn_y)\mathrm{sin}(k_zn_z)`$ (7)
$`\times \widehat{\mathrm{\Sigma }}^r(k_y,k_z)\mathrm{sin}(k_ym_y)\mathrm{sin}(k_zm_z),`$
where $`(𝐧,𝐦)`$ is the pair of sites on the surfaces inside the sample which are adjacent to the leads ($`L`$ or $`R`$). The self-energy $`\widehat{\mathrm{\Sigma }}^r(k_y,k_z)`$ in the channel $`(k_z,k_y)`$ is given by

$$\widehat{\mathrm{\Sigma }}^r(k_y,k_z)=\frac{t_\text{C}^2}{2t_\text{L}^2}\left(E_\mathrm{\Sigma }i\sqrt{4t_\text{L}^2E_\mathrm{\Sigma }^2}\right),$$
(8)
for $`|E_\mathrm{\Sigma }|<2t_\text{L}`$. We use the shorthand notation $`E_\mathrm{\Sigma }=E\epsilon (k_y,k_z)`$, where $`\epsilon (k_y,k_z)=2t_\text{L}[\mathrm{cos}(k_ya)+\mathrm{cos}(k_za)]`$ is the energy of quantized transverse levels in the lead. In the opposite case $`|E_\mathrm{\Sigma }|>2t_\text{L}`$ we have

$$\widehat{\mathrm{\Sigma }}^r(k_y,k_z)=\frac{t_\text{C}^2}{2t_\text{L}^2}\left(E_\mathrm{\Sigma }\text{sgn}E_\mathrm{\Sigma }\sqrt{E_\mathrm{\Sigma }^24t_\text{L}^2}\right).$$
(9)
In order to study the conductance as a function of two parameters $`t_\text{L}`$

and $`t_\text{C}`$ we change either one of them while holding the other fixed (at the unit of energy specified by $`t`$), or both at the same time. The first case is shown on Fig. 2 and Fig. 4 (upper panel), while the second one on Fig. 4 (lower panel). The conductance is depressed in all cases since these configurations of hopping parameters $`t_{\mathrm{𝐦𝐧}}`$ effectively act as a barriers. There is a reflection at the sample-lead interface due to the mismatch of the subbands in the lead and in the sample when $`t_\text{L}`$ differs from $`t`$. This demonstrates that adiabaticity is not necessary condition for CQ. In the general case, each set of channels which have the same energy subband is characterized by its own transmission function $`T_n(E_F)`$. When the coupling $`t_\text{C}=0.1`$ is small a double-barrier structure is obtained which has a resonant tunneling conductance. The electron tunnels from one lead to the other via discrete eigenstates. The transmission function is composed of peaks centered at $`E_r=2t[\mathrm{cos}(k_xa)+\mathrm{cos}(k_ya)+\mathrm{cos}(k_za)]`$, where $`k_x=k\pi /(N_x+1)a`$ is now quantized inside the sample, i.e. $`k`$ runs from $`1`$ to $`N_x`$. The magnitude and width of peaks is defined by the rate at which an electron placed between barriers leaks out into the lead. These rates are related to the level widths generated through the coupling to the leads. In our model they are energy (i.e. mode) dependent. For example at $`E_F=0`$ seven transmission eigenvalues are non-zero (in accordance with open channels on Fig. 3) and exactly at $`E_F=0`$ three of them have $`T=1`$ and four $`T=0.5`$. Upon decreasing $`t_\text{C}`$ further all conductance peaks, except the one at $`E_F=0`$, become negligible. Singular behavior of $`G(E_F)`$ at subband edges of the leads was observed before.

It is worth mentioning that the same results are obtained using a non-standard version of Kubo-Greenwood formula  for the volume averaged conductance

$`G`$ $`=`$ $`{\displaystyle \frac{4e^2}{h}}{\displaystyle \frac{1}{L_x^2}}\text{Tr}\left(\mathrm{}\widehat{v}_x\text{Im}\widehat{G}\mathrm{}\widehat{v}_x\text{Im}\widehat{G}\right),`$ (11)
$`\text{Im}\widehat{G}`$ $`=`$ $`{\displaystyle \frac{1}{2i}}(\widehat{G}^r\widehat{G}^a),`$ (12)
where $`v_x`$ is the $`x`$ component of the velocity operator. This was originally derived for an infinite system without any notion of leads and reservoirs. The crucial non-standard aspect is use of the Green function (5) in formula (4). This takes into account, through lead self-energy (7), the boundary conditions at the reservoirs. The reservoirs are necessary in both Landauer and Kubo formulations of linear transport for open finite systems. They provide thermalization and thereby steady state of the transport in the central region. Semi-infinite leads  are a convenient method to model the macroscopic reservoirs. When employing the Kubo formula (4) one can use current conservation and compute the trace only on two adjacent layers inside the sample. To get the correct results in this scheme $`L_x`$ in Eq. (4) should be replaced  by a lattice constant $`a`$.

In the quantum transport theory of disordered systems the influence of the leads on the conductance of the sample is understood as follows.  An isolated sample has a discrete energy spectrum. Attaching leads necessary for transport measurements will broaden energy levels. If the level width $`\mathrm{\Gamma }`$ due to the coupling to leads is larger than the Thouless energy $`E_{\text{Th}}=\mathrm{}/\tau _\text{D}\mathrm{}D/L^2`$, ($`D=v_F\mathrm{}/3`$ being the diffusion constant) the level discreteness is unimportant for transport. For our case of ballistic conduction, $`E_{\text{Th}}`$ is replaced by the inverse time of flight $`\mathrm{}v_F/L`$. In the disordered sample where $`\mathrm{\Gamma }E_{\text{Th}}`$, varying the strength of the coupling to the leads will not change the transport coefficients. In other words, the intrinsic resistance of the sample is much larger than the resistance of the lead-sample contact.  In the opposite case, discreteness of levels becomes important and the strength of the coupling defines the conductance. This is the realm of quantum dots  where weak enough coupling can make the charging energy $`e^2/2C`$ of a single electron important as well. Changing the properties of the dot-lead contact affects the conductance, i.e. the result of measurement depends on the measuring process. The decay width $`\mathrm{\Gamma }=\mathrm{}/\tau _{\text{dwell}}`$ of the electron emission into one of the leads is determined by transmission probabilities of channels through the contact and mean level spacing.  This means that mean dwell time $`\tau _{\text{dwell}}`$ inside our sample depends on both $`t_\text{C}`$ and $`t_\text{L}`$. Changing the hopping parameters will make $`\tau _{\text{dwell}}`$ greater than the time of flight $`\tau _f=L/v_F`$. Thus we find that ballistic conductance sensitively depends on the parameters of the dephasing environment (i.e. the leads).

In conclusion, we have studied the transport properties of a ballistic nanocrystal placed between two semi-infinite leads in the simplest strip geometry. We observe extreme sensitivity of the conductance to changes in the hopping parameter in the leads as well as the coupling between the leads and the sample. As can be easily anticipated, the conductance evolves from perfect quantization (as a result of an ideal adiabatic geometry) to resonant tunneling. Nevertheless, it is quite amusing that vastly different $`G(E_F)`$ are obtained between these two limits (e.g. Fig. 4). The results are of relevance for the analogous theoretical studies in disordered conductors as well as in the experiments using clean metal junctions with different effective electron mass throughout the circuit.

This work was supported in part by NSF grant no. DMR 9725037. We thank I. L. Aleiner for interesting discussions.