text:
# Assessment of exchange-correlation functionals for the calculation of dynamical properties of small clusters in TDDFT

## I Introduction

Density Functional Theory (DFT) became, in the two last decades, the method of election for the ab initio calculation of material properties. Systems with a couple thousand atoms are now routinely investigated, and the calculated energies, geometries, etc., very often agree spectacularly with the experimental data. This success, however, was not immediate: More than twenty years of research were necessary to obtain an exchange-correlation (xc) energy functional which was precise enough to satisfy the quantum-chemistry needs. The accurate and sophisticated Generalized Gradient Approximations (GGA) that we now have at our disposal are indeed the result of a long history of attempts, tests and failures.

Despite these remarkable achievements, there are some quantities which are beyond the reach of the conventional DFT theory. In this Article we will be concerned by one of such properties, namely electronic excitation energies and optical spectra. Several extensions of the original framework were put forward to obtain excited-state properties. From all different approaches, Time-Dependent DFT (TDDFT) emerged the last years as the main DFT formalism for these calculations. One of the main advantages of the TDDFT formulation is that it allows one to deal with problems beyond the perturbative regime as, for example, the response of atoms to ultra-short and intense laser beams. In this way, by solving the time-dependent Kohn-Sham (KS) equations, we can get, e.g., information concerning the harmonic spectrum or the ionization yields. This type of approach is very important in condensed matter science where femtosecond laser pulses are used to monitor the dynamics of electrons in a solid. Furthermore, with larger computational power, we may deal with non-adiabatic couplings between the electronic and the ionic degrees of freedom in presence of these high electromagnetic fields (we may observe phonon assisted structural transformation induced by these external fields).

Excitation energies can be obtained from TDDFT either from the position of the poles of the KS linear response function , or from the time-dependent density. In the second case, we will be using an approximated xc time-dependent potential, $`v_{\mathrm{xc}}(𝐫,t)`$, while in the former, the key quantity is the xc kernel, $`f_{\mathrm{xc}}(𝐫,𝐫^{},t,t^{})`$, defined by the functional derivative:

$$f_{\mathrm{xc}}(𝐫,𝐫^{},t,t^{})=\frac{\delta v_{\mathrm{xc}}(𝐫,t)}{\delta n(𝐫^{},t^{})}.$$
(1)
As expected, the xc kernel, due to its complicated structure, is much harder to model than the local time-dependent xc potential. Furthermore, both calculations usually start from the KS ground-state of the system, which depends on the approximation used for the static $`v_{\mathrm{xc}}(𝐫)`$.

One of the most widely used approximation to time-dependent phenomena is the adiabatic local density approximation (ALDA). Although this functional is constructed using the xc correlation energy of the (constant density) electron gas, it yields, like the LDA does in the static case, rather accurate results for systems with rapidly varying densities, as atoms, surfaces or clusters . For example, the photo-absorption spectrum of rare-atom gases was computed in Ref.  and the agreement with experiment is remarkably good. Results of similar quality have been achieved for the photo-response of small metallic and semiconducting clusters . For the metallic clusters it was shown that the inclusion of xc effects in the dielectric response was important to get the correct red-shift, as compared to experiments, while for the semiconducting (silicon) clusters it was found that the spectra of different isomers was sufficiently different to distinguish between them. Surface and confinement effects were responsible for the appearance of absorption in the optical range for the two silicon clusters, Si4 and Si6. Similar results were found for C20 clusters .

The purpose of this work is twofold: To address the impact of a good xc-potential in the optical spectrum of small clusters and to estimate how relevant the $`f_{\mathrm{xc}}`$ kernel is in this respect. We note that in Ref. the quality of different functionals was tested for light atoms, and in Ref. for the calculation of the correlation energy of the homogeneous electron gas. In the first of these works, which included calculations using the exact static xc potential, the correct description of the xc-potential used to obtain the ground state seemed to be the dominant factor, whereas the kernel played a marginal role. However this was not the case when looking at the total correlation energy of the homogeneous electron gas. In particular, the nonzero spatial range of $`f_{xc}(r,r^{},\omega )`$ could not be neglected, whereas the frequency dependence appeared to be less important. Our work clearly demonstrates that the optical spectrum of small systems is not only determined by the static xc potential, but also the kernel plays a role.

The rest of the Article is structured in the following way: In Section II, we give a brief overview of the xc functional zoo used in static KS calculations; We then proceed by explaining how these functionals can be changed in order to be used in TDDFT calculations; In Section IV we give some details on the techniques used in our calculations, and present our results for some selected sodium and silicon clusters; We finally conclude and give some remarks on the quality of the functionals tested.

## II XC Functionals for Static DFT

The first (and simplest) approximation for the xc functional, the Local Density Approximation (LDA) introduced by Kohn and Sham in 1965 , yielded remarkably good results for such a modest effort, but the method didn’t prove accurate enough to be used in theoretical chemistry. The LDA xc energy functional can be written as

$$E_{\mathrm{xc}}^{\mathrm{LDA}}=d^3rn(𝐫)ϵ_{\mathrm{xc}}^{\mathrm{hom}}(n(𝐫)),$$
(2)
where $`n(𝐫)`$ is the electronic density at point $`𝐫`$, and $`ϵ_{\mathrm{xc}}^{\mathrm{hom}}(n)`$ stands for the xc energy density of the homogeneous electron gas with density $`n`$. This function is calculated very accurately using Monte-Carlo techniques , and then fitted to some simple analytical function . The functional (2) has several shortcomings, among which we point out the following:

* It neglects non-local effects, i.e. the LDA xc energy density at point $`𝐫`$ only depends on the density at that point. We therefore should not expect this functional to work in cases where the density has very strong spatial variations.
* The exchange part of the functional does not cancel exactly the self-energy part of the Hartree term. This leads to a wrong asymptotic behavior of the xc potential for finite systems (it goes exponentially to zero, instead as $`e^2/r`$). Properties that strongly depend on this asymptotic behavior, like the ionization potential of atoms and molecules, come out with very large errors. Also, there are no Rydberg states within the LDA, and negative ions usually do not bind, which renders the calculation of electron-affinities impossible.
* LDA usually overbinds, giving too short bond-lengths, etc.

The next generation of functionals included the so-called Generalized Gradient Approximated (GGA) functionals . They can be written as:

$$E_{\mathrm{xc}}^{\mathrm{GGA}}=d^3rn(𝐫)ϵ_{\mathrm{xc}}^{\mathrm{GGA}}(n(𝐫),n(𝐫)),$$
(3)
where $`n(𝐫)`$ is the gradient of the density at the point $`𝐫`$. $`ϵ_{\mathrm{xc}}^{\mathrm{GGA}}`$ is usually some analytic function with some free parameters that are either fitted to experiment, or determined by some exact sum-rules. GGAs solve some of the problems present in the LDA, and in some cases yield results with high enough precision to be used as a tool in quantum chemistry calculations. Recently, a new class of functionals generalizing the GGAs has been proposed. These so-called Meta-GGAs , depend explicitly not only on the density and its gradient, but also on the kinetic energy density The MGGAs, due to the explicit dependence on the kinetic energy density are, in fact, not simple GGAs, but orbital functionals. . This extra dependence adds more flexibility, and allows better approximations to the exact xc functional to be built. Unfortunately, all the GGAs and MGGAs proposed till now suffer, to different degrees, from the same self-interaction problem as the original LDA. A way to circumvent it was proposed by van Leeuwen and Baerends in 1994 : They applied Becke’s construction  not to the derivation of the exchange energy functional (like in Becke’s original work), but to the modeling of the xc potential directly. By imposing the correct asymptotic behavior to the potential, they were able to get much better ionization potentials (and eigenvalues in general).

To complete the zoo of energy functionals, we still have to refer to orbital functionals, in the so-called Optimized Effective Potential (OEP), or Optimized Potential Method (OPM) . This third generation of functionals is written explicitly in terms of the Kohn-Sham orbitals (being nevertheless implicit functionals of the density):

$$E_{\mathrm{xc}}^{\mathrm{OEP}}=E_{\mathrm{xc}}^{\mathrm{OEP}}\left[\phi _1(𝐫)\mathrm{}\phi _N(𝐫)\right].$$
(4)
The xc potential is then calculated using twice the chain-rule for functional derivatives:

$`v_{\mathrm{xc}}^{\mathrm{OEP}}(𝐫)={\displaystyle \frac{\delta E_{\mathrm{xc}}^{\mathrm{OEP}}}{\delta n(𝐫)}}`$ (5)
$`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle }d^3r^{}d^3r^{\prime \prime }[{\displaystyle \frac{\delta E_{\mathrm{xc}}^{\mathrm{OEP}}}{\delta \phi _i(𝐫^{})}}{\displaystyle \frac{\delta \phi _i(𝐫^{})}{\delta v_{\mathrm{KS}}(𝐫^{\prime \prime })}}+\mathrm{c}.\mathrm{c}.]{\displaystyle \frac{\delta v_{\mathrm{KS}}(𝐫^{\prime \prime })}{\delta n(𝐫)}}.`$ (6)
The first term of the right, $`\frac{\delta E_{\mathrm{xc}}^{\mathrm{OEP}}}{\delta \phi _i(𝐫^{})}`$, can easily be obtained from the expression for $`E_{\mathrm{xc}}^{\mathrm{OEP}}\left[\phi _1(𝐫)\mathrm{}\phi _N(𝐫)\right]`$, while the second, $`\frac{\delta \phi _i(𝐫^{})}{\delta v_{\mathrm{KS}}(𝐫^{\prime \prime })}`$ can be calculated using first-order perturbation theory. Finally, the remaining functional derivative can be identified with the inverse of the response function for non-interacting electrons, $`\chi _{\mathrm{KS}}^1(𝐫,𝐫^{})`$. Rearranging the terms in Eq. (5), we obtain an integral equation for $`v_{\mathrm{xc}}^{\mathrm{OEP}}`$. The solution of this integral equation is numerically very involved, and has only been achieved for systems with very high symmetry. An alternative is to perform an approximation first proposed by Krieger, Lee and Iafrate (KLI), which transforms the hard task of solving an integral equation in three dimensional space into the simple one of solving a small set of linear equations.

Two examples of these orbital functionals are the exact exchange (EXX), and self-interaction corrected (SIC) LDA functionals . In the first, one employs the exact expression for the exchange energy:

$$E_\mathrm{x}^{\mathrm{EXX}}=\frac{1}{2}\underset{\sigma }{}\underset{i,j=1}{\overset{N}{}}d^3rd^3r^{}\frac{\phi _{j\sigma }^{}(𝐫)\phi _{i\sigma }^{}(𝐫^{})\phi _{i\sigma }(𝐫)\phi _{j\sigma }(𝐫^{})}{\left|𝐫𝐫^{}\right|}.$$
(7)
The SIC-LDA functional, originally proposed by Perdew and Zunger , can be written as:

$`E_{\mathrm{xc}}^{\mathrm{SIC}}`$ $`=`$ $`E_{\mathrm{xc}}^{\mathrm{LDA}}[n_{}(𝐫),n_{}(𝐫)]{\displaystyle \underset{\sigma }{}}{\displaystyle \underset{i=1}{\overset{N}{}}}E_{\mathrm{xc}}^{\mathrm{LDA}}[\left|\phi _i(𝐫)\right|^2,0]`$ (9)
$`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\sigma }{}}{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle d^3rd^3r^{}\frac{\left|\phi _{i\sigma }(𝐫)\right|^2\left|\phi _{i\sigma }(𝐫^{})\right|^2}{\left|𝐫𝐫^{}\right|}}.`$
It is clear that the SIC-LDA obeys two of the features of the exact functional: It exactly cancels the self-interaction part of the Hartree energy and it vanishes for one electron systems. The SIC-LDA functional is nevertheless ill-defined, for it is not invariant upon an unitary transformation of the Kohn-Sham wave-functions.

## III XC Functionals for TDDFT

The simplest, and most commonly used, approximation to the xc functional in TDDFT is the adiabatic LDA (ALDA), in which the static LDA functional is used for the dynamical properties, but evaluated at the time-dependent density:

$$v_{\mathrm{xc}}^{\mathrm{ALDA}}(𝐫,t)=v_{\mathrm{xc}}^{\mathrm{hom}}(n(𝐫,t)).$$
(10)
In the ALDA, the $`f_{\mathrm{xc}}`$ kernel is a contact function in time and space:

$$f_{\mathrm{xc}}^{\mathrm{ALDA}}(𝐫,t;𝐫^{},t^{})=\delta (tt^{})\delta (𝐫𝐫^{})\frac{dv_{\mathrm{xc}}^{\mathrm{hom}}(n)}{dn}|_{n=n(𝐫,t)}.$$
(11)
Following the same reasoning, it is straightforward to write adiabatic GGA potentials and $`f_{\mathrm{xc}}`$ kernels. Unfortunately, the thresholds for the onset of absorption calculated either with adiabatic LDA or GGA functionals are typically below the observed ones (by several eV in the case of atoms). This is an intrinsic drawback of the ALDA/GGA because, in principle, TDDFT should yield the correct thresholds. This problem is once more related to the wrong asymptotic behavior of the effective Kohn-Sham potential that goes exponentially to zero instead of $`e^2/r`$ for neutral systems. This is due, as already mentioned in the previous Section, to the insufficient correction of the self-interaction part of the Hartree potential. Also the xc kernel, $`f_{\mathrm{xc}}`$, suffers from a self-interaction error. A simple way to correct the asymptotic part of the adiabatic potential is using the Adiabatic LB94 approximation. However, the high lying excitation energies calculated with this functional are usually overestimated for small molecules, and the LB94 does not do quite as well as the ALDA for the low-lying states .

A generalization of the EXX potential to the time-domain is also possible . It starts from the perturbative expansion of the action functional and then uses the chain-rule for functional derivatives to derive an integral equation for the xc potential. A KLI-like approximation can then be used to simplify the task of solving this integral equation. Using the same procedure, one can derive an integral equation for $`f_{\mathrm{xc}}`$. However, for practical calculations, it is desirable to devise a simple analytical expression for the xc kernel. To this end, it is common practice to further simplify the KLI potential by neglecting one of its terms . In the static exchange-only case, this leads to the so-called Slater approximation . $`f_{\mathrm{xc}}`$ then reads:

$$f_{\mathrm{xc}\sigma \sigma ^{}}^{\mathrm{EXX}\mathrm{approx}.}(𝐫t,𝐫^{}t^{})=\delta (tt^{})\delta _{\sigma \sigma ^{}}\frac{\left|_{i=1}^N\phi _{i\sigma }(𝐫)\phi _{i\sigma }^{}(𝐫^{})\right|^2}{\left|𝐫𝐫^{}\right|n_\sigma (𝐫)n_\sigma (𝐫^{})}.$$
(12)
One should note that, due to the extra functional derivative in the definition of the xc kernel, it is much more complicated, in this case, to evaluate directly $`f_{\mathrm{xc}}`$ than to time-propagate the KS states with $`v_{\mathrm{xc}}`$.

Finally, an adiabatic SIC-LDA can be easily written in the spirit of Eq. (9).

Other functionals designed to retain features which are important in time-dependent calculations have appeared in Quantum-Chemistry literature over the past years. A brief account of these can be found in Ref. .

## IV Method and Results

We used a real-space, real-time approach to solve the TD Kohn-Sham equations for $`\mathrm{Na}_2`$, $`\mathrm{Na}_4`$, silane ($`\mathrm{SiH}_4`$) and disilane ($`\mathrm{Si}_2\mathrm{H}_6`$). This time-evolution method is adequate to be combined with real-space calculations of the ground-state by means of the finite-difference pseudopotential method, or with adaptative coordinates. In both cases, a real-space discretization of the kinetic energy operator leads to sparse Hamiltonian matrices which do not need to be stored in memory, and are easy to handle. Furthermore, in order to propagate the states in time, we do not need to compute the complicated $`f_{\mathrm{xc}}`$ kernel, but only the much simpler $`v_{\mathrm{xc}}`$ potential. In this way, we are able to do time-dependent EXX calculations without having to perform Slater’s approximation .

The space was uniformly discretized, with mesh-spacing ranging from $`0.4`$Å for the sodium clusters to $`0.25`$Å for the silicon ones, and the points were contained inside a sphere with radius $`10`$Å for sodium and $`7`$Å for silicon. The ionic potential was modeled by a soft-core Troullier-Martins pseudopotential. By using these settings, a convergence of better than $`0.1`$ eV was reached for the eigenvalues. The time-evolution was performed with a modified Krank-Nicholson rule, with a time-step of 0.0025 $`\mathrm{}`$/eV, for a total evolution time, $`T_{max}`$, of 50 $`\mathrm{}`$/eV for the sodium clusters and 25 $`\mathrm{}`$/eV for the silicon clusters. With this parameters we can resolve excitation energies within $`\mathrm{\Delta }E=2\pi \mathrm{}/T_{max}`$ that corresponds to tenths of eV. We added absorbing boundary conditions to improve the quality of the spectrum above the ionization threshold (avoiding artificial formation of standing waves in the spherical box used to confine the cluster).

To calculate the dipole strength function (which is simply proportional to the absorption cross-section), we prepare our system in the ground-state, and then excite it with a delta electric field, $`E_0\delta (t)`$. The dipole strength can then be related to the imaginary part of the dynamical polarizability:

$$S(\omega )=\frac{4\pi m_e}{h^2}\omega \mathrm{}\alpha (\omega ),$$
(13)
where $`h`$ is Planck’s constant, $`m_e`$ the electron’s mass, and the dynamical polarizability,

$$\alpha (\omega )=\frac{2}{E_0}𝑑rz\delta n(𝐫,\omega ).$$
(14)
In the last expression, $`\delta n(𝐫,\omega )`$ stands for the Fourier transform of $`n(𝐫,t)n(𝐫,t=0)`$. The calculation of the dipole strength involves two important approximations: The choice of the static xc potential used to obtain the ground-state, and the time-dependent xc potential (directly related to the $`f_{\mathrm{xc}}`$ kernel) used to propagate the state. The relative importance of both is still fairly obscure. Petersilka et al. linear-response results for the Helium and Beryllium atoms indicates that the choice of the $`f_{\mathrm{xc}}`$ kernel is less dramatic than the choice of a good static $`v_{\mathrm{xc}}`$ potential. However, this is no longer true for the lowest excitation energies of Beryllium and if we look at the singlet-triplet splittings, where the effects arising from $`v_{\mathrm{xc}}`$ cancel. The fairly accurate results of Lithium and Beryllium within the LDA were related to the quality of the ground state calculation for few electron systems as light atoms and the $`H_2`$ molecule.

By performing the time evolution, in the limit $`t\mathrm{}`$, we have all the information about the linear response of the system, and should therefore be able to study the xc kernel without having to perform explicitly the functional derivative of the xc potential. Thus, we can study the effect of different $`f_{\mathrm{xc}}`$ kernels for a fixed $`v_{\mathrm{xc}}`$ potential in the time-evolution method. To address this question, we propose the following method: First, we obtain the ground state of the system with the static functional $`v_{\mathrm{xc}}^{(1)}`$, and then perform the time evolution with the xc potential:

$$v_{\mathrm{xc}}(𝐫,t)=v_{\mathrm{xc}}^{(1)}(𝐫,t=0)+\left[v_{\mathrm{xc}}^{(2)}(𝐫,t)v_{\mathrm{xc}}^{(2)}(𝐫,t=0)\right],$$
(15)
where $`v_{\mathrm{xc}}^{(2)}`$ is some other xc functional. In this way, we probe the importance of the $`f_{\mathrm{xc}}^{(2)}`$ kernel, without having to spend more computational time in calculating this rather complex quantity. This method will be used for silane, for a varied selection of $`v_{\mathrm{xc}}^{(1)}`$ and $`v_{\mathrm{xc}}^{(2)}`$ to quantify the dependence of the dynamical spectrum with the xc-potential and kernel.

In the following, PZ stands for LDA, with $`ϵ_\mathrm{c}^{\mathrm{hom}}`$ taken from Monte-Carlo calculations  and then parameterized by Perdew and Zunger ; PBE is the Perdew, Becke and Ernzerhof GGA functional ; LB94 represents the van Leeuwen and Baerends potential ; SIC is the SIC-LDA functional ; and, finally, EXX stands for exact-exchange . SIC and EXX are both treated within the KLI approximation . “exp” will be used to denote experimental values. All energies are measured in eV unless otherwise stated.

### A Sodium

In Table I we show the ionization potentials (IP), obtained from the highest occupied KS eigenvalue, and HOMO-LUMO gap for the two sodium clusters. The IPs obtained either by using the LDA or the GGA functionals are much smaller than the experimental value. This is a well known problem related to the wrong asymptotic behavior of the xc potentials. The situation does improve enormously if we correct for this deficiency, as can be seen from the LB94, SIC and EXX values. Although the calculated IPs differ by more than $`2`$ eV, the HOMO-LUMO gap stays essentially constant, i.e, correcting for the asymptotic part of the exchange potential amounts to a constant shift in the eigenvalues, and not to an opening of the gap.

All the calculated optical spectra for Na2 are quite similar (see Fig. 1), regardless of the xc potential used, and exhibit three clear peaks in the 2-5 eV range (the third peak in the LB94 curve is almost completely smeared out). They all compare quite well with experiment, although the DFT peaks are all shifted towards higher energies. This shift can be understood in terms of the competition between the Coulomb repulsion in the electronic kernel and the electron/hole attraction from the xc part of the response (that is, $`\mathrm{\Delta }ϵ=ϵ_cϵ_v+<\phi _v\phi _c|\frac{1}{|r_1r_2|}|\phi _v\phi _c>+\mathrm{\Delta }_{xc}`$). In particular for a given valence-conduction transition $`vc`$, the exchange correlation kernel within the simple ALDA introduced only a local and static attractive electron-hole interaction:

$`\mathrm{\Delta }_{xc}`$ $`=`$ $`{\displaystyle d^3r_1d^3r_2\phi _v^{}(𝐫_1)\phi _c^{}(𝐫_1)\delta (𝐫_1𝐫_2)\frac{V_{xc}(𝐫_1)}{n}\phi _v(𝐫_2)\phi _c(𝐫_2)}`$ (16)
$`=`$ $`{\displaystyle d^3r\phi _v^{}(𝐫)\phi _v(𝐫)\frac{V_{xc}(𝐫)}{n}\phi _c^{}(𝐫)\phi _c(𝐫)}.`$ (17)
This expression is more complicated for the other kernels, but it is clear that they do introduce an effective attractive interaction that it is not complete. We might have to recall that we have neglected dynamical effects and they may need to be included in the xc kernel to recover this minor effect. Also, we should keep in mind that temperature effects (vibrational motion) of the molecule may introduce a broadening of the spectrum as well as a shift of the peaks to lower frequencies . The functional which yields the best results for the dimer is, by a small margin, the EXX, while the strongest depart from experiment is the LB94 curve. We note that the 3rd peak in the LDA and PBE spectra lies in the continuum of states (above the ionization threshold). Thus it is more a resonance than a well defined bound transition as observed in experiments. This deficiency is observed in all the clusters studied in the present work and it is, once more, mainly due to the wrong description of the asymptotic part of the exchange potential.

For Na4 (Fig. 2) all DFT calculations yield very similar spectra (and similar to the many-body results based on a $`GW`$ quasiparticle calculation and including the electron-hole interaction through the solution of the Bethe-Salpeter equation). The spectra consist of three well separated peaks in the 1.5-3.5 eV range, and a broader feature at around 4.5 eV. The comparison with the experimental values is quite good, although the peaks appear, once more, at higher energies. The best results were obtained using the LDA, the GGA and the EXX, while the LB94 spectrum showed the strongest deviation from experiment. Again, the error is larger for the high energy peaks, where transitions to conduction states close to the ionization threshold are involved. Due to the small overlap between these Rydberg-like states and the low-lying states, exchange processes are not relevant and the polarization contribution to correlation, although being generally weak becomes the dominant contribution to the renormalization of these single particle excitations.

Finally, in Table II, we present the static polarizabilities of Na2 and Na4, obtained either through a finite field (static) calculation, or from the Fourier transform of the time-dependent density. For Na2 all functionals perform equally well, being the results smaller than the experimental value by around $`10\%`$. This is consistent with other DFT calculations that give static polarizabilities in the range 33.1-38.2 Å3 for Na2 and 67.1-78.7 Å3 for Na4 (see Ref.  and results cited therein). In the case of Na4, the orbital functionals SIC and the EXX yield slightly better static polarizabilities than the LDA and the GGAs, but the results are still smaller than the experimental value. We also rejoice that the two methods used to calculate the static polarizability yield very similar results. The neglect of correlation as well as temperature effects is responsible for the obtained smaller polarizabilities values. We indicate that the simple argument that a potential with the correct asymptotic will lead to more localized charge and therefore lower polarizability does not hold for these small systems. These results give support to previous studies on simple metallic jellium spheres.

### B Silicon

We have calculated the two simplest hydrogen terminated silicon clusters: silane (SiH4) and disilane (Si2H6). These systems pose a much harder challenge than the simple sodium clusters (where even jellium calculations within the LDA yield reasonable results), not only because of the presence of the p-electrons, but also due to the hydrogen, which is very hard to describe by a reasonably soft-core pseudopotential. As expected, the five xc functionals we tried yielded quite dissimilar results.

Our calculations for the IP and HOMO-LUMO gaps for silane and disilane are summarized in Table III. Undoubtfully, EXX and LB94 gave the best IPs of all five functionals tested, almost at the level of the much more involved GW with exciton effects or Monte-Carlo calculations. SIC is slightly worse, and LDA and PBE yield, as usual, completely irreal IP. Although the IP changes by more than $`4`$ eV going from the LDA to the EXX, the HOMO-LUMO gap increases by just around $`0.5`$ eV: The main difference is, once more, a nearly rigid shift of the eigenvalues. However, this shift brings the relevant single particle transitions below the ionization threshold. The results are also consistent with the finding that EXX calculations provide larger gaps for semiconductors.

We will discuss separately the silane and disilane spectra. The silane spectrum (Fig. 3) consists of three peaks between 8 and 12 eV (the two peaks derive from a Jahn-Teller splitting of the triply degenerate $`2\mathrm{t}_24\mathrm{s}`$ transition), followed by a much broader feature at higher energies. The curves obtained with the two traditional functionals, LDA and PBE, are quite similar to each other, and the onset for absorption is underestimated by around 1 eV. The underestimation of the onset of absorption is a well known deficiency of LDA based functionals that is even more dramatic in the case of infinite bulk systems where excitonic effects and band-gap renormalization are not properly described by these simple functionals. The SIC spectrum is unphysically shifted to lower energies (fact that could be anticipated by looking at the low SIC HOMO-LUMO gap), and the first peak is split. This second effect is an artifact of SIC, which, in this implementation, spontaneously breaks the 3-fold degeneracy of the HOMO state. LB94 and EXX behave quite well: the onset for absorption is now correct, and the error in the position of the first three peaks is reduced by a factor of 2 from the LDA/PBE results (see Table V). van Leeuwen’s functional performs marginally better in this case than EXX, which overestimates slightly the excitation energies. The transitions close to the ionization threshold give rise to a peak in the spectrum that in experiments is usually broader than in the calculations, as new decaying channels are available that are not included in the present calculations. Although the functionals we used do not properly include the damping of the excitations, finite lifetimes can be simulated by convoluting the calculated spectrum with some Lorentzian function.

The disilane LDA, PBE, LB94 and SIC spectra are all very similar (see Fig. 4) and consist of five peaks in the 7-12 eV interval, followed by a broader feature at higher energies. All these curves suffer from the same deficiency: The separation between the second and third peak is too small, so the second peak appears as a shoulder of the third due to the limited resolution of the calculations (in Table V the second line for Si2H6 describes the position of this shoulder). The curve closest to experiment was calculated with the EXX functional (note its first peak at exactly the experimental value). SIC yielded a quite unreasonable spectrum, consisting of broad peaks at too low energies, and consolidated its position of the worst of the five functionals tested (as what concerns the calculation of optical absorption spectra, of course).

For these clusters we also provide the static polarizabilities in Table IV. For silane we find, in agreement with the results for sodium clusters, that all the methods lead to polarizabilities lower than the experimental values. In spite of the proper asymptotics of the EXX functional, LDA yields the best results, which indicates that a subtle error cancellation is taking place, and that it is important to have some sort of correlation in the functional. This result seems to hold true for these small systems and deserves a more detailed analysis.

### C Xc potential vs xc kernel

In order to study the sensitivity of the optical spectra to the xc functional used to calculated the ground-state, and to the time-dependent xc potential used to propagate in time, we plot, in Fig. 5, the spectrum of silane, for some choices of $`v_{\mathrm{xc}}^{(1)}`$ and $`v_{\mathrm{xc}}^{(2)}`$ (see Eq. (15)). In the first plot we take ground states calculated with several xc functionals (EXX, LB94 and LDA), and propagate these states with the LDA, while in the second we fix the ground state (calculated with the EXX functional), and use three different potentials to propagate it (LDA, LB94 and EXX). The results clearly indicate that the optical spectrum of SiH4 is much more sensitive to the choice of the starting state, than to the quality of $`v_{\mathrm{xc}}^{(2)}`$. This numerical evidence is in good agreement with the findings of Petersilka et al.. From the second plot of Fig. 5 we see that the ALDA provides the strongest attractive $`f_{\mathrm{xc}}`$, shifting the transitions to lower energies by as much as 0.4 eV, the effect being stronger for the high energy transitions. Nevertheless, the role of these functionals needs to be quantified for larger systems and, in particular, for extended systems where the form of the initial potential is expected to be less important and where correlations are supposed to be the dominant contribution to the excitation spectra. This different behavior of finite and infinite systems poses a theoretical problem when designing functionals as the two classes of systems are mainly sensitive either to the spatial non locality or to the frequency dependence of the xc functional.

## V Conclusion

We presented a detailed study of the impact of different xc potentials and kernels in describing the static and dynamical properties of small clusters. The ionization threshold as well as specific single-particle transitions are very sensitive to the specific form of the xc potential used to obtain the ground state of the system. The static polarizabilities turn out to be rather close to experiment but always smaller than the observed values. Surprisingly, the polarizabilities calculated with the LDA or the GGA are better than the ones calculated with LB94, SIC or EXX. In the case of the EXX, this is related to the absence of dynamical correlations, while the LB94 and the SIC functionals seem to break the subtle error cancellation between exchange and correlation always present in the LDA and the GGA functionals.

For the optical spectra a different conclusion is reached. EXX and LB94 perform remarkably well, even if the high lying excitations (below the ionization threshold) appear shifted to higher energies. This error is traced back to the weak effective attractive xc kernel. It seems clear that the right asymptotics are important but not enough to get a more quantitative agreement with experimental data. SIC performs badly, especially for the silicon clusters. The more traditional functionals, the LDA and the GGA, showed their well known problems: for example, the onset for optical absorption is underestimated by around 1.5 eV in SiH4.

Finally, our spectra calculated with different xc-functional for the ground state calculation and for the time evolution ($`v_{\mathrm{xc}}^{(1)}v_{\mathrm{xc}}^{(2)}`$), suggest that the choice xc potential used to calculated the ground state is more important than the choice of the xc kernel for the case of finite systems. This conclusion should not be valid when one is dealing with infinite periodic systems.

## VI Acknowledgments

We benefit from discussions with E. K. U. Gross, G. Bertsch, K. Yabana, L. Reining, and K. Capelle. We acknowledge support by DGES (PB98-0345), JCyL (VA28/99), RTN program of the European Union NANOPHASE (contract HPRN-CT-2000-00167). M. A. L. M. would like to thank the kind hospitality and encouraging discussions with F. Nogueira. A. R. acknowledges the support from the sabbatical program Salvador de Madariaga of the Spanish MEC (PR2000-0335) and the Ecole Polytechnique. A. C. acknowledges the support from the MEC under the graduate fellowship program. Part of the calculations present in this Article were performed at the Centro de Física Computacional of the Universidade de Coimbra.

## Table Captions

Table I: Calculated ionization potentials, (IP), obtained from the highest occupied KS eigenvalue, and HOMO-LUMO gaps for Na2 and Na4 for different xc functionals compared to available experimental data. All values are in eV.

Table II: Static polarizabilities for the small sodium clusters estimated using a finite electrical field (FF), or calculated from the $`\omega =0`$ Fourier transform of the time-dependent dipole moment (FT). The results are compared to experiments  and other calculations . All values are in Å3.

Table III: Calculated ionization potentials (IP) obtained from the highest occupied KS eigenvalue, and HOMO-LUMO gaps for SiH4 and Si2H6. The available experimental data as well as results from many body calculations are given for comparison. DMC stands for the diffusion quantum Monte-Carlo calculations of Ref. , GW for the GW with exciton effects results (GW-Bethe-Salpeter equation) of Ref. , and finally, QMC for quantum Monte-Carlo calculations . All values are in eV.

Table IV: Static dipolar polarizabilities for the hydrogenated silicon clusters obtained from the $`\omega =0`$ Fourier transform of the time-dependent dipole moment. All values are in Å3.

Table V: Lowest excitation energies for SiH4 and Si2H6 obtained with different xc-functionals compared to experiments and other calculations based on a many body formalism. The SIC results are absent due to the difficulty to identify the peaks in the spectra. DMC stands for the diffusion quantum Monte-Carlo calculations of Ref. , and the calculated value of 9.47 should be compared with the average of the two first, Jahn-Teller splitted, levels. GW are the GW with exciton effects results (GW-Bethe-Salpeter equation) of Ref. , and finally, QMC for Quantum Monte-Carlo calculations . All values are in eV.

## Figure Captions

Fig. 1: Averaged dipole strength for Na2, calculated for several xc functionals. The “exp” curve is from Ref. , who adapted the experimental results of Ref. . The experimental curve is plotted in arbitrary units.

Fig. 2: Averaged dipole strength for Na4, calculated for several xc functionals. “exp” stands for the experimental photodepletion cross sections of Ref. , while “GW” is the many-body calculation of Ref. , which includes self-energy and excitonic effects. These two curves are shown in arbitrary units for the sake of comparison.

Fig. 3: Averaged dipole strength for SiH4, calculated for several xc functionals. The “exp” curve is from Ref. .

Fig. 4: Averaged dipole strength for Si2H6, calculated for several xc functionals. The “exp” curve is from Ref. .

Fig. 5: Averaged dipole strength for SiH4, calculated using different xc-functionals fro the calculation of the ground state and for the time evolution ($`v_{\mathrm{xc}}^{(1)}v_{\mathrm{xc}}^{(2)}`$ ) (see Eq. 15). For each curve, the first functional mentioned corresponds to $`v_{\mathrm{xc}}^{(1)}`$ and the second to $`v_{\mathrm{xc}}^{(2)}`$


format: text

text:
# INTEGRAL capabilities for faint Gamma-Ray Bursts

## 1 Introduction

The study of gamma-ray bursts (GRBs) is just one of the many objectives of INTEGRAL. But developments in the GRB-field over the past few years have made it increasingly clear that INTEGRAL may make a very significant contribution to this fast developing field. We know today that most GRBs originate in the very distant universe. In fact, we believe that their intrinsic brightness allows us to detect these events at epochs corresponding to the formation of the earliest stellar populations. Thus they may be used as probes into the first stages of star formation and their spectra may reveal the early heavy-element enrichment of the interstellar medium. Lamb & Reichart (2000) claim that GRBs could be originated by primitive stellar populations up to redshift $`z50`$.

## 2 The needs of the GRB community

The discovery of the GRB afterglow, in X-rays, radio and the optical (Costa et al. 1997, van Paradijs et al. 1997) have finally provided us the long sought tool for associating a burst with a concrete object in the sky. Although the duration of the afterglow is measured in hours or days instead of the seconds which are characteristic of the bursts themselves, there are still two mandatory requirements for successful afterglow searches: accurate initial positions and rapid dissemination of the alerts. The successes achieved by the SAX team have been based on the positions accurate to maybe 10 square arcminutes and delays of a few hours (Boella et al. 1997). Lately, several afterglows have also been detected based on arcminute positions provided by the IPN with delays of the order of 24 hours (Hurley et al. 2000a, 2000b). No doubt, IPN positions will continue to be useful and warrant follow up for a number of years to come, but if we want to fully exploit the potential of high resolution spectroscopy of the intense phases of the afterglows - and if want to check in higher detail our models for the afterglow process itself - then we must provide arcminute positions with delays of only a few minutes or even less.

## 3 INTEGRAL capabilities

INTEGRAL will be the first gamma-ray spacecraft which combines imaging instruments of high precision and a continuous real time telemetry link. In Table 1 we compare in a simple minded way a number of different space missions with capabilities for GRB research. The missions are divided in three groups depending on their energy range. In the first group the INTEGRAL/IBIS sensitivity is normalized respect to Swift sensitivity, and in the second group the INTEGRAL/JEM-X and SAX sensitivity is given respect to the HETE-2 one. For consistency we have included R$`ø`$mer/WATCH and CGRO/BATSE in a third group.

For calculating a value of the instrumental sensitivity (hereafter named as $`P_{ins}`$) for each experiment, we assume that the sensitivity to a burst is proportional to the square-root of the detector area, and inversely proportional to the square-root of the background countrate. Therefore, this simple formula becomes; $`S\sqrt{A/\mathrm{\Omega }}`$. where $`S`$ is the sensitivity, $`A`$ is the detector area and $`\mathrm{\Omega }`$ is the sky coverage. The most reliable comparison can be done from Table 1 is between Swift and IBIS (the ISGRI part) because these are very comparable detector technologies (CdZnTl in Swift, CdTl in IBIS). Then assuming a sensitivity value of 0.04 ph s-1cm-2 for Swift (Gehrels 1999), we obtain a 0.013 ph s-1cm-2 sensitivity for IBIS. HETE-2 and JEM-X share the same energy range so one could try to obtain the sensitivity of JEM-X. So, considering a sensitivity of 0.2 ph s-1cm-2 for HETE-2 (Ricker 1998) and the relative sensitivity of JEM-X respect to HETE-2 (displayed in Table 1), we obtain a sensitivity of 0.012 ph s-1cm$`2`$ for JEM-X. However, their different detector technology make of this number a preliminary calculation of the JEM-X sensitivity. Besides the reduced number of bursts that JEM-X will detect ($`2`$ per year) does not make worthwhile to perform a specific calculation in order to study the JEM-X capabilities for faint bursts. It is evident from Table 1 that IBIS on INTEGRAL will be the most sensitive gamma-ray burst detector ever flown and not likely to be matched (excluding JEM-X which is detecting a few number of them), sensitivity wise, by any other mission within the coming decade. With this mission we have a unique opportunity!

## 4  Detectability of a faint population of GRBs

We have selected the most sensitive future missions (INTEGRAL/IBIS, Swift and HETE-2) to calculate their capabilities of detecting a high redshift population of bursts. For the estimate of the number of GRBs that these missions will detect, we assume:

* GRB spectra are power-laws; f$`{}_{\nu }{}^{}\nu ^\alpha `$
* The GRB rate is proportional to the star formation rate (SFR) in the universe. The SFR considered are the one given by Rowan-Robinson (1999) for z $`<`$ 5 and the one calculated by Gnedin & Ostriker (1997) for z $``$ 5. (See Fig.1).
* The GRB Luminosity function is given by;

$$S(L)=\{\begin{array}{cc}L^\beta \hfill & \text{ }L_{min}<L<L_{max}\hfill \\ 0\hfill & \text{ Otherwise }\hfill \end{array}$$

being $`L`$ the peak photon luminosity and $`\beta `$ the luminosity function index. $`L_{min}`$, $`L_{max}`$ determine the width of the luminosity function.
* Although the effect of several universe models have been tried, the cosmological parameters presented in this paper are $`\mathrm{\Omega }_m`$=0.3, $`\mathrm{\Omega }_\mathrm{\Lambda }`$=0.7, Ho = 65 km s-1 Mpc-1.

According to the former assumptions the differential GRB detection rate at a given photon peak flux $`P`$ at the detector (ph cm-2$`s^1`$) is given by the following convolution integral;

$$N_{GRB}(P)=C\mathrm{\Omega }e_0^{inf}R_{GRB}S(L)𝑑L$$
(1)
$`e`$ is the efficiency of the orbit, $`\mathrm{\Omega }`$ is the instrumental coverage of the sky and $`R_{GRB}`$ is the GRB detection rate if they were standard candles, i,e: $`R_{GRB}=\frac{SFR(z)}{(1+z)}\frac{dV(z)}{dz}\frac{dz}{dP}`$, being $`V`$ the comoving volume and $`SFR(z)`$ the star formation rate. The value of the proportionality constant $`C`$ is unknown. Fig 2. shows $`N_{GRB}(P)`$ as well as the detection thresholds of several instruments.

The relationship between $`L`$, $`z`$ and $`P`$ is given by the next expression;

$$P=\frac{L}{4\pi D(z)^2(1+z)^\alpha }$$

Where $`D(z)`$ is the comoving distance. In our calculations different values of $`\alpha `$, Lmin, Lmax and $`\beta `$ are considered. The values of $`\alpha `$, $`\beta `$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$ and $`\mathrm{\Omega }_m`$ do not change the final result qualitatively. Instead the values of $`L_{min}`$, $`L_{max}`$ are very relevant to the determination of the number of high redshift GRB detections. We consider the most pessimistic case where $`L_{min}=10^{57}`$ ph s-1 and $`L_{max}=10^{58}`$ ph s-1 (according to the GRB redshifts measured so far the GRB luminosity function seems to be wider). We can calculate the contribution to the integral (1) by the GRBs with redshift larger than $`z_{edge}`$;

$$N_{GRB}\left(P|z_{edge}\right)=C\mathrm{\Omega }e_0^{inf}H(z\left(L\right),z_{edge})R_{GRB}S\left(L\right)𝑑L$$

Where $`H(z(L),z_{edge})`$ is a step function that vanishes unless $`z(L)>z_{edge}`$. Obviously, $`N_{GRB}(P)=N_{GRB}(P|0)`$, and $`\frac{N_{GRB}(P|z_{edge})}{N_{GRB}(P)}1`$. Finally we can calculate the number of GRBs detected above a given instrumental photon flux threshold $`P_{ins}`$ that have redshifts larger than $`z_{edge}`$;

$$N_{GRB}(z_{edge}|P_{ins})=_{P_{ins}}^{inf}N_{GRB}(P|z_{edge})𝑑P$$

The ignorance of the proportionality constant $`C`$ prevents us to derive an absolute value for $`N_{GRB}(z_{edge}|P_{ins})`$. However, we can determine the relative quantity $`\frac{N_{GRB}(z_{edge}|P_{ins})}{N_{GRB}(0,P_{ins})}`$, which provide us the proportion of detections that have a redshift larger than $`z_{edge}`$ (see Fig. 3).

## 5 IBIS vs. Swift; Comparison of the number of the GRB detections

As it is shown in Fig. 3, $``$10% of the GRBs detected by IBIS will have a redshift larger than 8.4. For Swift the $`z>8.4`$ population will be just $``$ 4% of the total number of detections. HETE-2 is the less sensitivity detector, being constrained to detect GRBs with redshifts $`z<6`$. Therefore we will not consider HETE-2 for the further study aimed to calculating the relative number of detections as a function of the redshift. We will be centered in comparing IBIS and Swift capabilities. Besides, as we noted, the similar energy range and detector technologies of IBIS and Swift guarantee a reliable calculation of this fraction.

For determining the relative number of detections between two experiments, $`A`$ and $`B`$, the next expression has to be calculated:

$`f_{A/B}(z_{edge})=\frac{\mathrm{\Omega }_Ae_A_{P_A}^{inf}_0^{inf}H(z\left(L\right),z_{edge})R_{GRB}S\left(L\right)𝑑L𝑑P}{\mathrm{\Omega }_Be_B_{P_B}^{inf}_0^{inf}H(z\left(L\right),z_{edge})R_{GRB}S\left(L\right)𝑑L𝑑P}`$

This function will give the relative number of GRB detections with $`z>z_{edge}`$. We have applied the former expression to derive the fraction $`f_{IBIS/Swift}`$ as a function of the GRB redshift.

If we consider $`z_{edge}=0`$ the fraction $`f_{IBIS/Swift}`$ gives us the fraction of GRBs detected with $`z>0`$, i.e, all the detections independently of their redshifts are considered. We obtain a value of $`f_{IBIS/Swift}(0)=1/7.8`$ (see Fig. 4), which is consistent with the value of $`f_{IBIS/Swift}`$ derived from last column of Table 1. The large field of view (FOV) of Swift in comparison to IBIS makes that for $`z_{edge}<11.6`$ $`f_{IBIS/Swift}<1`$. Instead, for further redshifts than 11.6, IBIS sensitivity becomes the governing factor and $`f_{IBIS/Swift}>1`$.

## 6 Conclusion

Fig. 3 shows $`\frac{N(z_{edge},P_{ins})}{N(0,P_{ins})}`$, for HETE-2, IBIS, and Swift. For IBIS the tail of $`\frac{N(z_{edge},P_{ins})}{N(0,P_{ins})}`$ extends to redshifts $`z_{edge}>11`$. Swift and HETE-2 would detect a closer population of burst, specially HETE-2 would be constrained to redshifts $`z_{edge}<6`$.

If we consider all the GRBs ($`z_{edge}>0`$) the detection fraction $`f_{IBIS/Swift}`$=1/7.8. Thus, at low redshifts the large FOV of Swift in comparison to INTEGRAL instrumentation governs the number of detections. However, at high redshifts the better sensitivity of IBIS makes the fraction of detections $`f_{IBIS/Swift}>1`$ (for redshifts $`z_{edge}>11.6`$ $`f_{IBIS/Swift}>1`$, see Fig. 4). Although JEM-X FOV and sensitivity are less suitable than the one of IBIS to detect GRBs, the spectral peak of the high redshift GRBs (usually at 500–1200 keV) will be in the detection range of JEM-X. Therefore JEM-X will be also a very valuable tool to study the high redshift GRBs.

In conclusion, the capabilities of studying GRBs of JEM-X and IBIS on board INTEGRAL are complementary to the ones of missions like Swift and HETE-2 specially devoted to prompt localizations of GRBs. Whereas Swift and HETE-2 would detect more GRBs than INTEGRAL, JEM-X and IBIS instruments would detect very high redshift GRBs unreachable to the above mentioned missions. Therefore, INTEGRAL and specially IBIS will be a very valuable tool to trace the SFR rate in the early universe.

## Acknowledgments

J. Gorosabel acknowledges the receipt of a Marie Curie Research Grant from the European Commission.


text:
# Accretion in stellar clusters and the IMF

## 1 Introduction

One of the most fundamental unsolved problems in astronomy is the origin of the distribution of stellar masses, the so-called Initial Mass Function (IMF). Significant observational work has gone into establishing the exact form of the mass spectrum (c.f. Scalo 1986; Kroupa, Tout & Gilmore 1990; Kennicut 1998) and investigating it’s universality (Kennicut 1998; Leitherer 1998; Scalo 1998). Unfortunately this work has not led to an increase in our understanding of how the mass spectrum originates. This shortcoming is intrinsically linked to our lack of understanding of the star formation process. Thus, theories of the IMF have generally involved sufficient free-parameters to ensure that they can reproduce the general form of the mass spectrum but also, unfortunately, ensuring that they do not help us discern the physics involved.

Observations of star forming regions have clearly shown that most stars are formed in clusters (e.g. Lada et al. 1991; Lada, Strom and Myers 1993; Meyer & Lada 1999; Clarke, Bonnell & Hillenbrand 2000). Furthermore, the individual stellar masses found in these clusters span the observed range of stellar masses and are consistent with being drawn from the field star IMF (Hillenbrand 1997; Hillenbrand & Carpenter 2000). Thus, the basic physics which results in the observed mass spectrum is present in this fundamental building block of the star formation process: the young stellar cluster. Competitive accretion in gas-rich stellar clusters results in a spectrum of stellar masses (Bonnell et al. 1997, 2000) and in this paper we investigate whether this process can lead to the observed IMF.

There have been several theories advanced to explain the observed distribution of stellar masses (cf. Clarke 1998; Meyer et al. 2000). These theories generally involve either fragmentation or accretion plus feedback as the basic physics in setting the stellar masses. Fragmentation of a giant molecular cloud can set the stellar mass distribution if the GMC clump-mass distribution is, for some reason, of the same form as the IMF (Elmegreen 1997), or if feedback changes the basic fragment (Jeans) mass (Silk 1977). Alternatively, accretion onto a small isolated protostellar core can set the final masses but this requires a halting mechanism (eg feedback) that really determines the mass distribution (Adams & Fattuzzo 1996). The main criticisms of these theories are that a) they rely on additional processes (e.g. feedback or an initial mass distribution) and b) that they neglect the fact that the stars are not formed in isolation but in clusters and therefore can interact during their formation. Furthermore, they are generally unable to explain why clusters are initially mass-segregated (eg, Hillenbrand & Hartmann 1998; Bonnell & Davies 1998) which should be a requirement of any theory that attempts to explain the IMF. In this paper, we discuss how competitive accretion in clusters can set the distribution of stellar masses and result in a mass-segregated cluster.

## 2 Accretion in Clusters

Young stellar clusters are generally found to contain significant amounts of gas. In fact, the total gas mass is typically many times the mass in stars (Lada 1991). This gas is a remnant of the star formation process and can be understood in terms of the inefficiency of fragmentation where the majority of a system’s mass is not immediately included in the fragments but is accreted by them over longer time periods (Larson 1978; Boss 1986; Bonnell & Bastien 1992; Burkert & Bodenheimer 1993; Klessen, Burkert & Bate 1998). This occurs due to the non-homologous nature of gravitational collapse (Larson 1969) so that the amount of mass that has reached the region where and when fragmentation occurs, is small. Thus, the initial fragment mass should generally be small, $`M<<M_{}`$. Furthermore, in a cluster environment, the fragment mass, if it is proportional to the Jeans mass, should be smallest in the centre of the cluster where the density is highest, in direct contradiction to the observed mass segregation (Bonnell & Davies 1998; Bonnell et al. 1998). It is thus unlikely that the mass distribution is solely determined by the fragmentation process.

Accretion of the residual gas in clusters can play a crucial role in setting the final stellar masses (Zinnecker 1982; Larson 1992). As all the stars are embedded in the same envelope of gas, they compete for the gas as they move through the cluster. This competition results in unequal mass accretion rates and hence a mass spectrum

Simulations of accretion in small stellar clusters (Bonnell et al. 1997) have shown that even if the stars have equal initial masses, they accrete unevenly with those stars near the centre accreting more than those near the outer parts of the cluster. This dependency of the accretion rates on position in the cluster is even clearer in larger clusters ($`N=100`$; Bonnell et al. 2000) where accretion rates can vary by over an order of magnitude from the outside to the cluster centre when the cluster is centrally-condensed. The physics behind this differential accretion is the effect of the cluster potential in funneling material down to the centre plus the individual star’s ability to capture this material. Initially uniform clusters also accrete differentially as they evolve towards a centrally-condensed configuration, due to the local competition between stars. The final cluster displays a significant degree of mass segregation due to the differential accretion (Bonnell et al. 2000).

A star’s accretion rate depends on it’s cross section for accreting material, $`\pi R_{\mathrm{acc}}^2`$, the gas density through which it moves, $`\rho `$ and the velocity with which it moves relative to the gas, $`V_{\mathrm{rel}}`$:

$$\dot{M}_{}=\pi \rho V_{\mathrm{rel}}R_{\mathrm{acc}}^2.$$
(1)
The question is what is the accretion radius, $`R_{\mathrm{acc}}`$. An isolated star of mass $`M_{}`$, moving in a uniform medium with velocity $`V_{\mathrm{rel}}`$ relative to the gas with sound speed $`c_s`$, has an accretion radius given by Bondi-Hoyle accretion where

$$R_{\mathrm{BH}}=2GM_{}/(V_{\mathrm{rel}}^2+c_s^2)$$
(2)
(e.g. Bondi & Hoyle 1944; Ruffert 1996). On the other hand, in a cluster environment where cold gas dominates the potential, the accretion radius can be described by the star’s tidal-lobe radius.

$$R_{\mathrm{tidal}}0.5\left(\frac{M_{}}{M_{\mathrm{enc}}}\right)^{\frac{1}{3}}R,$$
(3)
where $`M_{}`$ is the star’s mass, $`M_{\mathrm{enc}}`$ is the cluster mass interior to the star’s position in the cluster $`R`$ (e.g. Paczynski 1971; Frank, King & Raine 1985). Material inside this tidal-lobe is more attracted to the star whereas outside this lobe the acceleration is dominated by the cluster potential.

Simulations of accretion in gas-dominated stellar clusters have shown that the accretion is well modelled by a tidal-lobe radius as long as the gas dominates the gravitational potential (Bonnell et al. 2000). When this is the case and the gas is unsupported, as expected in order for the cluster to have fragmented, both the gas and the stars collapse towards the cluster centre. The stars are impeded from virialising due to the changing potential, and the gas accretion. Therefore the gas and stellar velocities are locally correlated such that the relative gas velocity is small. This results in a large $`R_{\mathrm{BH}}`$ such that $`R_{\mathrm{BH}}>R_{\mathrm{tidal}}`$ and therefore the smaller $`R_{\mathrm{tidal}}`$ is the appropriate accretion radius.

In contrast, once the stars dominate the gravitational potential, which occurs first in the core of the cluster due to the higher accretion rates there and the sinking of the more massive stars, the stars virialise and their velocities relative to the gas become large. When this occurs the Bondi Hoyle radius becomes much smaller due to the large gas velocities such that $`R_{\mathrm{BH}}<R_{\mathrm{tidal}}`$. At this point, the Bondi-Hoyle radius is a better description of the accretion radius (Bonnell et al. 2000).

## 3 The Initial Mass Function

In order to use the above formulations of the accretion rate to determine the resultant mass spectrum, we need a basic model of a stellar cluster. The cluster is initially gas dominated in order to provide the mass reservoir for accretion. Observations of the youngest clusters estimate that up to 90 per cent of the mass is in gas (Lada 1991). Such systems are unlikely to be in equilibrium, and regardless of their initial configurations, will collapse towards a

$$\rho R^2$$
(4)
density distribution (Larson 1969; Hunter 1977). The stars will also collapse, even if initially virialised, due to the change in the potential. Thus we assume that the stars follow the same distribution:

$$nR^2$$
(5)
This profile is a good approximation to the stellar distribution of most clusters (Binney & Tremaine 1987) including the youngest clusters that have been studied such as the Orion Nebula cluster (Hillenbrand & Hartmann 1998). Furthermore, it is the expected distribution produced by a violent relaxation of a non-equilibrium purely stellar cluster (Lynden-Bell 1967).

The gas distribution will be modified by the accretion and hence removal of material. The accretion rates are largest in the centre and hence it is there that they will first affect the density profile. This gas removal destabilises the equilibrium of the isothermal sphere allowing gas to infall from further out in the cluster. In the case of infall onto a point-mass, the density profile then takes on the form (Hunter 1977; Shu 1977; Foster & Chevalier 1993)

$$\rho R^{\frac{3}{2}}.$$
(6)
We thus assume an initial density profile of an isothermal sphere (equation 4) which evolves towards a profile composed of the shallower, equation(6) profile in the central regions while the outside maintains the isothermal sphere profile. The transition between these two regimes occurs when gas accretion has increased the stellar masses sufficiently that the stars dominate the gravitational potentialin the cluster centre. We discuss the asymptotic limits to the accretion and mass spectrum in these two cases below.

### 3.1 Gas dominated potentials

In addition to the gas density we require estimates of the relative gas velocity and of the accretion radius. Where gas dominates the potential ($`\rho R^2`$), the accretion radius is given by the tidal-lobe radius (Bonnell et al. 2000)

$$R_{\mathrm{acc}}R_{\mathrm{tidal}}.$$
(7)
In order to evaluate the tidal-lobe radius we need the mass enclosed at a particular radius in the cluster, $`R`$,

$$M_{\mathrm{enc}}=_0^R4\pi \rho r^2𝑑r,$$
(8)
which in a $`\rho R^2`$ distribution is

$$M_{\mathrm{enc}}R.$$
(9)
Simulations of accretion in stellar clusters have shown that in gas-dominated clusters, the relative gas velocity is generally subsonic as both the stellar and gas velocities are dominated by the acceleration of the cluster potential (Bonnell et al. 2000). Alternatively, if the stars are virialised then the relative gas velocity is the velocity dispersion,

$$V_{\mathrm{rel}}\sqrt{\frac{GM_{\mathrm{enc}}}{R}},$$
(10)
Thus, if the gas distribution follows equation(4), then we have

$$V_{\mathrm{rel}}=constant.$$
(11)
In either case the relative gas velocity is independent of position in the cluster.

We can now use the above to consider how accretion in a stellar cluster results in a mass spectrum of the form

$$dNM_{}^\gamma dM_{}$$
(12)
where $`\gamma =2.35`$ is the Salpeter IMF. First of all, we note that from equation(5) the number of stars at a given radius is constant:

$$dN=n(R)4\pi R^2dR\mathrm{const}\times dR.$$
(13)
Starting from the general equation for the accretion rates,

$$\dot{M}_{}=\pi \rho V_{\mathrm{rel}}R_{\mathrm{acc}}^2,$$
(14)
and adapting it for the case where gas dominates the potential such that $`\rho R^2`$ and $`R_{\mathrm{acc}}R_{\mathrm{tidal}}`$, then we have

$$\begin{array}{cc}\dot{M}_{}\hfill & g(t)R^2\times \left(\left[\frac{M_{}}{M_{\mathrm{enc}}}\right]^{\frac{1}{3}}R\right)^2\hfill \\ & g(t)\left(\frac{M_{}}{M_{\mathrm{enc}}}\right)^{\frac{2}{3}}\hfill \\ & g(t)\left(\frac{M_{}}{R}\right)^{\frac{2}{3}}.\hfill \end{array}$$
(15)
The term $`g(t)`$ includes all time-dependence of the homologously evolving cluster and $`R`$ is the initial stellar position in the cluster.

Integrating equation(15) with respect to time, $`t`$, and in the limit where accretion dominates the final mass, we get

$$M_{}R^2h(t)^3,$$
(16)
where $`h(t)=g(t)𝑑t`$. Thus, at a time $`t_{\mathrm{final}}`$, when the accretion is halted, we have

$$RM_{}^{\frac{1}{2}}$$
(17)
and thus

$$dRM_{}^{\frac{3}{2}}dM_{}.$$
(18)
This results in a mass spectrum (equation 13)

$$dNM_{}^{\frac{3}{2}}dM_{}.$$
(19)
In general, the $`\rho R^2`$ profile will not extend right to the centre of the cluster. Instead, the central regions must have a flatter profile to stop the density from becoming infinite. A collapsing isothermal gas cloud has a $`\rho R^2`$ density profile in the outer regions and an ever decreasing central region with near uniform density (Larson 1969). Stars in this region will see a shallower density profile which will affect their accretion rates (equation 14). Adapting equation(15) for a density profile $`\rho R^\alpha `$ while maintaining a uniform $`V_{\mathrm{rel}}`$ gives

$$M_{}R^\alpha $$
(20)
and analogous to equation(18),

$$dRM_{}^{\frac{1}{\alpha }1}dM_{}.$$
(21)
Taking the stellar density to be proportional to the gas density, $`n(R)\rho (R)`$, equation(13) becomes

$$\begin{array}{cc}dN\hfill & =n(R)4\pi R^2dRR^{2\alpha }dR,\hfill \\ & M_{}^{\frac{2}{\alpha }}M_{}M_{}^{\frac{1}{\alpha }1}dM_{},\hfill \\ & M_{}^{\frac{3}{\alpha }}dM_{},\hfill \end{array}$$
(22)
and $`\gamma `$ is now $`\gamma =3/\alpha `$.

From equation(22) we see that when the density profile is shallower than $`R^2`$, the resulting mass spectrum is steeper than $`M_{}^{3/2}`$. Figure 1 plots the resulting IMF using equation(22) and a density profile that combines an outer region with a $`\rho R^2`$ profile with an inner Gaussian density profile (fitted onto the outer region where it has the same slope). The mass spectrum is plotted as a function of the logarithm of the mass such that the slope is $`\mathrm{\Gamma }=\gamma +1`$. The low-mass stars that form in the outer region have the $`\gamma =3/2`$ mass spectrum whereas the higher-mass stars formed in the inner regions have an ever steeper mass spectrum due to the central Gaussian density profile that converges to uniform density. It is also worth noting that if the collapse gives a density profile that is steeper than $`R^2`$ as was found in the outer regions of the collapse of a Bonnor-Ebert spheroid (Foster & Chevalier 1993), then this translates (equation) into a shallower IMF for the low-mass stars that are in this region.

### 3.2 Stellar dominated potentials

In contrast to a gas dominated system, when the cluster potential is dominated by the stars and they are virialised, the appropriate accretion radius is the Bondi-Hoyle radius (Bonnell et al. 2000). The Bondi-Hoyle radius models better the accretion radius as the gas and stellar velocities are no longer correlated. The gas is primarily infalling from further out in the cluster onto the stellar dominated core. When infall is onto a point mass, the resulting velocity profile is

$$V_{\mathrm{rel}}R^{1/2}$$
(23)
which assuming a constant mass infall rate with radius implies a density profile of the form of equation(6),

$$\rho R^{3/2}.$$
(24)
In the case of accretion onto a stellar cluster, the distributed mass and the presence of multiple accretion sources will modify these profiles. We can generalise them as

$$V_{\mathrm{rel}}R^\eta ,$$
(25)
and

$$\rho R^\zeta ,$$
(26)
under the constraint that the mass infall, $`4\pi R^2\rho V_{\mathrm{rel}}`$ decreases with decreasing radius, $`\eta +\zeta 2`$. we have

$$\begin{array}{cc}\dot{M}_{}\hfill & V_{\mathrm{rel}}\rho R_{\mathrm{BH}}^2\hfill \\ & V_{\mathrm{rel}}\rho \left(\frac{M_{}}{V_{\mathrm{rel}}^2}\right)^2\hfill \\ & \rho \frac{M_{}^2}{V_{\mathrm{rel}}^3}\hfill \\ & g(t)R^\zeta M_{}^2R^{3\eta }\hfill \\ & g(t)M_{}^2R^{3\eta \zeta }.\hfill \end{array}$$
(27)
As above, $`g(t)`$ includes the time-dependence of the homologously evolving cluster and $`R`$ is the initial stellar position in the cluster. In the case where $`3\eta \zeta =0`$ as occurs for accretion onto a point mass, the accretion rates are independent of position in the cluster.

Assuming that the cluster maintains the same form ($`3\eta \zeta `$ constant), solutions to equation(27) are

$$M_{}=\frac{M_{}(0)}{1\beta M_{}(0)h(t)}$$
(28)
where $`M_{}(0)`$ is the initial stellar mass, $`h(t)=g(t)𝑑t`$, and $`\beta `$ is basically the accretion clock which depends on the gas content and on $`R^{3\eta \zeta }`$ but not on time. The primary difference between this case and the preceding tidal-lobe accretion is the higher dependency on the stellar mass such that the initial mass is always important. This is the case for any accretion rate which depends on the stellar mass to a power greater than one.

The resulting mass spectrum from such an accretion rate was studied by Zinnecker (1982) and we follow his derivation here. As there is a direct correlation between initial and final mass, then there is a direct mapping from the initial mass spectrum to the final one,

$$F(M_{})dM_{}=F(M_{}(0))dM_{}(0).$$
(29)
Inverting equation(28) and inserting it into equation(29), as long as $`\beta `$ is independent of the initial masses, we have

$$\begin{array}{cc}F(M_{})\hfill & =F(M_{}(0))\left(\frac{dM_{}}{dM_{}(0)}\right)^1\hfill \\ & =F\left(M_{}(0)\right)\left(\frac{d}{dM_{}(0)}\frac{M_{}(0)}{1\beta M_{}(0)h(t)}\right)^1,\hfill \\ & =F\left(M_{}(0)\right)\left(1\beta M_{}(0)h(t)\right)^2,\hfill \\ & =F\left(\frac{M_{}}{1+\beta M_{}h(t)}\right)\left(\frac{M_{}(0)}{M_{}}\right)^2.\hfill \end{array}$$
(30)
In the limit $`\beta M_{}h(t)>>1`$ ($`M_{}>>M_{}(0)`$) this gives

$$F(M_{})M_{}^2,$$
(31)
and thus as long as there is an initial (small) range of stellar masses, the asymptotic limit of the mass spectrum is

$$dNM_{}^2dM_{}.$$
(32)
This is somewhat complicated if there is a correlation between the initial masses and radii in the cluster. In this case $`\beta `$ is dependent on $`M_{}(0)`$ and we cannot derive an analytical expression for the resulting mass spectrum. In the case where the incoming mass distribution is from the tidal-lobe accretion, then from §3.1 we have an expression for the initial mass-radius correlation (equation 17). Using this, we can generate a mass spectrum from accretion for a cluster in a $`nR^2`$ configuration. Numerical simulations of accretion in stellar dominated clusters have values of $`\eta =1/3`$ and $`\zeta =3/2`$ (Bonnell et al. 2000). Equation (27) then becomes

$$\dot{M}_{}M_{}^2R^{\frac{1}{2}}.$$
(33)
Figure 2 shows the resultant IMF for clusters accreting from an initial delta function mass distribution and those with an initially $`M_{}^{3/2}`$ mass spectrum that accretes according to equation (33). In the former case, the stars, in a $`nr^2`$ distribution, accrete with a spatially-independent accretion rate ($`\dot{M}_{}M_{}^2`$). These clusters develop an $`M_{}^2`$ mass spectrum. In contrast, with the initial mass segregation, the IMF develops into a steeper ($`\gamma 2.5`$) IMF as the stars accrete in the stellar dominated regime. This is somewhat steeper than the analytical $`M_{}^2`$ mass spectrum due to the initial correlation between stellar mass and position in the cluster.

We expect stellar dominated clusters to accrete towards a mass spectrum with asymptotic limits of $`2.5\mathrm{¡}\mathrm{}\gamma 2`$. This mass spectrum is expected for the higher mass stars as it is they that are located in the centre of the cluster and have attained a higher mass due to their previous accretion in the gas dominated regime. Clusters that have mass distributions that are uncorrelated with cluster position going into the stellar dominated phase result in mass spectra of $`\gamma 2`$ whereas those that are mass segregated initially will have slightly steeper mass spectra $`\gamma 2.5`$.

It is important to note that the above derivation relies upon all the stars stopping their accretion at the same time, as the final masses are a strong function of the time (equation 28). This should not be too strong a limitation as in the stellar dominated potential the accretion timescale must be much longer than the local crossing time as $`M_{\mathrm{gas}}<<M_{\mathrm{stars}}`$. In a rich young stellar cluster, gas removal is due to the presence of massive stars. These stars ionise the gas such that the sound speed ($`10`$ km/s) is comparable to or greater than the velocity dispersion. Thus the gas should be removed on a timescale shorter than or comparable to the crossing time, and all the stars in the core should stop accreting quasi-simultaneously. If this is not the case, then the resulting mass spectrum will exhibit a local maximum corresponding to stars that are just outside the cleared region.

In summary, we expect that the combination of a gas dominated and a stellar dominated regimes and the different accretion physics operating in each, results in a two power-law IMF. The lower-mass stars have a shallower ($`\gamma 1.5`$) slope as their mass accumulation is dominated by tidal-lobe accretion wheras the higher mass stars have a steeper ($`2\gamma \mathrm{¡}\mathrm{}2.5`$) slope as their mass accumulation is dominated by Bondi-Hoyle accretion in the stellar dominated core.

## 4 Simulated Cluster IMFs

In order to explore the relevance of the asymptotic limits for the mass spectrum, we performed a number of simulations of accretion onto a cluster of 1000 stars. The stars are initially of equal mass ($`0.1M_{}`$) and are embedded in gas which comprises 91 per cent of the total cluster mass The gas is cold such that it contains 1000 Jeans masses,

$$M_J=\left(\frac{5R_gT}{2G\mu }\right)^{3/2}\left(\frac{4}{3}\pi \rho \right)^{1/2},$$
(34)
where $`\rho `$ is the gas density, $`T`$ is the gas temperature, $`R_g`$ is the gas constant, $`G`$ is the gravitational constant, and $`\mu `$ is the mean molecular weight. The simulations were performed with a hybrid N-body SPH code (Bate, Bonnell & Price 1995) which uses standard SPH particles to model the gas and sink-particles to model the stars (for more details see also Bonnell et al. 2000). These sink-particles interact only gravitationally with the rest of the cluster and by accreting gas particles that come within their sink radius. This sink or accretion radius was taken to be the smaller of the tidal-lobe radius, the Bondi-Hoyle radius or the stellar separation. In addition, gas particles can only be accreted if they are bound to the star. The simulations were performed with either 9000 or 90000 SPH particles. The low-resolution simulations resulted in an excess of stars in the lowest (initial) mass bins as they were not able to resolve the low accretion rates involved. The higher-resolution simulations resolved the accretion onto all stars and both low and high-resolution simulations resulted in similar high-mass mass-spectra.

Both the stellar and gas distribution are initially uniform but collapse down towards a centrally condensed distribution. Gas is accreted by the stars and removed from the simulation. As the gas is accreted preferentially near the centre once a power-law density profile is established ($`\rho R^2`$), the stars soon dominate the potential there. The gas must then infall from further out in the cluster. The simulations result in mass spectra that are broadly consistent with our analytical expectations. The low-mass stars that accrete their mass during the gas dominated phase display shallow ($`\gamma 3/2`$) mass spectra whereas the high-mass stars that accrete the majority of their mass during the stellar dominated phase display steeper ($`2.5\mathrm{¡}\mathrm{}\gamma 2`$) mass spectra. As the clusters are initially uniform and cold, only the end of the simulations develop a power-law density profile and that just before entering the stellar dominated regime. Because of this, the expected low-mass mass spectrum does not extend over a large mass range.

In Figure 3, we plot the resultant mass function for one of the higher resolution simulations, (performed with 90000 SPH particles) at three different times during the evolution. The mass spectrum increases in breadth with time and develops into a power-law for the high-mass stars. The low-mass mass spectrum is broadly consistent with tidal-lobe accretion in a mostly uniform cloud, although only over a small range in mass. Due to the uniform initial conditions, only a fraction of the stars experience the power-law density profile before entering the stellar-dominated phase. This region would extend further in mass if the cluster was initially centrally condensed when accretion begins. At later times, the high-mass mass spectrum converges to a slope between $`\gamma =2`$ and $`\gamma =2.5`$ consistent with the accretion during the stellar dominated phase, and to the Salpeter slope of $`\gamma =2.35`$.

### 4.1 Accretion dynamics

In order ascertain whether we are correct in our interpretation of a two-power law IMF resulting from accretion in gas-dominated and stellar-dominated regimes, we evaluated how much of the eventual stellar mass was added in each regime. Figure 4 plots the amount of mass accumulated by each star during the stellar-dominated phase against the final stellar mass. High-mass stars accumulate the majority of their mass during this phase whereas most of the low-mass stars only accrete a small fraction of their mass during this phase. The break between the two regimes occurs at approximately the same final mass ($`1M_{}`$) where the mass function displays a break between the two slopes. This supports the assertion that the two different power-laws in the IMF derive from different physical regimes which affect how the stars accrete. The low-mass stars derive their mass from tidal-lobe accretion during the gas-dominated regime whereas the high-mass stars derive the majority of their mass from a Bondi-Hoyle type accretion that occurs in the inner parts of the cluster where the potential is dominated by the stars themselves.

## 5 Discussion

In addition to producing the two power-law mass spectrum, competitive accretion naturally results in a certain degree of mass segregation. This arises due to the accretion in the gas-dominated phase where there is a strong correlation between accretion rate, and thus the final mass, and position in the cluster (see equations(16) and (20). This direct correlation between the final mass and position in the cluster neglects variations in the initial masses and the relative movements of the stars due to their interactions. If the cluster is mass segregated entering the stellar dominated phase, then the $`\dot{M}_{}M_{}^2`$ implies that the mass segregation will persist. Simulations of accretion in clusters show that the mass segregation does result but that there is not a one to one correlation between mass and radius (Bonnell et al. 2000). In fact, low-mass stars are located throughout the cluster, including in the core, but the high-mass stars are predominantly located in the central regions as is found in young stellar clusters such as the ONC (Hillenbrand 1997).

It is also worth noting that although the models presented here are meant to consider accretion onto young stars, they are equally appropriate for the growth of clumps in a molecular cloud. As the clumps evolve towards gravitational instability, they will accrete from the surrounding gas and this accretion will be governed by the physics described here. Thus, for example, the clump mass-function found by Motte, André & Neri (1998) for the $`\rho `$ Oph molecular cloud could be due to the accretion by the pre-stellar clumps as the whole system collapses down to form a cluster. The observed $`\gamma =3/2`$ slope would imply that the whole system is in a $`\rho R^2`$ density configuration and is subvirial (dominated by the diffuse gas not in the clumps). The steeper slope found by Motte et al. (1998) at the high-mass end of the mass spectrum can be interpreted as arising from a region of near-uniform gas density. A test of such a possibility is to estimate the degree of mass segregation of the clumps in this pre-stellar cluster system.

Finally, it is possible that the mass spectrum for massive stars, $`M_{}\mathrm{¿}\mathrm{}10M_{}`$, is significantly different than that for lower mass stars if they do not form in a similar fashion. If massive stars cannot accrete above $`10M_{}`$ due to the effect of radiation pressure on the infalling dust (Yorke & Krügel 1977; Yorke 1993) but form through a merger process in a dense core (Bonnell et al. 1998) then the expected mass spectrum could be significantly different from that presented here.

## 6 Conclusions

Competitive accretion in young stellar clusters results in a two power-law mass spectrum with a slope of $`\gamma 3/2`$ for low-mass stars and a steeper slope of $`2\gamma \mathrm{¿}\mathrm{}2.5`$ for high-mass stars. The different slopes are due to whether the gas or the stars dominate the cluster potential. When gas dominates the cluster potential, the accretion rates are given by a tidal-lobe accretion radius and results in a $`dNM_{}^{3/2}dM_{}`$ mass spectrum in an isothermal sphere distribution. The shallower density profile in the central regions implies a steeper mass spectrum, $`3/2\gamma \mathrm{¿}\mathrm{}4`$. Once the stars have accreted enough to dominate the potential, which occurs from the inside out, the appropriate accretion radius is the Bondi-Hoyle radius. When this happens the strong dependence of the accretion rate on stellar mass $`\dot{M}_{}M_{}^2`$ results in a mass spectrum of $`dNM_{}^2dM_{}`$ to $`dNM_{}^{2.5}dM_{}`$.

Simulations of accretion in clusters containing 1000 stars result in mass spectra broadly consistent with these analytic estimates. The high-mass stars are confirmed to accumulate the majority of their mass in a stellar-dominated potential. Lastly, competitive accretion in clusters naturally results in mass segregation as the accretion rates are higher near the cluster centre during the gas dominated phase.

## 7 Acknowledgements

We thank Hans Zinnecker for useful discussions and continual enthusiasm. IAB acknowledges support from a PPARC advanced fellowship.


text:
# Bosonization Theory of Excitons in One-Dimensional Narrow Gap Semiconductors

## Abstract

Excitons in one-dimensional narrow gap semiconductors of anti-crossing quantum Hall edge states are investigated using a bosonization method. The excitonic states are studied by mapping the problem into a non-integrable sine-Gordon type model. We also find that many-body interactions lead to a strong enhancement of the band gap. We have estimated when an exciton instability may occur.

One-dimensional (1D) narrow gap semiconductors can be realized in anti-crossing quantum Hall edge states and nanotubes. The gaps in these systems are single particle gaps and not many body gaps. Theoretically, they provide the unusual condition that the bare band gap $`t`$ is much smaller than the characteristic Coulomb energy scale $`E_c`$. Moreover, strong quantum fluctuations are present in these systems, reflecting the 1D character. The ground state of these systems may be unstable against the spontaneous formation of excitons if the exciton binding energy exceeds the band gap. In three dimensional semiconductors excitons can be treated successfully by solving the Bethe-Salpeter equation . However, in the strong coupling regime of 1D systems, the perturbative approaches are not expected to be reliable due to large quantum fluctuations. If the Coulomb scale is much larger than the gap one might naively expect that exciton instability would occur. However, this simple picture neglects screening which is expected to be large due to the smallness of the gap. It is unclear whether a bound state of an ”electron” and a ”hole” can exists in the presence of strong quantum fluctuations. Not much is known about the physics of excitons in 1D narrow gap semiconductors.

This problem can be addressed within a bosonization approach, which is applicable as long as the the characteristic energy scale of the problem is smaller than the band width $`W`$. Thus bosonization provides a natural framework for studying 1D excitons. In this work we will consider anti-crossing of quantum Hall edge states in the barrier region of between two 2D electron gases. Assume that the applied magnetic field is sufficiently strong that the system is spin-polarized with the filling factor $`\nu =1`$ and that the Fermi level is in the gap of anti-crossing edge states. We find that the problem can be mapped into a sine-Gordon (sG) model with an extra term representing the long range Coulomb interaction. Due to this extra term the model becomes non-integrable and an exact solution is unavailable. However, the non-linear cosine term of this sG model can be expanded provided that $`\delta _{ex}\frac{v/v_0}{\mathrm{ln}W/t}<1`$, where $`v`$ is the Fermi velocity, $`v_0=e^2/ϵ\mathrm{}`$, $`W=\mathrm{}v/a`$, and $`a`$ is the smallest length scale in the problem. In this regime a perturbative theory may be applied to calculate excitonic energies. Our work also indicates that a large enhancement of band gap occurs. Our perturbative approach suggests an approximate estimate for when an exciton instability may occur: $`\delta _{ex}>1`$ . (We cannot exclude the possibility that the higher order corrections neglected in our perturbative calculation may work against such an instability).

The system is modeled by the following Hamiltonian

$`H`$ $`=`$ $`H_0+H_{\mathrm{coul}}+H_t,`$ (1)
$`H_0`$ $`=`$ $`v{\displaystyle 𝑑x\left[i\psi _R^{}_x\psi _R+i\psi _L^{}_x\psi _L\right]}`$ (2)
$`=`$ $`\pi v{\displaystyle 𝑑x\left[\rho _R^2+\rho _L^2\right]}`$ (3)
$`H_{\mathrm{coul}}`$ $`=`$ $`{\displaystyle }dxdy[{\displaystyle \frac{V(xy)}{2}}(\rho _R(x)\rho _R(y)+\rho _L(x)\rho _L(y))`$ (4)
$`+`$ $`V(xy)\left(\rho _L(x)\rho _R(y)\right)]`$ (5)
$`H_t`$ $`=`$ $`t{\displaystyle }dx[\psi _R^{}(x)\psi _L(x)+\mathrm{H}.\mathrm{C}],`$ (6)
The operators $`\psi _R(\psi _L)`$ is the right-moving (left-moving) edge electron operators. $`\rho _R=:\psi _R^{}\psi _R:`$ is the (normal-ordered) right-moving edge electron density operator ($`\rho _L`$ is similarly defined). In practice, we can take $`V(x)=\frac{e^2}{ϵ}\frac{1}{\sqrt{x^2+a^2}}`$ and the Coulomb matrix element is $`V(k)=\frac{2e^2}{ϵ}K_0(a|k|)`$. The tunneling between the right-moving and left-moving is modeled by $`H_t`$. The tunneling amplitude $`t`$ is assumed to be significant only near the Fermi wave vector $`k_F`$(The Fermi wavevector is set to zero from now on). The cut-off value of the wave vector is $`1/a`$. Note that a single particle gap opens up near the Fermi points due to this tunneling term. The commutation relations of density operators in momentum space are given by  $`[\rho _R(q),\rho _R(q^{})]=\frac{q}{2\pi }\delta _{q+q^{}},[\rho _L(q),\rho _L(q^{})]=+\frac{q}{2\pi }\delta _{q+q^{}}`$ The interacting electron systems can be bosonized in a standard way. The explicit relations between the electron operators and the bosonic variable $`\varphi _R`$, $`\varphi _L`$ are $`\rho _R=\frac{1}{2\pi }_x\varphi _R,\rho _L=\frac{1}{2\pi }_x\varphi _L,\psi _R=\frac{1}{\sqrt{2\pi a}}e^{i\varphi _R},\psi _L=\frac{1}{\sqrt{2\pi a}}e^{i\varphi _L}`$. It is convenient to introduce the conjugate fields $`\varphi _\pm `$: $`\varphi _\pm \varphi _R\pm \varphi _L`$. The effective bosonized action in imaginary time reads

$`S`$ $`=`$ $`{\displaystyle 𝑑x𝑑\tau \left[\frac{v}{8\pi }\left((_x\varphi _+)^2+(_x\varphi _{})^2\right)+\frac{i}{4\pi }_\tau \varphi _+_x\varphi _{}\right]}`$ (7)
$`+`$ $`{\displaystyle \frac{1}{8\pi ^2}}{\displaystyle 𝑑x𝑑y𝑑\tau V(xy)_x\varphi _+(x)_y\varphi _+(y)}`$ (8)
$``$ $`{\displaystyle \frac{t}{2\pi a}}{\displaystyle }dxd\tau [e^{i\varphi _+(x)}+\mathrm{H}.\mathrm{C}].`$ (9)
The $`\varphi _{}`$ can be integrated out, leaving us with

$`S`$ $`=`$ $`{\displaystyle \frac{1}{8\pi }}{\displaystyle \frac{dkd\omega }{(2\pi )^2}\left[vk^2(1+\frac{V(k)}{\pi v})+\frac{\omega ^2}{v}\right]}`$ (10)
$`\times `$ $`\varphi _+(\omega ,k)\varphi _+(\omega ,k)`$ (11)
$``$ $`{\displaystyle \frac{t}{\pi a}}{\displaystyle 𝑑x𝑑\tau \mathrm{cos}(\varphi _+(x,\tau ))}.`$ (12)
The above action looks very similar to the sG model, except for the momentum-dependent Coulomb interaction $`V(k)`$. If $`V(k)`$ were momentum independent (local interaction in real space), the action would be exactly that of sG model. The Euclidean action of the sine-Gordon model is given by

$$A_{SG}=d^2x\left[\frac{1}{16\pi }(_\mu \varphi )^22\mu \mathrm{cos}(\beta \varphi )\right],$$
(13)
where the speed of ”light” has been set to unity. It is useful to show the equivalence between the above sG action and our action (10) for the short range interaction ($`V(k)=V`$ is constant). This is achieved by the relations

$`\varphi `$ $`=`$ $`\varphi _+[4(1+{\displaystyle \frac{V}{\pi v}})]^{1/4},\beta =[4(1+{\displaystyle \frac{V}{\pi v}})]^{1/4},`$ (14)
$`\mu `$ $`=`$ $`{\displaystyle \frac{t}{\pi av(1+V/\pi v)^{1/2}}}.`$ (15)
From the above relations it is clear that the strong coupling regime of the original electron system (large $`V`$) is mapped to the weak coupling regime (small $`\beta `$)of sG model.

Before we investigate the long-range case it is instructive to review some known the physical properties of sG model, which is exactly solvable. The spectrum of sG model consist of the breathers ($`B_n`$), the soliton, and the anti-solition. The breathers are the bound states of the soliton and anti-soliton . The lightest bound state $`B_1`$ coincides with the fundamental boson $`\varphi `$ of the sG model (13), and it is associated with the perturbative treatment. The number of breather is given by $`n=1,2,3,\mathrm{}<1/\xi `$, where $`\xi =\frac{\beta ^2}{1\beta ^2}`$. Note that even for infinitesimally small short range repulsion $`V`$, at least one breather exists. The breather energy (excitation energy) is given by

$$m=2M\mathrm{sin}\frac{\pi \xi }{2},$$
(16)
where $`M`$ is the energy of solition. The parameter $`\mu `$ and the soltion energy $`M`$ are related through

$$\mu =\frac{\mathrm{\Gamma }(\beta ^2)}{\pi \mathrm{\Gamma }(1\beta ^2)}\left[M\frac{\sqrt{\pi }\mathrm{\Gamma }(\frac{1+\xi }{2})}{2\mathrm{\Gamma }(\frac{\xi }{2})}\right]^{22\beta ^2}.$$
(17)
Note that in the $`\beta 0`$ limit, $`Mv_\rho ^2\mu /\beta ^2\mathrm{}`$ and $`mv_\rho ^2(\mu /\beta ^2)\beta ^4v_\rho ^2\mu \beta ^2t/va`$, where ”the speed of light” $`v_\rho `$ has been reinstated for clarity . In other words, in the large $`V`$ limit the soliton becomes very massive and leads to a large enhancement of the band gap. In contrast, the lightest breather energy approaches a constant value given by the coefficient of the cosine term of Eq.(10).

In the absence of the Coulomb interaction the exact value of the single particle gap should be $`t`$. However, the bosonized action (10) contain the factor $`a`$ , and it is unclear how this factor disappears in the final result for the gap. Let us try to understand how this happens. The crucial fact is that the cosine term of the sG action (13) is normal ordered and it gives the dimension to the cosine operator. The exact results Eq.(16, 17) were obtained with

$$\mu ^2:\mathrm{cos}(\beta \varphi (x))::\mathrm{cos}(\beta \varphi (y)):\frac{1}{2}\frac{\mu ^2}{|xy|^{4\beta ^2}},$$
(18)
as $`|xy|0`$. In the absence of Coulomb interaction, corresponding to $`\beta ^2=1/2`$ , the same correlation functions (two-point function of tunneling term) can be computed exactly when expressed in terms of electron operators

$$[t\psi _R^{}(x)\psi _L(x)][t\psi _L^{}(y)\psi _R(y)]=\frac{t^2}{|xy|^2}.$$
(19)
The formal application of bosonization formula to this electron tunneling term gives $`\frac{t}{\pi a}\mathrm{cos}(\beta \varphi )`$ as in our action (10). As it stands it is not normal ordered, and consequently the short-distance singularity must be regularized in the calculation of correlation function.

$`\left({\displaystyle \frac{t}{2\pi a}}\right)^2e^{i\beta \varphi (x)}e^{i\beta \varphi (y)}=`$ (20)
$`\left({\displaystyle \frac{t}{2\pi a}}\right)^2\mathrm{exp}\left[\beta ^2(\varphi (0)\varphi (0)\varphi (x)\varphi (y))\right].`$ (21)
We have to regularize $`\varphi (0)\varphi (0)`$ as $`\varphi (0)\varphi (a)=2\mathrm{ln}\frac{L^2}{a^2}`$, where $`L`$ is the system size. $`L`$ always appears in the definition of propagator but they are canceled for the physical correlation functions (see below). Since $`\varphi (x)\varphi (y)=2\mathrm{ln}\frac{L^2}{|xy|^2}`$ we obtain

$$\left(\frac{t}{2\pi a}\right)^2e^{i\beta \varphi (x)}e^{i\beta \varphi (y)}=\left(\frac{t}{2\pi a}\right)^2\frac{1}{(|xy|/a)^{4\beta ^2}}.$$
(22)
Only at $`\beta ^2=1/2`$, the length scale $`a`$ drops out. Comparing Eq.(18) with Eq.(22), the $`t`$ can identified with $`\mu `$. Then by applying formula Eq.(17) at $`\beta ^2=1/2`$ , we find that $`t=M`$, which is the expected result. The above argument checks the consistency of the bosonization formulation of our problem. We emphasize again, the above results holds only if the Coulomb interaction is absent.

When the long range Coulomb interaction is present the model becomes non-integrable and an exact solution is unavailable. But, under certain conditions, the cosine term of the action (10) can be expanded and excitons may be studied perturbatively. (See the perturbative calculation below). Expanding the cosine term of (10) up to the 4th order, we get an effective action of exciton states.

$`S`$ $``$ $`{\displaystyle \frac{1}{8\pi }}{\displaystyle \frac{dkd\omega }{(2\pi )^2}\left[vk^2(1+\frac{V(k)}{\pi v})+\frac{\omega ^2}{v}+\frac{4t}{\mathrm{}a}\right]}`$ (23)
$`\times `$ $`\varphi _+(\omega ,k)\varphi _+(\omega ,k)`$ (24)
$``$ $`{\displaystyle \frac{1}{4!}}{\displaystyle \frac{t}{\pi \mathrm{}a}}{\displaystyle 𝑑x𝑑\tau }:\varphi _+^4:.`$ (25)
At this point, we note that the parameter $`\frac{t}{\pi a}`$ in action (10) should be understood as a renormalized quantity. The application of bosonization formula to the fermion bilinear mass term assumes the implicit normal ordering . Such a normal ordering sums the tadpole diagrams (see Fig.1) and the effect of the normal ordering is contained in the cut-off parameter $`\frac{1}{2\pi a}`$. Therefore, in the perturbative expansion of cosine term we have to exclude all the tadpole type diagrams, and then all other terms of perturbative expansion are finite and well-defined. In analogy with the sG model, we will associate the solition-like (quantum) solution of (10) with the renormalized electrons, and the breather-like $`\varphi _+`$ mode with the excitonic electron-hole bound state.

If the quantum correction due to $`\varphi _+^4`$ is neglected, the zeroth order dispersion relation of exciton state can be read off from the quadratic part of (23).

$$E^{(0)}(k)=\mathrm{}\left[v^2k^2(1+\frac{v_0}{2v}\mathrm{ln}\frac{1}{a|k|})+\frac{4tv}{\mathrm{}a}\right]^{1/2}.$$
(26)
The energy gap is given by

$$\mathrm{\Delta }_{ex}^{(0)}=\sqrt{\frac{4\mathrm{}tv}{a}}=2\sqrt{tW}.$$
(27)
Next, we compute the correction to the exciton excitation energy due to the quantum fluctuation of $`\varphi _+^4`$ term. The first such correction appears in the second order perturbation of $`\varphi _+^4`$. Because the $`\varphi _+^4`$term is normal ordered there are no tadpole diagrams. Only the ”sunset” diagram (see Fig.1) contributes, which is both ultraviolet and infrared convergent.

It contributes a negative correction to the gap as can be checked by a simple calculation $`e^{(\mathrm{\Delta }_{ex}^{(0)})^2\varphi ^2}(1+\lambda ^2\varphi ^2\varphi ^3\varphi ^3)e^{((\mathrm{\Delta }_{ex}^{(0)})^2\lambda ^2\varphi ^3\varphi ^3)\varphi ^2},`$ where $`\lambda t/\pi \mathrm{}a`$. The symmetry factor of the sun-set diagram is $`\frac{1}{4!}\frac{1}{4!}4432=1/6`$. In terms of self-energy the sunset diagram gives

$`\mathrm{\Sigma }(i\omega ,k)`$ $`=`$ $`{\displaystyle \frac{1}{6}}\left({\displaystyle \frac{t}{\pi \mathrm{}a}}\right)^2{\displaystyle 𝑑x𝑑\tau [D(x,\tau )]^3e^{ikx+i\omega \tau }},`$ (28)
$`D(x,\tau )`$ $`=`$ $`{\displaystyle \frac{d\omega dk}{(2\pi )^2}\frac{ve^{ikx+i\omega \tau }}{\omega ^2+v^2k^2(1+\frac{\alpha }{2}\mathrm{ln}\frac{1}{|k|a})+\frac{4tv}{\mathrm{}a}}}.`$ (29)
The excitation energy which includes the above quantum correction is given by

$$\mathrm{\Delta }_{ex}=\left[\frac{4\mathrm{}tv}{a}8\pi \mathrm{}^2v\mathrm{\Sigma }(0,0)\right]^{1/2}.$$
(30)
Write the above equation in the form

$$\mathrm{\Delta }_{ex}^2=(\mathrm{\Delta }_{ex}^{(0)})^2\left[1\frac{\mathrm{\Sigma }(0,0)}{t/(2\pi \mathrm{}a)}\right].$$
(31)
The second term in the bracket of Eq.(31) is the quantum correction in the dimensionless form. Explicitly, ($`y_1,y_2,\eta _1,\eta _2`$ are dimensionless).

$`\delta _{ex}={\displaystyle \frac{\mathrm{\Sigma }(0,0)}{t/(2\pi \mathrm{}a)}}={\displaystyle \frac{1}{4\pi }}{\displaystyle 𝑑y_1𝑑y_2\left[D(y_1,y_2)\right]^3},`$ (32)
$`D(y_1,y_2)={\displaystyle \frac{d\eta _1d\eta _2}{(2\pi )^2}\frac{e^{i\eta _1y_1+i\eta _2y_2}}{4+\eta _1^2+\eta _2^2f(|\eta _2|)}},`$ (33)
where $`f(|\eta _2|)=1+\frac{v_0}{2v}\mathrm{ln}[\frac{1}{|\eta _2|}\sqrt{\frac{W}{t}}]`$. We may extract the effective expansion parameter by rescaling. It is easy to see that

$$\delta _{ex}\frac{\delta \mathrm{\Sigma }(0,0)}{t/(2\pi \mathrm{}a)}\frac{v/v_0}{\mathrm{ln}W/t}.$$
(34)
We note that, in the perturbative regime $`\delta _{ex}0`$ the excitation energy approaches the value of the coefficient of the cosine term. This same feature also appeared in the exactly solvable sine-Gordon model. Our perturbative treatment for the excition states is thus valid when $`\delta _{ex}<1`$. This result suggests the possibility that when $`\delta _{ex}>1`$ is satisfied exciton instability may occur. However, we cannot exclude the possibility that the higher order corrections neglected in our perturbative calculation may prohibit such an instability.

The action Eq.(10) supports topological soliton excitations due to the topologically inequivalent vacuum of the cosine potential . Thus, we can not use the perturbative expansion of the cosine term in the study of the solition excitation. The original sine-Gordon model can be solved exactly, and the energy gaps of soliton and breather are known exactly (see Eq.(16), Eq.(17)) . Our action (10) is not exactly solvable, and here we will give only an estimate of the gap of soliton excitation. In estimating the gap of soliton excitation, it is sufficient to consider the static solition. In our perturbative regime we may neglect $`1`$ compared to $`\frac{V(k)}{\pi v}`$. Then, in real space, the energy functional for static solution reads

$`E[\varphi _+(x)]`$ $`=`$ $`{\displaystyle \frac{1}{8\pi }}{\displaystyle 𝑑x𝑑y\frac{e^2}{ϵ|xy|}_x\varphi _+(x)_y\varphi _+(y)}`$ (35)
$`+`$ $`{\displaystyle \frac{t}{\pi a}}{\displaystyle }dx(1\mathrm{cos}(\varphi _+(x)),`$ (36)
(A constant term is added to make the expression positive definite). We define a dimensionless variable $`\overline{x}`$, $`\overline{x}=x\left(\frac{t}{e^2ϵa}\right)^{1/2}=x\left(\frac{t}{E_c}\right)^{1/2}\frac{1}{a}`$. Then, the energy functional becomes

$`E[\varphi _+(x)]`$ $`=`$ $`\sqrt{{\displaystyle \frac{e^2t}{ϵa}}}[{\displaystyle \frac{1}{8\pi }}{\displaystyle }d\overline{x}d\overline{y}{\displaystyle \frac{1}{|\overline{x}\overline{y}|}}_{\overline{x}}\varphi _+(\overline{x})_{\overline{y}}\varphi _+(\overline{y})`$ (37)
$`+`$ $`{\displaystyle \frac{1}{\pi }}{\displaystyle }d\overline{x}(1\mathrm{cos}(\varphi _+(\overline{x}))].`$ (38)
The coefficient in front of bracket has the dimension of energy and it gives a characteristic gap scale of solitions. The ground state, which sits at a minimum of the cosine potential, has zero classical energy. The soliton excitation connects two adjacent classical ground states, and according the above estimate the gap of solition $`E_G`$ is of the order $`\sqrt{\frac{e^2t}{ϵa}}=\sqrt{E_ct}`$. Since we expect $`E_c>>t`$ the Coulomb interaction enhances the value of band gap significantly.

In summary, we have studied excitons in 1D narrow gap semiconductors of two anticrossing quantum Hall edges. According to our perturbative approach the exciton state may lie in the gap when $`\delta _{ex}<1`$. Our study indicates that many body interactions enhances the value of band gap significantly. Our result suggests that an exciton instability may occur when $`\delta _{ex}1`$. However, it is desirable to calculate how the higher order corrections neglected in our perturbative calculation may change this condition. The actual values of $`v`$, and $`t`$ in anticrossing edges states are not known well and it is difficult to estimate precisely the actual value of $`\delta _{ex}`$. Experimental observation of excitons and the investigation of exciton instability would be most interesting.

We are grateful to Al. B. Zamolodchikov and M.P.A. Fisher for useful comments. H. C. L. was supported by the Korea Science and Engineering Foundation (KOSEF) through the grant No. 1999-2-11400-005-5, and by the Ministry of Education through Brain Korea 21 SNU-SKKU Program. S.R.E.Y was supported by Korea Research Foundation Grant KRF-2000- 015-DP0125.