format: text

text:
# Chemical and spectrophotometric evolution of Low Surface Brightness galaxies

## 1 Introduction

Galaxies with a central disc surface brightness (in the $`B`$ band) well below the Freeman value of $`\mu _{B,0}`$=21.65 mag arcsec$`{}_{}{}^{2})`$ \- typical of the average previously catalogued “classical” High Surface Brightness galaxies or HSBs - are classified as LSBs. In the last decade, in particular, a considerable body of observational data has been collected on those objects (e.g. de Blok, van der Hulst & Bothun 1995 ; O’Neil, Bothun & Cornell 1997a ; O’Neil et al. 1997b ; Beijersbergen, de Blok & van der Hulst 1999 ; O’Neil, Bothun & Schombert 2000a ; van den Hoek et al. 2000 ; Matthews, van Driel & Monnier Ragaigne 2001; Galaz et al. 2002). Although there is no unambiguous definition of an LSB galaxy, in the following we will adopt as the limit between HSBs and LSBs a central disc surface brightness value of $`\mu _{B,0}`$=22 mag arcsec-2, a commonly used criterion.

LSB galaxies may turn out to be crucial for studies of galaxy formation and evolution and of the ‘cosmic’ chemical evolution if they constitute the majority of galaxies (as suggested by O’Neil & Bothun 2000). They are also of interest in the field of high redshift quasar absorbers: some Damped Lyman Alpha systems (DLAs) are identified with LSBs, and DLAs could be due as much to LSBs as to ‘classical’ spirals (e.g. Boissier, Péroux & Pettini 2003 and references therein).

Observations have shown a large variety in the properties of of LSB galaxies, which range from dwarfs to massive systems (e.g. Bothun, Impey & McGaugh 1997) and from the “blue” to the “red” edge of galactic colours (O’Neil et al. 2000a). This variety suggests that they form a heterogeneous family, as argued by Bell et al. (2000).

Nevertheless, their potentially interesting implications for galactic astrophysics motivated various theoretical studies, aiming towards an understanding of their nature. For instance, Gerritsen & de Blok (1999) used N-body simulations to study their star formation, while van den Hoek et al. (2000) compiled observational characteristics of 24 LSBs and tried to derive their star formation histories by comparison with simplified chemical evolution models. On the other hand, Jimenez et al. (1998) modelled the chemical and spectro-photometric evolution of LSBs by assuming that they are discs with larger angular momentum than HSBs; this idea had been suggested through the study of the formation of disc galaxies by Dalcanton, Spergel & Summers (1997), who considered a gravitationally self-consistent model for the formation of both LSB and HSB discs.

In a recent series of papers (Boissier & Prantzos 2000, Prantzos and Boisier 2000, Boissier and Prantzos 2001, Boissier et al. 2001) the chemical and spectrophotometric evolution of spiral galaxies of various masses and spin parameters was computed and compared quite successfully to a large number of observational data: Tully-Fisher relation, colours, integrated spectra, star formation efficiency, gas fraction etc. Those works showed that for relatively large values of the spin parameter $`\lambda `$ the properties of discs ressembled closely those of LSBs (in particular the central surface brightness).

Following the ideas of Dalcanton et al. (1997) and Jimenez et al. (1998) in the present paper we extend the same models to larger values of the spin parameter, keeping the remaining galactic physics (stellar Initial Mass Function, prescriptions for Star Formation Rate and Infall rate etc.) the same. We compare then our results to an extended set of observed properties of LSBs (H i mass, colours, abundances). Our main motivation is to answer the question: can HSB spirals and LSBs be modelled in the same framework, where LSBs merely have a larger angular momentum, or is it necessary to invoke specific assumptions in order to explain their properties?

In section 2 we present the compilation of data on LSBs properties that we adopted for our study. In section 3 we explain how the simplest (1st order) models of LSBs are constructed on the basis of our - quite successful - HSBs models, by merely increasing the spin parameter. These simple models are then compared to the observations in section 4. In section 5 bursts and truncations are added to the smooth star formation history of the “simple models”, in order to obtain better fits to the ensemble of the observational data. Section 6 gives a summary of the obtained results.

## 2 Low Surface Brightness galaxies: the data sample

Low Surface Brightness galaxies (LSBs) display a very broad range of sizes, masses and spectral line widths. Their typical observed 21 cm H i line width is about 100 km/s, though it can be as high as 600 km/s. “Giant” LSBs have also very large scalelengths, the most striking case being Malin-1. In some sense, their large velocities and sizes make them similar to HSB spirals, as argued in Bothun et al. (1997).

In this paper, we explore possible links in the physical properties/histories between HSB and LSB discs. Since we have no a priori idea of which of the LSBs could actually be linked to HSBs, we consider observations for the whole range of LSBs, from dwarfs to giants. Our sample consists of the following data sets, available in the literature.

O’Neil et al. LSBs sample

The observational data on LSBs from O’Neil et al. (1997a,b) provides information on the surface brightness, scalelength, absolute magnitude of the disc and colours for 127 galaxies. 80% of them are well fitted by an exponential profile, which justifies the use of a disc model for studying LSBs. The data suggest that a large fraction of LSBs are redder than previously thought. H i data for some of these galaxies are given in O’Neil, Bothun, & Schombert (2000) and Chung et al. (2002). Note that some of these objects are members of galaxy clusters (Pegasus and Cancer).

de Blok et al. LSBs sample

For 21 late-type field LSBs from the sample of Schombert & Bothun (1988) and the UGC, de Blok et al. (1995) published disc scalelengths and absolute magnitudes and colours; for subsets of that sample H i data are given in de Blok et al. (1996), abundances in de Blok & van der Hulst (1998a), CO line upper limits in de Blok & van der Hulst (1998b) and H$`\alpha `$ rotation curves in McGaugh, Rubin & de Blok (2001). The gas content of those galaxies is also available from the work of van den Hoek et al. (2000).

Bulge-dominated LSBs sample

Although most of the LSBs seem to be late-type objects without any apparent bulge, some of them clearly do have bulges, and we do not exclude these from our analysis, to avoid introducing a bias against the most evolved LSB discs. For this reason, we also use the data on surface brightness, scalelength, and colours of 20 galaxies from Beijersbergen et al. (1999). As our models concern only the disc component of LSBs, we used the “area-weighted” colours given by Beijersbergen et al., which should be representative of the discs.

Giant LSBs sample

We also use the data of Matthews et al. (2001), which provides surface photometry and H i data for 16 giant LSB galaxies, with high luminosities ($``$1010 L) and large scalelengths ($`>6`$ kpc). Though LSB Giants are relatively rare, they may play an important role in our investigation of the link between HSBs and LSBs because of their similar sizes/velocities (note that there are many more HSBs with large velocity widths than with small ones, like for LSBs).

Infra-red selected LSBs sample

A new sample of LSBs selected from the 2MASS near-infrared survey (Monnier Ragaigne et al., 2002, MR2002) is used as well. Though these 4,000 objects were selected on their low $`K_s`$-band central disc surface brightness ($`\mu _{K,0}`$$`<`$18 mag arcsec-2) they turned out to be bluer than expected; the $`B`$-band central disc surface brightnesses of the 20 objects for which we have optical surface photometry actually lie in the range of $``$21-22.5 mag arcsec-2, straddling the 22 mag arcsec-2 $`B`$-band selection criterion for LSBs we adopted in the present paper. Here we use the subset of 20 galaxies for which we presently have available $`BVRI`$ surface photometry, H i spectra and $`JHK`$ photometry. This sample allows the exploration of the link between optically HSB and LSB galaxies.

Homogenisation of the data

Because the data sets have different origins, they need homogenisation before a proper comparison with our disc models. We adopted the following procedure:

$``$ All distance-dependent quantities are reduced to distances based on a Hubble constant of 65 km/s/Mpc.

$``$ The models presented in the next section were computed under the hypothesis that the discs are observed face-on, and the data were corrected accordingly to an inclination angle i=0. This affects magnitudes and surface brightness by a quantity equal to $`2.5C\mathrm{log}(\mathrm{cos}i)`$, where $`C`$ ranges from 0 for an optically thick disc to 1 for an optically thin one. Since LSBs are usually not much evolved chemically (they have low metallicities), it is likely that they are dust-free and we assume $`C`$=1. This assumption is justified by the fact that the colours of LSBs seem to be independent of inclination (e.g. O’Neil et al. 1997a).

$``$ The absolute total disc magnitude used in this paper, $`M_{B,d}`$, is deduced from the fit with an exponential profile (e.g. O’Neil et al. 1997a), i.e. $`M_{B,d}=\mu _02.5\mathrm{log}(2\pi R_d^2)`$, where $`\mu _0`$ and $`R_d`$ are the central surface brightness and the scalelength, in mag arcsec-2 and arcsec, respectively.

$``$ The total gas mass is obtained by correcting the H i mass for the helium fraction ($`M_G`$/$`M_{HI}`$=1.4). The molecular component in LSBs is generally negligible, though a number of LSBs have now been detected in CO (Matthews & Gao 2001; O’Neil, Hofner & Schinnerer 2000, and references therein).

## 3 Models: from “classical” HSB spirals to LSBs

Dalcanton et al. (1997) suggested that LSBs could be the large angular momentum equivalents of “classical” spirals. That idea was adopted in Jimenez et al. (1998) who constructed models relatively similar to those used in the present study; they showed that this assumption can lead to reasonably good results, at least for the limited observational sample they used (which did not include objects like the red LSBs of the O’Neil et al. sample).

In this paper, we have adopted the same assumption, i.e. that LSBs can be explained with models similar to those of HSB spirals, but with a larger angular momentum. The strength of this work is that it relies on detailed models of the chemical and spectro-photometric evolution of disc galaxies, already shown to be successful in reproducing the main properties of the Milky Way and of nearby HSB spirals: see Boissier & Prantzos 1999 (BP99) for the Milky Way model, Boissier & Prantzos 2000 (BP2000) for the extension to other spirals and Boissier et al. (2001) on the agreement of the models with the observed gas fraction, star formation efficiency and ages of nearby spirals.

Notice that Schombert, McGaugh & Eder (2001) determined the gas fraction of LSB dwarf galaxies, and argued that their data were compatible with the BP2000 models of discs with the largest angular momenta. As an alternative possibility, they also suggested that LSBs could be on average younger than HSBs. That conclusion is also in agreement with our models, since large angular momentum galaxies (LSBs) form the bulk of their stars later, on average, and thus they appear younger than HSBs. We shall note that their galaxies are dwarf-LSBs, corresponding to the lower masses found in the samples we collected (section 2).

Since the LSBs models presented here are a simple extrapolation of those presented previously in this series, we will recall here only briefly their main features (section 3.1); we discuss, in particular, the role of the spin parameter and of its distribution (section 3.2), since we assume that it is at the origin of the differences between HSBs and LSBs.

### 3.1 General features

In BP99, the disc of our Galaxy is simulated as an ensemble of concentric rings gradually built up by infall of gas of primordial composition. The adopted Star Formation Rate (SFR), as suggested by the density wave theory, is

$$\psi =\alpha \mathrm{\Sigma }_{gas}^{1.5}\frac{V(R)}{R}$$
(1)
where V(R) is the disc rotation velocity as function of radius, $`\mathrm{\Sigma }_{gas}`$ the gas density (in M pc-2), $`\psi `$ the star formation rate (in M pc-2 Gyr-1), and $`\alpha `$ is an efficiency parameter adjusted as to reproduce the properties of the Solar Neighborhood. That form of the SFR is in agreement with the observed current SFR profile of the Milky Way (BP1999) and of some nearby spiral galaxies (Boissier et al. 2003, in preparation). We adopt the Initial Mass Function (IMF) of Kroupa et al. (1993), characterised by a flattening at low masses, and consider it as universal, i.e. not evolving in time or space. The infall rate is declining exponentially, with a timescale $`\tau `$ = 7 Gyr in the local disc in order to reproduce the observed G-dwarf metallicity distribution. Inner galactic zones are formed earlier than the outer disc (“inside-out” disc formation). The evolution is computed independently at each radius and the global properties (masses, luminosities, colours) are obtained by integration or by fitting the obtained profiles to recover scalelengths and central surface brightnesses. For further details on the ingredients and the model of the Milky Way, see BP99.

In BP2000, the mass and scalelength of the discs of HSB spirals are computed with similar models, using simple “scaling laws”, as suggested by Mo, Mao & White (1998) in the framework of Cold Dark Matter scenarios of galaxy formation. Within this simplified approach, discs are characterised by two parameters: their circular velocity $`V_C`$ (determined by the mass of the dark halo), and their spin parameter $`\lambda =J|E|^{1/2}G^1M^{5/2}`$ where $`J`$, $`E`$, $`M`$ are, respectively, the angular momentum, total energy and mass of the dark halo, whereas $`\lambda `$ is a dimensionless quantity measuring the specific angular momentum of the dark halo. The disc scalelength can then be expressed as:

$$R_d=R_{d,MW}\times \frac{V_C}{V_{C,MW}}\times \frac{\lambda }{\lambda _{MW}}$$
(2)
and the central surface density as

$$\mathrm{\Sigma }_0=\mathrm{\Sigma }_{0,MW}\times \frac{V_C}{V_{C,MW}}\times \left(\frac{\lambda }{\lambda _{MW}}\right)^2$$
(3)
where the index $`MW`$ refers to the corresponding value in the Milky Way (see BP2000 for details). The disc mass varies with $`V_C^3`$. From equation 3 it is obvious that the central surface density depends strongly on the spin parameter and that large values of it are likely to represent LSB galaxies. Actually, Dekel and Woo (2002) present some observational evidences (based on the Sloan Digital Sky Survey data) that the surface brightness decreases with the spin parameter $`\lambda `$, even if they consider $`\lambda `$ as a “second parameter” for the LSB class of galaxies, the first parameter being the stellar mass. This is not surprising since they include all dwarf galaxies in their study, and have then a much larger range of stellar mass than we consider. Note that our surface brightnesses also depends on the stellar mass (what can be deduced from figure 7).

On the other hand, since the same star formation rate prescription (equation 1) is applied for all galaxies, a lower star formation rate in LSBs will result naturally from the dependence of the SFR on radius and on gas surface density.

In BP2000 it was argued that observational data on nearby spirals suggest that the infall time-scale decreases with galaxy mass, i.e. that less massive discs were formed later, on average; this dependence of infall time-scales on the circular velocity $`V_C`$ is also included here. Infall rates depend also on the local surface density (in order to produce an inside-out formation of the galactic discs) and this results in a slower formation of LSB galaxies.

Because LSBs display a large range of H i line widths, we computed a new set of models with six values for the circular velocity: 40, 80, 150, 220, 290 and 360 km/s (only the first one was not considered in BP2000).

### 3.2 Spin parameter distribution

In BP2000, models of HSB spirals were obtained for 80 $`<V_C(\mathrm{km}/\mathrm{s})<`$ 360 and $`0.02\stackrel{<}{}\lambda \stackrel{<}{}0.07`$. It was noted that for $`0.07\stackrel{<}{}\lambda \stackrel{<}{}0.09`$, the surface densities were so low that the discs were in fact at the limit of LSB galaxies ($`\mu _{B,0}`$ between 22.5 and 23 mag arcsec-2).

The models of BP2000 with $`\lambda \stackrel{>}{}0.06`$ have a central surface brightness close to 22 mag arcsec-2. In this work we consider larger values of the spin parameter $`\lambda `$ that produce “naturally” LSB galaxies.

Mo et al. (1998) provide a largely used spin parameter distribution (figure 1) which gives a reasonable fit to the results of numerical simulations of the formation of dark matter halos; it can be expressed as follows:

$$P_\lambda (\lambda )d\lambda =\frac{1}{\sqrt{2\pi }\sigma _\lambda }exp\left[\frac{ln^2(\lambda /\overline{\lambda })}{2\sigma _\lambda ^2}\right]\frac{d\lambda }{\lambda }$$
(4)
where the average spin parameter $`\overline{\lambda }`$ = 0.05 and $`\sigma _\lambda `$ = 0.5.

Values of $`\lambda `$ larger than 0.06 will produce discs with a central blue surface brightness lower than 22. mag arcsec-2. According to the distribution shown in figure 1, 35% of the galaxies will have spin parameters in the interval \[0.06,0.21\], while the 60% in the interval \[0.02,0.06\] are likely to be HSB spirals (see BP2000). Note that the limits of the intervals are a bit arbitrary since the surface brightness depends slightly on the velocity $`V_C`$ too. Especially, galaxies with the lower rotational velocities can be LSBs even with lower spin parameters. The limits adopted here represent however a conservative estimate for the number of LSB galaxies.

These figures are not fully consistent with the claim by O’Neil & Bothun (2000) that “the majority of the galaxies” are LSBs. However they suggest that a significant part of the population of LSB galaxies could indeed be discs of large spins. They are also at odds with the claim of Dekel and Woo (2002) that 95 % of galaxies are LSBs, but this is mainly due to differences in definitions. Dekel and Woo (2002) include in their “LSB” class all dwarf galaxies, while we limit ourselves to discs with rotation velocities larger than 40 km/s. This choice excludes of course a large number of galaxies because of the shape of the luminosity function and is justified by the fact that we are indeed interested in the behaviours of (relatively) massive LSB disc galaxies, comparing them to similar HSBs.

The LSB disc models constructed for this paper were therefore computed for 4 values of $`\lambda `$: 0.07, 0.09, 0.15 and 0.21 (the first two values were already used in BP2000). All other ingredients of the models are identical to those used in the HSB disc models, as discussed in section 3.1.

The present work can not be considered as a proof that LSB galaxies are large-spin parameter discs since we did considerate other possibilities as for instance that LSB galaxies are settling in halos with larger baryonic fractions or concentration parameters. We prefer to consider the large spin parameter hypothesis because discs with large spin parameters are expected to exist since the spin parameter distribution of figure 1 does present an important tail at large values, and that these discs with large spin parameter values were not examined in the grid of models previously published (BP2000).

## 4 Model results and comparison to observations

### 4.1 Rotation curves

In order to compute the Star Formation Rate (equation 1), we need to adopt a rotation curve $`V(R)`$ in our models. As in BP2000, it is computed as the sum of the contributions of the baryonic matter (exponential disc) and of the dark matter halo (non-singular isothermal sphere). The use of an exponential disc is justified by the fact that 80 % of the galaxies of the O’Neil sample (the larger) are well fitted by an exponential disc (see also Bell et al., 2000). The shape of the halo used to compute the rotation curve is that of a non-singular isothermal sphere. We take the core radius to be a constant fraction of the scalelength of the disk obtained by equation 2. The CDM halo density profiles are known to be too ”cuspy” with respect to observations (e.g. de Blok et al. 2001). It is thus difficult to decide which shape to use in the models. Our aim is not to include the most realistic profile of dark matter haloes, but to end up with rotation curves in rough agreement with the observed ones (this is actually all that matters for our simple models through the equation 1 used to compute the SFR). We make such a comparison in the following.

The shape and amplitude of the resulting curves depend on the spin parameter and the circular velocity, as a result of the scaling relations.

It is often argued that LSBs are dominated by dark matter (e.g. McGaugh et al. 2001, and references therein), and, indeed, in our models the contribution of the baryonic disc to the total rotation curve is lower for LSBs than for HSBs, as a result of the lower densities in the disc.

To check that our rotation curve are realistic enough for our purpose, we compare in figure 2 rotation curves computed with spin parameters larger than 0.07 to a number of typical rotation curves constructed by de Blok, McGaugh & Rubin (2001) for 26 LSB galaxies with available high quality H$`\alpha `$ and H i data.

For most galaxies a reasonable agreement is obtained between model curves and observations. We note that in 2 out of the 26 objects (U6614 and U11748), a peak appears close to the centre, which may be due to a pronounced bulge component. Since our models concern discs only, we cannot reproduce such features and we clearly underestimate the rotation velocity (and hence the SFR) in the inner $``$10 kpc of these two objects. Still, these are large and massive galaxies and only their inner parts are affected.

We consider this comparison as satisfactory and we conclude that our simple rotation curves are realistic enough to compute the SFR with equation 1.

### 4.2 Model results

Figure 3 to 5 display some of our results concerning respectively the chemical and the spectro-photometric evolution of LSBs. On the left side of each figure, results are given for a typical value of the spin parameter $`\lambda `$=0.15 and for five values of the rotational velocity $`V_C`$ of our grid of models. On the right side results are given for a typical value of the rotational velocity $`V_C`$=150 km/s and for three values of the spin parameter $`\lambda `$ of our grid of LSB models.

In Figure 3, we present from top to bottom the SFR, the Star Formation Efficiency (SFE) and the gas fraction. The trends obtained in the case of HSBs (BP2000) are also reproduced here: more massive discs formed their stars earlier, had higher star formation efficiencies in the past and have smaller gas fractions than their lower mass counterparts; these features characterise also discs with low angular momentum. LSBs are however different from HSBs: the models have lower SFE because of the lower gas densities. As a result, LSBs at the current epoch are less evolved, with large gas fractions and young stellar populations.

In Figure 4, the magnitude (top), B-K colour (middle) and central surface brightness (bottom) are presented. An important property is again similar to what was found in HSBs: the more massive galaxies (larger rotational velocities) are redder than low mass ones. This figure shows also that LSB galaxies were generally fainter in the past, and that their surface brightness has stayed low during their evolution, according to the models.

The mass accretion history of our models depends on the mass of the galaxy. This dependence was constrained by observations in spirals (e.g. Boissier et al., 2001): it is purely empirical and is not derived from a cosmological context. In figure 5, we compare the obtained histories with the predictions for dark matter halos in our hierarchical universe with the cosmology $`\mathrm{\Omega }_M`$=0.3, $`\mathrm{\Omega }_\mathrm{\Lambda }`$=0.7, h=0.65). We assume for each disc that the cold dark matter halo is ten times as massive as the baryonic disc (at redshift 0), and that the halo follows the “universal mass accretion history” of van den Bosch (2002), depending mainly on the mass of the halo.

For intermediate circular velocities, the mass history empirically established follows a similar trend as the one of the dark matter halo. For lower masses, we clearly observe that the baryon accretion in our models occurs later than the dark matter accretion. This might be due to feedback issues in low-mass galaxies (e.g. Dekel and Woo, 2002). In the simple models, no feedback is included explicitly, and the mass accretion rates we were conducted to adopt mimic its effect.

When we move towards more massive galaxies, our baryonic discs build up more and more rapidly, and for the most massive galaxies even more rapidly than the dark matter halo. It has been already noticed that the properties of massive spirals (red colours, high metallicities, low gas fractions) are well reproduced in the simple models by assuming very high redshift formation. Some more realistic cosmological models might be useful to study this paradoxal situation, but it is beyond the scope of our paper. Note that very massive galaxies represent a small fraction of the galaxies we are studying in this paper.

### 4.3 Tully-Fisher relation

The Tully-Fisher relation, the linear relation between the absolute magnitude and the logarithm of the inclination-corrected H i line width, is one of the best-known properties of spiral galaxies. Observationally it is still difficult to establish if LSBs follow the same Tully-Fisher relation as HSB spirals or not. Some studies conclude that LSBs follow the same relation as classical spirals (Zwaan et al. 1995; Sprayberry et al. 1995; Chung et al. 2002), while in others LSB are found to be underluminous for a given line width (Persic & Salucci, 1991; Matthews et al. 1998). Chung et al. (2002) suggest that some of these differences may be due to selection effects.

Although O’Neil et al. (2000a) find that a few of their LSB galaxies show large departures from the standard Tully-Fisher relation, subsequent VLA H i imaging of four of them and a new analysis of the single-dish data (Chung et al. 2002) show that all (or most) of their original H i profiles may have been contaminated by H i from a nearby galaxy.

In figure 6 we plot the absolute blue magnitude as a function of the H i line-width, corrected for inclination, for the LSBs samples of De Blok et al. (crosses), O’Neil et al. (triangles), Giant LSBs (triangles) and the infrared-selected galaxies (filled dots). The data for the O’Neil et al. galaxies were taken from Chung et al. (2002), after cleaning from known contaminations. The combination of the different sample shows the existence of a Tully-Fisher relation, though with a considerable scatter. Our LSB models are overplotted, with each curve corresponding to a different spin parameter. Models with a lower spin parameter are also shown (stars) for reference. It is assumed that the circular velocity of our models, $`V_C`$, corresponds to half the inclination-corrected H i line-width, 0.5W50/sin(i).

The models predict the existence of a Tully-Fisher relation for both HSB and LSB disc galaxies. However, the adopted scaling relations (section 3.1) lead to a difference between the HSB and LSB TF relations: at a given rotational velocity the luminosity of LSBs is lower than for HSB spirals because of differences in their star formation history and gas fraction. The amplitude of the effect depending on the velocity $`V_C`$, the slope of the relation changes with the spin parameter.

Although the modeled and observed LSB Tully-Fisher relations show relatively similar slopes, the model values lie on average about 2 mag below the mean observed absolute magnitudes. This shift could be due to systematic differences between the measured and modeled circular velocities and to uncertainties in both of them.

### 4.4 Scalelength and central surface brightness

Figure 7 displays the disc scalelength and central surface brightness in the B-band as a function of $`M_{B,d}`$, the absolute magnitude of the disc component (as defined in section 2). ¿From the adopted scaling relations (section 3.1), we expect that the scalelength will increase with luminosity, while the central surface brightness will depend more strongly on the spin parameter. This is indeed the case, as seen in figure 7.

We note that the lack of observed galaxies in some parts of these diagrams (specifically, those with large luminosities but very faint central surface brightnesses, the so-called “Malin 1 cousins”) might result from selection effects. On the other hand, in view of the distribution of the spin parameter (figure 1) and of the circular velocity function (Gonzalez et al. 2000) it is expected that only a small number of LSB galaxies with large rotational velocities should exist. Therefore, one would expect such galaxies to be quite rare compared to LSBs of considerably smaller size or less extreme surface brightness.

Our models produce LSB disc galaxies similar in size and surface brightness to those of the presently available observational samples. However, within the framework of our models, giant discs with lower surface brightness could in principle be produced by systems with even larger $`\lambda `$ values than those adopted in our present calculations.

Finally, the scaling relationships indicate that LSBs should have larger scalelengths than HSBs for the same mass. The models for HSBs (stars) in Figure 7 are indeed located at smaller scalelenths than the models for LSBs. The same result is obtained in the observations since the MR2002 sample is characterised by smaller scalelengths since the other data, and correspondingly higher surfaces brightnesses.

### 4.5 Gas content

Our models predict a correlation between the absolute disc magnitude and the gaseous mass (bottom panel of figure 8). This results mainly from the mass-luminosity relation: more massive galaxies are more luminous *and* contain more gas (as expected if the gas-to-total mass ratio does not depend strongly on luminosity). The data follow the same trend, but the observational scatter is larger than what our models can account for, especially in the faint LSB range (where most of the data come from the O’Neil et al. sample).

Note that in a recent study, Galaz et al. (2002) presented a clear correlation between the HI mass and the Ks magnitude for a sample of galaxies they studied with near-infrared imagery, confirming the trend we observe in Figure 8. As a direct consequence of the scaling relationships, using the words of Galaz et al. (2002), “bigger galaxies have more of everything”.

The gas mass-to-light ratio ($`M_{gas}/L_B`$) is an interesting intensive quantity, more related to the “chemical state” of the galaxy than the absolute gas mass. Our models predict larger values of $`M_{gas}/L_B`$ for LSBs than for HSB spirals, as a result of the lower star formation efficiency in the former. In figure 8 it appears that our models predict a slightly declining trend with absolute magnitude: more massive galaxies have proportionally smaller gas fractions, due to the smaller time-scale of infall rates. The data show a steeper trend (arising mainly from the giant LSBs) with a broader dispersion than allowed by our models. For the LSBs data of O’Neil et al., two values of this ratio are shown, giving an indication of the corresponding uncertainty: one value is computed with a luminosity derived from the absolute magnitude $`M_{B,d}`$ of the disc only and one is based on the total absolute magnitude of the galaxy (the two filled square linked by a vertical line).

Note that the HSBs models have lower gas contents than the LSBs models, what is in agreement with the observations: the MR2002 sample is characterized by smaller values of the gas to luminosity ratio than the other samples.

A way to obtain a stronger decrease of $`M_{gas}/L_B`$ with $`M_{B,d}`$ in our models would be to introduce a different parameterisation of the key physical ingredients in order to consume more gas in massive discs, either by either making the gas available earlier (more rapid infall) or by changing the prescription for star formation (enhancing the SFR efficiency). Still, given the large observational uncertainties and the crudeness of our “1st order” models, we consider that the simple “large-spin” models presented here reproduce satisfactorily the observed $`M_{gas}/L_B`$ values. We also note that departures from these values could be due to star formation events that occur in addition to the smooth history of the simple models, as will be discussed in section 5.

### 4.6 The luminosity-metallicity relation

The “luminosity-metallicity” relation is an important property of (HSB) spiral galaxies (e.g. Zaritsky, Kennicutt & Huchra 1994): massive galaxies contain proportionally more metals than low-mass ones. The metallicity is estimated from the oxygen abundance at a given radius in the disc of the galaxy (see Zaritsky et al. 1994 for details).

This relation is very important for studies of galactic evolution since it suggests that massive galaxies have processed more material into stars, or that low mass galaxies have lost part of their metals or a combination of both occurred. The former hypothesis was adopted in Boissier et al. (2001), on the basis of observed star formation efficiencies and gas fractions, for spirals with rotational velocities as low as 80 km/s. Recently, Garnett (2002) used a large data sample of spiral and irregular galaxies with rotational velocities as low as 5 km/s and convincingly argued that below $``$120 km/s the effective yield of galaxies is reduced, most probably as a result of mass loss. We note that this effect, by itself, is not sufficient to explain the observed colour-luminosity relation of spirals, with the most massive ones being redder than their low mass counterparts; a younger age for the latter has to be assumed in order to account for the data (see Prantzos 2002 for a recent review).

Because our LSBs models are an extrapolation of the HSB ones, we expect them to produce a luminosity-metallicity relation in the case of LSBs also. However, McGaugh (1994) compiled metallicity measurements in LSBs and concluded that they do not display such a relation \[his data are shown in figure 9 as a function of $`M_{B,d}`$\]. It should be noted that McGaugh (1994) excluded giant LSBs from his analysis, while we consider a continuous range of galaxies, from dwarfs to giants; if we add the giants to the sample, a small trend is observed (see figure 9). On the other hand, the abundance measurement in LSBs concern only a few H ii regions per object, located at various galactocentric radii. Therefore, it is impossible to determine an abundance at a “characteristic” radius and the situation is very different from the work of Zaritsky et al. (1994) concerning HSB spirals. This difficulty introduces a large scatter in the luminosity-metallicity relation (if there exists one), compared to the case of the HSBs.

In figure 9 we compare the data of McGaugh (1994) with the metallicities of our LSB models. Since the observations were done at various radii, for a comparison with the models we computed the distribution expected in the absolute magnitude - log(O/H) plane, by weighting the results with the distributions of the spin parameter (figure 1), circular velocity (we adopt the distribution of Gonzalez et al. 2000), and Star Formation Rate. The reason for introducing the latter is that abundances are measured in H ii regions, implying the existence of massive stars and, hence, important star formation activity.

As can be seen from figure 9, the abundances measured in LSBs are in satisfactory agreement with the values expected from our simple models, at least in a statistical sense.

Galaz et al. (2002) argue that the J-Ks colour index they measured in their near-infrared galaxy sample is a good indicator of the metallicity. On this base, they deduce the existence of a luminosity-metallicity relationship, the redder colours corresponding to the galaxies with the largest stellar mass. From the variation of the J-Ks index between the less luminous and the more luminous galaxies, they conclude that the metallicity increases by a factor between 20 and 100 between them. In a logarithmic scale, this is 1.3 to 2 dex, what compares well with the slope we can estimate from Figure 9.

### 4.7 Colours

Three colour indices are presented as a function of absolute disc magnitude in figure 10. The models are characterized by blue colours and an extremely small colour dispersion, especially for low-mass discs. Both features stem from the fact that most of the model galaxies have formed the bulk of their stars toward the end of their history (see figure 3). The differences in their early star formation histories play a minor role (especially in low-mass discs) because the colours are always dominated by the youngest population of stars.

Although these model results are in broad agreement with the “blue” LSBs of the observational sample, they do not reproduce the colours observed in the large number of “red” LSBs. It is worth noting that the red LSBs are mainly found in the O’Neil et al. (2000) sample, part of which consists of cluster galaxies, though it should also be noted that O’Neil et al. (2000a) find no difference between the properties of their LSBs in clusters and in the field.

We also point out that colour uncertainties may be larger than usually quoted ($`0.1`$ mag), since O’Neil et al. (2000a) and Bell et al. (2000) find for the same object (P1-7) $`BV0.9`$ and $`BV=0.6\pm 0.1`$, respectively. Notwitstanding this difference and the importance of the uncertainties in LSBs’ colours, the existence of “red” LSB discs is undisputable and that property is clearly not reproduced by the simple models presented here. In section 5, we suggest that another “ingredient” is necessary to reconciliate the simple models with the existence of such red galaxies.

### 4.8 Summary of the “simple” model calculations

The results of our simulations can be summarised as follows:

$``$ Simple models (i.e. with smooth SF histories), differing by HSBs only by their larger spin parameters ($`\lambda `$), are in reasonable agreement with several observed properties of LSBs, e.g. the relations between absolute disc magnitude and central surface brightness, disc scalelength, metallicity, as well as total gas content. The models also predict a Tully-Fisher relation, in broad agreement with observations.

$``$ Some important discrepancies exist between the observations and the “simple” models: the scatter in many of the observed relations (Tully-Fisher, gas mass vs. luminosity, colours) is much larger than predicted by the models. In particular, the “red” LSB galaxies are not reproduced by any of our models, since “by construction” they are dominated by the young stellar population.

The discrepancies between models and observations are particularly important. In fact, the models of BP2000 (on which this work relies) were developed mainly for field galaxies with smooth star formation histories. If the environment of LSB is different or if they are more sensitive to perturbations (because of their lower densities), some of their properties may be due to departures from the ideal smooth star formation histories. This idea is tested in the next section, where we try to reconcile the large-spin models and the discordant observations by adding “events” (bursts and truncations) on otherwise smooth star formation histories. Indeed, O’Neil et al. (2000a) suggested that the properties of both blue and red LSB galaxies can be reproduced by population synthesis models, assuming “on” and “off” states of the star formation for different stellar populations. Van den Hoek et al. (2000) already suggested that small amplitude bursts may play a role in determining the colours of LSBs. In their N-body simulations, Gerritsen & de Blok (1999) argued that Star Formation Rate fluctuations are responsible for the large range of colours observed in LSBs and predicted that less than 20% of gas-rich LSBs should be “red”.

## 5 Starbursts and truncated star formation histories

In this section, we investigate a way to produce red LSB galaxies while introducing only small changes to the “simple” models presented in section 3. We will refer to the star formation histories derived with the simple models as “standard” histories, which represent a smoothed average of the true evolution of the LSBs, to which we will now add star formation “events”.

First, we study the effect of a starburst added to the “standard” SF history of one of our models, and then we consider the effect of a truncation of the “standard” star formation at a recent epoch.

There are reasons to believe that the star formation rate behaves in such a way in LSB galaxies.

Martin and Kennicutt (2001) show that in nearby spirals, beyond the “threshold radius”, where the density of gas becomes smaller than the critical density (according to the Toomre instability criteria) a few HII regions (and then recently formed stars) are still found. They argue that local perturbations and density waves can temporarily and locally enhance the density to values high enough to allow star formation. This way of forming stars is obviously less smooth than the one described in the simple averaged models. Despite a large scatter, we can derive from the figure 9 of Kennicutt (1989) that the critical gas density at which the threshold is found range between 3 and 10 $`M_{}pc^2`$ for most of the galaxies. This range of values is indicated in the Figure 11 as a grey area. The gas surface density in some of our LSB models is indicated as a function of radius at the present epoch (top panel). In the bottom panel, each curve presents the evolution with time of the gas surface density at a radius equal to one scalelength, for the same models as in the top panel.

Obviously, the gas surface density has been lower than (or close to) the threshold value observed in spiral galaxies, in most of our model galaxies and most of the time. Thus, the star formation is more likely to proceed through a succession of quiescent phases (when the density is lower than the threshold), and relatively active phases (when enough gas has been accumulated by infall and condensed by local perturbations or spiral modes to reach the threshold). This scenario is indeed suggested by Galaz et al. (2002) in order to explain the scatter in optical data by a sporadic past star formation history.

Another possibility is that the star formation history is affected by feedback from supernovae. Feedback effects may be crucial to the history of LSB galaxies (see for instance the study of Dekel and Woo, 2002). We can speculate that the effect of supernovae could be similar to the empirical threshold discussed above: stopping the star formation by blowing the gas away and delaying the birth of the next generation of stars until enough gas has settled again in the disk.

It is beyond the capacities of the simple models to test the details of these ideas concerning the way the star formation rate is affected in low density environment by threshold effects, spiral modes, perturbation, feedback. This task could be undertaken only with much more sophisticated hydrodynamical models (e.g. Samland and Gerhard, 2003). We compute instead two simple examples for illustrative purpose.

Of course, the star formation history of a galaxy may consist of a series of “bursts” and “drops” in its star formation rate. However, the “chemical” properties of a 13 Gyr old galaxy with such a jittery SFR would be virtually identical to those resulting from an averaged star formation history, while its colours may be affected mainly by the last of such events. Therefore, our adopted “smooth plus last event” model, while simplistic, may be a good enough approximation for the present study of LSBs.

### 5.1 Adding a star burst

Figure 12 shows the $`BI`$ colour and the gas mass-to-light ratio in the $`B`$-band for the galaxies from the various samples and the grid of “standard” models.

It would be difficult to get a clear insight by adding bursts of various ages, intensities, metallicities, etc. to all of the “simple” models. Instead, we prefer to study one illustrative example, by adding a single burst of star formation to the “standard” star formation history of one of our models ($`\lambda `$=0.07, $`V_C`$=80 km/s). This burst creates 5 106 M of stars with a metallicity equal to one third of solar and a “normal” IMF (from Kroupa et al. 1993 in all our models).

When the burst is 0.005 Gyr old, the model results correspond to the position indicated by the tip of the arrow labelled “B” in figure 12; for older bursts the results are located along the curve, as indicated by the corresponding numbers, and they are progressively closer to those of the “standard” (smooth SF) models.

For very recent bursts, it seems likely that the brightest part of the galaxy would be the burst itself, in which case the galaxy would probably not be classified as an LSB, but rather as a “blue compact” (if the burst is centrally located). Legrand et al. (2000) proposed such a scenario for i Zw 18, a blue compact galaxy of very low metallicity, which could result from a star burst occurring in an LSB galaxy after a long period of low and constant star formation rate.

For burst ages under 0.5 Gyr, the colours of the galaxy are bluer than without the burst; beyond that age the old stars of the burst dominate the stellar populations of the galaxy, which becomes redder than in the “standard” model. However, even for a very old burst, we do not obtain a very red galaxy. In fact, after 12.5 Gyrs the colours of the galaxy are quite close to those of the smooth “standard” model.

Assuming that all the gas taking part in the burst is removed from the galaxy, we can also compute a total gas mass-to-light ratio. For young bursts, $`M_{gas}/L_B`$ is lower than without the burst (because the luminosity increases significantly and the gas amount is reduced). For older bursts, $`M_{gas}/L_B`$ proceeds towards “standard” values, as the luminosity of the burst (which dominates the stellar population) decreases with age.

The results obtained here suggest that adding star bursts to the smooth SF history of the “standard” model can increase the scatter in the gas mass-to-light ratio (see figure 12); it will also increase the scatter around the Tully-Fisher relation, especially for luminous LSBs with a relatively low gas mass-to-light ratio.

### 5.2 Truncating the star formation

To illustrate the effect of truncations of the star formation history, we added to the same model ($`\lambda `$=0.07, $`V_C`$=80 km/s) the condition that 1 Gyr before the present epoch the star formation stopped abruptly. As a result, the stellar population reddens passively, since no young and blue stars are created. The resulting stellar population is different from that of a “normal” galaxy, without the need to invoke a special IMF.

The effect of truncating the star formation history in this way is shown in figures 12 and 13 by arrows labelled with a “T”. The luminosity decreases by about 1.5 magnitudes due to the lack of newly formed stars, while the colours become quite red. This model corresponds rather well to the faint red LSB galaxies of the O’Neil et al. sample.

Because of the decrease in luminosity, galaxies with larger than average gas mass-to-light ratio are obtained. Notice, however, that the scatter in the $`M_{gas}/L_B`$ vs $`M_B`$ relation that is obtained by adding truncations to a smooth SF history is smaller than that due to the diversity in the values of the spin parameter $`\lambda `$. Another consequence is that truncations of the SF will increase the scatter in the Tully-Fisher relation, especially for low-luminosity systems (as is the case with starbursts)

### 5.3 The effects of bursts and truncations

The two simple examples presented in the previous sections illustrate that a starburst, or a truncation of the star formation, may strongly affect the colours, luminosity, and gas mass-to-light ratio of LSB galaxies, with respect to the “simple” models with smooth histories discussed in section 3.

If the simple models reproduce the average history of LSB galaxies, then the occurrence of star bursts of various ages, masses and metallicities, and of truncations of star formation at various galactic ages, may explain the existence of the “red” LSBs as well as the large scatter found in several observational properties of LSBs, like their Tully-Fisher relation, gas mass-to-light ratio, and colours.

If some LSBs have had recent truncations in their star formation history their colours may be a poor indicator of their true nature; indeed, a massive LSB disc (having an old stellar population according to our models) may have the same $`BI`$ colour as a less massive one with a truncated star formation history, evolving passively for e.g. a Gyr. This is in agreement with the finding of Bell et al. (2000) that “red” LSB galaxies form an heterogeneous group.

Finally, we note that the events we invoked (bursts and truncations) should affect more strongly the optical than the infra-red properties of LSBs, as the latter are more sensitive to the old, average, stellar population. We present in figure 13 the equivalent of figure 12 for the near infra-red $`K`$-band, showing the effects of the burst and the truncation. A comparison between figures 12 and 13 shows that infra-red properties are indeed less affected, and that we obtain less extreme colours in the infra-red, indicating that adding star formation events will introduce a relatively much smaller scatter in the colours, Tully-Fisher relation and gas mass-to-light ratio in the infra-red than in the optical.

In their recent study of a near-infrared sample of galaxies, Galaz et al. (2002) foundd that the near-infrared data and colours are less scattered than the optical ones. This is in full agreement with our previous deduction. They also found a decrease of the HI mass to light ratio with the Ks magnitude, together with an increase of the value of the J-K colour index, both corresponding qualitatively to the results of figure 13.

## 6 Conclusions

We have compiled data from various observational samples of Low Surface Brightness galaxies (LSBs) in order to cover the whole presently known range of their properties, from “dwarfs” to “giants”, from “blue” to “red” galaxies and from those at the limits of detectability in central disc surface brightness to objects straddling the adopted boundary between LSBs and HSBs (High Surface Brightness galaxies). We have not included in our study very faint dwarfs, but only gas-rich LSB discs, with rotational velocities larger than 40 km/s or so. We compared them to the results of models of their chemical and photometric evolution, obtained under the assumption that these LSBs are similar to HSBs except for a larger angular momentum. This was done on the base of successful models of HSBs that were presented in a recent series of papers.

The comparison is relatively satisfactory: observed disc scalelengths, central surface brightnesses and total gas masses are compatible with the model predictions. In their recent work, Schombert et al. (2001) also found that the gas fraction of their LSB dwarfs is compatible with our main assumption. However, the star formation history of some of the LSBs appears to be different from this simple picture, as indicated by the large scatter found in several observational relations (Tully-Fisher, total gas mass-to-luminosity) as well as by the existence of red galaxies, which cannot be not explained within the framework of the simplest models (i.e. large angular momentum discs with a smooth star formation history).

In order to account for the discrepancies between the simple models and the observations we suggest that star formation “events”, bursts and truncations, should be included to the simple model. We show that these ingredients may easily reconcile the observations with the models, since they reproduce the observed large scatter in certain relations; besides, truncations in the SF history may account for the existence of the red LSBs found by O’Neil et al. (1997b). Such “events” may result from the fact that the gas surface density in most LSBs is probably lower than the threshold for star formation found by Kennicutt (1989), or from supernovae feedback (Dekel and Woo, 2002). Such events are however unlikely to alter as much the near infrared data as the optical ones, quite in agreement with the observations of Galaz et al. (2002)

We conclude that LSBs may be the equivalent of HSBs with larger angular momentum (spin parameter $`\lambda `$), but that their observed properties could be largely affected by recent events in their star formation history; in that respect, they present similarities to the case of blue compact galaxies.


text:
# Doping- and size-dependent suppression of tunneling in carbon nanotubes

\[

## Abstract

We study the effect of doping in the suppression of tunneling observed in multi-walled nanotubes, incorporating as well the influence of the finite dimensions of the system. A scaling approach allows us to encompass the different values of the critical exponent $`\alpha `$ measured for the tunneling density of states in carbon nanotubes. We predict that further reduction of $`\alpha `$ should be observed in multi-walled nanotubes with a sizeable amount of doping. In the case of nanotubes with a very large radius, we find a pronounced crossover between a high-energy regime with persistent quasiparticles and a low-energy regime with the properties of a one-dimensional conductor.

71.10.Pm, 71.20.Tx, 72.80.Rj

\]

During recent years there has been much interest in the investigation of the electronic properties of carbon nanotubes (CN). The reduced dimensionality of these systems leads to the appearance of unconventional effects, such as a suppression of the tunneling conductance at low energy scales. This has been interpreted as a signature of the so-called Luttinger liquid behavior, characterized by the absence of electron quasiparticles in the spectrum. Evidence of a power-law behavior in the tunneling density of states at low energies has been obtained from measurements in ropes, individual single-walled nanotubes (SWNT) and multi-walled nanotubes (MWNT).

In the Luttinger liquid picture, critical exponents of observables like the density of states are not universal and depend on the interaction strength. In the case of carbon nanotubes, this refers to the ubiquitous Coulomb interaction. However, a precise determination of the interaction strength is precluded by the fact that the electron-electron interaction is actually long-ranged, and it can be assimilated to a coupling constant, only after introducing a suitable infrared cutoff in the singular expression of the one-dimensional (1D) Coulomb potential.

Furthermore, the shape of CN becomes relevant in setting the strength of the electron correlations. In the case of MWNT, measurements of the conductance refer usually to the outer layer, whose electronic properties are influenced by the interaction with inner metallic cylinders. Also, MWNT use to be significantly doped, what leads to the presence of a large number of subbands at the Fermi level. The contribution of a large number of modes at low energies has then an appreciable impact in the enhancement of observables like the tunneling density of states.

The purpose of the present paper is to study the combined effect of the finite length of the electron system, on the one hand, and of the number of subbands at low energies, on the other hand, in the Luttinger liquid description of CN. With this aim, we adopt a renormalization group (RG) approach, which is well-suited to obtain the energy dependence of quantities like the quasiparticle weight. In order to cure the infrared singularities arising from the long-range interaction, we implement a dimensional regularization of the theory, writing formally all bare quantities slightly away from dimension $`D=1`$. This formal artifact has a real physical meaning, since the deviation $`D1`$ yields a measure of the finite size of the system, with the limit $`D1`$ corresponding to the case of increasingly long nanotubes.

The Coulomb potential $`1/|𝐫|`$ can be represented in three spatial dimensions as the Fourier transform of the propagator $`1/𝐤^2`$. If the interaction is projected onto one spatial dimension, by integrating for instance the modes in the two transverse dimensions, then the Fourier transform has the usual logarithmic dependence on the momentum, $`V_C(k)=(1/2\pi )\mathrm{ln}(k_c/k)`$. In this approach, $`k_c`$ is the memory that the system keeps of the finite size of the transverse dimensions. We choose instead to integrate formally a number $`3D`$ of dimensions, so that the long-range potential gets the representation

$$\frac{1}{|x|}=\frac{d^Dk}{(2\pi )^D}e^{ikx}\frac{c(D)}{|𝐤|^{D1}}$$
(1)
where $`c(D)=\mathrm{\Gamma }((D1)/2)/(2\sqrt{\pi })^{3D}`$. In the limit $`D1`$, the expression of the potential in momentum space becomes

$$\frac{c(D)}{|𝐤|^{D1}}\frac{1}{2\pi }\left(\frac{1}{D1}\mathrm{ln}(|𝐤|)+\mathrm{}\right)$$
(2)
This shows that $`1/(D1)`$ corresponds actually to the logarithm of the high-energy cutoff $`k_c`$. In the RG framework, where all quantities are scaled down to low momenta, the dimensionless quantity $`1/(D1)`$ has to be traded for the logarithm of the length of the system, measured in units of the nanotube finite radius.

The great advantage of dealing with the above representation is that the Coulomb interaction gives rise classically to a scale invariant theory at any dimension above $`D=1`$, which allows to proceed with the RG program. It is only at $`D=1`$ that the mentioned logarithmic dependence on the momentum leads to an imperfect scaling behavior, making the use of the RG approach not quite appropriate. Then, we write the hamiltonian for the linear branches of CN in the form

$`H`$ $`=`$ $`v_F{\displaystyle \underset{\alpha \sigma }{}}{\displaystyle _0^\mathrm{\Lambda }}𝑑p|𝐩|^{D1}{\displaystyle \frac{d\mathrm{\Omega }}{(2\pi )^D}\mathrm{\Psi }_{\alpha \sigma }^+(𝐩)𝝈\mathbf{}𝐩\mathrm{\Psi }_{\alpha \sigma }(𝐩)}`$ (3)
$`+e^2{\displaystyle _0^\mathrm{\Lambda }}𝑑p|𝐩|^{D1}{\displaystyle \frac{d\mathrm{\Omega }}{(2\pi )^D}\rho (𝐩)\frac{c(D)}{|𝐩|^{D1}}\rho (𝐩)}`$ (4)
where the $`\sigma _i`$ matrices are defined formally by $`\{\sigma _i,\sigma _j\}=2\delta _{ij}`$ and $`\rho (𝐩)`$ are density operators made of the electron fields $`\mathrm{\Psi }_{\alpha \sigma }(𝐩)`$, with $`\alpha `$ labeling the Fermi point and $`\sigma `$ the spin projection. The sum in Eq. (4) runs then over the usual four modes present in a single-walled nanotube, but it may include also the contribution from a large number of subbands in the case of a doped multi-walled nanotube.

We focus on the scaling properties of the model as the cutoff $`\mathrm{\Lambda }`$ is lowered, when a large number $`N`$ of subbands contribute to the electronic properties down to the Fermi level. Each subband is labeled with a different quantum number, which corresponds to the momentum in the dimension around the nanotube. For this reason, the dominant processes are those where each scattered electron remains in the same linear branch. The main effect of the interaction is to dress the bare electron propagator with the polarization of the $`N`$ different subbands given by

$$\mathrm{\Pi }(𝐤,\omega _k)=2Nb(D)\frac{v_F^{2D}𝐤^2}{|v_F^2𝐤^2\omega _k^2|^{(3D)/2}}$$
(5)
where $`b(D)=\frac{2}{\sqrt{\pi }}\frac{\mathrm{\Gamma }((D+1)/2)^2\mathrm{\Gamma }((3D)/2)}{(2\sqrt{\pi })^D\mathrm{\Gamma }(D+1)}`$ . The polarization (5) is the analytic continuation of the known result for two linear branches with opposite chirality, which we take away from $`D=1`$ in order to carry out a consistent regularization of the Coulomb interaction. After dressing the interaction with the polarization (5), the electron self-energy is given by the expression

$`\mathrm{\Sigma }(𝐤,\omega _k)`$ $`=`$ $`e^2{\displaystyle _0^\mathrm{\Lambda }}𝑑p|𝐩|^{D1}{\displaystyle \frac{d\mathrm{\Omega }}{(2\pi )^D}\frac{d\omega _p}{2\pi }}`$ (6)
$`G(𝐤𝐩,\omega _k\omega _p){\displaystyle \frac{i}{\frac{|𝐩|^{D1}}{c(D)}+e^2\mathrm{\Pi }(𝐩,\omega _p)}}`$ (7)
The low-energy properties of the theory are investigated by taking the limit $`\mathrm{\Lambda }0`$, where the self-energy $`\mathrm{\Sigma }`$ turns out to have terms linear in $`\omega _k`$ and $`𝐤`$ that depend logarithmically on the cutoff. This is the signal that the scale of the electron wavefunction $`Z^{1/2}`$ and the Fermi velocity $`v_F`$ are renormalized at low energies. The divergent contributions to the electron propagator read

$`{\displaystyle \frac{1}{G}}`$ $`=`$ $`{\displaystyle \frac{1}{G_0}}\mathrm{\Sigma }Z^1(\omega _kv_F𝝈\mathbf{}𝐤)`$ (8)
$`Z^1{\displaystyle \frac{f(D)}{2N}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}(1)^ng^{n+1}({\displaystyle \frac{n(3D)}{n(3D)+2}}\omega _k`$ (9)
$`+(1{\displaystyle \frac{2}{D}}{\displaystyle \frac{n(3D)+1}{n(3D)+2}})v_F𝝈\mathbf{}𝐤)h_n(D)\mathrm{log}(\mathrm{\Lambda })`$ (10)
where $`g=2Nb(D)c(D)e^2/v_F`$, $`h_n(D)=\frac{\mathrm{\Gamma }(n(3D)/2+1/2)}{\mathrm{\Gamma }(n(3D)/2+1)}`$, and $`f(D)=\frac{1}{2^D\pi ^{(D+1)/2}\mathrm{\Gamma }(D/2)b(D)}`$ .

The usual RG argument is that the renormalized propagator $`G`$ must be a finite quantity, so that the divergent dependences on the cutoff $`\mathrm{\Lambda }`$ have to be reabsorbed in the scale of the wavefunction $`Z^{1/2}`$ and the Fermi velocity $`v_F`$ . Under a differential variation of $`\mathrm{\Lambda }`$, $`Z^{1/2}`$ is renormalized according to the equation

$`\mathrm{\Lambda }{\displaystyle \frac{d}{d\mathrm{\Lambda }}}\mathrm{log}Z(\mathrm{\Lambda })`$ $`=`$ $`{\displaystyle \frac{f(D)}{2N}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}(1)^ng^{n+1}`$ (11)
$`{\displaystyle \frac{n(3D)}{n(3D)+2}}h_n(D)`$ (12)
The renormalization of $`v_F`$ can be translated into that of the effective coupling $`g=2Nb(D)c(D)e^2/v_F`$, since the electron charge $`e`$ is not renormalized in our model. The RG equation for $`g`$ becomes

$`\mathrm{\Lambda }{\displaystyle \frac{d}{d\mathrm{\Lambda }}}g(\mathrm{\Lambda })`$ $`=`$ $`{\displaystyle \frac{f(D)}{2N}}{\displaystyle \frac{2(D1)}{D}}g^2{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}(g)^n`$ (13)
$`\left({\displaystyle \frac{(3D)n+1}{(3D)n+2}}\right)h_n(D)`$ (14)
We are now in a position to study the influence of the long-range Coulomb interaction in the limit $`D1`$. For this purpose, we start by considering the RG equation (14). The series showing in the r.h.s. can be summed up at $`D=1`$, the flow equation taking then the following form in the neighborhood of that point:

$$\mathrm{\Lambda }\frac{d}{d\mathrm{\Lambda }}g(\mathrm{\Lambda })\frac{1}{2N}(D1)g\left(1\frac{1}{\sqrt{1+g}}\right)$$
(15)
It follows from Eq. (15) that, at $`D=1`$, there is formally a line of fixed-points covering all values of the interaction strength. However, we still have to take into account that the effective coupling $`g(D)`$ becomes singular in the limit $`D1`$. We can shuffle this divergence into the initial value of the coupling, $`g_0(D)`$, which turns out to have near $`D=1`$ the asymptotic behavior

$$g_0(D)N\frac{e^2}{\pi ^2v_F}\frac{1}{D1}$$
(16)
Then, by matching the behavior of $`g_0(D)`$ with that of Eq. (15), we observe that the $`D1`$ factor in the r.h.s. of the RG equation is not completely canceled out in the limit $`D1`$ . This shows in a rigorous way that the bare 1D long-range interaction is at a RG fixed-point for arbitrary values of the interaction strength.

We consider next the RG equation (12) for the electron wavefunction scale in the limit $`D1`$. The series in the r.h.s. can be also summed up at $`D=1`$, with the result that the scaling equation in that limit reads

$$\mathrm{\Lambda }\frac{d}{d\mathrm{\Lambda }}\mathrm{log}Z(\mathrm{\Lambda })\frac{1}{4N}\left(\sqrt{1+g}+\frac{1}{\sqrt{1+g}}2\right)$$
(17)
The function in the r.h.s. of Eq. (17) coincides with the anomalous electron dimension found in the exact solution of the Luttinger model. This ensures that the RG approach is a sensible way to obtain the low-energy properties of the model. We have to bear in mind that our description introduces however two important differences with respect to the usual treatment of 1D interacting electrons. On the one hand, the slight deviation of the model from $`D=1`$ allows to control the approach to the bare long-range interaction as the length of the system is increased. On the other hand, we have also incorporated the effect of the number $`N`$ of subbands that contribute at low energies, in order to account for the influence of doping in MWNT.

With the RG approach we can face two different experimental conditions, depending on the magnitude of the typical energy scale involved in the measurement process, as compared to the spacing between subbands in the nanotube sample. When the latter is larger than the temperature or bias voltage applied to the sample, we are in a situation where the number of subbands $`N`$ can be taken as constant along the RG flow. Otherwise, for large enough temperature or bias voltage, the number of subbands that contribute at the scale of the high-energy cutoff is a decreasing variable as $`\mathrm{\Lambda }0`$. We deal separately with the two instances in what follows.

RG approach with a constant number $`N`$ of subbands.— In transport experiments, the typical scale of temperature or bias voltage lies usually below the scale of the spacing between subbands. This has been so, even in the measurements made in MWNT, where the spacing corresponding to a typical diameter $`d17\mathrm{nm}`$ is as small as $`\mathrm{}v_F/d29\mathrm{meV}`$. In these conditions, the only subbands that contribute to the properties measured experimentally are those crossing the Fermi level.

We apply then Eqs. (16) and (17) to confront the experimental results on the tunneling density of states gathered from different nanotube samples. Starting with the measurements made in SWNT, we take a number of subbands $`N=2`$ in the equations and adjust the deviation from $`D=1`$ in accordance to the length of the experimental sample. Following the argument below Eq. (2), we use the correspondence $`1/(D1)\mathrm{ln}(L/d)`$, $`L`$ being the nanotube length and $`d`$ the nanotube diameter. A suitable choice corresponding to the experiments reported in Ref. is $`L/d10^3`$, which gives in turn $`D1.14`$.

The measurements of the tunneling density of states in SWNT have shown a power-law dependence on energy, with values of the critical exponent $`\alpha `$ accumulating around $`0.35`$ . We have represented in Fig. 1 the estimates obtained from the r.h.s. of Eq. (17), for small values of $`D1`$. The values of $`\alpha `$ have a smooth dependence on the length $`L`$ of the system and fall around $`\alpha 0.35`$ for $`N=2`$ and $`D1.14`$, with a reasonable choice of the coupling $`e^2/\pi ^2v_F1.5`$.

Our results show also an overall agreement with the exponents measured in MWNT. It has been noticed that such systems are significantly doped, so that a large number of subbands are found at the Fermi level. The experimental conditions refer to a situation where $`N510`$ (without taking into account the spin degeneracy). It has been reported that the values of the critical exponent $`\alpha `$ measured in 11 different samples range from 0.24 to 0.37 . This variation can be accounted for within our RG approach by assuming that the number of subbands used in the renormalization may shift from $`N=2`$ to $`N=10`$. Part of the drift observed in the critical exponent may be also due to the smaller aspect ratio of MWNT, although this fact is difficult to assess given the lack of information about the total length of the experimental samples.

We show then that the suppression of tunneling in MWNT can be softened by increasing the doping level. Our results may be also relevant in general for nanotubes of large radius where there are a large number of subbands crossing the Fermi level. These instances can be considered as midway in the process of making contact with the physical properties of a graphene sheet. This requires taking systems with larger transverse size, what in turn may lead to a situation where the spacing between subbands is smaller than the typical energy scale in the experimental measurements. One has then to change slightly the computational scheme, as discussed below.

RG approach with a cutoff-dependent number of subbands $`N(\mathrm{\Lambda })`$.— In samples of a very large diameter, one may envisage conditions where the temperature or the bias voltage are much larger than the subband spacing in the nanotube. The RG approach can be still implemented, but taking into account that the number of subbands contributing in the partial integration of modes at energy $`\mathrm{\Lambda }`$ depends on the value of the high-energy cutoff. Let us suppose for simplicity that the system is at half-filling, with the typical structure of two subbands crossing at the Fermi level. For not too large energies, the number of subbands crossing the energy level $`\epsilon `$ has then a linear dependence on $`\epsilon `$, $`N(\epsilon )=N_0+n_p\epsilon `$, where $`N_0=2`$ is the number of subbands at the Fermi level.

We obtain the scaling of the quasiparticle weight $`Z`$ and the different observables by imposing a dependence of $`N`$ on the cutoff $`\mathrm{\Lambda }`$ according to the above formula. This description yields a sensible prediction for experiments where the average effect of a large number of subbands is measured. The dependence of the quasiparticle weight $`Z`$ on energy as the cutoff is sent towards the Fermi level is represented in Fig. 2, for different values of $`n_p`$. In the present instance, the dominant contribution to the power-law behavior of the tunneling density of states comes from the number of subbands varying with the energy scale. This sets the value of $`\alpha `$ close to 1 from the start, what is the natural way of recovering the characteristic linear density of states of a graphite layer from our 1D point of view.

We observe from the results in Fig. 2 that the quasiparticle weight $`Z`$ tends to have a flat behavior at high energies for large values of the density of subbands $`n_p`$. This is in contrast to the rapid decrease signaling the typical power-law behavior for small values of $`n_p`$. In the curves for $`n_p=100`$ and $`n_p=1000`$, we see the existence of a crossover between a regime with persistent quasiparticles and another characteristic of the Luttinger liquid behavior. The physical interpretation is that, for high energies above the crossover scale, the system has similar properties to the 2D graphene, while one has to look at sufficiently small energy scales (or large length scales) to measure the properties of the 1D wire.

We have obtained several results that may be checked against future measurements carried out in MWNT and nanotubes of very large radius. We predict that the exponent giving the power-law behavior of the tunneling density of states may suffer a significant reduction upon doping those systems, with the possibility of reaching values as small as 0.1 for $`N40`$. We have also shown that, when dealing with nanotubes of very large radius, there is a high-energy regime with persistent electron quasiparticles which has properties closer to 2D graphene than to the Luttinger liquid. We believe such features may be of interest when developing carbon-based devices made of graphene and nanotube structures with different shapes.

This work is partly supported by the Italian Research Ministry, National Interest Program.


text:
# References

Multiscale Algorithms for Eigenvalue Problems

Nimal Wijesekera1, Guogang Feng2, Thomas L. Beck2∗

1Department of Physics

2Department of Chemistry

University of Cincinnati

Cincinnati, OH 45221-0172

thomas.beck@uc.edu

\*corresponding author

Abstract

Iterative multiscale methods for electronic structure calculations offer several advantages for large-scale problems. Here we examine a nonlinear full approximation scheme (FAS) multigrid method for solving fixed potential and self-consistent eigenvalue problems. In principle, the expensive orthogonalization and Ritz projection operations can be moved to coarse levels, thus substantially reducing the overall computational expense. Results of the nonlinear multiscale approach are presented for simple fixed potential problems and for self-consistent pseudopotential calculations on large molecules. It is shown that, while excellent efficiencies can be obtained for problems with small numbers of states or well-defined eigenvalue cluster structure, the algorithm in its original form stalls for large-molecule problems with tens of occupied levels. Work is in progress to attempt to alleviate those difficulties.

Introduction

Electronic structure methods for large-scale problems can be divided into three general categories: plane-wave , traditional basis set , and real-space methods . Real-space methods result in a banded Hamiltonian, through finite difference , finite element , or wavelet  representations of the Laplacian operator. Gaussian basis sets lead to a smaller total matrix size of the Hamiltonian (relative to real-space methods), but the matrix is less banded. The plane-wave basis set, on the other hand, is completely delocalized in real space. The bandedness of the Hamiltonian in real-space methods is advantageous in several respects. First, it leads to ease of developing parallel codes. Second, methods developed recently which scale linearly with the system size generally rely on localization of the orbitals and real-space methods ‘mesh’ well with those algorithms . Third, efficient multiscale methods accelerate convergence by decimating errors over a wide range of length scales. Fourth, finite clusters or periodic systems can be treated with equal effort. Finally, local mesh refinements can be incorporated without degrading the efficiency of the solver . Several groups have developed real-space solvers for the Kohn-Sham equations of density functional theory (DFT) in the last decade . Multigrid (MG) methods have been employed extensively to accelerate the convergence rate .

A central feature of the Kohn-Sham problem is its nonlinearity. It is nonlinear in two respects. First, the eigenvalue problem itself is nonlinear since one solves for both the eigenvalues and eigenfunctions. Second, the self-consistent potential depends nonlinearly on the charge density obtained from the squares of the eigenstates. Therefore, it can be expected that nonlinear multigrid methods will lead to increased efficiencies. This has indeed been observed previously in studies which compare a full approximation scheme (FAS) nonlinear approach to linearized MG methods .

The most costly operation in the nonlinear eigenvalue approach of Brandt et al.  is the Ritz projection (preceeded by Gram-Schmidt orthogonalization) on the fine level. If $`q`$ orbitals span the whole physical domain, with $`N_g^h`$ grid points on the fine level (labelled by $`h`$), this step of the algorithm scales as $`q^2N_g^h`$. In the present paper, we carry the nonlinear FAS scheme further by implementing an extension of the Brandt et al. algorithm proposed by Costiner and Ta’asan . In their approach, the orthogonalization and Ritz projection are moved to coarse levels within the FAS strategy. This results in an 8-fold decrease in computational cost for the Gram-Schmidt/Ritz operation for each coarser level in three dimensions. We discuss relevant details of their method, and then present numerical results on fixed potential and self-consistent eigenvalue problems related to atomic and molecular structure.

Algorithms

Multigrid methods accelerate the convergence rate of iterative relaxation solvers for partial differential equations by decimating errors on a wide range of length scales ; it is the long wavelength modes of the error which degrade single-level relaxation efficiency. Nonlinear multigrid methods incorporate a full representation of the problem on coarse levels. Modifications of the fine grid matrix equation are necessary on the coarse grids to obtain zero correction at convergence. These modifications are termed defect corrections. In this paper, we examine the FAS multigrid method of Costiner and Ta’asan . They presented a detailed account of the algorithm in two papers: the first concerns fixed potential problems and the second addresses self-consistency. A brief summary will be given here. We will follow the notation from their work.

FAS Eigenvalue Method

Suppose we wish to solve an eigenvalue problem represented in real space with finite differences. This leads to the matrix equation

$$AU=U\mathrm{\Lambda }.$$
(1)
The matrix $`A`$ is the $`N_g^h\times N_g^h`$ (banded) Hamiltonian, $`U`$ is the $`q\times N_g^h`$ matrix of the eigenvectors, and $`\mathrm{\Lambda }`$ is the diagonal $`q\times q`$ matrix of eigenvalues.

In the FAS approach, we express the coarse grid (level $`i`$) problem as

$$F_iU_i=T_i.$$
(2)
Here $`F_iU_i=A_iU_iU_i\mathrm{\Lambda }_i`$ and $`T_i`$ is the defect correction. On the finest grid $`T_i=0`$. On grids $`j`$ coarser than $`i`$,

$$T_j=I_i^j(T_iF_iU_i)+F_jI_i^jU_i,$$
(3)
where the operator $`I_i^j`$ is the restriction operator. We use full weighting restriction throughout, which just involves a local trapezoid-rule integration of the function values from the fine grid. In the above equations, if the exact numerical fine grid solution is inserted into the coarse grid equations, identities are obtained. This is equivalent to the zero correction at convergence condition.

An initial approximation is first obtained on the fine grid. We obtain this approximation by implementing a full multigrid (FMG) cycle , beginning on the coarsest level, interpolating to the next finer grid, performing MG cycles there, and so on until the finest grid is reached. In this way, a good initial approximation is obtained for very little numerical effort. Following relaxation (typically 2-5 Gauss-Seidel or successive over-relaxation/SOR steps), the fine grid approximation is passed to the coarse level by restricting the eigenfunctions, the potential, and the defect correction from the finer grid. Relaxations (and generalized Ritz projections, see below) are performed on the current coarse grid, and the problem is then passed again to the next coarser level. This process is repeated until the coarsest grid is reached. We typically utilize three grid levels when the finest grid is a $`65^3`$ mesh. Once relaxation is done on the coarsest level, a correction step for the next finer level is performed:

$$U_i^{new}=U_i^{old}+I_j^i(U_jI_i^jU_i^{old}).$$
(4)
$`I_j^i`$ is the interpolation operator. Linear interpolation by lines is used throughout except during the FMG process when passing to the next finer level, where cubic interpolation by lines is used (to obtain a better initial guess on the fine grid). The correction steps are continued until the finest grid functions are corrected followed by relaxation steps there. The full cycle through all the levels is termed a V-cycle. One may then repeat the MG V-cycles until a desired convergence is obtained.

Generalized Ritz Projection, GRP

Consider a new eigenvalue relation, in which the matrix $`V`$ results from a linear combination of the current approximation $`U`$:

$$AV=V\mathrm{\Lambda },$$
(5)
where $`V=UE`$ and $`E`$ is a $`q\times q`$ matrix of normalized vectors, the columns of which determine the coefficients for the linear combination of old vectors from $`U`$. Then on the fine grid we have the relation

$$AUE=UE\mathrm{\Lambda }.$$
(6)
When this problem is passed to the coarser levels, the proper FAS transfer is

$$AUE=UE\mathrm{\Lambda }+TE.$$
(7)
If we multiply on the left by $`U^T`$, we obtain the following generalized (nonsymmetric) eigenvalue problem on the coarse grid:

$$U^T(AUT)E=(U^TU)E\mathrm{\Lambda }.$$
(8)
Notice that the eigenfunctions are no longer orthonormal when passed to the coarse grid. Also, it is easy to show that, at convergence, the eigenvalues are the same on all grid levels; in principle, there is no need to compute them on the finest grid. We solve this $`q\times q`$ eigenproblem with standard linear algebra packages. Once solved, we obtain new eigenvalues, eigenfunctions (linear combinations of the current approximation), and defect corrections (also linear combinations of the old defect corrections). On the finest grid, where the defect correction $`T`$ is zero, the GRP reduces to the usual Ritz projection employed in the Brandt et al. algorithm .

Backrotation, BR

A subtle aspect of the correction scheme outlined above is that the coarse-grid eigenfunctions must properly match their fine-grid counterparts for the correction step. Therefore, Costiner and Ta’asan  introduced a process called backrotation in their solver. This operation is designed to prevent rotations in degenerate or near-degenerate subspaces, and to prevent sign changes, rescalings, and permutations of the eigenvectors. In the backrotation, the $`E`$ matrix is modified towards the ends listed above. If this step is not employed in the algorithm of Ref. , the solver typically does not fully converge. As a simple example, imagine that the sign of one of the eigenfunctions changes during GRP on a coarse level. If the correction is then interpolated to the fine grid, the corrected function will be severely distorted. At convergence, the $`E`$ matrix should approach the unit matrix. During processing (prior to backrotation), it tends to have block diagonal form, where the blocks correspond to degenerate or near-degenerate subspaces. The dimensions of the blocks determine the eigenvalue cluster sizes. An extensive discussion of the backrotation is given in the original paper.

Relaxation

A major feature of MG methods is that relatively simple relaxation strategies can be employed so long as they decimate errors with wavelengths comparable to the grid spacing on a given level . Gauss-Seidel is the most common one utilized. We have investigated several relaxation strategies for smoothing on each level. Originally, we used the Gauss-Seidel with shift form given in Ref. . Recently, we have extended this relaxation method to an SOR form, and find improved convergence. On the coarsest level, we have employed Gauss-Seidel directly, Gauss-Seidel with constraints designed to maintain eigenfunction orthonormality on the fine grid , and Kaczmarz  relaxation. Kaczmarz relaxation is guaranteed to converge, but it exhibits slower convergence relative to Gauss-Seidel or SOR (this is not a significant issue on coarse levels). It will be specified below which method was used for each application. Further details of relaxation methods will be presented in an extensive account of our algorithm .

Self-Consistent Problems

Some of the applications presented below concern self-consistent solution of the Kohn-Sham equations. In the work presented here, we update the eigenfunctions and self-consistent potential sequentially. That is, given an initial approximation to the effective potential, an MG V-cycle is performed to update the eigenfunctions. From the new eigenfunctions, a new charge density is computed, from which a new effective potential is obtained. The main computational step for updating the effective potential is solution of a Poisson problem. This equation is also solved with MG V-cycles. The total time to solve the Poisson equation is less than that for updating a single eigenfunction, and this operation scales linearly with system size. We are currently exploring approaches to update the effective potential on coarse levels simultaneously with the eigenfunctions . This can be expected to accelerate the self-consistent convergence rate. We note that in our calculations so far, we have found no need for potential or charge density mixing of old and new solutions; we believe this is due to long-wavelength stabilization of the charge density during the FMG preconditioning process.

Pseudopotentials

For our calculations on atoms and molecules, we have incorporated the separable dual-space Gaussian pseudopotentials developed by Goedecker et al. . These pseudopotentials have analytic forms with only several parameters per atom, and they exhibit optimal decay properties in both real and reciprocal space. We have implemented the real-space relativistic version of these pseudopotentials for the present calculations. For calculations on the coarse grids, we simply compute the function values just as we do on the fine level. Of course fewer grid points are necessary to sample the pseudopotential on coarse grids due to the decay properties of the projectors. It has been shown that application of pseudopotentials in real space is more efficient than in reciprocal space .

Numerical Results

Fixed Potential Problems

We first present results of the FAS algorithm on fixed potential problems. As a benchmark, we utilized the original algorithm of Ref.  and solved the same two dimensional eigenvalue problem addressed in that paper. Their kinetic energy operator is scaled by a factor of two. A second order finite difference approximation was assumed for the Laplacian. Half of their potential is

$`v=5y\mathrm{sin}(3\pi x).`$ (9)
The total domain size was taken as one, and four grid levels were utilized in the FAS process. Gauss-Seidel relaxation (with a shift parameter of zero) was employed on all four levels. One relaxation step was performed on the coarsest level while enforcing constraints designed to maintain eigenfunction orthonormality on the finest level. On all other levels, two relaxation steps were performed. In total, 7 FAS V-cycles were implemented on the finest level. The computed eigenvalues and residuals are shown in Table 1. The residual for each eigenfunction is defined as

$`r=\sqrt{{\displaystyle \frac{\left|H\psi E\psi \right|^2}{N_g^h}}},`$ (10)
where the sum is over all the fine grid points and $`N_g^h`$ is the total number of grid points. The algorithm with Ritz projection performed on the finest level converges nicely to the numerically exact eigenfunction/eigenvalue pairs.

We next compare the convergence rates of three FAS algorithms: 1) Ritz on fine grid (Ref. , 2) generalized Ritz with backrotation (GRBR) on coarse grids with no orthonormalization on the fine level , and 3) GRBR on coarse grids with periodic (every 5 V-cycles) Gram-Schmidt orthonormalization on the fine grid. These approaches were tested on simple fixed potential problems, namely the three-dimensional harmonic oscillator and the hydrogen atom. Both of these physical problems have degenerate subspaces leading to eigenvalue clusters which must be handled in the backrotation step. For these problems, Gauss-Seidel relaxation was employed on all levels, except Kaczmarz relaxation was utilized for the coarsest grid relaxation for the hydrogen atom problem (this led to increased efficiency).

The convergence results are presented in Figs. 1 and 2. For the harmonic oscillator problem, we solved for ten states and used a 12th-order finite difference representation for the Laplacian. The total domain size is ten, and the fine grid is a $`65^3`$ mesh; three grid levels were utilized. Atomic units are used throughout. A shift parameter equal to the current eigenvalue was employed in the Gauss-Seidel relaxations (5 steps per level). The GRBR step was performed on the coarsest level. A total of 15 V-cycles were conducted. The compute times per V-cycle for the three algorithms listed above were 50.9, 34.6, and 36.6 sec., respectively, on an 800 MHz machine. The $`q`$ values are small enough that the orthogonalization and Ritz projections do not yet dominate the compute time. For the hydrogen atom case, we generated the fixed potential by numerical solution of the Poisson equation for a fixed central unit charge. Again, 12th-order finite differences were employed, and we solved for 14 states on a domain with a total side length of 28. The fine grid is a $`65^3`$ mesh. For this case, it was found that a shift parameter equal to the potential was more efficient. A total of 20 V-cycles were conducted. For both cases, the fastest (lowest eigenvalue) and slowest (highest eigenvalue) converging cases are shown. The GRBR (without orthonormalization) operation was conducted on the middle grid level, while for the third algorithm with periodic Gram-Schmidt operations on the fine grid, the GRBR step was implemented on the coarsest level. The compute times per V-cycle were 81.5, 59.7, and 52.8 sec. for the three algorithms.

For both physical problems, performing the orthonormalization and Ritz projection on the fine level leads to the most efficient convergence. The harmonic oscillator potential is smooth, and the GRBR algorithm exhibits good convergence behavior. However, the convergence rate is slightly slower than when orthonormalization and Ritz projection are performed on the finest level. The GRBR convergence rate is the same whether or not orthonormalization is periodically done on the finest level, indicating good separation can be obtained without processes on the fine level. For the hydrogen atom case, the overall convergence rate is reduced relative to the harmonic oscillator, presumably due to the singular nature of the potential. Also, periodic Gram-Schmidt operations on the finest level increase the convergence rate slightly.

As part of the backrotation step, the degenerate subspaces (eigenvalue clusters) must be identified so as to prevent rotations within those clusters. Both the harmonic oscillator and the hydrogen atom problems possess clear eigenvalue structure which can be directly incorporated or determined during the solution process. This ensures healthy convergence of the GRBR algorithm, behavior which was also observed in the original work of Costiner and Ta’asan . We will see below that, for large self-consistent molecular cases without such well-defined eigenvalue cluster structure, the algorithm may stall due to mixing during the backrotation step.

Self-Consistent Pseudopotential Calculations

Computational results for self-consistent pseudopotential Schrödinger-Poisson eigenvalue problems are presented. Just as for the fixed potential problems presented above, the algorithms used were Ritz projection on the fine level coupled with Gram-Schmidt orthogonalisation and the GRBR algorithm without and with periodic fine grid Gram-Schmidt operations. In both GRBR algorithms, the eigenfunctions were normalized on the fine level to ensure charge conservation. Three examples, Ne, CO, and the benzene dithiol molecule were used in the study. All three cases are three dimensional and were treated with 12th order finite difference representations. Three grid levels were utilized comprising $`17^3,33^3`$ and $`65^3`$ total points . The first example, the Ne atom, possesses 4 occupied states. Similarly 5 and 21 eigenvectors were required for the CO and benzene dithiol molecules, respectively. The Ne and CO examples possess well-defined eigenvalue cluster structure (triply degenerate p states for Ne and a doubly degenerate $`\pi `$ bonding orbital for CO). Convergence results are presented in Figs. 3-5.

Choosing the optimal parameters for the relaxation scheme was important. The shift parameter for Gauss-Seidel relaxation was taken as $`\lambda +v`$, where $`\lambda `$ is the eigenvalue and $`v`$ is symbolically the effective potential (nonlocal in the case of the pseudopotential). Gauss-Seidel relaxation ($`\omega =1`$, where $`\omega `$ is the overrelaxation parameter) was used on all levels except on the finest grid where SOR relaxation was employed ($`\omega =1.7`$). These near-optimal relaxation parameters were determined by numerical experiments. In the final V-cycle involving the three grid levels, three relaxation steps were done on each level.

All three algorithms converged for the Ne and CO cases. The fine grid Ritz plus Gram-Schmidt and coarse-grid GRBR with periodic fine level Gram-Schmidt algorithms both exhibited excellent convergence rates. Implementing GRBR with no fine grid separation (orthogonalization) led to slower convergence; the total energy convergence slows after a few self-consistency iterations. In the GRP, the fine level separation of wavefunctions comes as a result of GRBR done on the middle level. But complete separation may not always be achieved on the fine level (Table 2). Performing only a few fine grid Gram-Schmidt operations during the entire solution process restores the convergence rate and leads to acceptable fine grid eigenfunction orthogonality at the end. Clearly the form of the potential (self-consistent pseudopotential in this case) and resulting eigenvalue structure have impacts on the convergence behavior of the GRBR algorithm.

In the case of benzene dithiol, the fine grid separation algorithm converged nicely as for the smaller molecules. However, neither variant of the GRBR algorithm converged properly; the solver stalled with only modest energy convergence. To probe for the reason for this lack of convergence, we first converged the system using fine-grid separation for 15 V-cycles, and then used the resulting potential in a fixed potential calculation. We observed that the $`E`$ matrix determined from GRBR gradually begins to destablize rather than to converge to the unit matrix as expected. This suggests again that the form of the potential and the resulting eigenvalue cluster structure affect the convergence of the GRBR algorithms. Benzene dithiol is a relatively large molecule with 14 atoms (C,O,S and H) and 21 valence eigenstates. We believe that the difficulties arise from the ambiguous cluster structure of the eigenvalues which leads to mixing during the backrotation operation. Similar calculations were performed on the benzene molecule, which possesses clear symmetry and cluster structure, and the algorithm converged. Since determination of the cluster structure is crucial for convergence, and this must be performed automatically in the algorithm in order to treat general systems, this issue must be addressed for the GRBR approach to provide a generally convergent scheme. Work is in progress investigating these difficulties. One possible solution is suggested below.

Discussion

The objective of this paper has been to provide a test of the nonlinear FAS multigrid eigenvalue method of Refs.  for solving the Kohn-Sham equations. This method is promising because it removes the expensive orthogonalization and Ritz projection operations to coarse levels (with a corresponding 8-fold reduction in cost per level in three dimensions). The necessary fine grid work only involves relaxation, normalization of the eigenfunctions, and perhaps orthogonalization within degenerate clusters. The model problems treated in Refs.  possess relatively smooth potentials and well-defined eigenvalue cluster structure. Similar to the results of Costiner and Ta’asan, we find good convergence of the GRBR approach for fixed potential and self-consistent eigenvalue problems with well-defined eigenvalue clusters. However, we found that for a larger molecular case with complicated eigenvalue structure, the GRBR approach did not converge properly. We linked these difficulties to the backrotation step which appears to be highly sensitive to the determination of the clusters. Since determination of the clusters must be done numerically during the course of the solution process, this issue must be addressed to develop a generally convergent solver for large systems.

We have recently begun investigating one possible approach to deal with these difficulties. Notice that the $`E`$ matrix of Eqs. 6 and 7 is formally the same on all levels at convergence (just as are the eigenvalues). In the algorithm of Ref. , the backrotation involves a modification of the $`E`$ matrix designed to prevent rotations in degenerate or near-degenerate subspaces, sign changes, rescalings, and permutations of eigenvectors. Utilizing the fact that the $`E`$ matrix is the same on the coarse and fine levels (at convergence), we can circumvent the backrotation by first applying the $`E`$ matrix to the current fine-grid approximation to the eigenfunctions prior to the correction step. This approach is in a sense a hybrid of the fine and coarse grid Ritz projections; we use the coarse grid to generate the new eigenvalues and the $`E`$ matrix, but we use that $`E`$ matrix to alter the fine-grid occupied subspace. Thus the expensive step of constructing the Ritz matrix has been moved to the coarser level. The use of the $`E`$ matrix to update the fine-grid function is a relatively cheap operation. Formally, it does scale as $`q^2N_g^h`$, but the $`E`$ matrix is of block diagonal form, with the blocks of dimension of the corresponding degenerate cluster. These clusters are generally very small. Therefore realistically the update operation scales as $`qN_g^h`$ if the eigenfunctions span the whole domain. Some discussion along these lines was already given in Ref. . We are currently exploring the use of this idea in our nonlinear FAS multigrid solver. In preliminary results, we have found it to converge properly for all of the physical problems examined in this paper.

Acknowledgments

We gratefully acknowledge the support of the National Science Foundation (CHE-0112322) for this research. We also thank Achi Brandt and Shlomo Ta’asan for many helpful discussions.