text:
# The abundance of dark galaxies

## 1 Introduction

There is now growing evidence that the observed number of dwarf galaxies surrounding the Milky Way or M31 does not agree with the predictions of high–resolution N-body simulations for the currently favored Cold Dark Matter cosmological model (flat, low matter-density Universe – LCDM) (Klypin et al., 1999; Moore et al., 1999). Given the successes of the current LCDM paradigm in reproducing many of the large-scale observables (e.g., Jaffe et al (2001); Peacock et al (2001); Verde et al. (2001); Lahav et al. (2001)), it seems natural to try to make changes to the existing paradigm without modifying its large-scale properties. There are two routes in this direction: one is to prevent the formation of small dark matter halos, the second is to hide them by suppressing star formation. Some modifications belonging to the first group consist of changing the nature of the dark matter itself to have different physical properties: self–interacting (Spergel & Steinhardt, 2000), warm (Colombi, Dodelson & Widrow, 1996; Bode, Ostriker & Turok, 2001) or other variants (e.g., Goodman (2000); Hu, Barkana & Gruzinov (2000); Cen (2001)). It has also been proposed that the solution to the problem may lie within the nature of the initial conditions and that some modification to the shape of the inflaton potential may suppress small scale structure (Kamionkowski & Liddle, 2000), although this solution faces some problems in reproducing the observed amount of structure in the Lyman-$`\alpha `$ forest. On the other hand, the solution of suppressing star formation via supernova feedback was suggested a long time ago (e.g., White & Rees (1978)). This process expels the gas from the halo and therefore suppresses star formation. The accretion and cooling of gas can also be suppressed in the presence of a strong photo–ionizing background. This route has also been known for a number of years now (Doroshkevich, Zel’Dovich & Novikov, 1967) and has been investigated recently in the context of semi–analytic models of galaxy formation (Bullock et al., 2001; Chiu et al., 2001; Somerville, 2001; Benson et al., 2001a).

A closely related problem is the fact that the observed faint end of the luminosity function has a shallower slope (Blanton et al., 2001) than predicted from high-resolution LCDM simulations, if one simply assumes a constant mass-to-light ratio (e.g., Jenkins et al. (2001)). Stellar feedback is again often used in semi–analytic models to suppress star formation in halos with shallow potential wells (e.g., Cole et al. (2000)).

All of the above models seek to systematically increase the mass-to-light ratio in small halos, by depleting the cold gas fraction available for star formation. In this paper we show that such drastic expulsion of gas is, perhaps, unnecessary: provided that the gas conserves angular momentum during collapse, a large fraction of low mass halos will be Toomre stable (Toomre, 1964; Kennicutt, 1989), i.e., $`Q=\mathrm{\Sigma }_c/\mathrm{\Sigma }>1`$ (where $`\mathrm{\Sigma }`$ is the disk surface mass density and $`\mathrm{\Sigma }_c`$ is the critical surface density to trigger gravitational instability and thus start formation). If the baryon fraction ($`f_d`$) in disks is $`f_d0.5\mathrm{\Omega }_b/\mathrm{\Omega }_m`$, all disks in halos with masses $`M<10^9M_{}`$ will be Toomre stable and fail to form stars. The “missing” dwarf galaxies predicted to be around the Milky Way and M31 may fall in this category. We also show that, if a small fraction of baryons initially present in the dark matter halo settle in a disk, as some observational evidence suggests (Jimenez, Verde & Oh, 2002; Guzik & Seljak, 2001; van den Bosch, Burkert & Swaters, 2001; Burkert, 2000; Fukugita, Hogan & Peebles, 1998), even more massive halos have some probability of remaining dark, provided they have high spin. This could potentially reconcile the observed luminosity function with the mass function predicted in high–resolution LCDM N-body simulations. We show that the predictions of this model are consistent with our previous analysis of the rotation curves of galaxies (Jimenez, Verde & Oh (2002), hereafter JVO02): none of the galaxies lie in the region of parameter space forbidden by the Toomre instability criterion.

Throughout this paper we assume a LCDM cosmology given by:$`(\mathrm{\Omega }_m,\mathrm{\Omega }_\mathrm{\Lambda },\mathrm{\Omega }_b,h,\sigma _8h^1)=(0.3,0.7,0.039,0.7,1.0)`$.

## 2 The mass distribution of dark galaxies

We model spiral galaxies as exponential disks embedded in a cold dark matter halo. Following Dalcanton et al. (1997); Jimenez et al. (1997); Mo et al. (1998) we assume a NFW (Navarro, Frenk & White, 1997) profile $`\rho [(r/r_c)(1+r/r_c)^2]^1`$ where $`r_c`$ is the break radius and the surface mass density of the disk is given by $`\mathrm{\Sigma }=\mathrm{\Sigma }_0\mathrm{exp}[r/Rd]`$, where $`R_d`$ is the scale length of the disk.

For a baryonic disk to be locally gravitationally unstable, despite the stabilizing effects of tidal shears and pressure forces, we require the Toomre parameter $`Q<1`$, where:

$$Q=\frac{c_\mathrm{s}\kappa }{\pi G\mathrm{\Sigma }},$$
(1)
$`\mathrm{\Sigma }`$ is the disk surface mass density, $`\kappa =1.41(V/r)(1+\mathrm{d}\mathrm{ln}V/\mathrm{d}\mathrm{ln}r)^{1/2}`$ is the epicyclic frequency and $`c_s`$ is the gas sound speed.

Regions in local disk galaxies where $`Q>1`$ are observationally associated with very little star formation, indicating that the Toomre criterion is obeyed remarkably well, if a gas sound speed of $`c_s6\mathrm{km}\mathrm{s}^1`$ is assumed (Kennicutt, 1989). Ferguson et al. (1996) noted that for one galaxy (NGC 6946) the agreement with the Toomre criterion was not perfect: there was no truncation of star formation for $`Q>1.6`$ if the observed vertical velocity dispersion was adopted. However, agreement with the Toomre criterion could be obtained if the velocity dispersion was assumed to be constant with radius. Wong & Blitz (2001) have measured the value of $`Q`$ in the inner part of rotation curves for a sample of 8 galaxies. For half of the sample the predicted value from Toomre theory for the abundance of gas is in excellent agreement with observations. In the other cases the agreement is not as remarkable, but still there is no evidence for large deviations from the Toomre stability predictions.

Jimenez et al. (1997) developed a disk model to study the dependence of $`Q`$ on mass, radius and spin parameter of the dark matter halo and found that for a certain range of the above parameters, some disks would be dark. An important ingredient in their study was the assumption of a constant baryon fraction, namely they adopted the nucleosynthesis value from Walker et al. (1991). Using this constant disk fraction for all galaxies, they concluded that for an Einstein–de Sitter Universe the minimum dark matter halo mass for a visible disk galaxy is about $`10^9`$ M. On the other hand, for a low matter density universe (e.g., LCDM), this minimum mass drops substantially, by about two orders of magnitude, and thus the fraction of dark halos also falls drastically. The reason for this is that for a low density Universe with baryon density fixed by nucleosynthesis, the baryons constitute a higher fraction of the halo mass.

It seems unlikely that all baryons will be able to cool and form stars. Recent studies from galaxy–galaxy lensing from the Sloan Digital Sky Survey (SDSS) show that for $`L^{}`$ galaxies the baryon fraction is about half of the nucleosynthesis value (Guzik & Seljak, 2001). On the other hand, theoretical considerations suggest that the baryon fraction for galaxies, either less massive or more massive than $`L_{}`$, is likely to be much smaller due to feedback and inefficient cooling respectively (e.g. van den Bosch (2001); Benson et al. (2001b)). This implies that a large fraction of disks should be Toomre stable. Below we quantify how this affects the predicted number density of observable galaxies.

### 2.1 Luminosity function

CDM models typically predict many more galaxies at the faint end of the luminosity function than are actually observed. This can most easily be seen by comparing the galaxy circular velocity function with theoretical predictions. The galaxy velocity function is obtained by combining the observed luminosity function with empirically determined luminosity–velocity relations (Tully-Fisher relation for spirals, and Faber-Jackson for ellipticals). Gonzalez et al. (2000) find that treating the entire galaxy population as spirals does not significantly alter the derived velocity function. In fig. 1, we plot a representative velocity function from their paper, derived from the SSRS2 (Marzke et al., 1998) B band survey and the Yasuda et al. (1997) Tully-Fisher relation. We also construct a new velocity function from the Sloan Digital Sky Survey (SDSS), which has published luminosity functions in 5 bands (Blanton et al., 2001). We choose to use the R band, as the R band Tully-Fisher relation is tighter and has less scatter than the Tully-Fisher relation at shorter wavelengths. Blanton et al. (2001) perform the color transformations from SDSS magnitudes to LCRS (Shectman et al., 1996) R-band magnitudes (see their Table 3); we perform the color transformations and corrections for internal extinction as in Gonzalez et al. (2000) to the $`R_{\mathrm{courteau}}`$ magnitudes used by Courteau (1997). We then use the Courteau (1997) $`R_{\mathrm{courteau}}`$ band Tully-Fisher relation to obtain the velocity function. Note that the Sloan velocity function has a significantly higher normalization than the SSRS2 velocity function. This is due to the use of Petrosian magnitudes for SDSS, which allows detection of galaxies of significantly lower surface brightness. The luminosity density of galaxies in the SDSS is therefore significantly higher than other redshift surveys. The SDSS luminosity function is very similar to the luminosity functions derived by other surveys such as the Las Campanas Redshift Survey (LCRS) and Two Degree Field Galaxy Redshift Survey (2dFGRS) if they re-analyze their data using the shallower isophotal limits for galaxy magnitude employed by these surveys; see Blanton et al. (2001) for discussion.

The theoretical velocity function is given by:

$$\mathrm{\Psi }(V_c)=\frac{dN}{d\mathrm{log}V_c}=\frac{dN}{dM}\frac{dM}{d\mathrm{log}(V_c)}f_{\mathrm{bright}},$$
(2)
where $`V_c`$ is the circular velocity in the flat part of the rotation curve and $`f_{\mathrm{bright}}`$ is the fraction of halos which form stars. For the mass function $`dN/dM`$, we use the Sheth & Tormen (1999) mass function, which gives a good fit to the mass function seen in the Hubble Volume simulations (see Jenkins et al. (2001)). For evaluating the derivative $`\frac{dM}{d\mathrm{log}(V_c)}`$, we use the fitting formula of Mo, Mao & White (1998), which gives the circular velocity $`V_c`$ as a function of the virial mass of the dark halo ($`M_{200}`$), the concentration parameter ($`c`$), the spin of the dark halo ($`\lambda `$) and $`f_d`$. As in Navarro, Frenk & White (1997), for a given $`M_{200}`$ we evaluate the concentration parameter from the typical collapse redshift. Below we discuss our choice of $`f_d,\lambda `$ further and evaluate $`f_{\mathrm{bright}}`$ using Toomre instability arguments.

Let us start by assuming $`f_{\mathrm{bright}}=1`$, a constant disk mass fraction $`f_d=0.06`$ (i.e., about half the nucleosynthesis value (Fukugita et al., 1998)), and the disk spin parameter to be the mean value found in N-body simulations $`\lambda _{\mathrm{mean}}=0.042`$ (Bullock et al., 2001). In figure 1 we plot the theoretical velocity function (upper solid line labeled by $`f_d=0.06`$). At low velocities, the predicted number density of galaxies is much larger than observed, which is the discrepancy we seek to resolve.

This cannot be resolved by making $`f_d`$ a function of $`M_{200}`$, unless the $`f_d(M_{200})`$ dependence is unrealistically drastic. To illustrate this we assume a different functional form for $`f_d(M_{200})`$, as shown in fig. 2. Our ansatz for the $`f_d(M_{200})`$ relation models a fairly mild factor of two decline in the value of $`f_d`$ from high mass to low mass halos, which fits the median value of $`f_d`$ at high masses (Mo, Mao & White, 1998) and low masses (van den Bosch & Swaters, 2001). If this was the only process at play it would produce only a mild change in the theoretical velocity function, as shown by the thin solid line labeled by NDG (no dark galaxies) in figure 1.

To include the effects of the Toomre instability criterion we proceed as follows. We define a critical spin parameter $`\lambda _{\mathrm{crit}}`$ such that the disc is marginally stable: $`Q[f_d(M_{200}),\lambda _{\mathrm{crit}},M_{200},c(M_{200})]=1.5`$, where $`Q`$ is evaluated at one disk scale-length, and we assume a gas sound speed $`c_s=6\mathrm{km}\mathrm{s}^1`$ (the gas will in general be photoionized by the meta-galactic ionizing background to $`T10^4`$K). The fraction of bright galaxies is given by $`f_{\mathrm{bright}}=_0^{\lambda _{\mathrm{crit}}}p(\lambda )𝑑\lambda `$, where $`p(\lambda )`$ is the probability distribution of $`\lambda `$ derived from numerical simulations (Bullock et al., 2001). Note that since halos with $`\lambda >\lambda _{\mathrm{crit}}`$ are assumed not to host visible galaxies, the probability distribution of $`\lambda `$ in observed galaxies is skewed. We evaluate the circular velocity assuming a median value $`\lambda _{\mathrm{median}}`$, defined by the implicit equation $`_0^{\lambda _{\mathrm{median}}}p(\lambda )𝑑\lambda /f_{\mathrm{bright}}=0.5`$. The velocity function we obtain from computing $`f_{\mathrm{bright}}`$ in this manner is shown in figure 1 by the thick solid line labeled by DG (dark galaxies). Thus, once we include the suppression of high-spin galaxies by explicitly computing $`f_{\mathrm{bright}}`$, the theoretical predictions match the observations well.

The suppression of accretion after reionization (Bullock et al., 2001; Chiu et al., 2001; Somerville, 2001; Benson et al., 2001a), is an attractive mechanism for regulating $`f_d`$ since (unlike for instance, supernova explosions) it affects dark galaxies as well. We model this effect by using the fitting formula (Gnedin, 2000; Somerville, 2001):

$$f_d^{\mathrm{photo}}(M_{200},z_f)=\frac{\mathrm{\Omega }_b}{\mathrm{\Omega }_m}\frac{1}{[1+0.26\frac{M_\mathrm{C}}{M_{200}}]^3}$$
(3)
where we approximate $`M_C`$ as the mass corresponding to a halo with virial velocity $`V_c=50\mathrm{km}\mathrm{s}^1`$ (Somerville, 2001) at the redshift of formation $`z_\mathrm{f}`$. For a given halo mass $`M_{200}`$, extended Press-Schechter theory (Lacey & Cole, 1994) gives the probability distribution of collapse redshifts $`p(z_f)`$; we can then use equation (3) to obtain the probability distribution $`p(f_d)`$. Note that $`f_d=\mathrm{\Omega }_b/\mathrm{\Omega }_m`$ for $`z>z_{\mathrm{reion}}`$, where $`z_{\mathrm{reion}}`$ is the redshift of reionization; we assume $`z_{\mathrm{reion}}7`$ (the results depend only very weakly on $`z_{\mathrm{reion}}`$).

In the photoionization model, the fraction of dark galaxies can be computed via equation 4, with the probability distribution $`p(f_d)`$ given above. Significant suppression of the number density of low circular velocity galaxies only takes place off the scales plotted in fig 1. Photoionization only successfully suppresses gas accretion in very small halos; it succeeds in alleviating the dwarf satellite problem, but leaves unaffected halos of $`V_c>75`$ Km s<sup>-1</sup>. By itself it therefore cannot produce the low values of $`f_d`$ required to reconcile the theoretical and observed velocity function.

We do not attempt to model the $`f_d(M_{200})`$ relation in more detail (which is the task of hydrodynamic galaxy formation modelling), but merely note that it is consistent with observations and theoretical prejudice that gas expulsion and suppression of accretion are more important for shallow potential wells. Possible mechanisms include suppression of accretion due to entropy injection in the intergalactic medium by supernovae and/or AGN jets, or ram pressure stripping as a smaller halo falls into a larger halo. We note merely that by invoking the Toomre instability criterion, we do not require the very strong variation of $`f_d`$ with $`M_{200}`$ required in most semi-analytic models of galaxy formation. In particular, we do not require a very sharp drop in $`f_d`$ at the faint end.

### 2.2 The fraction of dark galaxies

In figure 3 we plot the fraction of dark galaxies $`f_{\mathrm{dark}}`$ as a function of halo mass.

The dashed line is $`f_{\mathrm{dark}}=1f_{\mathrm{bright}}`$, where $`f_{\mathrm{bright}}`$ has been computed in the previous section. The dotted line shows the resulting $`f_{\mathrm{dark}}`$ if photoionization was the only process at play. As already noted this mechanism only suppresses accretion in halos of low circular velocity and cannot reconcile the observed and theoretical velocity function. Assuming $`f_d`$ as a deterministic function of $`M_{200}`$ might be an oversimplification: in reality $`f_d`$ is likely to be highly stochastic (e.g., van den Bosch (2001); Benson et al. (2001b)). Also, analysis of rotation curves (JVO02) confirms that there is a significant scatter in the relation (see §3). We thus compute also the expected fraction of dark galaxies by neglecting any possible mass dependence and using the distribution of $`f_d`$ ($`p(f_d)`$) empirically obtained by fitting the rotation curves for the NFW profile (see JVO02 and §3):

$$f_{\mathrm{dark}}=_{\lambda _{crit}}^{\mathrm{}}𝑑\lambda _0^{f_{d\mathrm{crit}}}𝑑f_dp(f_d)p(\lambda )$$
(4)
where $`\lambda _{\mathrm{crit}}`$ and $`f_{d\mathrm{crit}}`$ are defined by $`Q(M_{200},f_{d\mathrm{crit}},\lambda _{\mathrm{crit}},c)=1`$. The result is shown in figure 3 as a solid line. The two approaches (the one described in §2.1 or using the empirically obtained $`f_d`$ distribution) yield very similar values for $`f_{\mathrm{dark}}`$, indication of the fact that the main mechanism to create dark galaxies is the Toomre criterion.

### 2.3 Dwarf satellites

From fig. 3 it is clear that star formation within dark matter halos with masses below $`10^{10}`$ M will be almost completely suppressed (not visible because they will be Toomre stable). Thus this provides a mechanism to hide these low mass halos. To make a more quantitative estimate of the viability of this effect, in figure 4 we compare the observed cumulative mass function of dwarf satellite with predictions from N-body simulations. The solid line is from Moore et al. (1999), where the velocity has been converted in mass. The gray area is an estimate of the error from fig. 2 of Moore et al. (1999). The solid and dashed lines are the predicted cumulative number of satellites once the probability that some galaxies will remain dark has be taken into account. In particular we have used the predictions of the solid and dashed lines of figure 3 respectively. This prediction agrees nicely with observations (triangles). We have not shown results for masses below $`10^9`$ M because the uncertainties in the procedure to obtain the corrected line become too important; only a very small fraction ($`<1`$%) of dark matter halos will be visible thus making this fraction very dependent on the exact value of the predicted number of dark galaxies.

## 3 Rotation Curves

In this section we present some indication of the existence of dark galaxies from analysis of a large sample of rotation curves.

We use the same set of rotation curves as in JVO02 where we determined the best fit disk parameters within the context of both a disk within a NFW profile halo and within a pseudo-isothermal profile halo (see JVO02 for details).

The two sets of observed rotation curves are:

* Courteau catalog (Courteau, 1997) which consists of optical ($`H_\alpha `$) long–slit rotation curves for 300 Sb–Sc UGC galaxies.
* 64 spiral galaxies (Sa–Sd) observed in $`H_\alpha `$ by Palunas & Williams (2000)

In total we use 364 optical rotation curves (to avoid any problems with the smearing of the beam in radio observations (van den Bosch & Swaters, 2001)) corrected for inclination and with error bars. For these galaxies the scale length of the exponential disk ($`R_d`$) is also measured. The two models have three free parameters: $`M_{200}`$, $`r_c`$, and $`f_d`$. In principle there is another disk parameter: the modified spin parameter $`\lambda ^{}`$. If $`\lambda `$ is the spin parameter of the dark halo, then $`\lambda ^{}=j_d/f_d\lambda `$ where $`j_dJ_d/J`$ and $`J_d,J`$ are the total angular momentum of the disk and halo respectively. $`\lambda ^{}`$ can be readily obtained using the measured value of $`R_d`$ and the best fit parameters (see JVO02). If angular momentum is conserved, then $`\lambda =\lambda ^{}`$, i.e. $`j_d/f_d=1`$.

In JVO02 we used the sample of 364 galaxies for which $`R_d`$ is measured, we found the best fitting parameters ($`M_{200}`$, $`f_d`$, and $`c`$) for the above two models and recovered $`\lambda ^{}`$. Within both models for the majority of the galaxies, a meaningful set of fitting parameters was recovered. For the remaining galaxies we obtained values for $`f_d`$ that where either implausibly low ($`f_d0.001`$) or may be too high ($`f_d>\mathrm{\hspace{0.33em}0.2}`$).

### 3.1 Fit to the disk models

As described in detail in JVO02, for each galaxy rotation curve and for both models, the best fitting parameters were obtained using a standard $`\chi ^2`$ minimization and exploring the whole likelihood surface. The sample of galaxies studied is very non-uniform: some curves have much smaller error-bars than others, some curves do not extend to large enough radii to show the flattening of the rotation curve, while others show strong evidence of spiral arms and bars in the rotation pattern, moreover we have not attempted to model any bulge nor bar component. Thus it is is important to keep in mind that galaxies might be more complicated than our model for galaxy rotation curves. However, for most of these galaxies, the high-quality of the rotation curves measurements allow degeneracies among the disk parameters to be lifted (JVO02). Here, we are not interested in a detailed modeling of individual galaxy dynamics but rather the general statistical trends of the recovered disk parameters from the whole sample.

Fig. 5 shows the distribution of the best fitting disk parameters in the $`f_d`$-$`\lambda ^{}`$ plane and $`M_{200}`$-$`f_d`$ plane, for the NFW profile (top) and pseudoisothermal profile (bottom); diamonds correspond to galaxies from the Courteau (1997) sample and triangles are from Palunas & Williams (2000). In JVO02 we showed that low surface brightness galaxies (LSBs) follow the same trend as the rest of the sample, so any observational bias against LSBs does not seem to affect this result: the only difference is that LSBs have systematically higher values of $`\lambda ^{}`$ for a fixed value of $`f_d`$. In both models there is a correlation between $`\lambda ^{}`$ and $`f_d`$, albeit with a scatter, and a weak anti-correlation between $`M_{200}`$ and $`f_d`$.

For $`f_d<0.02`$, $`\lambda ^{}`$ values are between 0.01 and 0.03, while for $`f_d>0.02`$, $`\lambda ^{}`$ takes all values between 0.02 and 0.4.

The origin of this correlation is controversial: van den Bosch (2001) argues that when feedback by star formation is included in models of galaxy formation, such a correlation between $`\lambda ^{}`$ and $`f_d`$ arises naturally for low mass galaxies, where feedback is most efficient; Burkert (2000) showed that it can be an effect of disk growth happening inside-out; or disks might preferentially lose baryons with high angular momentum. On the other hand, Burkert, van den Bosch & Swaters (2002) argue that the correlation may be entirely due to intrinsic degeneracies due to the fact that from rotation curves analysis the total dark halo mass is not well constrained. In JVO02 (figure 14) we find that this is the case for rotation curves that are not measured at large enough radii, but, for most of the galaxies in the sample (70%), the high-quality rotation curves allow model degeneracies to be lifted. This indicates that the random error in recovered parameters is small. On the other hand, the systematic errors (especially due to the uncertainty in the dark matter profile) may be larger, possibly causing the correlations.

We argue here that the most significant feature of the left panels of fig. 1 is not the correlation itself, but the fact that spirals avoid the upper left portion of the $`\lambda ^{}`$$`f_d`$ plane and that this is in agreement with the Toomre instability criterion. In passing we also note that disk instability might make galaxies with disk parameters in the bottom left corner of the plot to drop out of the sample: for these parameter values discs might be too concentrated and form a spheroidal instead of a disc or galaxies might become bulge dominated.

The right panels show an interesting, although weak, correlation between $`M_{200}`$ and $`f_d`$: low mass spirals ($`M_{200}<10^{11}M_{}`$) tend to have $`f_d>0.1`$ since no galaxies are found in the lower left portion of the plot. There can be several different effects that module the dependence of $`f_d`$ on the halo mass. Linear theory evolution from inflation predicts a weak anti correlation between spin parameter and peak height (Heavens & Peacock, 1988), and hence mass of the halo. The expected anti correlation is, however, very weak. On the other hand, cooling is less efficient in massive halos ($`M_{200}>\mathrm{\hspace{0.33em}10}^{12}M_{}`$), thus $`f_d`$ is expected to decrease strongly with increasing $`M_{200}`$ (e.g., van den Bosch (2001)) at the high mass end. However, in small mass halos, reheating of cooled gas through feedback and photoionization processes are much more efficient than in massive systems, thus implying that $`f_d`$ should increase with $`M_{200}`$ (e.g., Benson et al. (2001b)) at the low mass end. The correlation we find from rotation curves analysis is in rough agreement with the high-mass end trend (see fig. 5 of van den Bosch (2001)). For his no feedback model, the baryon fraction decreases for increasing masses. For the feedback model this is only true for $`M_{200}`$ bigger than $`5\times 10^{11}`$ M. For lower masses this trend reverses and we do not find evidence for it from our analysis of rotation curves. This may be due to the fact that for these low masses, only those galaxies with highest $`f_d`$ are visible.

Here we show that Toomre instability criterion predictions are in agreement with the zone of avoidance for combinations of disk parameters that involve low $`M_{200}`$ and low $`f_d`$ or low $`\lambda `$.

### 3.2 Stability of discs

The void areas found in fig. 5 can be explained by resorting to the Toomre instability criterion: for $`Q>1`$ the gas disc is stable and thus star formation is strongly suppressed. Expressions for $`Q`$ and $`\mathrm{\Sigma }`$ for the NFW profile have been presented in Jimenez et al. (1997), here we derive the corresponding expressions for the pseudoisothermal profile.

Recall that in the pseudoisothermal profile the density $`\rho `$ as a function of the radius $`r`$ is $`\rho (r)=\rho _0[1+(r/rc)^2]^1`$, where $`\rho _0`$ denotes the central finite density. The disk surface density as a function of the radius is $`\mathrm{\Sigma }(r)=\mathrm{\Sigma }_0\mathrm{exp}(r/R_d)`$, where $`\mathrm{\Sigma }_0=f_dM_{200}/(2\pi R_d^2)`$. Appendix I of JVO02 gives a relation between $`\rho _0`$ and the concentration parameter $`c`$ and a fitting formula for $`R_d`$ as a function of the disk parameters: $`M_{200}`$, $`\rho _0`$, $`r_c`$ and $`\lambda ^{}`$.

The shaded areas in figure 5 represent the void regions predicted by the Toomre criterion assuming $`c_s=6`$ km s<sup>-1</sup> for all values of masses and radii. The dashed line shows how the avoidance region would change if $`c_s=2`$ km s<sup>-1</sup>. Such a low value is based on the findings of Corbelli & Salpeter (1995), who argue that when the quasar UV ionizing background is the only source of ionization of the HI gas, as should be the case for dark galaxies with no supernova, then the sound speed of the gas cannot be higher than 2 km s<sup>-1</sup>. Measured values of the sound speed in the Milky Way and in nearby galaxies (see e.g., Kennicutt (1989)) range between 3 and 10 km s<sup>-1</sup>. Note that using a low value for $`c_\mathrm{s}`$ is a conservative assumption, since for small values of the velocity dispersion the formation of dark galaxy is more difficult.

More specifically, for a grid in $`M_{200}`$, $`\lambda ^{}`$, $`f_d`$, we compute the surface density of the exponential disc $`\mathrm{\Sigma }`$ and whether $`Q>1`$. The shaded area in the left panels of fig. 5 are obtained assuming a fixed mass of $`5\times 10^{11}`$ M, approximately in the middle of the mass range considered. The shaded area in the right panels is obtained as follows: we assume a log-normal distribution of spin parameters with mean 0.042 and dispersion 0.5 (Bullock et al., 2001), thus 68% of galaxies are expected to have $`\lambda ^{}>\lambda _{68}^{}0.033`$. We then assume that a combination of $`f_d`$ and $`M_{200}`$ parameters correspond to a dark galaxy when for all $`\lambda ^{}>\lambda _{68}^{}`$ the disc is Toomre stable. The Toomre criterion seems to work better for the NFW profile than for the pseudoisothermal.

However, for the pseudoisothermal profile model there is no theoretical prediction for the expected $`\lambda ^{}`$ distribution and the dependence of the concentration parameter $`cR_{200}/r_c`$ on the halo mass. To produce the gray areas in fig. 5 we assume the same log-normal distribution for $`\lambda `$ as predicted from CDM simulations. This assumption should be valid since JVO02 find that the empirically recovered $`\lambda ^{}`$ distribution is well approximated by the CDM prediction. For setting the concentration parameter we consider that, when a galaxy rotation curve is fitted by the pseudoisothermal profile, $`r_c`$ obtained is about 1/10 of that obtained when the curve is fitted by the NFW profile (fig. 5 of JVO02). Thus, for a given mass, we compute the concentration parameter as expected in the NFW profile, and convert it to the pseudoisothermal one. The scatter around this relation is, however, quite large, thus introducing some uncertainty in the grey areas so obtained.

The Toomre criterion is a local instability criterion; global instabilities can in principle trigger star formation. For example discs in which the self gravity is dominant are likely to be unstable to the formation of a bar. Here we use the results of Efstathiou, Lake & Negroponte (1982) and Mo, Mao & White (1998) to conclude that the onset of these instabilities does not affect the grey areas.

We conclude that the Toomre criteria at least in the context of a NFW profile seems to predict correctly the region of the parameter space that galaxies tend to avoid (i.e., these galaxies should be dark).

## 4 Discussion and Conclusions

In this paper we propose that gas in a large fraction of low mass dark matter halos may form Toomre stable disks. Such halos would be stable to star formation and therefore remain dark. This may potentially explain the discrepancy between the predicted and observed number of dwarf satellites in the Local Group, as well as the deviation between the predicted and the observed faint end slope of the luminosity function. We show that model fits to rotation curves are consistent with this hypothesis: none of the observed galaxies lie in the region of parameter space forbidden by the Toomre stability criterion. Such Toomre stable disks may be the origin of damped Ly$`\alpha `$ systems seen at high redshift (Kauffmann, 1996; Jimenez et al., 1998; Mo et al., 1998).

Most semi-analytic models of galaxy formation achieve a reconciliation between the observed and predicted abundance of low luminosity galaxies by drastically decreasing $`f_d`$ for faint galaxies. At present, there is no evidence from rotation curve modelling that low circular velocity disks are dark matter dominated (which would be the case if $`f_d`$ were very small). Our model only requires a very gentle decline in $`f_d`$ toward low masses. It magnifies the effect of previously proposed mechanisms (e.g. supernovae feedback, suppression of accretion), since $`Qf_d^1`$ and low $`f_d`$ disks are more likely to be Toomre stable; such mechanisms can therefore be ’tuned down’ to lower levels. This may be relevant to the problem of trying to simultaneously fit the luminosity function and the Tully-Fisher relation. In our model, the mass-to-light ratio is much more stable as a function of halo mass (except in dark halos, where it is infinite). Indeed, by construction we fit both the Tully-Fisher relation and the luminosity function, since both were used to obtain the velocity function.

Jimenez et al. (1997) proposed that disks in high spin halos may not be low surface brightness galaxies but instead be completely dark. This was at variance with the proposal that a large fraction of the ’missing galaxies’ lie in low surface brightness galaxies which formed preferentially in high spin halos (Dalcanton, Spergel & Summers, 1997). Here, we have extended the Jimenez et al. (1997) model to a LCDM scenario and show that it is supported from an analysis of rotation curves. If it is true that disks in such high spin halos may not be low surface brightness but instead be completely dark, then there should be a cutoff in the surface brightness distribution of galaxies; as surveys probe increasingly low surface brightness galaxies the luminosity density will not continue to increase but instead will plateau.

It is important to consider independent indications for the existence of this dark galaxy population and devise possible observational tests. For example, Augusto & Wilkinson (2001) and Trentham, Möller & Ramirez-Ruiz (2001) have investigated the detectability of dark objects.

These dark galaxies have very low surface densities of HI (below 10<sup>20</sup> cm<sup>-2</sup> at one scale length of the disk), still high enough to be detected by current surveys (Zwaan, 2001) if the HI was in its neutral phase. Current observations that are sensitive enough to detect such amount of neutral HI have been carried out recently by Charlton et al. (2000); Zwaan & Briggs (2000); Zwaan (2001). In particular, Zwaan (2001) found no significant detection of HI, down to a detection limit of $`7\times 10^6`$ M, not associated with visible galaxies. Therefore, it seems not likely that large masses of neutral HI are harbored in dark galaxies. This should not come as a surprise, since, given the low surface densities, the hydrogen might be ionized by the extragalactic background. At densities below 10<sup>19</sup> cm<sup>-2</sup>, HI is not able to self-shield from the external radiation field and most of it will be in its ionized phase. Another possibility is that, due to pressure ram stripping from the hot gas of the hosting large galaxy, most of the HI gas could be removed from the dark galaxy (Quilis, Moore & Bower, 2000).

A more promising route to detect the presence of dark galaxies is through their gravitational lensing properties. Recently Metcalf & Zhao (2001); Dalal & Kochanek (2001); Keeton (2001); Bradac et al. (2001) have discussed the need of substructure to explain the relative fluxes of multiple lensed systems. They show that the flux ratios observed in each lens can only be explained by the presence of substructure within a large smooth halo. This population would be similar to our dark galaxies.

Dark galaxies with dark halo masses above $`10^{10}`$ M should be directly visible via strong gravitational lensing. They will produce splitting angles $``$$`<`$ 0.1 arc second, which is just below detection with current technology, but will be easily observable with future radio interferometers (e.g., VLBI)

Using the metallicity distributions of globular clusters, Cote, West & Marzke (2001) conclude that the mass spectrum of proto-galactic fragments for the galaxies in their sample has a power law with index $`2`$, indistinguishable from that predicted from N-body simulations. They argue that the missing satellite population must therefore belong to a class of dark galaxies, similar to the ones we consider here.

The biggest caveat in our model is the assumption of conservation of angular momentum. In principle angular momentum loss is possible due to torquing by substructure in the dark matter halo; indeed, it is invariably seen in SPH simulations of disk galaxy formation. This results in more compact disks which are much more susceptible to gravitational instability. However, both fits to the scale lengths of disks (Mo, Mao & White, 1998) and modeling of the rotation curves (Jimenez, Verde & Oh, 2002; van den Bosch, Burkert & Swaters, 2001) recovers spin parameters consistent with the log-normal distribution seen in dissipationless N-body simulations, indicating little or no loss of angular momentum. In addition, we have ignored the effect of mergers; as halos merge their disks are likely to be disrupted and transformed into spheroidal systems, which again are more likely to become self-gravitating. The merger history of halos is well handled in semi-analytic models of galaxy formation which employ Press-Schechter based Monte-Carlo merger trees. It would be very interesting to incorporate the additional physics of Toomre stability into pre-existing semi-analytic (or hydrodynamic) galaxy formation models to compare its importance against various other proposed schemes for suppressing star formation in low mass halos.

## acknowledgments

We thank David Spergel for insightful comments and stimulating discussions. LV is supported in part by NASA grant NAG5-7154. SPO is supported by NSF grant AST-0096023. LV and RJ thank the TAPIR group at Caltech for hospitality.


format: text

text:
# Can radiative correction cause large neutrino mixing?

## Abstract

We investigate whether radiative corrections can be responsible for the generation of large mixing in a pair of degenerate neutrino with same CP parity. We find that this mechanism is fine tuned and doesn’t work for arbitrary mixing at the high scale.

The data from the solar neutrino and the atmospheric neutrino experiments can be explained through flavor oscillations of massive neutrinos. The favoured solution of both these problems are a pair of neutrinos with very small mass squared difference and large mixing. The mass square difference required for the solar neutrino oscillation is less than $`10^5eV^2`$ while for the atmospheric neutrinos it is $`10^3eV^2`$. If neutrinos form the hot Dark Matter then the masses of all the neutrinos have to be of the order of 1eV. This, along with the very small mass square differences suggests that neutrino mass spectrum is almost degenerate. If neutrinos have Majorana mass, then one can have a pair of degenerate neutrinos with opposite CP parities. These are called Pseudo Dirac neutrinos. The mixing between such pairs can be shown to be maximal . The smallness of neutrino mass can be obtained through the seesaw mechanism in which neutrinos can have both Dirac and Majorana masses. Even in this mechanism Pseudo Dirac neutrinos with maximal mixing are possible for a suitable choice of the structure of right handed Majorana mass matrix. While these models can generate only one maximal mixing in 3 generation scenario, the solar neutrino problem and the atmospheric neutrino problem requires two pairs with maximal mixing. Another problem is that for degenerate Majorana neutrinos, large mixing may not be stable under small radiative corrections. A pair of Majorana neutrinos with opposite CP parities is stable but for a pair with same CP parity large mixing collapse to 0 due to radiative correction. A detailed analysis on the stability of the neutrino mass matrix and mixing for various type of mass hierarchy is done in . Some authors  suggested that the same instability can cause large mixing at the low scale from an arbitrary mixing at the high scale. This mechanism is similar to the MSW mechanism where the mixing angle gets enhanced at the resonance. As this requires the present scale of the universe to be at the ‘resonance’ as far as neutrino mixing is concerned, this is a fine tuned solution. This is because 0 mixing is a fixed point of the renormalisation group evolution of neutrino masses and mixing whereas maximal mixing is not . The fact that for degenerate mass structure, large mixing is unstable have been discussed by other authors  We investigate whether such a mechanism can be at work to generate the required large mixing considering the range of mass square differences from the solar and atmospheric neutrino experiments.

Let $`m_1`$ and $`m_2`$ be the mass eigenvalues of two neutrinos at a high scale $`\mathrm{\Lambda }`$. Let $`\theta `$ be the mixing angle at this scale. Then in the flavour basis the mass matrix will be given as

$`M_F(\mathrm{\Lambda })=\left(\begin{array}{cc}\hfill m_1\mathrm{cos}^2\theta +m_2\mathrm{sin}^2\theta & \hfill \frac{m_2m_1}{2}\mathrm{sin}(2\theta )\\ \hfill \frac{m_2m_1}{2}\mathrm{sin}(2\theta )& \hfill m_1\mathrm{sin}^2\theta +m_2\mathrm{cos}^2\theta \end{array}\right)`$

Due to radiative correction this mass matrix gets modified in the low scale $`M_Z`$ to

$`M_F(M_Z)`$ $`=`$ $`\left(\begin{array}{cc}1+\delta _\alpha & 0\\ 0& 1+\delta _\beta \end{array}\right)M_F(\mathrm{\Lambda })\left(\begin{array}{cc}1+\delta _\alpha & 0\\ 0& 1+\delta _\beta \end{array}\right)`$ (5)
$`=`$ $`\left(\begin{array}{cc}\hfill (1+2\delta _\alpha )(m_1\mathrm{cos}^2\theta +m_2\mathrm{sin}^2\theta )& \hfill (1+\delta _\alpha +\delta _\beta )\frac{m_2m_1}{2}\mathrm{sin}(2\theta )\\ \hfill (1+\delta _\alpha +\delta _\beta )\frac{m_2m_1}{2}\mathrm{sin}(2\theta )& \hfill (1+2\delta _\beta )(m_1\mathrm{sin}^2\theta +m_2\mathrm{cos}^2\theta )\end{array}\right)`$ (8)
$``$ $`\left(\begin{array}{cc}A& B\\ B& C\end{array}\right)`$ (11)
Here $`\delta _\alpha `$ and $`\delta _\beta `$ are the radiative corrections to the flavours $`\alpha `$ and $`\beta `$.

The mixing angle at the low scale $`\theta ^{}`$ in the limit of near degeneracy $`m_1m_2m`$ is given as

$`\mathrm{tan}2\theta ^{}`$ $`=`$ $`{\displaystyle \frac{2B}{CA}}`$ (12)
$`=`$ $`{\displaystyle \frac{(1+\delta _\alpha +\delta _\beta )(m_2m_1)\mathrm{sin}(2\theta )}{(m_2m_1)\mathrm{cos}(2\theta )+2(\delta _\beta \delta _\alpha )m}}`$ (13)
If

$$2(\delta _\beta \delta _\alpha )=\frac{m_1m_2}{m}\mathrm{cos}(2\theta )$$
(14)
then $`\theta ^{}45^o`$

Let $`ϵ=2(\delta _\beta \delta _\alpha )`$. Without loss of generality let $`m_1>m_2`$. Whatever mixing angle $`\theta `$ one takes at high scale $`\mathrm{\Lambda }`$, when one goes to low scale the mixing angle $`\theta ^{}`$ will be maximal if condition (14) is satisfied. We call this value of $`ϵ`$ given by eq.(14) as $`ϵ_r`$ ($`r`$ referring to resonance). However it is crucial to investigate whether this condition can be satisfied at an arbitrary mixing angle $`\theta `$ at the high scale. This is because the masses $`m_1`$ and $`m_2`$ in eq.(14) are the masses at the high scale $`\mathrm{\Lambda }`$. Though the masses themselves don’t have significant evolution through radiative correction the same cannot be said about the mass differences near degeneracy. To see this, consider the mass square difference at the low scale

$`\mathrm{\Delta }m^2(M_Z)=(C+A)\sqrt{(CA)^2+4B^2}`$

From eq.(12) this becomes

$`\mathrm{\Delta }m^2(M_Z)\mathrm{sin}(2\theta ^{})=2B(C+A)`$

From eq.(11), in the limit of near degeneracy we get

$$\mathrm{\Delta }m^2(M_Z)\mathrm{sin}(2\theta ^{})=\mathrm{\Delta }m^2(\mathrm{\Lambda })\mathrm{sin}(2\theta )$$
(15)
Here we have used the near degeneracy approximation in the following

$`{\displaystyle \frac{m_2m_1}{m}}={\displaystyle \frac{(m_2m_1)(m_2+m_1)}{m(m_2+m_1)}}{\displaystyle \frac{\mathrm{\Delta }m^2}{2m^2}}`$

This is as expected in a MSW like scenario, where due to significant change in mass square difference, the mixing angle undergoes drastic enhancement in the resonance region. Eq.(15) says that once $`\mathrm{\Delta }m^2(M_Z)`$ and $`\mathrm{sin}(2\theta ^{})`$ are determined from oscillation experiments, the mass square difference and the mixing angle at the high scale $`\mathrm{\Lambda }`$ are related. So if the mixing angle $`\theta `$ is small at the high scale then the mass square difference $`\mathrm{\Delta }m^2(\mathrm{\Lambda })`$ should be large. This would need a large $`ϵ`$ to satisfy the resonance condition (14). We will see this in detail in the numerical plots below.

When $`\theta ^{}`$ is maximal, we have from (15)

$$\mathrm{\Delta }m^2(M_Z)=\mathrm{\Delta }m^2(\mathrm{\Lambda })\mathrm{sin}(2\theta )$$
(16)
We see here that once the evolution of $`\mathrm{\Delta }m^2`$ is fixed from the scale $`\mathrm{\Lambda }`$ to the scale $`M_Z`$, the mixing angle at the high scale $`\mathrm{\Lambda }`$ is fixed and not arbitrary.

Eq.(14) is a relation between $`ϵ_r`$ and the mass difference and mixing angles at the high scale. It would be appropriate to get $`ϵ`$ and $`ϵ_r`$ in terms of the quantities at the low scale. Rearranging eq.(12) and using eq.(15) we get

$$ϵ=\frac{\mathrm{\Delta }m^2(M_Z)}{2m^2}\mathrm{sin}(2\theta ^{})\left(\frac{1}{\mathrm{tan}(2\theta ^{})}\frac{1}{\mathrm{tan}(2\theta )}\right)$$
(17)
Since the combination $`\mathrm{\Delta }m^2(M_Z)\mathrm{sin}(2\theta ^{})`$ does not change under renormalisation group evolution as evident from eq.(15), we see from eq.(17) that given an $`ϵ`$ if a small $`\theta `$ is enhanced to a large $`\theta ^{}`$ around $`45^o`$ then a large $`\theta `$ will end up at a $`\theta ^{}`$ around $`90^o`$ which is small mixing but with inverted hierarchy. This observation is consistent with the explicit evolution of the neutrino mixing in SM, MSSM and two Higgs doublet model done in ref.  and .

$`ϵ_r`$ is given by $`\theta ^{}=45^o`$ in eq.(17)

$$ϵ_r=\frac{\mathrm{\Delta }m^2(M_Z)}{2m^2}\left(\frac{1}{\mathrm{tan}(2\theta )}\right)$$
(18)
From eq.(17) and (18) we get the range of $`ϵ`$ that generates large mixing around the maximal. This is given as

$$\delta ϵ=ϵ_rϵ\frac{\mathrm{\Delta }m^2(M_Z)}{2m^2}\left(\frac{\mathrm{sin}(2\theta ^{})1}{\mathrm{tan}(2\theta )}\right)$$
(19)
where $`\mathrm{tan}(2\theta ^{})>>\mathrm{tan}(2\theta )`$ around $`\theta ^{}=\pi /4`$ and small $`\theta `$. So if we want the mixing at the low scale to be anything in the range of $`30^o`$ to $`45^o`$ the range of $`ϵ`$ can be obtained by putting $`\theta ^{}=30^o`$ in eq.(19). This will be about $`15\%`$ of $`ϵ_r`$. If we also consider the reversed hierarchy after maximal mixing, this range becomes $`30\%`$. With the current possibility of atmospheric neutrino mixing angle to be very near maximal $`(4050)^o`$, this range gets restricted to about $`4\%`$ around $`ϵ_r`$.

While it may be possible for this mechanism to enhance some mixing angle from high scale to low scale with a few $`\%`$ range of $`ϵ`$ around $`ϵ_r`$, another question we ask is, given an $`ϵ`$ what is the range of mixing angle at the high scale that can be enhanced at the low scale. Figure (1) shows a plot of $`ϵ`$ verses the high scale mixing angle $`\theta `$ at various values of the low scale mixing $`\theta ^{}`$ as obtained from eq.(17). We have taken $`\mathrm{\Delta }m^2(M_Z)=10^3eV^2`$ and $`m=1eV`$. So this is the case for atmospheric neutrino solution. It is clear from the figure that for very small $`\theta `$ ($`<0.5^o`$) the $`ϵ`$ needed for the enhancement to occur is as large as 0.1. As $`ϵ`$ is expected to be small, this mechanism cannot produce a large mixing from such small mixing at high scales. If $`ϵ0.01`$ only $`\theta `$ in a small range $`1.5^o2^o`$ can be enhanced to near maximal mixing ($`30^o`$ to $`45^o`$) at the low scale. But if $`\theta `$ is over $`2.5^o`$ then this mechanism becomes an overkill as the low scale mixing then again becomes small but with an inverted hierarchy. This is as expected in MSW mechanism when the neutrinos pass through the resonance and then go far away from it, the mixing is no longer maximal and the hierarchy gets inverted. If $`ϵ0.001`$, then only $`\theta `$ between $`10^o`$ to $`15^o`$ get enhanced. While for $`ϵ=0.0001`$ high scale mixing upto about $`25^o`$ don’t get any significant enhancement. If the same analysis is done for the solution of solar neutrinos then the corresponding $`ϵ`$ should be lowered by two orders of magnitude as $`\mathrm{\Delta }m^2(M_Z)`$ is $`10^5eV^2`$ instead of $`10^3eV^2`$ in eq.(18). So we see that there is no range of $`ϵ`$ that can enhance a large range of mixing angles at the high scale. Thus in order to generate maximal mixing at the low scale through radiative correction of the masses at the high scale, one has to start with specific mixing angle at the high scale depending upon the particular model one considers. Similar conclusion have been obtained in specific models  through renormalisation group evolution of neutrino mass matrix.

To summarise, we see that radiative correction cannot give a general enhancement of any mixing angles at the high scales to maximal mixing at the low scale. If it works due to fine tuning of parameters within a few percent for low mixing angles, it would kill the large mixing at high scales to small mixing at low scale. This is somewhat to be expected as large mixing between neutrinos of same CP parity is unstable near mass degeneracy . Radiative magnification doesn’t give the freedom to start with arbitrary mixing angle at high scale in building models of neutrino mass and mixing structures.


text:
# Symmetry breaking in the Hubbard model at weak coupling

## Abstract

The phase diagram of the Hubbard model is studied at weak coupling in two and three spatial dimensions. It is shown that the Néel temperature and the order parameter in $`d=3`$ are smaller than the Hartree-Fock predictions by a factor of $`q=0.2599`$. For $`d=2`$ we show that the self-consistent (sc) perturbation series bears no relevance to the behavior of the exact solution of the Hubbard model in the symmetry-broken phase. We also investigate an anisotropic model and show that the coupling between planes is essential for the validity of mean-field-type order parameters.

The Hubbard model is one of the most important and most prominent theoretical models in modern condensed matter physics. Originally introduced in order to describe magnetism in transition metals, the Hubbard model has been most intensively investigated since Anderson’s proposal that the model should capture the essential physics of cuprate superconductors. The Hubbard-Hamiltonian

$$=t\underset{\mathrm{𝐢𝐣}\sigma }{}c_{𝐢\sigma }^{}c_{𝐣\sigma }^{}+U\underset{𝐢}{}n_𝐢n_𝐢$$
(1)
describes itinerant electrons with spin $`\sigma `$ on a lattice with nearest neighbor sites $`\mathrm{𝐢𝐣}`$, interacting through short-ranged Coulomb repulsion $`U`$. The success of this relatively simple lattice model for mobile interacting electrons is based on its ability to explain a number of important phenomena in condensed matter physics. Among these are the (Mott-Hubbard) metal-insulator transition, antiferromagnetism, ferromagnetism, incommensurate phases, phase separation, and normal-state properties of high-$`T_c`$ materials.

In spite of the prominence of the Hubbard model in condensed matter theory and the deceptively simple two-parameter form of its Hamiltonian, comparatively little is known exactly (or even accurately) about its solution.

In $`d=1`$ the exact solution has been determined by the Bethe Ansatz technique. Even in the extreme weak coupling regime the ground state of the one-dimensional Hubbard model is nonperturbative, possessing a discontinuity at $`U=0^+`$ where the Mott-Hubbard gap opens.

Much has recently been learned about the solution in high spatial dimensions ($`d=\mathrm{}`$), but even here the exact analytical solution is unknown, so that one has to take recourse to perturbative or numerical methods. There is no doubt that the half-filled Hubbard model exhibits long-ranged antiferromagnetism in the high-dimensional limit, whereas the situation away from half-filling is less clear.

In $`d=2`$ much is known from exact diagonalization or Monte Carlo simulation of finite systems. In addition, much effort is spent on the analytical and numerical analysis of perturbative and nonperturbative results. While there seems to be now a consensus that at half-filling the ground state has long-range antiferromagnetic order, there are still controversies when the system is doped and – even at half-filling – concerning the question of whether there is or is not a precursor pseudogap.

In so far as numerical work on the three-dimensional Hubbard model has been carried out, the extrapolation of the results to the thermodynamic limit is clearly made difficult by the small linear system size. Much effort has been invested in the analysis of the phase diagram using perturbative techniques. In the strong coupling regime ($`U\mathrm{}`$), where the half-filled Hubbard model reduces to an effective Heisenberg model, high-temperature series expansions and $`1/S`$-expansions can reliably be used to estimate the Néel temperature and the ground state energy, respectively. However the analysis of the intermediate and weak coupling regime is much more difficult: For this regime a variety of approaches has been proposed, e.g., by Kakehashi, Logan, Cyrot,  and Daré. These approaches all yield estimates for the Néel temperature which reduce to the Hartree-Fock result at weak coupling. However, it is well-known from studies of the Hubbard model, based on the $`1/d`$-expansion or the local approximation, that the Hartree-Fock approximation overestimates Néel temperature and order parameter in $`d=3`$ by a factor of the order of four, even in the extreme weak-coupling limit. The precise value of this renormalization factor in $`d=3`$ is as yet unknown.

The goal of this paper is, first, to present exact results for the broken-symmetry phase of the three-dimensional Hubbard model at weak coupling. In particular, we present asymptotically exact formulas for the Néel temperature and order parameter. It will be seen that this asymptotic formula also yields a useful approximation formula for the Néel temperature for all $`UD`$, where $`D=12t`$ is the band width of a three-dimensional hypercubic lattice. Second, we address sc perturbation theory for the broken-symmetry phase in low dimensions ($`d=2`$) in the ground state. Our main result here is that sc perturbation theory breaks down altogether, and that low-order sc perturbation theory bears no relevance to the behavior of the exact solution of the Hubbard model in $`d=2`$. We conclude that the antiferromagnetic order parameter cannot have a mean-field form. This result has obvious relevance for theories of high-$`T_c`$ superconductivity, many of which are based on perturbation theory in strictly two-dimensional systems. Third, we investigate the order parameter of the anisotropic Hubbard model, where two-dimensional planes are weakly coupled in $`c`$-direction by a small hopping amplitude $`t_{}t_{}`$. We demonstrate that the sc perturbation series converges – however sluggishly – as long as the interplane hopping $`t_{}`$ is finite.

Since the behavior of the Hubbard model even in the weak-coupling limit is nonperturbative one has to apply self-consistent theories. There are several ways of imposing self-consistency. The method at fixed order parameter  and some of its results are described in Refs.  for the special case of the Hubbard model in high spatial dimensions ($`d\mathrm{}`$). Here we extend these investigations to all finite dimensions. As described in Refs.  the order parameter $`\mathrm{\Delta }`$ and the Néel temperature $`T_N`$ are determined by the roots of the optimization equation

$$\frac{df}{d\mathrm{\Delta }}=0,$$
(2)
where $`f`$ is the free energy density. It then follows that one has to determine the order parameter from the selfconsistency condition

$$\mathrm{\Delta }=2h_o_0^{\mathrm{}}𝑑ϵN_d(ϵ)\frac{\mathrm{tanh}\left[\frac{1}{2}\beta \eta (ϵ)\right]}{\eta (ϵ)}.$$
(3)
In Eq. (3) $`h_0`$ denotes the symmetry-breaking field, $`\eta (ϵ)=\text{sgn}(ϵ)\sqrt{ϵ^2+h_0^2}`$ the dispersion, and $`\beta =1/T`$ the inverse temperature. $`N_d(y)`$ being the density of states in $`d`$ dimensions. It is important to note that in each order $`h_0`$ and $`\mathrm{\Delta }`$ have to be determined self-consistently from the Hartree-Fock contribution and the fluctuations together to obtain systematic corrections to mean-field theory.

First we consider sc perturbation theory in the broken-symmetry phase in dimensions $`d3`$. Fortunately, although the evaluation of the various higher order diagrams in this approach at finite values of the interaction ($`U>0`$) is difficult, if possible at all, the results at weak coupling (i.e., for $`U0`$) are simple: For all $`d3`$ one finds that the exact value of the Néel temperature $`T_N`$ and the exact order parameter $`\mathrm{\Delta }`$, can be expressed in terms of their Hartree equivalents and a scaling factor $`q_d`$. As is well-known, the Hartree Néel temperature is exponentially small for $`U0`$:

$`T_N^He^{I_d1/UN_d(0)},`$

where $`I_d`$ can be expressed in terms of an integral,

$`I_d={\displaystyle _0^{\mathrm{}}}𝑑y{\displaystyle \frac{1}{y}}\left(\mathrm{tanh}y1+{\displaystyle \frac{N_d(y)}{N_d(0)}}\right)\mathrm{ln}2,`$

Now, we find that the exact expressions for $`T_N`$ and $`\mathrm{\Delta }(T)`$ differ from the mean-field results by a scaling factor $`q_d`$:

$$\begin{array}{ccc}\hfill T_N& & q_d^{}T_N^H\hfill \\ \hfill \mathrm{\Delta }(T)& & q_d^{}\mathrm{\Delta }^H(T/q_d^{})\hfill \end{array}(U0).$$
(4)
For the special case $`d=\mathrm{}`$ we know from Ref.  that $`q_{\mathrm{}}=\mathrm{exp}(\gamma )`$, where $`\gamma =\sqrt{2}\text{arcoth}(\sqrt{2})`$, so that $`q_{\mathrm{}}0.2875`$. The renormalization factor $`q_d^{}`$ for all $`d3`$ is determined by the second order diagrams only and can be conveniently written in the form of an exponential of a lattice sum:

$$q_d=e^{S_d},$$
(5)
where $`S_d`$ is given by

$`S_d={\displaystyle \frac{1}{t}}{\displaystyle \underset{|𝐣|\text{even}}{}}{\displaystyle \frac{F_𝐣^2}{F_\mathrm{𝟎}^3}}{\displaystyle _0^{\mathrm{}}}𝑑\tau \left[G_𝐣(\tau )\right]^2`$

and $`F_𝐣`$ and $`G_𝐣(\tau )`$ are integrals over Bessel functions:

$`F_𝐣`$ $`=`$ $`{\displaystyle \frac{1}{2\pi t}}{\displaystyle _0^{\mathrm{}}}𝑑z{\displaystyle \underset{\mathrm{}=1}{\overset{d}{}}}\left[J_{|j_{\mathrm{}}|}(z)\right]`$
$`G_𝐣(\tau )`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle _0^{\mathrm{}}}𝑑z{\displaystyle \frac{\tau }{\tau ^2+z^2}}{\displaystyle \underset{\mathrm{}=1}{\overset{d}{}}}\left[J_{|j_{\mathrm{}}|}(z)\right].`$
This expression is formally exact for all $`d>2`$. We have numerically evaluated the renormalization factor $`q_d`$ for $`d=3`$ and find the following result:

$`q_3=0.2599.`$

Thus, the exact Néel temperature and the exact order parameter are smaller than the Hartree-Fock predictions by almost a factor of four.

Several remarks are in order: ($`i`$) The term labeled by $`𝐣`$ in the lattice sum $`S_d`$ stems from second order diagrams containing two Hubbard vertices a distance $`|𝐣|`$ apart. ($`ii`$) Since we know from previous calculations in high dimensions that sites at a relative distance $`|𝐣|`$ give a contribution of order $`d^{|𝐣|}`$ to the free energy, one expects that the lattice sum $`S_d`$ converges extremely rapidly if $`d`$ is large. ($`iii`$) Keeping only the $`𝐣=\mathrm{𝟎}`$ term in $`S_d`$ corresponds to the local approximation of Ref. . In $`d=3`$ we find that $`q_3^{\text{loc}}=0.2673`$. Comparison of $`q_3^{\text{loc}}0.2673`$ with the exact result $`q_3^{}0.2599`$ shows that the local approximation works very well even in $`d=3`$. ($`iv`$) Note that odd lattice sites (i.e., interactions between the two sublattices) do not contribute to $`q_d^{}`$ to leading order in $`U`$ for $`U0`$. This is due to a different symmetry of the Green functions for odd and even values of $`|𝐣|`$. ($`v`$) Also note that our result (5) for the renormalization factor breaks down for $`d=2`$ since in this case the integrals $`F_𝐣`$ diverge logarithmically for all even $`|𝐣|`$.

The asymptotic result for the Néel temperature in Eq. (4) has a relevance beyond the pure weak coupling limit. This is particularly clear in infinite dimensions, where accurate quantum Monte Carlo data are available. A comparison of the Monte Carlo results with the appropriate asymptotic formula for $`d=\mathrm{}`$  reveals excellent agreement for all $`UD`$. Therefore one expects that (4) yields an equally good approximation for $`UD`$ in $`d=3`$. The phase diagram of the three-dimensional Hubbard model was calulated with a quantum Monte Carlo technique by Hirsch and Scalettar et al. A comparison of the Monte Carlo data from Ref.  with the expected behavior (4) reveals a significant discrepancy, suggesting (as was also pointed out in Refs. ) that the existing Monte Carlo estimates for $`T_N`$ are too high. Improved Monte Carlo simulations show a significantly reduced Néel temperature. Yet it is still more than 15% higher than the asymptotic value (see Fig. 1). Our weak-coupling approximation formula (4) and the strong-coupling approximation formula $`T_N^{\text{sc}}=3.83t^2/U`$ from Ref.  can hence serve as benchmarks for future Monte Carlo simulations.

Since the upper critical dimension of the Hubbard model is presumably $`d_u=4`$, one expects that the critical behavior of the order parameter in $`d=3`$ is characterized by a nontrivial critical exponent for all $`U>0`$. Yet the exact asymptotic behavior of $`\mathrm{\Delta }(T)`$ in (4) displays mean-field behavior near $`T_N`$, with a critical exponent $`\beta =1/2`$. These two observations can easily be reconciled by calculating the size $`\mathrm{\Delta }T`$ of the Ginzburg region in $`d=3`$, which requires as input the correlation length $`\xi `$ and the jump in the specific heat $`\mathrm{\Delta }C`$ at $`T_N`$. An explicit calculation shows that $`\xi `$ is exponentially large at weak coupling, $`\xi \left[T_N(T_NT)\right]^{1/2}`$, while $`\mathrm{\Delta }C`$ is exponentially small, $`\mathrm{\Delta }CT_N`$. Combination gives for the relative size of the Ginzburg region: $`\mathrm{\Delta }T/T_NT_N^4`$, which is exceedingly small at weak-coupling and vanishes for $`U0`$.

Now we address sc perturbation theory for the broken-symmetry phase in $`d=2`$. Since $`T_N=0`$ in $`d=2`$, we focus on the renormalization of the ground state order parameter: $`\mathrm{\Delta }=q_2\mathrm{\Delta }^H`$. Calculating only the second order contribution to Eq. (2) we find $`df_2/d\mathrm{\Delta }h_0x^2/6`$ in the limit $`U0`$ with $`x=\frac{U}{4\pi ^2t}\left[\mathrm{ln}\left(\frac{4t}{h_0}\right)\right]^2`$. The solution of the self-consistency equation (3) shows an interaction-dependent renormalization:

$$\mathrm{\Delta }=\frac{8t}{U}\left(3\sqrt{3}\right)\mathrm{exp}\left(\sqrt{\frac{4\pi ^2t(3\sqrt{3})}{U}}\right).$$
(6)
In contrast to the result in $`d3`$, where $`q_d`$ is constant, the renormalization factor in $`d=2`$ vanishes exponentially as $`U0`$. This result already indicates that the sc perturbation series in two dimensions has no small expansion parameter at weak coupling. Hence higher order fluctuation terms are important to decide whether the mean-field prediction for the form of the order parameter is at least qualitatively correct or does totally fail in $`d2`$ in the symmetry-broken phase.

To gain some insight into sc perturbation theory in higher orders we calculated bubble and ladder diagrams. We find

$$\frac{df_{\text{BL}}}{d\mathrm{\Delta }}h_0F_{\text{BL}}(x)$$
(7)
with the function

$`F_{\text{BL}}(x)`$ $`=`$ $`6+x{\displaystyle \frac{x^2}{3}}{\displaystyle \frac{1}{2}}\mathrm{ln}\left[(1+x)^{9\left(1+\frac{1}{x}\right)}(1x)^{\left(1\frac{1}{x}\right)}\right]`$
for $`0<x<1`$. The important result is that the optimization condition Eq. (2) has no roots for $`0<x<1`$, apart from the high temperature solution $`\mathrm{\Delta }=0`$.

We comment on the results. ($`i`$) Our analysis shows that Hartree-Fock theory bears no relevance to the behavior of the Hubbard model in $`d=2`$ at half-filling. The mean-field result is completely destroyed by quantum fluctuations. Hence the antiferromagnetic order parameter must have a completely different form: $`\mathrm{ln}(4t/h_0)\sqrt{t/U}`$. ($`ii`$) A calculation of bubble and ladder diagrams in $`d=\mathrm{}`$ shows that all contributions $`f_n`$ for $`n>2`$ yield small corrections to $`q`$ of the form $`q=q_{\mathrm{}}\mathrm{exp}\left[2N_{\mathrm{}}(0)\gamma ^2U+𝒪(U^2)\right]`$, where $`\gamma `$ has been defined above. ($`iii`$) Concerning the behavior of sc perturbation theory we conclude that only in dimensions $`d>2`$ low-order sc perturbation theory is capable of describing the physics of the Hubbard model correctly. In low dimensions sc perturbation theory in the symmetry broken phase diverges. ($`i\nu `$) Remarkably, in $`d=1`$ the structure of sc perturbation theory is exactly the same as in $`d=2`$, now with $`x=\frac{U}{2\pi t}\mathrm{ln}\left(\frac{2t}{h_0}\right)`$ in Eq. (7).

The observation that sc second order perturbation theory gives even quantitatively exact results in $`d=3`$ and diverges in $`d=2`$ indicates that between two and three spatial dimensions must be a transition where sc perturbation theory becomes inadequate. We introduced a weak hopping amplitude $`t_{}`$ in $`c`$-direction. The anisotropy is then given by the dimensionless quantity $`\lambda =t_{}/t_{}1`$ where $`t_{}`$ denotes the hopping in the planes. In this anisotropic three-dimensional model the renormalization factor can again be expressed as a lattice sum according to Eq. (5). However, the smaller $`\lambda `$ is the more sites have to be summed over. In the $`j_1,j_2`$-plane all points within a region $`1/\lambda ^2`$ around the origin contribute with $`\mathrm{ln}(1/\lambda )`$ to the renormalization factor, yielding

$`q_\lambda =\mathrm{exp}\left({\displaystyle \frac{\text{const.}}{\lambda ^2\mathrm{ln}(1/\lambda )}}\right).`$

The renormalization factor is thus strongly reduced in less than three dimensions and zero in $`d=2`$. This result demonstrates that the coupling to the third dimension is essential for mean-field theories to be qualitatively valid.

In summary, we have investigated the half-filled Hubbard model in two and three dimensions at weak coupling calculating systematic corrections to Hartree-Fock mean-field theory. We have shown that the mean-field antiferromagnetic Hartree-Fock solution does not yield exact results in the weak coupling regime. In particular, we showed that quantum fluctuation induced corrections make sc perturbation theory diverge in $`d2`$. In $`d>2`$ the exact order parameter and Néel temperture differ from the mean-field result by a scaling factor.

We acknowledge helpful discussions with F. Gebhard, M. Jarrell, A. Muramatsu, and T.M. Rice.


text:
# Scaling of nano-Schottky-diodes

## Abstract

A generally applicable model is presented to describe the potential barrier shape in ultra small Schottky diodes. It is shown that for diodes smaller than a characteristic length $`l_\mathrm{c}`$ (associated with the semiconductor doping level) the conventional description no longer holds. For such small diodes the Schottky barrier thickness decreases with decreasing diode size. As a consequence, the resistance of the diode is strongly reduced, due to enhanced tunneling. Without the necessity of assuming a reduced (non-bulk) Schottky barrier height, this effect provides an explanation for several experimental observations of enhanced conduction in small Schottky diodes.

The effect of downscaling the dimensions of a device on its electrical transport properties is an important topic today. Extremely small diodes have been experimentally realized and characterized in various systems, e.g. carbon nanotube heterojunctions yao99 , junctions between $`p`$-type and $`n`$-type Si nanowires cui01  or junctions between the metallic tip of a scanning tunneling microscope and a semiconductor surface avouris93 ; hasunuma98 . These experiments showed several deviations from conventional diode behavior. Despite some modelling in truly one-dimensional systems leonard99 ; odintsov00 , little work has been done on modelling the effects of downscaling a conventional diode, in the regime where quantum confinement does not play a role.

In this paper we present a simple model (based on the Poisson equation) describing the barrier shape in a diode, that is readily applicable to arbitrarily shaped small junctions. It is related to descriptions of inhomogeneities in the Schottky barrier height (SBH) in large diodes tung92 , barrier shapes in small semiconducting grains malagu02  and charge transfer to supported metal particles ioannides96 . Although we restrict ourselves to metal-semiconductor junctions, the model can easily be adapted for e.g. $`p`$-$`n`$-junctions. The main result is that if the size of the metal-semiconductor interface is smaller than a characteristic length $`l_\mathrm{c}`$, the thickness of the barrier is no longer determined by the doping level or the free carrier concentration, but instead by the size and shape of the diode. The resulting thin barrier in small diodes will give rise to enhanced tunneling, qualitatively explaining measurements of enhanced conduction avouris93 ; hasunuma98 ; smit02 , without the necessity of assuming a reduced SBH. Moreover, experimentally observed scaling behavior and deviating $`IV`$-curve shapes smit02  can be explained.

The transport properties of a Schottky diode are governed by the potential landscape which has to be traversed by the charge carriers. First, we study an easily scalable and highly symmetrical model system, namely a metallic sphere embedded in semiconductor (see Figure 1, upper left inset). The radius $`a`$ of the metallic sphere is a measure for the interface size: for large $`a`$, we expect to find the well known results for a conventional diode, while decreasing $`a`$ gives the opportunity to study finite size effects.

We only model the barrier shape in the semiconductor; the SBH $`\phi _\mathrm{B}`$ is accounted for in boundary conditions and is considered as a given quantity. For simplicity, the depletion approximation sze81  is adopted, which is valid for a wide range of realistic parameters. Moreover, the space charge region is assumed to be homogeneously charged, an assumption that will be discussed later. Solving the Poisson equation in $`n`$-type silicon with the boundary condition that the charge on the sphere cancels the total charge in the space charge region, we find for $`0xw`$

$$\begin{array}{c}\frac{e}{kT}V(x)=\hfill \\ \\ \frac{1}{2L_\mathrm{D}^2}\left[(a+w)^2\frac{2(a+w)^3}{3(a+x)}\frac{(a+x)^2}{3}\right],\hfill \end{array}$$
(1)
where $`x`$ is the radial distance from the interface, $`w`$ the depletion width and $`L_\mathrm{D}=\sqrt{\epsilon _\mathrm{s}kT/(e^2N_\mathrm{d})}`$ the Debye length. The zero-point of the potential is chosen in the semiconductor bulk. The value of $`w`$ is fixed by the second boundary condition $`V(0)=V_\mathrm{s}`$, where $`V_\mathrm{s}`$ is the total potential drop over the space charge region and satisfies $`V_\mathrm{s}=(\phi _\mathrm{B}\phi _\mathrm{s})/eV`$ (with $`\phi _\mathrm{s}=E_\mathrm{c}E_\mathrm{f}`$). Eq. (1) is valid for small bias voltage $`V`$. The limited validity of the depletion approximation at finite temperatures only affects the tail of the barrier (where $`|V(x)|kT`$), which is unimportant for the transport properties. From the equation, it can be seen that the characteristic length scale of this system is

$$l_\mathrm{c}\stackrel{\mathrm{def}}{=}L_\mathrm{D}\sqrt{2eV_\mathrm{s}/kT}=\sqrt{\frac{2\epsilon _\mathrm{s}V_\mathrm{s}}{eN_\mathrm{d}}}.$$

By comparing the diode size $`a`$ to $`l_\mathrm{c}`$ we can decide whether the diode is ‘small’ or ‘large’. In the lower right inset of Figure 1 the value of $`l_\mathrm{c}`$ is plotted versus doping concentration $`N_\mathrm{d}`$.

An important quantity for electrical transport is the Schottky barrier thickness. In Figure 1, the barrier full width at half maximum (FWHM, $`x_{1/2}`$) calculated from Eq. (1) is plotted as a function of diode size $`a`$. From the figure it is clear that for $`al_\mathrm{c}`$ the value of $`x_{1/2}`$ approaches a constant, which was expected for a large diode. Indeed, for $`al_\mathrm{c}`$, Eq. (1) reduces to $`V(x)=\frac{eN_\mathrm{d}}{2\epsilon _\mathrm{s}}(xw)^2`$, which is the well-known textbook sze81  result for band bending in the depletion approximation for an infinitely large diode. Both the depletion width $`w=\sqrt{(2\epsilon _\mathrm{s}/eN_\mathrm{d})V_\mathrm{s}}`$ and $`x_{1/2}`$ are in that regime independent of $`a`$.

Figure 1 shows that for $`al_\mathrm{c}`$ the value of $`x_{1/2}`$ is no longer constant, but decreases with decreasing $`a`$. For $`al_\mathrm{c}`$ it approaches $`x_{1/2}=a`$, i.e. the barrier thickness equals the diode size. This also follows from Eq. (1), which reduces to $`V(x)=V_\mathrm{s}a/(a+x)`$ for $`al_\mathrm{c}`$ and $`xw`$ (that is, close to the interface). Note that this is exactly the potential due to the charged sphere only. In this regime, the effect of the semiconductor space charge on the barrier shape and thickness can be neglected. This can be understood from the fact that the screening due to the space charge region takes place on a length scale $`l_\mathrm{c}`$, as in conventional (large) diodes. However, from Gauss’s law it follows that any charged object of typical size $`d<\mathrm{}`$ in a dielectric medium gives rise to a potential that behaves roughly as $`V(r)d/r`$. This Coulomb potential can be further screened by the formation of a space charge layer of opposite sign, but that additional screening can be neglected if $`dl_\mathrm{c}`$. This observation holds for *any* interface with typical dimensions much smaller than $`l_\mathrm{c}`$.

In a geometry that can actually be fabricated, the Poisson equation must be solved numerically. We have done this for $`n`$-doped silicon ($`N_\mathrm{d}=10^{15}\mathrm{cm}^3`$) in contact with metallic circular disks of various radii. In all further calculations $`\phi _\mathrm{B}=0.67`$eV was used, which is the barrier height of the CoSi2/Si(111)-interface tung92b . Figure 2 shows the FWHM-contours of the barriers as resulting from these calculations. Also shown are the FWHM-contours of the barrier due to the metallic contacts only, illustrating the negligible effect of the space charge region on the barrier thickness in very small diodes footnote1 .

To study the effect of the reduced barrier width on the transport properties of a small Schottky diode, a transmission coefficient $`T(E,V)`$ was obtained for the barrier shape from Eq. (1). This was done in a one-dimensional fully quantum mechanical calculation walker94 . Note that $`T(E,V)`$ is implicitly dependent on temperature and doping level, because these quantities influence the position of the Fermi-level in the bulk semiconductor. The current density is then given by

$$J(V)_0^{\mathrm{}}T(E,V)[f(\phi _s+E)f(\phi _s+E+V)]dE,$$

from which it follows that the zero bias differential conductance satisfies

$$\frac{\mathrm{d}J}{\mathrm{d}V}|_{V=0}_0^{\mathrm{}}T(E,V)f^{}(\phi _s+E)dE.$$

Here, $`f`$ is the Fermi-Dirac distribution function and $`E`$ the energy above the semiconductor conduction band edge. Transport due to electrons at energies below the barrier maximum ($`E<V_\mathrm{s}`$) is regarded as tunneling, while for $`E>V_\mathrm{s}`$ we speak of thermionic emission. Obviously, the contribution of thermionic emission is almost independent of the barrier thickness, while tunneling is strongly dependent on the barrier thickness.

In Figure 3, the calculated zero bias differential conduction is plotted as a function of diode size $`a`$ for several values of $`N_\mathrm{d}`$. For $`al_\mathrm{c}`$ this quantity is independent of $`a`$. For smaller values, the tunnel current starts to increase rapidly, eventually leading to a strong increase of the total conduction.

Moreover, the shape of the $`IV`$-curves changes with decreasing diode size. Our calculations (Fig. 4) show that for large diodes the $`IV`$-curve has exactly its expected exponential shape ($`I[\mathrm{exp}(eV/kT)1]`$). Apart from the total current increase, in small diodes the relative contribution of the reverse current starts to increase and eventually—in extremely small diodes—the reverse current exceeds the forward current, thus reversing the rectifying behavior of the diode.

Note that the image charge effect sze81  was neglected so far. However, inclusion of this effect would only enhance the phenomenon mentioned above, as it reduces the effective barrier height and width even further, especially in barriers which are narrow already.

One more issue that needs to be discussed is that of discrete random dopants. In our analysis, the dopants played a role in determining the Fermi-level position in the semiconductor bulk and were considered to provide a homogeneous space charge region. However, for the realistic parameters $`N_\mathrm{d}=10^{15}`$cm-3 and $`a=30`$nm (so $`al_\mathrm{c}`$) the volume in which the potential drops to half its initial value contains approximately one doping atom. Discrete energy levels of such a doping atom cannot be resolved at room temperature. More importantly, the potential well due to an ionized single dopant will locally distort the barrier shape. This effect complicates the potential landscape, but it can only significantly increase the conduction of the diode, when the dopant resides close to the interface footnote2 .

In conclusion, we have shown by means of a simple electrostatic argument that the Schottky barrier thickness becomes a function of the diode size for small diodes (e.g. smaller than $`l_\mathrm{c}80`$nm for $`N_\mathrm{d}=10^{17}`$cm-3). Consequently, the contribution of tunneling to the total conductance is greatly enhanced in small diodes. This effect explains several experimental results avouris93 ; hasunuma98 , without the assumption of a reduced SBH. Moreover, small diodes show $`IV`$-curve shapes that qualitatively differ from those of conventional diodes.

We wish to thank J. Caro and H.W.M. Salemink for detailed discussions concerning this work. This work is part of the research programme of the ’Stichting voor Fundamenteel Onderzoek der Materie (FOM)’, which is financially supported by the ’Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO)’. One of us, S.R., wishes to acknowledge fellowship support from the Royal Netherlands Academy of Arts and Sciences.