text:
# Influence of mass polydispersity on dynamics of simple liquids and colloids

## I Introduction

The dynamical behavior of liquids is an area of intense current interest. Much of this interest has been motivated by the desire to understand the progressively slower and more complex dynamics of dense, supercooled liquids as they are cooled toward the glass transition review . In the last few decades, numerous direct insights on dynamical motion in liquids have been obtained using molecular dynamics (MD) computer simulations, in which the spatial coordinates of particles as a function of time are calculated Allen . More recently, experimental studies of colloids have used confocal microscopy to track individual particles Kasper ; Marcus ; Weeks ; Kegel , thus generating the same type of data on microscopic particle motions as is obtained from simulations. For studying the glass transition, simulations and colloid experiments therefore serve as important model systems in which the implications of theory can be directly tested.

In both simulations and colloid experiments, fluids have been studied in which the particle size is polydisperse. In simulations, size polydispersity is often introduced to prevent crystallization of the deeply supercooled liquid (see e.g. KobAndersen ). In colloid experiments, some degree of polydispersity is always present, and depends on the process by which colloid particles are produced. To characterize the polydispersity of colloids, the distribution $`\theta (\sigma )`$ of particle diameters $`\sigma `$ is commonly found (or assumed) to be Gaussian:

$$\theta _G(\sigma )=\frac{1}{\delta \sqrt{2\pi }}\mathrm{exp}\left[\frac{1}{2}\left(\frac{\sigma \sigma _0}{\delta }\right)^2\right],$$
(1)
where $`\sigma _0`$ is the average particle diameter, and $`\delta `$ characterizes the width of the distribution Hunter . Polydispersity can then be quantified by the value of the dimensionless parameter $`c=\delta /\sigma _0`$.

In most experimental studies of colloids (see e.g. Hunter ) a system is regarded as effectively monodisperse if $`c<0.05`$. For many properties, such as the average liquid structure, this is a good approximation. However, dynamical properties, especially at the microscopic level, may depend sensitively on the nature of microscopic structural fluctuations, and so may be affected by even small polydispersities. In addition, size polydispersity in real colloids leads inevitably to a polydispersity of mass. However, some models of polydisperse liquids and colloids consist of systems in which particle size varies, but not particle mass KobAndersen ; DoliwaHeuer .

In this paper we seek to isolate and quantify the role of size and mass polydispersity on the dynamics of a simple liquid system, in particular to assess the need to incorporate mass polydispersity when simulating the dynamics of realistic systems having size polydispersity. To do so, we conduct MD simulations of a system of particles interacting via the Lennard-Jones (LJ) potential, polydisperse with respect to both mass and size, as a function of $`c`$. Our results show that a range of dynamical properties (the diffusion coefficient, the velocity autocorrelation function, and the non-Gaussian parameter) of a polydisperse fluid are systematically shifted from the corresponding monodisperse case. We discuss the implications of these results for observations of “dynamical heterogeneity” in simulations DoliwaHeuer ; kobdonati  and in colloid experiments Weeks ; Kegel .

## II Polydisperse Lennard-Jones Liquid

Since our aim is to study generic effects of polydispersity on liquid dynamics, we chose the well-studied LJ potential to model interparticle interactions. The LJ potential is popular for simulations of simple liquids, and there exist many studies to which to compare our results.

In simulations of colloids, the colloidal particles are often modeled as hard spheres, and for many cases this a good approximation. However, interactions among colloidal particles can take other forms, and can be explicitly controlled, for example, by attaching soluble polymer chains by one end to the particle surface to generate repulsion, or by adding nonadsorbing soluble polymers to the suspension to produce attraction GastRussel . Though the present work is motivated by the recent experiments studying the dynamics of colloidal particles, we do not address the question of how the behavior of a colloidal system depends on the shape of the interparticle interaction potential. We also do not take into account the influence of a solvent.

We perform equilibrium molecular dynamics simulations in three dimensions of a system of $`N=4000`$ particles interacting via the shifted-force LJ potential, a modification of the standard LJ potential,

$$V_{ij}(r)=4\epsilon \left[\left(\frac{\sigma _{ij}}{r}\right)^{12}\left(\frac{\sigma _{ij}}{r}\right)^6\right].$$
(2)
Here, $`V_{ij}`$ is the potential of interaction of two particles $`i`$ and $`j`$, separated by a distance $`r`$. $`\epsilon `$ characterizes the strength of the pair interaction and is constant for all particle pairs. In the shifted-force LJ interaction, the LJ potential and force are modified so as to go to zero continuously at $`r=2.5\sigma _0`$, and interactions beyond $`2.5\sigma _0`$ are ignored Allen .

Polydispersity is introduced through the particle size: $`\sigma _{ij}=(\sigma _i+\sigma _j)/2`$ where $`\sigma _i`$ ($`\sigma _j`$) characterizes the diameter of a particle $`i`$ ($`j`$). Particles are assigned $`\sigma `$ values by random sampling from the Gaussian distribution in Eq. 1. We also impose a mass polydispersity appropriate for the given size polydispersity. The mass of a particle $`i`$ is $`m_i=m_0(\sigma _i/\sigma _0)^3`$, where $`m_0`$ is the mass of a particle of size $`\sigma _0`$. Particle trajectories are evaluated using the leap-frog Verlet algorithm Allen , using the appropriate value of $`m_i`$ in the equation of motion of each particle.

Throughout this work we use reduced units. Energy is expressed in units of $`\epsilon `$, length in units of $`\sigma _0`$, the number density of particles $`\rho `$ in units of $`\sigma _0^3`$, and temperature $`T`$ in units of $`\epsilon /k`$, where $`k`$ is Boltzmann’s constant. Time $`t`$ is expressed in units of $`\sqrt{m_0\sigma _0^2/\epsilon }`$. In these units the time step used for integrating the particle equations of motion is $`0.01`$.

After equilibration, all quantities are evaluated in the microcanonical ensemble. We present data for $`\rho =0.75`$ and $`T=0.66`$, a state not far from the triple point of the monodisperse LJ fluid ($`\rho =0.85`$, $`T=0.76`$nic ; BoonYip . We chose this state point so as to avoid the dense, deeply supercooled liquid region of the phase diagram of the monodisperse LJ system, where spontaneous crystallization could interfere with the evaluation of equilibrium properties. We conduct separate simulations for $`c=0`$, $`0.05`$ and $`0.1`$.

## III Pair distribution function, diffusion coefficient, and velocity autocorrelation function

The pair distribution function $`g(r)`$ that characterizes the average liquid structure HansenMcDonald  is shown in Fig. 1 for each $`c`$ studied. The effect of increasing polydispersity is to reduce the height of, and broaden the peaks associated with the successive neighbor shells around each particle. However the mean position of each neighbor shell does not change noticeably.

We test for a dependence on $`c`$ of the bulk transport properties by evaluating the diffusion coefficient $`D`$. We obtain $`D`$ from $`r^2(t)`$ using the Einstein relation,

$$D=\underset{t\mathrm{}}{lim}\frac{r^2(t)}{6t}.$$
(3)
Fig. 2 shows the dependence of $`D`$ on $`c`$. We find that at fixed $`\rho `$ and $`T`$, $`D`$ decreases systematically by about 10% as $`c`$ increases from zero to $`0.1`$.

Fig. 3 shows the dependence on $`c`$ of the velocity autocorrelation function $`\psi (t)`$ BoonYip :

$$\psi (t)=\frac{𝐯(0)𝐯(t)}{|𝐯|^2},$$
(4)
where $`𝐯(t)`$ is the velocity of a particle at time $`t`$. $`D`$ is related to the integral of $`\psi (t)`$, and consistent with the decrease of $`D`$, the negative part of $`\psi (t)`$ becomes larger in magnitude with increasing $`c`$. This trend reflects an increase with $`c`$ of the strength with which single particle dynamical properties of the system are retained on a time scale comparable to the collision time.

## IV Non-Gaussian parameter

The general non-Gaussian parameter $`\alpha _n(t)`$ is defined for integers $`n1`$ as BoonYip ,

$$\alpha _n(t)=\frac{r^{2n}(t)}{c_nr^2(t)^n}1,$$
(5)
where $`c_n=[135\mathrm{}(2n+1)]/3^n`$. $`r^{2n}(t)`$ is the ensemble average of the $`2n^{\mathrm{th}}`$ power of the particle displacements after a time $`t`$ Rahman :

$$r^{2n}(t)=\frac{1}{N}\underset{i=1}{\overset{N}{}}|𝐫_i(t)𝐫_i(0)|^{2n}.$$
(6)
Here, $`𝐫_i(t)`$ denotes the position of particle $`i`$ after a time $`t`$ following a reference time $`t=0`$ in equilibrium. $`N`$ is the total number of particles in the system.

The functions $`r^{2n}(t)`$ also represent the even moments of $`G_s(𝐫,t)`$, the self-part of the van Hove correlation function HansenMcDonald . For an isotropic fluid made up of particles with spherically symmetric interactions, we can restrict our attention to $`G_s(r,t)`$, the probability density that a particle located at the origin at time $`t=0`$ will be found within $`dr`$ of a distance $`r`$ at time $`t`$ vanHove :

$$G_s(r,t)=\frac{1}{N}\underset{i=1}{\overset{N}{}}\delta (r|𝐫_i(t)𝐫_i(0)|).$$
(7)
In terms of $`G_s(r,t)`$, $`r^{2n}(t)`$ can be written,

$$r^{2n}(t)=4\pi _0^{\mathrm{}}r^{2n}G_s(r,t)r^2𝑑r.$$
(8)
In the case of an ideal gas of non-interacting particles having a Maxwell-Boltzmann velocity distribution HansenMcDonald , $`G_s(r,t)`$ is a Gaussian function of $`r`$:

$$G_s(r,t)=\left(\frac{\beta m}{2\pi t^2}\right)^{3/2}\mathrm{exp}\left(\frac{\beta mr^2}{2t^2}\right),$$
(9)
where $`\beta =1/kT`$. In this case, it is readily shown that $`r^{2n}(t)=c_nr^2(t)^n`$ and so $`\alpha _n(t)=0`$. For systems in which correlations of particle motions are prominent, $`G_s(r,t)`$ is not Gaussian, and the deviation of $`\alpha _n(t)`$ from zero serves to quantify the deviation of $`G_s(r,t)`$ from the Gaussian form.

In the present work we present results for $`\alpha _2(t)`$, the most commonly calculated non-Gaussian parameter:

$$\alpha _2(t)=\frac{3r^4(t)}{5r^2(t)^2}1.$$
(10)
In Fig. 4 we plot $`\alpha _2(t)`$ for three different values of polydispersity $`c=0`$, $`0.05`$ and $`0.1`$. Qualitatively, there are two effects induced by increasing polydispersity: (i) the characteristic, intermediate-time peak of $`\alpha _2`$, at approximately $`t=1`$, increases in magnitude as $`c`$ increases; and, (ii) the value of $`\alpha _2(t)`$ does not start from zero in the limit $`t0`$ when $`c0`$. We clarify each of these effects in turn below.

To distinguish the influence of mass and size polydispersity separately, we conduct two new simulations, one (“size-only”) for a system in which the size polydispersity is $`c=0.1`$, but in which all the particle masses are set equal to $`m_0`$; and another (“mass-only”) for which the mass polydispersity is taken from our previous “mass-and-size” $`c=0.1`$ case, but with all the particle sizes then set equal to $`\sigma _0`$. We compare in Fig. 5 the resulting behavior of $`\alpha _2`$ as a function of $`t`$ with the behavior found for the monodisperse $`c=0`$ case; and with the case where both size and mass are polydisperse with $`c=0.1`$.

First we focus on the behavior observed near the maximum of $`\alpha _2(t)`$ at approximately $`t=1`$. Although the mass-only curve in Fig. 5 lies above that of the monodisperse system, it still lies well below that corresponding to polydispersity of both mass and size. Hence, mass polydispersity is not solely responsible for the increase of the maximum of $`\alpha _2`$ with $`c`$. Interestingly, the size-only curve is also well below that corresponding to polydispersity of both mass and size. Even the sum of the deviations from the monodisperse case of the mass-only and size-only curves, is insufficient to account for the height of the curve for the system with polydispersity of both mass and size. That is, polydispersity of both mass and size together has a greater impact on dynamical properties at intermediate times than can be obtained from polydispersity of mass or size alone.

Next we turn our attention to the behavior of $`\alpha _2`$ as $`t0`$. (For the remainder of this work we will denote the limit as $`t0`$ of $`\alpha _2`$ as $`\alpha _2^{}`$.) In MD simulations of a one-component LJ system Rahman ; Desai  and of a binary LJ mixture KobAndersen ; kobdonati , $`\alpha _2^{}=0`$. However, in the binary LJ mixture studied, the two species differ in size only and have the same mass, in contrast to our system in which the masses of particles differ in accordance with the polydispersity of their sizes. Two of the curves in Fig. 5 ($`c=0`$ and “size-only”) correspond to systems with no mass polydispersity, and in both cases $`\alpha _2^{}=0`$. For the other two curves ($`c=0.1`$ and “mass-only”), $`\alpha _2^{}`$ adopts the same nonzero value. It is clear that the mass polydispersity is solely responsible for the behavior of $`\alpha _2^{}`$.

As $`t0`$, the atomic motions in the fluid correspond to those of free particles, and the distribution of velocities is the Gaussian function given in Eq. 9. For a monodisperse system, this means that $`\alpha _2^{}=0`$. For a system with polydisperse masses, each particle of a given mass also samples the Gaussian velocity distribution given in Eq. 9. However, the Gaussian distributions sampled will have different widths for particles of different $`m`$. Consequently, the form of the total $`G_s(r,t)`$ function is not in general Gaussian because it is a superposition of individual Gaussians of different width. The result is a non-zero value of $`\alpha _2^{}`$, as found in our simulations.

Since the limit $`t0`$ corresponds to the free-particle limit for atomic motion, we can calculate the non-Gaussian parameter for a polydisperse system as $`t0`$ exactly. Consider a system of $`N`$ particles in which there are $`M`$ species (labeled by index $`j`$) each having $`N_j`$ particles of mass $`m_j`$. The moments $`r^{2n}(t)`$ can be found using a modified form of Eq. 6 appropriate for an $`M`$-component system:

$$r^{2n}(t)=\frac{1}{N}\underset{j=1}{\overset{M}{}}\underset{i=1}{\overset{N_j}{}}|𝐫_i(t)𝐫_i(0)|^{2n}.$$
(11)
Eq. 11 can be rewritten as,

$$r^{2n}(t)=\underset{j=1}{\overset{M}{}}f_jr^{2n}(t)_j,$$
(12)
where $`f_j=N_j/N`$ is the fraction of particles of species $`j`$ and

$$r^{2n}(t)_j=\frac{1}{N_j}\underset{i=1}{\overset{N_j}{}}|𝐫_i(t)𝐫_i(0)|^{2n},$$
(13)
where the sum is over particles only of species $`j`$. For each species, the atomic motion as $`t0`$ is also described by Eq. 9 with the appropriate value of $`m=m_j`$, and hence, the moments $`r^2(t)_j`$ and $`r^4(t)_j`$ can be found in the limit $`t0`$ by substituting Eq. 9 into Eq. 8 for each $`j`$. The result is,

$$r^2(t)_j=\frac{3t^2}{\beta m_j},$$
(14)
and,

$$r^4(t)_j=\frac{15t^4}{\beta ^2m_j^2}.$$
(15)
The value of $`\alpha _2^{}`$ for the multicomponent system can then be found by using Eqs. 12, 14 and 15 in Eq. 10:

$$\alpha _2^{}=\frac{_{j=1}^Mm_j^2f_j}{\left(_{j=1}^Mm_j^1f_j\right)^2}1$$
(16)
This expression highlights that $`\alpha _2^{}`$ may not equal zero for a system with polydisperse masses.

If the polydispersity is expressed as a continuous distribution of masses $`\varphi (m)`$, Eq. 16 generalizes to,

$$\alpha _2^{}=\frac{_0^{\mathrm{}}m^2\varphi (m)𝑑m}{\left(_0^{\mathrm{}}m^1\varphi (m)𝑑m\right)^2}1.$$
(17)
Experimental studies of colloids typically characterize polydispersity not in terms of the masses, but in terms of the particle diameters, described by $`\theta (\sigma )`$. Assuming that the particle mass $`m`$ is proportional to $`\sigma ^3`$, and that $`\varphi (m)dm=\theta (\sigma )d\sigma `$, Eq. 17 becomes,

$$\alpha _2^{}=\frac{_0^{\mathrm{}}\sigma ^6\theta (\sigma )𝑑\sigma }{\left(_0^{\mathrm{}}\sigma ^3\theta (\sigma )𝑑\sigma \right)^2}1.$$
(18)
Note that the value of $`\alpha _2^{}`$ therefore depends only on the shape of the mass distribution function, and is otherwise constant for all choices of $`\rho `$, $`T`$ and interparticle interaction.

We apply the above result to the case of the Gaussian distribution of particle diameters given in Eq. 1. We substitute $`\theta =\theta _G`$ with $`\sigma _0=1`$ in Eq. 18 and calculate $`\alpha _2^{}`$ as a function of $`c`$ (Fig. 6). We evaluate the integrals in Eq. 18 numerically, replacing the limits of integration $`(0,\mathrm{})`$ with $`[0.01\sigma _0,2\sigma _0]`$. This avoids the divergence of the integrands at $`\sigma =0`$, and in any case is more physical, since a real distribution of particle sizes would have a non-zero lower bound, and a finite upper bound. As seen in Fig. 6, $`\alpha _2^{}c^2`$ for small $`c`$, but increases more rapidly than this for $`c>0.1`$ The predictions of Eq. 6 are in agreement with our simulation results. For $`c=0.05`$ Eq. 18 gives $`\alpha _2^{}=0.02344`$, while our simulation gives $`\alpha _2^{}=0.027`$; for $`c=0.1`$ we obtain $`\alpha _2^{}=0.10781`$ and $`0.113`$, respectively.

We can also use Eq. 18 to calculate $`\alpha _2^{}`$ for a system in which the distribution of sizes is not Gaussian. For example, Sear Sears  simulated a system of hard spheres, polydisperse in both size and mass, using a “hat” function of width $`w`$: $`\theta (\sigma )=(w\sigma _0)^1`$ for $`\sigma _0(1w/2)\sigma \sigma _0(1+w/2)`$, and $`\theta (\sigma )=0`$ otherwise, and $`m\sigma ^3`$. For this case, we are able to solve Eq. 18 exactly, giving,

$$\alpha _2^{}=\frac{(1+\omega /2)^5(1\omega /2)^5}{5\omega (1\omega ^2/4)}1.$$
(19)
For $`\omega =0.3`$, Eq. 19 gives $`\alpha _2=0.06916`$, while for $`\omega =0.7`$, $`\alpha _2=0.4222`$, in agreement with the simulation results in Fig. 6 of Ref. Sears .

## V Discussion

When $`\alpha _2`$ has been extracted via confocal microscopy in colloid experiments, values as high as $`\alpha _2^{}1.5`$ have been observed Weeks ; Marcus ; Kasper . In these studies the polydispersity ranged from $`c=0.01`$ to $`0.1`$. Kasper, et al. Kasper  observed that $`\alpha _2^{}`$ is not zero for all values of the volume fraction occupied by the colloidal particles. Marcus, et al. Marcus  found that $`\alpha _2^{}`$ is nonzero and increases with increasing volume fraction. Weeks, et al. Weeks  found that of $`\alpha _2^{}`$ is approximately constant for small volume fractions but grows for higher volume fractions. In general, $`\alpha _2^{}`$ was found to increase with the volume fraction occupied by the colloid particles, in contrast to the absence of any density dependence in Eq. 18. In the case of real colloids, the behavior of $`\alpha _2`$ as $`t0`$ is complicated by the fact that solvent-induced hydrodynamic forces among particles potentially introduce strong, short-time-scale correlations of particle velocities, invalidating the free-particle assumption that is the basis of Eq. 9. The large difference between the behavior of $`\alpha _2^{}`$ found for these systems, and that predicted by Eq. 18 demonstrates that polydispersity alone cannot account for the observed values of $`\alpha _2^{}`$, and that hydrodynamic effects indeed dominate the short-time dynamical behavior of real colloids.

At intermediate times, we find that the peak value of $`\alpha _2`$ increases with $`c`$. The maximum value of $`\alpha _2`$ has been shown kobdonati  to correlate to the degree of dynamical heterogeneity present in the system: that is, transient, spatially correlated groups of particles whose characteristic structural relaxation time differs from the mean. Our results therefore suggest that dynamical heterogeneity, prominently observed at lower $`T`$ and higher $`\rho `$ than studied here, may be enhanced as polydispersity increases. One source of this enhancement may be the influence of mass polydispersity on spontaneously occurring density fluctuations, that in turn control the development of dynamical heterogeneities. In general, the occurrence and quantification of dynamical heterogeneity in a polydisperse system is likely to be more complicated than in a monodisperse (or even bi-disperse) system. At the same time, since we observe a slowing of the dynamics with increasing polydispersity, the dynamical heterogeneity may be more prominent and longer-lived in a polydisperse system, and so may facilitate the study of these complex structures.

In summary, we have illustrated that a system of particles with polydispersity of both mass and size is a more realistic model (compared to models without mass polydispersity) for studying the dynamics of colloids in MD simulations. Our results show that typical polydispersities found in real systems can induce an influence of the order of 10% on dynamical properties. This is a small effect for studies, such as those near a glass transition, where relaxation times may vary by several orders of magnitude. At the same time, knowledge of the amount and direction of the impact of polydispersity on dynamics is required because polydispersity is so commonly found in systems studied both in simulations and experiments. This knowledge is also crucial for precise tests of theories, particularly those formulated for perfectly monodisperse systems. We also note that our results can be tested experimentally in colloids by deliberately varying the polydispersity of the studied colloidal system. Indeed, by varying $`c`$ alone it might be possible to study an “isothermal-isochoric glass transition” (i.e. a glass transition at both fixed $`T`$ and $`\rho `$) by setting up an appropriate series of colloidal systems where the polydispersity is progressively increased.

## VI Acknowledgments

NK thanks Sergiy Kuzmin for valuable discussions. NK and PHP acknowledge financial support from NSERC (Canada).


format: text

text:
# Potential formation sites of super star clusters in ultra-luminous infrared galaxies

## 1 Introduction

Recent high resolution imaging by $`Hubble`$ $`Space`$ $`Telescope`$ ($`HST`$) has revealed super star clusters (hereafter referred to as SSCs) with masses of $`10^510^9`$ $`M_{}`$, very young ages and large dust extinction in ultra-luminous infrared galaxies (ULIGs), most of which are ongoing mergers (e.g., Shaya et al. 1994 Surace et al. 1998; Surace & Sanders 1999; Scoville et al. 2000). These SSCs are frequently observed not only in merging galaxies with and without ultra-luminous infrared activity (e.g, in the “Antenna” by Whitmore et al. 1999) but also in starburst galaxies such as M82 (e.g., de Grijs et al. 2001), which implies that the presence of SSCs is a generic feature in intensively starbursting or merging galaxies such as ULIGs (Lutz 1991; Ashman & Zepf 1992; Surace & Sanders 1999). These discoveries of SSCs in ULIGs have raised the following three questions (e.g., Shaya et al. 1994; Forbes, Brodie, & Grillmair 1997; McLaughlin 1999; Larsen & Richtler 2000): (1) What is the dominant mechanism of SSC formation ? (2) How can the global dynamics of merging galaxies be associated with the SSC formation ? (3) Is the formation efficiency of SSCs different between merging galaxies and normal spirals ?

The purpose of this Letter is to investigate the formation and dynamical evolution of SSCs in ULIGs formed by gas-rich major mergers and thereby to try and provide some plausible answers for the above three questions. In particular, we adopt the plausible assumption that the rather high pressure of the warm interstellar gas ($`P_\mathrm{g}`$ $`>`$ $`10^5`$ $`k_\mathrm{B}`$: $`k_\mathrm{B}`$ is Boltzmann’s constant) can induce global collapse of giant molecular clouds to form massive compact star clusters corresponding to SSCs (Jog & Solomon 1992; Harris & Pudritz 1994; Elmegreen & Efremov 1997). Accordingly we investigate numerically the time evolution of gaseous pressure in gas-rich mergers with different merger parameters. Previous theoretical results have already demonstrated that if a compact stellar system has a very large number of stellar-mass black holes left after the normal processes of stellar evolution (in massive stars), these black holes can evolve into a single massive black hole (e.g., Quinlan & Shapiro 1989; Lee 1995). Based on simple analytic arguments, we suggest that SSCs formed in an ULIG can be transformed into a compact self-gravitating stellar system (i.e., a considerably large SSC) either with a few tens massive black holes (MBHs) with masses ranging from $`10^4`$ to $`10^6`$ $`M_{}`$ or with $`10^6`$$`10^7`$ stellar-mass black holes. Thus the present numerical results combined with analytical ones suggest that the dynamical evolution of SSCs in the central regions of ULIGs can be closely associated with the formation of super-massive black holes.

## 2 Model

Our numerical methods for modeling chemodynamical evolution of dusty starbursts associated with galaxy mergers have already been described by Bekki & Shioya (1998) Furthermore the details of the adopted TREESPH, which is based on that used in early numerical works (e.g., Heller & Shlosman 1994), are given by Bekki (1995). We accordingly give these numerical techniques only a brief review here. We construct models of galaxy mergers between gas-rich disks by using the model of Fall-Efstathiou (1980). The total mass and the size of a progenitor exponential disk are $`M_\mathrm{d}`$ and $`R_\mathrm{d}`$, respectively. From now on, all the masses and lengths are measured in units of $`M_\mathrm{d}`$ and $`R_\mathrm{d}`$, respectively, unless otherwise specified. Velocity and time are measured in units of $`v`$ = $`(GM_\mathrm{d}/R_\mathrm{d})^{1/2}`$ and $`t_{\mathrm{dyn}}`$ = $`(R_\mathrm{d}^3/GM_\mathrm{d})^{1/2}`$, respectively, where $`G`$ is the gravitational constant and assumed to be 1.0 in the present study. If we adopt $`M_\mathrm{d}`$ = 6.0 $`\times `$ $`10^{10}`$ $`\mathrm{M}_{}`$ and $`R_\mathrm{d}`$ = 17.5 kpc as fiducial values, then $`v`$ = 1.21 $`\times `$ $`10^2`$ km $`\mathrm{s}^1`$ and $`t_{\mathrm{dyn}}`$ = 1.41 $`\times `$ $`10^8`$ yr, respectively. The dark-to-disk halo mass ratio and the star-to-gas mass ratio are set equal to 4.0 and 9.0, respectively. Bulge component is not included in the present study. An isothermal equation of state is used for the gas with a temperature of $`7.3\times 10^3`$ K (corresponding to a sound speed of 10 km $`\mathrm{s}^1`$). The initial metallicity, $`Z_{}`$, for each gaseous particle in a given galactic radius $`R`$ (kpc) from the center of a disk is given according to the observed relation $`Z_{}=0.06\times 10^{0.197\times (R/3.5)}`$ of typical late-type disk galaxies (e.g., Zaritsky, Kennicutt, & Huchra 1994). We present the results mostly for merger models with nearly prograde-prograde orbital configurations, parabolic encounters, and a pericentric distance of 0.5 $`R_\mathrm{d}`$. We change the mass ratio ($`m_2`$) of the two merging disks and thereby investigate the dependence of SSC formation processes on $`m_2`$ over the range 0.1 (corresponding to minor merging) $``$ $`m_2`$ $``$ 1.0 (major). For comparison, we also show the results of an isolated disk model and a tidally interacting one (in which two galaxies do not merge with each other owing to the large pericentric radius of 2 $`R_\mathrm{d}`$ adopted for the parabolic orbit). Total particle number in each simulation is 23178 for the collisionless components and 20000 for the collisional ones.

We adopt the following two different star formation laws and thereby convert the gaseous components into collisionless new stellar particles. In the first law, we adopt the formation model for SSCs (globular clusters) of Harris & Pudritz (1994) in which interstellar gaseous pressure ($`P_g`$) in star forming regions of a galaxy can drive pressure confined, magnetized self-gravitating proto-cluster molecular clouds to collapse to form SSCs if $`P_g`$ is larger than the surface pressure ($`P_\mathrm{s}`$) of the clouds:

$$P_\mathrm{g}P_\mathrm{s}2.0\times 10^5k_\mathrm{B}.$$
(1)
In the simulations, a gas particle becomes $`one`$ new SSC if the particle’s pressure is larger than $`P_\mathrm{s}`$ = $`2.0\times 10^5k_\mathrm{B}`$. Although we cannot investigate the detailed physical processes of SSC formation in the present $`global`$ (from $``$ 100 pc to 10 kpc scale) simulation, the adopted phenomenological approach of SSC formation enables us to trace at least the most promising formation sites of SSCs in an admittedly plausible way. The other is the Schmidt law (Schmidt 1959) with the exponent of 1.5 (Kennicutt 1989) and newly formed stars due to the Schmidt law are referred to as “normal stars” (represented by NS from now on) for convenience. If both conditions are satisfied above, only SSC is assumed to be formed. The effects of radiative cooling and those of supernovae feedback effects on SSC formation will be explored in our future papers. Furthermore, if $`P_\mathrm{s}`$ = $`2.0\times 10^6k_\mathrm{B}`$, the formation efficiency is $``$ 6 times smaller.

## 3 Results

Figure 1 clearly demonstrates how the gaseous pressure of the interstellar medium becomes dramatically higher ($`10^4`$-$`10^5`$ $`k_\mathrm{B}`$) in some regions of the major merger model with $`m_2`$ = 1.0 at $`T`$ = 0.28 Gyr (starburst phase with the star formation rate of $``$ 75 $`M_{}`$ $`\mathrm{yr}^1`$) in comparison with the initial disks. Furthermore, gaseous pressure for some fraction of the gas particles (in particular, in the central regions) exceeds the threshold value of $`2.0\times 10^5k_\mathrm{B}`$ for SSC formation, essentially because rapid transfer of gas to the central region of the merger and efficient gaseous shock dissipation increase the gaseous density and pressure greatly. These results clearly explain the reason why SSCs can be formed in ULIGs and imply that the formation efficiency of SSCs are much higher in merging galaxies with massive starbursts than in “normal” disk galaxies. Figure 2 shows that SSCs are formed not only in the central region, where gas fueling is more efficient, but also in the outer tidal tails and “bridges” among the two cores, where the tidal force of the merger forms thin and high-density gaseous layers owing to gaseous shock dissipation. This result is qualitatively consistent with recent $`HST`$ observational results that SSCs are located both in the central regions (e.g., Arp 220 by Shaya et al. 1994) and in between the two merging cores (e.g., VV114 by Scoville et al. 2000) of ultra-luminous infrared galaxy mergers.

As is shown in Figure 3, the peak metallicity in the metallicity distribution of SSCs is fairly high (0.05 corresponding to 2.5 solar value and a factor of 2 higher than metal-rich globular clusters in ellipticals), firstly because the initial mean metallicity of the merger precursor disk is assumed to be also high (1.4 solar), and secondly because chemical enrichment due to starbursts makes the metallicity of SSCs considerably higher. We stress here that the derived metallicity distribution depends strongly on the initial metallicity distributions of the gaseous disks: The peak metallicity of SSCs is smaller for a merger with low initial gaseous metallicity. Figure 3 also demonstrates that the age distribution is rather narrow (compared with the distribution in metallicity), which reflects the fact that SSC formation is efficient only when two disks interact with each other so violently that gaseous pressure becomes drastically high ($`>2.0\times 10^5k_\mathrm{B}`$).

Figure 4 provides the following three predictions on the formation of SSCs. Firstly, the formation efficiency of SSCs (i.e., the ratio of SSC mass to initial gas mass) is more likely to be higher in a merger with larger $`m_2`$. Total SSC number (i.e., formation frequency) is larger for larger $`m_2`$ and the SSC formation is more widely spread spatially for larger $`m_2`$. Assuming the same orbital and spin geometry and orientation, tidal disturbance that takes places in the galaxy merger is stronger and the maximum star formation rate is higher for a merger with larger $`m_2`$ (Bekki & Shioya 2000). Therefore these results imply that a merger with a more strongly perturbed morphology and a more massive starburst population can show a higher efficiency of SSC formation. Secondly, the ratio of the total SSC mass to that of the new stellar components (both SSCs and NS) is more likely to be larger for a merger with larger $`m_2`$. This implies that the luminosity fraction of SSCs, which is observationally estimated to be less than 20 % for ULIGs (e.g., Shioya et al. 2001 for the case of Arp 220), can be higher for a merger with larger $`m_2`$. Thirdly, the higher SSC formation efficiency and the larger relative SSC mass ratio in major mergers than in tidally interacting galaxies imply that major mergers, which more dramatically transform galactic morphologies than tidally interacting ones (e.g. M82), provide the most promising formation sites of SSCs.

## 4 Discussion and conclusion

One of the most longstanding and remarkable problems in the formation and evolution of ULIGs is whether there is an evolutionary link between the two different types of nuclear activities, starburst and AGN, (e.g., Sanders & Mirabel 1996). Previous theoretical studies have attempted to understand either physical mechanisms for the efficient radial transfer of gas to the nuclear starburst regions in major mergers (e.g., Barnes & Hernquist 1992; Mihos & Hernquist 1996) or the growth process of the already existing seed massive black holes (with masses of $``$ $`10^6`$ $`M_{}`$; Norman & Scoville 1988). We discuss here this problem in the context of how a cluster of SSCs in the circumnuclear region of an ULIG evolves dynamically into a seed MBH. We suggest that the following time scale of dynamical friction (Binney & Tremaine 1987) is important for discussing this problem of SSC evolution:

$$t_{\mathrm{fric}}=3.9\times 10^8(\frac{r_i}{500\mathrm{p}\mathrm{c}})^2(\frac{V_c}{400\mathrm{k}\mathrm{m}\mathrm{s}^1})(\frac{10^7M_{}}{M_{\mathrm{cl}}})\mathrm{yr}.$$
(2)
Here we neglect the ln($`\mathrm{\Lambda }`$) term ($``$ 6.6 for a plausible set of parameters), and $`r_\mathrm{i}`$, $`V_\mathrm{c}`$, and $`M_{\mathrm{cl}}`$ are the initial radius (from the center of an ULIG), circular velocity, and mass of the SSC, respectively. This time scale argument implies that SSCs formed by dusty starbursts in major mergers can sink into the central $`1030`$ pc region and consequently merge with one another to form a single massive nuclear SSC within a few dynamical time scales of the ULIG (i.e., $``$ $`10^8`$ yr). We here point out that the fate of SSCs depends strongly on whether or not these are tidally disrupted in merging galaxies: It is also possible that SSCs are disrupted by the central strong tidal field of ULIGs before they reach the central few pc in ULIGs. However, recent observational (e.g., Okuda et al. 1990; Nagata et al. 1995) and theoretical (e.g., Zwart et al. 2000) results demonstrated that the central region of the Galaxy can easily harbor compact star clusters (corresponding to SSCs). Although dynamical conditions could be different between the central region of the Galaxy and those of mergers, these results on the Galactic SSCs imply that SSCs in ULIGs can also survive from tidal disruption owing to their compact nature and consequently reach the central few pc via dynamical friction.

A high-density compact stellar system has been demonstrated to form a single MBH as a result of the coalescence of stellar-mass black holes or neutron stars left after the normal processes of stellar evolution, within the dynamical relaxation time scale of $`t_{\mathrm{rel}}`$ (e.g. Lee 1995). Therefore, if the central region of an ULIG has $``$ $`10^2`$ SSCs with typical mass $`10^6`$ $`M_{}`$ and if $`t_{\mathrm{merge}}`$ $`>`$ $`t_{\mathrm{rel}}`$, a single very massive nuclear cluster with $``$ $`10^2`$ seed MBHs (each of mass $`10^5`$ $`M_{}`$) can be formed after the merging of the SSCs. Furthermore, these MBHs in the cluster may well sink into the central sub-pc region as a result of dynamical friction and consequently form a self-gravitating MBH cluster. Although it could be possible that multiple merging of these MBHs through gravitational radiation loss in this MBH cluster can form a single MBH, the evolution of such multiple MBH interaction is obviously beyond the scope of this paper (see more detailed discussions in Valtonen et al. 1994 and Xu & Ostriker 1994). On the other hand, if $`t_{\mathrm{merge}}`$ $`<`$ $`t_{\mathrm{rel}}`$, then a single massive cluster with huge number of ($``$ $`10^7`$) stellar-mass black holes can be formed. In this case also, successive merging of these black holes in a cluster can drive the formation of a single super massive black hole within the Hubble time if both the density and the velocity dispersion of the cluster are rather high (Quinlan & Shapiro 1989). Although these discussions are highly speculative, in each case, it seems to be inevitable that a central supermassive black hole forms as a natural result of dynamical evolution of SSCs in the nuclear region of an ULIG.

Recent dynamical studies based on $`HST`$ photometry and ground-based spectroscopy of the central bulge regions of 36 nearby galaxies demonstrated that the nuclear MBH (or massive dark object) to bulge mass ratio ($`f_{\mathrm{BH}}`$) is $``$ 0.006 (Magorrian et al 1998). If we adopt the assumption that most giant elliptical galaxies and bulges are formed by gas-rich major mergers (e.g., Toomre 1977) and successors of ULIGs (e.g., Sanders et al. 1988), we can discuss whether the MBH formation due to merging of SSCs in ULIGs is consistent with the observed $`f_{\mathrm{BH}}`$ value. For a spheroidal galaxy formed by merging,

$$f_{\mathrm{BH}}=f_\mathrm{g}f_{\mathrm{sf}}f_{\mathrm{ssc}}f_{sbh},$$
(3)
where $`f_\mathrm{g}`$, $`f_{\mathrm{sf}}`$, $`f_{\mathrm{ssc}}`$, and $`f_{sbh}`$ are the gas mass fraction, the mass ratio of new stars formed during merging (including SSCs) to gas mass, that of SSCs to the total mass of new stars, and that of the total mass of stellar-mass black holes to the total mass of a SSC. If we adopt a set of plausible values, $`f_\mathrm{g}`$ = 0.2 (from the typical value observed for late-type spirals), $`f_{\mathrm{sf}}`$ = 0.8 (from the present study), $`f_{\mathrm{ssc}}`$ = 0.2 (Shioya et al. 2001), and $`f_{sbh}`$ = 0.1 (Norman & Scoville 1988), then $`f_{\mathrm{BH}}`$ = 0.0032, which is similar to the value observed. The most important point in the above estimation is that the origin of $`f_{\mathrm{BH}}`$ can be closely associated with the formation efficiency of SSCs in ULIGs ($`f_{\mathrm{ssc}}`$), which has not been pointed out by previous studies. Finally, we suggest that the origin of the physical relationship between the masses of MBHs and those of galactic spheroids (e.g., Magorrian et al. 1998) can be understood in terms of the formation and evolution of SSCs.

We are grateful to the anonymous referee for valuable comments, which contribute to improve the present paper.


text:
# On surface energy of graphene and carbon nanoclusters

## I Introduction

Recently, attention of physicists and chemists was recalled to the problem of reconstruction of edges of planar (and curvilinear) graphite. Considerable amount of articles on nanocarbons, surmounted with recent results on nanotubes and graphite polyhedral crystals, shows that this question is far from being clear, despite intensive experimental and theoretical research has been conducted.

The aim of this paper is to get theoretical understanding of the formation of facets and edges of graphite, graphene (monolayer of graphite) and small carbon sheet fragments. I will use simple analytical approach implementing the continuum energetics model assumptions. This approach bases on the count of dangling bonds along the free edge of graphene. That gives an energy of the facetHenceforth I will use the term facet for the 2D planar face of 3D graphite crystal and 1D straight edge of 2D graphene sheet. That will not lead to misunderstanding because the latter is a particular case of the former when the number of layers equals one. as a function of its orientation. The hexagonal lattice of the graphene has a three–fold symmetry and allows me to consider only a sector of the full plane between 0 and $`2\pi /6`$.

The function of specific facet energy, $`F(\omega )=E(\omega )/L`$, is a periodic function of the orientation angle, $`\omega `$. At the fixed $`\omega `$ it has a weak dependence on the facet length, $`L`$, (in 3D case it could depend on the length of the sheet edge or on the free surface area, which is defined by the specific type of misorientation). This additional dependence is fully neglected if considering the continuous edge, owing to full neglect of the atomic structure. A microscopic view is that at the scale of interatomic distance, $`b`$, one can not choose an arbitrary orientation of the edge. The natural approach to the description of free misoriented facet is to consider a number of one–atom kinks distributed along the edge so that their density depends on the (continuous) orientation angle.

The orientational dependence of the facet energy appears to be important for the edge zipping model of the single–wall nanotube nucleation. In order to understand the experimental fact that synthesized SWNTs had indices \[n,n$`\pm \delta `$n\] close to armchair helicity To mark a SWNT I follow the conventional index agreement proposed in Ref.. the zipping of the graphite edges with different orientations has to be considered.

## II Chirality dependent formation energy

The surface energy of free graphene facet depends on its orientation via the energy of dangling bonds (DBs). In this context, the term free facet means that there is no special reconstruction of the edge bonds. Of course, the reconstruction can change the facet energy substantially. I consider some type of the 3D facet reconstruction in Ref.. The reconstruction of the 2D facet is somewhat more complicated as the system possesses less degree of freedom. It will be discussed at length elsewhere.

The microscopic picture of the free edge is presented in Fig.1 (upper). Along the graphene edge the dangling bonds are schematically drawn. The orientation of the main edge is chosen along \[$`2\overline{1}0`$\] zigzag direction §§§There are multiple possible choice of the unit vector diads; the second (tilde) set, shown in the Figure 1 on the right, gives the edge symbol \[$`\overline{1}\overline{1}0`$\]. Both choices are equivalent to , which is defined respectively to the basis lattice unit vectors $`\stackrel{}{c}_1`$ and $`\stackrel{}{c}_2`$ (the third unit vector $`\stackrel{}{c}_3`$ is normal to the graphite plane, so the corresponding lattice number is always equal to 0 for graphene). The density of DBs increases at such discontinuities of the edge as the kink shown in the Figure1.

The total energy of DBs constitutes the free surface energy:

$$E_{\mathrm{edge}}=N_{\mathrm{layers}}\nu (\omega )LE_b,$$
(1)
here the edge density of DBs, $`\nu (\omega )`$, the edge length, $`L`$, and the number of layers, $`N_{\mathrm{layers}}`$, give the total number of DBs. $`E_b2.355`$ eV is the single DB energy. I do not consider here any passivation/contamination of the graphite surface, which could change $`E_b`$.

The less the edge density, the more stable the surface. It is the zigzag face normal to  (or symmetry equivalent) direction. The edge, being scrolled and zipped (as shown in Fig.2), will form a cylindrical sleeve (SWNT nucleus). Its axis is along direction \[$`\overline{1}10`$\] in the graphene plane. This corresponds exactly to armchair nanotube conformation. I speculate that it gives a plausible explanation of the preferable formation of armchair SWNTs at certain conditions. The direct observation of graphene scrolling was published in Ref.. It gives also an experimental support that the stable graphene edge is the zigzag one.

The edges of zigzag and armchair types have the minimum and maximum densities of the DBs. It is possible to expand the DB density around extremum points and write an interpolation formula in between. Surprising cuspA similar nonmonothonic behavior is known for the function of the surface energy of perfect crystals for a century at least (G.S.Wulff, Z.Krist.Mineral. v. 34, 449, 1901). is found for the minimum, where the density changes linearly in $`\omega `$, the chirality angle. The latter is a standard measure of how far does the SWNT axis deviate from the high symmetry direction, for example, the zigzag direction. Because my edge structures have 1D geometry close to tubular one (Fig. 2), it is natural to characterize their misorientation with the same chirality angle.

Let me define the edge density through the nanotube–like indices n,m which are the numbers of unit vectors constituting the vector $`\stackrel{}{A}`$ joining the left and right ends of the considered edge structure \[n m 0\] (Fig.1). The DB density is the number of boundary sites divided by the length of the boundary:

$$|A|=\sqrt{3}b\sqrt{n^2+nm+m^2},$$
(2)
where $`b=|c|/\sqrt{3}1.41`$ Å is the carbon bond length. The infinitesimal small shift of the vector $`A=`$\[n,n\] (or \[n,0\]) defining the armchair (or zigzag) basal edge results in the change of the edge length and the number of DBs. I expand both in series on the small parameter $`1/n\omega `$. For an armchair edge it reads as:

$$\nu _A=\frac{2n+1}{3bn\sqrt{1+1/n+1/3n^2}}\frac{2}{3b}\left(1\frac{1}{2n}+\frac{1}{2n}\frac{1}{6n^2}+\frac{3}{8n^2}\frac{1}{4n^2}+\mathrm{}\right)=\frac{2}{3b}\left(1\frac{1}{36}\omega ^2+\mathrm{}\right).$$
(3)
Naturally, the first order terms cancel while the second order terms give the square dependence in the chirality angle, $`\omega \frac{1}{n}`$, around the maximum. The coefficient is $`k_A=\frac{1}{36}`$. The behavior of a nearly zigzag edge density is different. The non–zero expansion terms appear in the first order:

$$\nu _Z=\frac{n+1}{\sqrt{3}bn\sqrt{1+1/n+1/n^2}}\frac{1}{\sqrt{3}b}\left(1+\frac{1}{n}\frac{1}{2n}\mathrm{}\right)=\frac{1}{\sqrt{3}b}\left(1+\frac{1}{2\sqrt{3}}\omega +\mathrm{}\right).$$
(4)
Hence, the linear dependence coefficient is $`s_Z=+\frac{1}{2\sqrt{3}}`$ and the density of DBs has a cusp near the point of the minimum.

An edge, which has the intermediate crystallography index \[n m 0\] is nothing but a series of kinks (single atom steps) along the basal edge of a definite type (Fig.1).

Let me consider the (approximately) circular fragment of graphene, which is naively expected to deliver the minimum to the surface energy in a continuum approximation. The actual surface energy will oscillate along the perimeter with the period of $`\pi /3`$ between maxima corresponding to equivalent armchair edges (Figure 3).

At zero temperature, the purely zigzag facing of a 2D–crystal is expected. At the thermodynamic equilibrium, the zigzag–like edges are formed and the energetic cost for the vicinal edge is proportional to $`\mathrm{exp}(\nu LE_b/T)`$, where $`\nu `$ depends on the orientation (facing angle) linearly, see Eq.(4).

Let me quantify this. The energy gain for the nucleation of the vicinal zigzag edge \[n, $`\delta `$n, 0\] is as follows:

$$\delta E(\omega )=N_{\mathrm{layers}}\delta \nu (\omega )E_bLN_{\mathrm{layers}}L\nu _ZE_bs_Z\omega \frac{N_{\mathrm{layers}}L\omega }{6b}E_b.$$
(5)
Here I subtract the DB energy of the optimum zigzag edge, $`N_{\mathrm{layers}}L\nu _ZE_b`$. The chirality equals $`\omega =\delta n/n`$.

For the zipping model of SWNT nucleation one has $`N_{\mathrm{layers}}=2`$, it yields the following exponentially small probability of the nucleation of chiral SWNT with indices \[n, n$`\pm \delta `$n\] at the temperature $`T`$:

$$ww_o\mathrm{exp}[\frac{1}{T}\frac{\omega L}{3b}E_b]=w_oe^{{\displaystyle \frac{\omega }{\omega _c}}},$$
(6)
where $`w_o`$ is the probability of nucleation of armchair \[n,n\] SWNT, $`\omega _c={\displaystyle \frac{3b}{L}}{\displaystyle \frac{T}{E_b}}`$ is a maximum chirality at given $`L`$ and $`T`$. The higher the temperature, the shorter the nucleus length. Thus the $`\omega _c`$ has a maximum as a function of temperature. Let me estimate $`L50`$Å at $`T2000`$ K. Then the $`\omega _c0.006`$ and the probability to nucleate  tube is $`10^7`$ smaller than for  SWNT

Let me consider the edge closest to zigzag one. The smallest edge defect is a couple of kinks (per total nucleus length) with the energy equals $`E_b`$. This minimum energy gain corresponds to the relative nucleation probability $`w/w_oe^{E_b/T}10^6`$ at $`T=2000`$K.

The calculated probability of chiral SWNT nucleation is too small to describe chirality scattering of the experiment. The real growing edge shape is likely far from equilibrium (ideal straight line). Therefore, the kinetics of graphene edge growth defines both the length and the chirality of the nucleus. I did not consider here the dangling bond passivation. Whence the kink energy is $`E_b/2`$, the passivation increases the maximum allowed chirality that could also explain the experimental scattering of SWNT angles.

## III Summary

The microscopic model describing the energy of the edge formation for different types of graphitic structures is developed. The most stable edge has a zigzag orientation owing to the smallest density of the dangling bonds along the facet. The DB density is calculated depending on the misorientation angle for the facet of the arbitrary direction.

In respect to proposed earlier SWNT nucleation mechanism via formation of a cylindrical sleeve along the graphene edge, this paper allows to speculate that the prevalent tubes have almost the armchair type which is complementary to the most stable zigzag graphene edge. The prevalent radius of such a sleeve nucleus is close to the optimum radius $`R_m`$ of nano–arch structure considered in . My estimation gives a lower limit of the uniformity of the SWNT synthesis yield. In ideal conditions, the novel nucleation mechanism predicts formation of  SWNT (and no chirality scattering is expected) which is best approximation of continual energetics model to the most abundant  SWNT observed experimentally.

Note, that in this paper I did not consider any kinetics of the formation, that might change the result.


text:
# 1 Introduction

## 1 Introduction

### 1.1 Status of the charm quark mass calculations

During the last few years a great effort has been devoted to the precise determination of the quark masses. Recent reviews about the present situation concerning the computation of the light quark masses can be found in refs. . As for the heavies, most of the studies done so far were focused on determining the $`b`$-quark mass. The charmed quark escapes the precision computation mainly because it is too heavy for the chiral perturbation theory to apply, and yet too light for the heavy quark expansion to set in.

A complete account of the presently available estimates of the charm quark mass value is given in the PDG review . They estimate the charm quark mass to be in the range 111Ref.  also contains a complete list of references.

$`1.15\mathrm{GeV}m_c^{\overline{\mathrm{MS}}}(m_c)1.35\mathrm{GeV}.`$ (1)
Very recently, two new QCD sum rule computations of this quantity appeared . After improving the calculation of the moment QCD sum rules for the charmonium states, the new result of ref.  is $`m_c^{\overline{\mathrm{MS}}}(m_c)=1.23(9)`$ GeV. Adopting quite a different QCD sum rule technique, in ref.  $`m_c^{\overline{\mathrm{MS}}}(m_c)=1.37(9)`$ GeV was obtained.

On the side of the lattice QCD computations there were a few attempts to compute the charm quark.

* By combining the QCD sum rule methodology with the lattice QCD computation of the moments of the heavy-heavy vector current correlation function, the authors of ref.  obtained $`m_c^{\overline{\mathrm{MS}}}(m_c)=1.22(5)`$ GeV. More complete discussion of the systematic uncertainties involved in their computation has been made in ref. . Unfortunately, this (elegant) method has not been followed by the other lattice groups.
* From the lattice QCD with Wilson quarks and by using the Ward identities, in ref.  a much larger value has been obtained, namely $`m_c^{\overline{\mathrm{MS}}}(m_c)=1.71(3)(20)`$ GeV. If we use their recent re-estimate of the mass renormalization constant , the above value gets down to $`m_c^{\overline{\mathrm{MS}}}(m_c)=1.50(3)(18)`$ GeV.
* Finally, in ref.  the charm quark mass has been computed by adopting the Fermilab way of interpreting the lattice QCD beyond the lattice cut-off. The quoted result is $`m_c^{\overline{\mathrm{MS}}}(m_c)=1.33(8)`$ GeV.

Since the actual situation with the lattice results is not clear, we decided to take advantage of the presently available non-perturbatively determined parameters, which are necessary for the complete elimination of the discretization effects that are linear in the lattice spacing (i.e. of $`𝒪(a)`$), and to make an estimate of the charm quark mass from the data that we obtained by working with the (relatively) fine grained lattice ($`a0.07`$ fm). Our result is

$`m_c^{\overline{\mathrm{MS}}}(m_c)=1.26(3)(12)\mathrm{GeV}.`$ (2)
### 1.2 Computation of the charm quark mass on the lattice

To compute the charm quark mass, we rely on the standard two strategies:

* Vector Ward identity (VWI)

$`m_Q^{(\mathrm{VWI})}(\mu )=Z_m(\mu )m_Q^{(\mathrm{VWI})}=Z_m(\mu ){\displaystyle \frac{1}{2}}\left({\displaystyle \frac{1}{\kappa _Q}}{\displaystyle \frac{1}{\kappa _{cr}}}\right),`$ (3)
where the critical value of the Wilson hopping parameter ($`\kappa _{cr}`$) corresponds to the chiral limit. 222Unless the physical units are explicitely displayed, all the quark masses mentioned in this section are assumed to be in lattice units. Since the lattice cutoff is finite and the charm quark mass is not negligible, it is highly important to improve the renormalization constant out of the chiral limit:

$`Z_m(\mu )={\displaystyle \frac{1}{Z_S^{(0)}(\mu )}}(1+b_mm_Q^{(\mathrm{VWI})}),`$ (4)
where the value of $`b_m=b_S/2`$, has been determined non-perturbatively in ref. , and $`Z_S^{(0)}(\mu )`$ is the renormalization constant of the scalar density computed in the chiral limit, which will be given in the next section.
* Axial-vector Ward identity (AWI)

$`m_Q^{(\mathrm{AWI})}(\mu )+m_q^{(\mathrm{AWI})}(\mu )=\overline{Z}_m(\mu ){\displaystyle \frac{{\displaystyle \underset{\stackrel{}{x}}{}}_4A_4^\mathrm{I}(x)P^{}(0)}{{\displaystyle \underset{\stackrel{}{x}}{}}P(x)P^{}(0)}},`$ (5)
where $`A_\mu =\overline{Q}\gamma _\mu \gamma _5q`$ and $`P=\overline{Q}\gamma _5q`$ are the axial vector current and the pseudoscalar density, respectively. For the full elimination of $`𝒪(a)`$ effects, the bare lattice axial current is improved in the chiral limit as $`A_\mu ^\mathrm{I}=A_\mu +c_A_\mu P`$, with $`c_A`$ already determined non-perturbatively . We used the symmetric definition of the derivative, e.g. $`_4f=(f(t+1)f(t1))/2`$. The mass renormalization constant is equal to $`\overline{Z}_m(\mu )=Z_A/Z_P(\mu )`$ and it is improved as

$`\overline{Z}_m(\mu )={\displaystyle \frac{Z_A^{(0)}}{Z_P^{(0)}(\mu )}}\left(1+(b_Ab_P){\displaystyle \frac{m_Q^{(\mathrm{VWI})}+m_q^{(\mathrm{VWI})}}{2}}\right).`$ (6)
Again, the value of $`(b_Ab_P)`$ has been computed non-perturbatively , whereas the value of $`\overline{Z}_m^{(0)}(\mu )=Z_A^{(0)}/Z_P^{(0)}(\mu )`$, will be provided in the next section.

## 2 Lattice details and results

### 2.1 A few details on the lattice computation

We work with two sets of lattice data, generated at $`\beta =6.2`$ by using the non-perturbatively improved Wilson action ($`c_{SW}=1.614`$ ). Each set contains 200 independent $`SU(3)`$ gauge field configurations. The value of the critical parameter, $`\kappa _{cr}`$, is fixed by requiring the bare $`m_q^{(\mathrm{AWI})}0`$.

The values of $`\kappa _{cr}`$, along with the values of the bare quark masses for the degenerate light quark combinations, as well as for the non-degenerate (heavy-light) ones, are given in tab. 1. Each $`m^{(\mathrm{AWI})}`$ and $`m^{(\mathrm{VWI})}`$ are obtained as

$`m^{(\mathrm{AWI})}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{{\displaystyle \underset{\stackrel{}{x}}{}}_4A_4^\mathrm{I}(x)P^{}(0)}{{\displaystyle \underset{\stackrel{}{x}}{}}P(x)P^{}(0)}},`$ (7)
$`m^{(\mathrm{VWI})}`$ $`=`$ $`{\displaystyle \frac{1}{4}}\left({\displaystyle \frac{1}{\kappa _Q}}+{\displaystyle \frac{1}{\kappa _q}}{\displaystyle \frac{2}{\kappa _{cr}}}\right),`$ (9)
where in the improved axial current, $`A_4^\mathrm{I}=\overline{Q}\gamma _4\gamma _5q+c_A_4\overline{Q}\gamma _5q`$, we used $`c_A=0.04(1)`$ . The statistical errors quoted in this work are obtained by using the standard jackknife procedure (with 5 configurations per jack). The quark masses listed in tab. 1 are the bare lattice ones, that we now need to renormalize. We next discuss the computation of the mass renormalization constants.

### 2.2 Mass renormalization constants

To evaluate the mass renormalization constants in the chiral limit we use the non-perturbative method  which allows one to compute these constants in the continuum $`\mathrm{RI}/\mathrm{MOM}`$ scheme. The use of this method out of the chiral limit, however, requires the improvement of the off-shell quantities for which several new counterterms appear, each with a coefficient that is to be fixed non-perturbatively (for more details, see ref.). For that reason, in this paper, we will employ the method to compute the renormalization constants in the chiral limit only . To improve the renormalization constants out of the chiral limit, we will rely on the results of ref.  in which the corresponding coefficients were computed non-perturbatively with the results, $`b_m=0.69(1)`$, $`b_Ab_P=0.04`$, at $`\beta =6.2`$.

We first focus on the mass renormalization constant, $`Z_m^{(0)}(\mu )=1/Z_S^{(0)}(\mu )`$. To that end, we combine the quark propagators of 4 different values of $`\kappa _q`$ in 10 different combinations (degenerate and non-degenerate in light quark mass), compute the amputated vertex function $`\mathrm{\Gamma }_S(\kappa _i,\kappa _j;a\mu )`$ and impose the standard RI/MOM renormalization condition 333Besides the values of the Wilson hopping parameters corresponding to the light quark masses that we already listed in tab 1, for this computation we also use the quark propagator with $`\kappa _q=0.1333`$.

$`{\displaystyle \frac{1}{Z_S^{(0)}(\mu )}}=\underset{\kappa _{1,2}\kappa _{cr}}{lim}{\displaystyle \frac{1}{Z_S(\kappa _1,\kappa _2;\mu )}}=\underset{\kappa _{1,2}\kappa _{cr}}{lim}{\displaystyle \frac{\mathrm{\Gamma }_S(\kappa _i,\kappa _j;p)}{Z_q^{1/2}(\kappa _i;p)Z_q^{1/2}(\kappa _j;p)}}|_{p^2=\mu ^2},`$ (10)
where $`Z_q(\mu )`$ is the quark field renormalization constant which is easily extracted by imposing the vector Ward identity on the quark propagator. Such an obtained $`Z_m^{\mathrm{RI}/\mathrm{MOM}}(\mu )=1/Z_S^{\mathrm{RI}/\mathrm{MOM}}(\mu )`$ is then converted to the renormalization group invariant (RGI) form, $`Z_m^{\mathrm{rgi}}=Z_m^{\mathrm{RI}/\mathrm{MOM}}(\mu )/c^{\mathrm{RI}/\mathrm{MOM}}(\mu )`$, by using the available mass anomalous dimension coefficients up to 4-loops encoded in the function $`c^{\mathrm{RI}/\mathrm{MOM}}(\mu )`$, computed in the same $`\mathrm{RI}/\mathrm{MOM}`$ scheme 444The explicit form of the function $`c^{\mathrm{RI}/\mathrm{MOM}}(\mu )`$ is also given in Appendix of the present letter. For consistency, we also use the 4-loop expression for the running coupling , and set $`n_F=0`$ with (quenched) $`\mathrm{\Lambda }_{\mathrm{QCD}}^{(n_F=0)}=0.25`$ GeV . A typical situation is shown in fig. 1(left). For every mass combination, we then extract $`Z_m^{\mathrm{rgi}}(\mu )`$ by fitting to a constant on the plateau $`1.0(a\mu )^21.8`$. After extrapolating in $`a\overline{m}_q=(m_1^{(\mathrm{VWI})}+m_2^{(\mathrm{VWI})})/2`$ to the chiral limit, we get

$`Z_m^{(0)\mathrm{rgi}}=3.391(22).`$ (11)
For an easier comparison with the other determinations of this renormalization constant we express this result also in the $`\overline{\mathrm{MS}}`$ scheme:

$`Z_m^{(0)\overline{\mathrm{MS}}}(2\mathrm{GeV})=1.332(9).`$ (12)
As it can be observed from fig. 1(right), the extrapolation to the chiral limit in this case is very smooth.

Next, we discuss the computation of the second mass renormalization constant, i.e. the one needed to compute the quark mass by using the axial Ward identity ($`\overline{Z}_m^{(0)}(\mu )=Z_A^{(0)}/Z_P^{(0)}(\mu )`$). We use the proposal of ref.  which, by judiciously combining the Ward identities, allows one to alleviate the problem of the contamination by the Goldstone boson . The renormalization condition, by which this task is achieved for the mass renormalization constant, can be simply written as

$`{\displaystyle \frac{Z_A^{(0)}}{Z_P^{(0)}}}(\mu )=\underset{\kappa _{1,2}\kappa _{cr}}{lim}{\displaystyle \frac{Z_A}{Z_P}}(\kappa _1,\kappa _2;\mu )=\underset{\kappa _{1,2}\kappa _{cr}}{lim}{\displaystyle \frac{m_1^{(\mathrm{VWI})}\mathrm{\Gamma }_P(\kappa _1;p)m_2^{(\mathrm{VWI})}\mathrm{\Gamma }_P(\kappa _2;p)}{(m_1^{(\mathrm{VWI})}m_2^{(\mathrm{VWI})})\mathrm{\Gamma }_A(\kappa _1,\kappa _2;p)}}|_{p^2=\mu ^2}.`$ (13)
Obviously, we can use only the non-degenerate quark mass combinations which, for our 4 values of $`\kappa _q`$, means 6 combinations. As in the previous case, we convert our result from the $`\mathrm{RI}/\mathrm{MOM}`$ to the renormalization group invariant constant (at 4-loop level) and fit in the same window as before, $`1.0(\mu a)^21.8`$ (see fig. 1(left)). With such extracted values for $`\overline{Z}_m^{\mathrm{rgi}}`$, for each combination of the Wilson hopping parameters ($`\kappa _1,\kappa _2`$), we then extrapolate to the chiral limit. This is also illustrated in fig. 1(right). Contrary to the first case, the mass dependence of the $`\overline{Z}_m^{\mathrm{rgi}}`$ is more pronounced. We extrapolate to the chiral limit linearly (filled square in fig. 1) to get our central value. In addition, we perform the quadratic extrapolation (the result of which is depicted by an empty square in fig. 1), and the difference between this and the central value is incorporated in the systematic uncertainty. Our result is

$`\overline{Z}_m^{(0)\mathrm{rgi}}=3.303(26)_{0.051}^{+0.000},`$ (14)
which in the $`\overline{\mathrm{MS}}`$ scheme reads

$`\overline{Z}_m^{(0)\overline{\mathrm{MS}}}(2\mathrm{GeV})=1.297(10)_{0.020}^{+0.000}.`$ (15)
This result agrees well with the one of ref. , $`\overline{Z}_m^{(0)\overline{\mathrm{MS}}}(2\mathrm{GeV})=1.316(14)(17)`$.

### 2.3 Putting it all together

Now we combine all the results from tab. 1 with the renormalization constants discussed in the previous section to get the renormalization group invariant quark masses. At this stage (after including the renormalization constants) one can identify the results of the axial Ward identity as the sum of the heavy and the light quark mass. Since we also computed the light quark mass separately, we can now simply subtract it from the sum and work with the heavies only. The results are presented in tab. 2. In the same table we give the values of the heavy-light mesons for which the light quark mass has been interpolated to the light $`s`$-quark in a usual way (see e.g. ref. ). The mass of the $`D_s`$ meson in lattice units is $`M_{D_s}0.73(3)`$, where we use $`a^1(m_K^{})=2.7(1)`$ GeV. Thus the charm quark mass is to be found through an interpolation between the results for $`\kappa _Q=0.125`$ and $`\kappa _Q=0.122`$.

Notice that in tab. 2 we give also the results for a larger quark mass corresponding to $`\kappa _Q=0.119`$. Although this value is not necessary for our final result for the charm quark mass, it will be helpful in assessing the amount of the systematic uncertainties which will be discussed in the next subsection.

Now, to get the value of the charm quark mass, we need to interpolate in the heavy meson mass to $`M_{D_s}`$ and then simply read off the charm quark mass. To do so we need to choose an interpolating formula. We consider the following ones:

* $`M_{P_s}=a_0+a_1m_Q+a_2m_Q^2`$;
* $`M_{P_s}=b_0+b_1/m_Q+b_2/m_Q^2`$;
* $`M_{P_s}=c_0+c_1m_Q+c_2/m_Q`$.

The first (i) is the naive linear interpolation ($`a_2=0`$), the second (ii) comes from the heavy quark expansion and the third (iii) is the hybrid of the two. For (iii) we obviously need at least three points, and thus we have to use also the heaviest of our quarks from tab. 2. The complete situation is presented in tab. 3. Since the heavier quark is more prone to the $`𝒪\left((am)^2\right)`$ artifacts, we first concentrate our discussion on the results of the first two interpolations in which we set $`a_2=b_2=0`$.

The results of the two interpolation formulae are totally consistent with each other, and we choose to quote the first one as our central number.

To get the final result, that can be confronted to the results of other approaches, we need to convert our values to the $`\overline{\mathrm{MS}}`$ scheme and express it in the physical units by using $`a^1(m_K^{})=2.7(1)`$ GeV. After recalling that $`m_c^{\mathrm{rgi}}c^{\overline{\mathrm{MS}}}(\mu )=m_c^{\overline{\mathrm{MS}}}(\mu )`$, where the function $`c^{\overline{\mathrm{MS}}}(\mu )`$ is known to 4-loop accuracy  (see App. of the present letter), one can easily solve that equation to obtain the standard value $`m_c^{\overline{\mathrm{MS}}}(m_c)`$. Our results are:

$`\mathrm{Lattice}\mathrm{I}`$ $`m_c^{\overline{\mathrm{MS}}}(m_c)_{\mathrm{VWI}}=1.144(3)\mathrm{GeV},`$ (16)
$`m_c^{\overline{\mathrm{MS}}}(m_c)_{\mathrm{AWI}}=1.373(34)\mathrm{GeV},`$
$`\mathrm{Lattice}\mathrm{II}`$ $`m_c^{\overline{\mathrm{MS}}}(m_c)_{\mathrm{VWI}}=1.132(3)\mathrm{GeV},`$ (17)
$`m_c^{\overline{\mathrm{MS}}}(m_c)_{\mathrm{AWI}}=1.325(43)\mathrm{GeV}.`$
The above results are obtained by using $`n_F=0`$ and $`\mathrm{\Lambda }_{\mathrm{QCD}}^{n_F=0}=0.25`$ GeV. We see that for both sets of our lattice data the two equivalent methods (VWI and AWI) give different results. The reason for that discrepancy most probably comes from the lattice artifacts which are $`(am_Q)^n`$ ($`n2`$). One way of seeing that is to include the higher order effects at tree-level by employing the so-called EKLM factors  (see also the discussion in ). The leading effect of the EKLM factors to our result is $`(am_Q)^2`$ and it modifies the RGI charm quark mass as follows

$`m_c^{\mathrm{rgi}}[\mathrm{VWI}]\left[1+{\displaystyle \frac{(m_c^{(\mathrm{VWI})})^2}{12}}\right]m_c^{\mathrm{rgi}}[\mathrm{VWI}],`$
$`m_c^{\mathrm{rgi}}[\mathrm{AWI}]\left[1{\displaystyle \frac{(m_c^{(\mathrm{VWI})})^2}{6}}\right]m_c^{\mathrm{rgi}}[\mathrm{AWI}],`$ (18)
where $`m_c^{(\mathrm{VWI})}`$ stands for the bare lattice charm quark mass whose value is $`m_c^{(\mathrm{VWI})}=0.375(25)`$. When the above modification is included and we pass onto $`m_c^{\overline{\mathrm{MS}}}(m_c)`$, from our Lattice-I simulation, we get

$`m_c^{\overline{\mathrm{MS}}}(m_c)_{\mathrm{VWI}}=1.161(3)\mathrm{GeV},`$
$`m_c^{\overline{\mathrm{MS}}}(m_c)_{\mathrm{AWI}}=1.342(33)\mathrm{GeV}.`$ (19)
After comparing these to the results (16), we see that the inclusion of the tree level $`𝒪\left((am_Q)^2\right)`$ effects makes our two results getting closer to each other, although it is not sufficient to remove the bulk of $`𝒪\left(a^2\right)`$ corrections.

### 2.4 Systematic uncertainties and our final result

We will now briefly summarize the sources of systematic uncertainties and comment each one of them.

* The most important source of the systematic error are the lattice artifacts of $`𝒪((am)^n)`$ $`(n2)`$. It is therefore highly important to repeat our calculation on the lattice with a smaller lattice spacing. As our central result we will quote the average of the two methods (VWI and AWI), for both our sets of data as given in eqs. (16,17). The larger statistical error will be attributed to our final result while the difference between any of the two methods and the averaged one is included in the systematic uncertainty. In other words we have:

$`\mathrm{Lattice}\mathrm{I}`$ $`m_c^{\overline{\mathrm{MS}}}(m_c)=1.26(3)(11)\mathrm{GeV},`$
$`\mathrm{Lattice}\mathrm{II}`$ $`m_c^{\overline{\mathrm{MS}}}(m_c)=1.23(4)(10)\mathrm{GeV}.`$ (20)
The results of our two simulations are very close to each other and for our final estimate of the charm quark mass we will quote the value obtained from the first lattice, which has the larger temporal extension.
* In our analysis, we used the value of the lattice spacing $`a^1=2.71(11)`$ GeV, as obtained from the comparison of the physical $`m_K^{}`$ and the one we computed on the lattice. We checked that if, instead of the above value, we use $`a^1=2.83(15)`$ as obtained from the $`f_K`$ decay constant, the final value of the charm quark gets larger by $`2\%`$.
* To be conservative, when taking the average of the masses obtained by using the two equivalent methods, we quoted the larger statistical error. However, the final result does not contain $`1.5\%`$ of the error on the renormalization constant $`\overline{Z}_m^{\mathrm{rgi}}`$ which we discussed in the text (see eq. (14)). Although the value of the renormalization constant $`Z_m^{\mathrm{rgi}}`$ does not suffer from the same uncertainty, to be on the safe side, we will add $`1.5\%`$ of error to our final result.
* Whenever needed, we used $`\mathrm{\Lambda }_{\mathrm{QCD}}^{n_F=0}=0.25`$ GeV. Varying this quantity by $`10\%`$ allows one to cover all the presently available lattice estimates for $`\mathrm{\Lambda }_{\mathrm{QCD}}^{n_F=0}`$ . The impact of that variation on the final charm quark mass value is $`\pm 1.6\%`$.
* In ref. , the authors also estimate the discretization errors on the improvement coefficient $`b_Ab_P`$. We have varied this coefficient by $`\pm 0.03`$ which introduces the change in the central value for the charm quark mass by less that $`1\%`$ (more precisely by $`\pm 0.6\%`$).
* As mentioned in the previous subsection, the result obtained from the naive (linear) interpolation practically coincides with the one that we get when employing the interpolation motivated by the heavy quark expansion. Since we have the third (heavier) quark mass, we checked that the quadratic terms to both formulae (i,ii), or the use of the “hybrid” formula (iii) affects the value of the charm quark mass by no more than $`1\%`$ (see tab. 3).
* For determination of the charm quark mass we used the masses of the heavy-light pseudoscalar mesons with the light quark mass (linearly) interpolated to the strange one. We checked that the quadratic interpolation in the light quark mass to get $`M_{P_s}`$, does not make any influence on our final charm quark mass value. In addition, we verified that our results remain unchanged if we use the vector heavy-light mesons instead of the pseudoscalars.

We sum the above sources of errors quadratically and obtain

$`m_c^{\overline{\mathrm{MS}}}(m_c)=1.26(3)(12)\mathrm{GeV},`$ (21)
which is our final result. We note also that in passing from the RGI to the $`\overline{\mathrm{MS}}`$ value, we might have as well used $`n_F=4`$ and the value $`\alpha _s(m_\tau )=0.334(22)`$ . This, however, does not modify the above result, namely we obtain $`m_c^{\overline{\mathrm{MS}}}(m_c)=1.27(3)(11)\mathrm{GeV}`$.

## 3 Summary and Perspectives

In this letter we have computed the charm quark mass on the lattice by taking the advantage of the full improvement of the Wilson action and operators by which all the lattice artifacts linear in lattice spacing are absent. Therefore, the computation of the charm quark mass on a reasonably fine grained lattice is expected to lead to the result close to the continuum limit. Results of our simulations, obtained at single value of the lattice spacing, indicate that the lattice artifacts $`𝒪(a^2)`$ are sensible. As a consequence the two equivalent methods to compute the quark mass on the lattice yield different values. This fact largely dominates the systematic uncertainty of our calculation, and for the precise determination of the charm quark mass it is therefore important to go to ever smaller lattice spacings. Finner lattices, in turn, require larger lattices (to keep the same physical volume) and thus more powerful computing resources. By using APE1000 we plan to make such a study in the near future. We stress also that our computation has been performed in the quenched approximation. The (unknown) uncertainty introduced by quenching is not included in our systematic error estimate. A naïve guess points towards a decrease of the quenched value by $`5\%`$555 This “guesstimate” is based on the present observations according to which the unquenching of the light quark masses reduces their values by $`10\%`$, whereas $`m_b`$ gets smaller by less than $`2\%`$ .  However we prefer to quote our result as the quenched one and wait for the partially unquenched computations to assess the amount of this source of systematic errors.

## Appendix

The RGI quark mass is defined as

$`m^{\mathrm{rgi}}={\displaystyle \frac{m^{\overline{\mathrm{MS}}}(\mu )}{c^{\overline{\mathrm{MS}}}(\mu )}}={\displaystyle \frac{m^{\mathrm{RI}/\mathrm{MOM}}(\mu )}{c^{\mathrm{RI}/\mathrm{MOM}}(\mu )}},`$ (22)
where the functions $`c^{\overline{\mathrm{MS}}}(\mu )`$ and $`c^{\mathrm{RI}/\mathrm{MOM}}(\mu )`$ are known to 4-loop accuracy . For $`n_F=0`$, they read

$`c^{\overline{\mathrm{MS}}}(\mu )=a_s(\mu )^{4/11}\left(1+0.68733a_s(\mu )+1.51211a_s(\mu )^2+4.05787a_s(\mu )^3\right),`$ (23)
$`c^{\mathrm{RI}/\mathrm{MOM}}(\mu )=a_s(\mu )^{4/11}\left(1+2.02067a_s(\mu )+14.21925a_s(\mu )^2+138.30689a_s(\mu )^3\right),`$ (24)
where, for short, we write $`a_s(\mu )\alpha _s(\mu )/\pi `$.


text:
# Untitled Document

Final version: 27 April 2001

A giant stream of metal-rich stars in the halo of the galaxy M31

Rodrigo Ibata, Michael Irwin, Geraint Lewis, Annette Ferguson and Nial Tanvir

Observatoire de Strasbourg, 11, rue de l’Universite, F-67000 Strasbourg, France

Institute of Astronomy, Madingley Road, Cambridge, CB3 0HA, UK

Anglo-Australian Observatory, PO Box 296, Epping, NSW 1710, Australia

Kapteyn Astronomical Institute, Postbus 800, 9700 AV Groningen, The Netherlands

Physical Sciences, Univ. of Hertfordshire, Hatfield, AL10 9AB, UK

Recent observations have revealed streams of gas and stars in the halo of the Milky Way1-3 that are the debris from interactions between our Galaxy and some of its dwarf companion galaxies; the Sagittarius dwarf galaxy and the Magellanic clouds. Analysis of the material has shown that much of the halo is made up of cannibalized satellite galaxies2,4, and that dark matter is distributed nearly spherically in the Milky Way. It remains unclear, however, whether cannibalized substructures are as common in the haloes of galaxies as predicted by galaxy-formation theory5. Here we report the discovery of a giant stream of metal-rich stars within the halo of the nearest large galaxy, M31 (the Andromeda galaxy). The source of this stream could be the dwarf galaxies M32 and NGC205, which are close companions of M31 and which may have lost a substantial number of stars owing to tidal interactions. The results demonstrate that the epoch of galaxy building still continues, albeit at a modest rate, and that tidal streams may be a generic feature of galaxy haloes.

Within the framework of hierarchical structure formation, large spiral galaxies like the Milky Way or Andromeda arose from the merger of many small galaxies and protogalaxies6. Later in their evolution, spiral galaxies become the dominant component in such mergers, cannibalizing smaller systems that fall within their sphere of influence. The complete destruction of the victim is usually progressive, and may take several orbits. However, the stellar debris from the destroyed dwarf galaxy follows a similar orbital trajectory to the progenitor, which is likely to have started life far away from the place of its final demise, and so the tidally disrupted matter tends to be deposited over a broad range in distance from the larger galaxy. Over time, with the accumulation of many such mergers, large galaxies develop an extensive stellar and dark-matter ’halo’, the latter being by far the most massive component of the galaxy. Meanwhile, part of the (dissipative) gas component of the smaller galaxies feeds the growth of the disk of the larger galaxy. This is seen in numerical simulations of galaxy formation, which result in galactic haloes comprising clumps of dark matter5. If this prediction is correct, then haloes should possess significant substructure-in contrast to previous suggestions7, which predict the dark and luminous components of haloes to be distributed smoothly.

Most of the halo of the Milky Way is metal-poor and, to first order, smoothly distributed; but recent studies have shown that the halo contains non-negligible stellar substructure. Evidence for the phase-space clumping of halo stars has even been found in the solar neighbourhood, where about 10single ancient accretion event4. The discovery that the Milky Way is surrounded by a giant rosette-like stream originating from the Sagittarius dwarf galaxy2,3 shows that-on the largest scales-the structure of the halo is substantially ’streamy’: approximately half of the intermediate-age stars at distances beyond about 15 kpc belong to the Sagittarius stream. This also implies2 that the last large accretion was that of the Sagittarius dwarf, and that the Milky Way has not cannibalized many other small galaxies for about 7 Gyr. Unless there has been a continual accretion of highly dark-matter-dominated ’galaxies’, or of galaxies containing exclusively old stars, the formation of the Galactic halo must have been essentially complete at a point in time less than half the age of the Universe.

We now consider whether the behaviour of the Milky Way is unusual. The only external galaxy in which halo substructure has been reported is NGC5907, which possesses a gigantic extraplanar stellar stream8, substantially brighter than similar structures in the Milky Way. However, this galaxy has long been known to be peculiar9, having a red, and luminous, flattened ’halo’, probably the result of the strong interaction that deposited the stream. To understand whether stream-like substructure is the generic morphology of galaxy haloes, it is necessary to investigate other ’normal’ galaxies like the Milky Way. The prime target for such a study is the Milky Way’s ’sister’ galaxy, M31-the Andromeda nebula, which at a distance of $`780`$ kpc is the closest large galaxy10.

The Wide Field camera11 on the 2.5-m Isaac Newton Telescope (INT WFC) is a four-chip charge-coupled device (CCD) mosaic camera imaging about 0.3 degree2 per exposure. On the nights of 3-9 September 2000, this instrument was used to survey the southeastern half of the halo of the Andromeda galaxy. To tile this region out to a distance of 4 ($`55`$ kpc in projection) from the centre of M31 required 58 contiguous fields. Images were taken in the equivalent of Johnson visual V and Gunn i bands under good atmospheric conditions, with 85% of fields taken in photometric conditions with seeing better than 1.2 arcsec.

Previous imaging studies of the M31 halo and outer disk have either sampled the outer parts of M31 at only a few discrete locations12-14, or have taken a panoramic, but much shallower, view15,16. In contrast, our deep survey allows us to make an uninterrupted study of the spatial variations in stellar density as a function of magnitude and colour over a large fraction of the halo. The continuous nature of this survey enables us to distinguish local density enhancements from both the large-scale structure of the halo of M31 and the underlying foreground Galactic distribution of stars. The exposure time of 800-1,000 s per passband per field reaches i-band magnitude i = 23.5 and V-band magnitude V = 24.5 (with a signal-to-noise ratio of $`5`$) and allows detection of individual red-giant branch (RGB) stars to an absolute V-band magnitude $`M_V=0`$ and main-sequence stars to $`M_V=1`$ in the halo of M31.

The WFC pipeline processing provides internal cross-calibration for the four INT WFC CCDs (each of 4,096 $`\times `$ 2,048 pixels) at a level of about 1% within each pointing. Field-to-field variations in photometric zero points were calibrated and cross-checked using a combination of multiple nightly photometric standard sequence observations and the overlap regions between adjacent WFC pointings. The overall derived photometric zero-points for the whole survey are good to the level of 1-2% in both bands. (Details of the data processing and calibration will be presented elsewhere.) Objects were classified as noise artefacts, galaxies or stars according to their morphological structure on all the images.

Figure 1 presents a summary of the results of our survey. Visual inspection of the surface density of sources classified as star-like on the i-band images and with magnitudes and colours consistent with the known properties of RGB stars at the distance of M31, shows the presence of a stream-shaped over-density of sources in the halo close to, but distinct from, the minor axis of M31. The RGB stellar density in the halo increases on average by a factor of two in the on-stream regions, and is statistically significant at $``$50-sigma. For stars brighter than the tip of the RGB in M31, the spatial density of sources is smooth and shows no sign of the feature: the stream therefore cannot be a foreground Galactic population. The on-stream stars follow a similar colour-magnitude sequence to the RGB stars in the remainder of the halo of M31, but with evidence of an enhanced metallicity relative to the ’normal’ M31 halo population owing to their redder colours. The average V-band surface brightness of the stream is $`\mathrm{\Sigma }_\mathrm{V}30\pm 0.5`$ mag arcsec-2, and the stream extends out to the current limit of our survey, at a projected distance of about $`40\mathrm{kpc}`$.

The metallicity of the stream, deduced from the colour-magnitude diagrams displayed in Fig. 2, covers a broad range with a mean slightly more metal-rich than the Galactic globular cluster 47 Tucanae (whose metallicity, that is, whose ratio of iron to hydrogen compared to that of the Sun, is $`[Fe/H]=0.7`$). We note that stars of near-solar metallicity are concentrated in the stream. Our survey reveals that such high-metallicity stars are also sparsely distributed throughout large parts of the halo of M31, consistent with the results of previous studies12,17. This high metallicity, with a mean approximately 10 times greater than that of the Milky Way’s halo, and the overly large stellar density of the halo-a factor of 10 greater14 than that of the Galactic halo-have until now been a puzzle.

What is this stream-like feature? The outer regions of galactic disks are metal-poor18, so the metallicity of the feature argues against a direct association with the outer disk of Andromeda. Furthermore, if it were a disk structure, it would be located at an impossible de-projected radial distance of about 140 kpc from the centre of the galaxy. This distance constraint is weakened if the feature is part of a pronounced disk warp, but that possibility is highly unlikely given the observed elongated structure. The only plausible explanation is that it is part of a large stellar stream within the halo. This stream lies along a line connecting the Andromeda satellites M32 and NGC205, and is aligned with the direction of elongation of the outer isophotes of NGC205, suggesting a relationship between the Andromeda stream and these two dwarf galaxies. The total absolute magnitude of the stellar stream is $`M_V14`$, which we estimate from the differential luminosity function between equivalent on- and off-stream fields (plausible assumptions were made to allow for fainter stream stars, the incompleteness of the survey and the uniformity of the stream). This luminosity is a factor of about 10 lower than that of either M32 and NGC205, consistent with the possibility that the stream is the debris stripped from one (or both) of these two dwarf galaxies during a recent interaction with M31.

M32 and NGC205 are both unusual dwarf elliptical galaxies which lie at projected distances of, respectively, only 5 kpc and 9 kpc from the centre of M31. NGC205 is still active in forming stars, and has had a varied and complex star-formation history19. The central regions are also known to contain dust, H I gas and molecular gas20, and there is substantial morphological21 and kinematic22 evidence that indicates that NGC205 is being tidally distorted and potentially disrupted. The H I distribution shows a well defined velocity gradient (unlike the stars), and is less extended than the overall stellar content of the galaxy20,23. The clear inconsistency between the gas and stars suggests that the gas may have recently been captured, possibly from the disk of M31.

M32, on the other hand, does not appear to be substantially distorted, yet it does seem to have a significant population of intermediate-age stars in addition to a classical old stellar component24,25. There is also clear evidence of a large metallicity spread in the giant branch, with a mean just below solar, supporting a prolonged star-forming epoch26.

The broad agreement of the metallicity distributions of the stream stars and these two dwarf satellites, together with their alignment and physical proximity to M31, point to a common origin. Furthermore, the unusual properties of most of the rest of the halo of M31 is also consistent with its relatively recent tidal origin. It seems quite likely that the apparently peculiar properties of M31’s halo are simply the result of a prolonged, aggressive bout of tidal stripping from either one, or both, of its two nearest neighbour satellite galaxies.

If this interpretation is correct, the stream (and possibly most of the stellar halo) has to be the result of previous interactions with M31, as there is otherwise not enough time to spatially separate if from either of the two dwarf galaxies. Furthermore, by comparison with numerical models of the Sagittarius stream27,28, there will be a similar stream (either leading or trailing the progenitor dwarf(s), depending on the sense of their orbit) on the opposite side of M31.

During disk-crossing episodes, gas in the disk of M31 and in these dwarf galaxies is shocked, leading to episodic bursts of star formation, and gas exchange. In the dwarf galaxies, this gas may either be recycled from earlier generations of stars or accreted from the M31 disk. These complex interactions may explain the unusual star-formation history of the dwarf galaxies and the Andromeda stream. In turn, the regular impacts perturb the disk of M31, possibly enough to induce its warped shape29 (it is plausible that a similar process is responsible for the warp in the disk of the Milky Way30).

The discovery of the Andromeda stream in the first deep, panoramic survey of the Milky Way’s nearest large companion suggests that halo substructure in the form of giant tidal streams may be a generic property of large spiral galaxies, and that the formation of galaxies continues at a moderate pace up to the present day. Although the dwarf galaxies M32 and NGC205 may be the source of this material, we cannot rule out the alternative possibility that the stream may be the fossil remnant of a third cannibalized system, found at a stage intermediate to those seen within the Milky Way2,4.

To understand more fully the origin and evolution of the Andromeda stream, a number of follow-up observational programmes are required, including an extension of the panoramic survey not only to the northern regions of the halo of M31, but also out to larger radii, as it is evident in the present data that the stream extends beyond the 40-kpc limits of the survey. The proximity of M31 also provides us with an opportunity to undertake a spectroscopic survey of individual stars within the stream, allowing us to map its kinematic and chemical properties. As with the tidal material detected in the halo of our own Galaxy, studies of the Andromeda stream would allow us to map the distribution of dark matter within the halo of our nearest neighbouring galaxy, as well as furthering our understanding of the process of galaxy formation.

REFERENCES

1. Putman, M. E. et al. Tidal disruption of the Magellanic Clouds by the Milky Way. Nature 394, 752-754 (1998).
2. Ibata, R., Lewis, G. F., Irwin, M., Totten, E. & Quinn, T. Great circle tidal streams: evidence for a nearly spherical massive dark halo around the Milky Way. Astrophys. J. 551, 294-311 (2001).
3. Ibata, R., Irwin, M., Lewis, G. F. & Stolte, A. Galactic halo substructure in the Sloan Digital Sky Survey: The ancient tidal stream from the Sagittarius dwarf galaxy. Astrophys. J. 547, L133-L136 (2001).
4. Helmi, A., White, S. D. M., de Zeeuw, P. T. & Zhao, H. Debris streams in the solar neighbourhood as relicts from the formation of the Milky Way. Nature 402, 53-55 (1999).
5. Klypin, A., Gottl ber, S., Kravtsov, A. V. & Khokhlov, A. M. Galaxies in N-body simulations: Overcoming the overmerging problem. Astrophys. J. 516, 530-551 (1999).
6. Cole, S., Aragon-Salamanca, A., Frenk, C. S., Navarro, J. F. & Zepf, S. E. A recipe for galaxy formation. Mon. Not. R. Astron. Soc. 271, 781-806 (1994).
7. Eggen, O. J., Lynden-Bell, D. & Sandage, A. R. Evidence from the motions of old stars that the Galaxy collapsed. Astrophys. J. 136, 748-766 (1962).
8. Shang, Z. et al. Ring structure and warp of NGC 5907: Interaction with dwarf galaxies. Astrophys. J. 504, L23-L26 (1998).
9. Sackett, P. D., Morrison, H. L., Harding, P. & Boroson, T. A. A faint luminous halo that may trace the dark matter around spiral galaxy NGC 5907. Nature 370, 441-443 (1994).
10. Stanek, K. Z. & Garnavich, P. M. Distance to M31 with the Hubble Space Telescope and HIPPARCOS red clump stars. Astrophys. J. 503, L131-L134 (1998).
11. Irwin, M. & Lewis, J. INT WFS pipeline processing. New Astron. Rev. 45, 105-110 (2001).
12. Holland, S., Fahlman, G. G. & Richer, H. B. Deep HST V- and I-band observations of the halo of M31: Evidence for multiple stellar populations. Astron. J. 112, 1035-1045 (1996).
13. Rich, R. M., Mighell, K. J., Freedman, W. L. & Neill, J. D. Local Group populations with the Hubble Space Telescope. I. The M31 globular cluster G1=Mayall II. Astron. J. 111, 768-776 (1996).
14. Reitzel, D. B., Guhathakurta, P. & Gould, A. Isolating red giant stars in M31’s elusive outer spheroid. Astron. J. 116, 707-722 (1998).
15. Walterbos, R. A. M. & Kennicutt, R. C. Multi-color photographic surface photometry of the Andromeda galaxy. Astron. Astrophys. Suppl. 69, 311-332 (1987).
16. Innanen, K. A., Kamper, K. W., van den Bergh, S. & Papp, K. A. The optical warp of M31. Astrophys. J. 254, 515-516 (1982).
17. Durrell, P. R., Harris, W. E., Pritchet, C. J. & Davidge, T. Photometry of the outer halo of M31. Am. Astron. Soc. Meeting 195, 4.04 (1999).
18. Ferguson, A. M. N., Gallagher, J. S. & Wyse, R. F. G. The extreme outer regions of disk galaxies. I. Chemical abundances of H II regions. Astron. J. 116, 673-690 (1998).
19. Lee, M. G. Stellar population is the central region of the dwarf elliptical galaxy NGC 205. Astron. J. 112, 1438-1449 (1996).
20. Young, L. M. & Lo, K. Y. The neutral interstellar medium in nearby dwarf galaxies. II. NGC 185, NGC 205, and NGC 147. Astrophys. J. 476, 127-143 (1997).
21. Hodge, P. W. The structure and content of NGC 205. Astrophys. J. 182, 671-696 (1973).
22. Bender, R., Paquet, A. & Nieto, J. Internal stellar kinematics of three dwarf ellipticals in the Local Group. Astron. Astrophys. 246, 349-353 (1991).
23. Welch, G. A., Sage, L. J. & Mitchell, G. F. The puzzling features of the interstellar medium in NGC 205. Astrophys. J. 499, 209-220 (1998).
24. Freedman, W. L. Stellar content of nearby galaxies. II–the Local Group dwarf elliptical galaxy M32. Astron. J. 98, 1285-1304 (1989).
25. Davidge, T. J. The evolved red stellar content of M32. Publ. Astron. Soc. Pacif. 112, 1177-1187 (2000).
26. Grillmair, C. J. et al. Hubble Space Telescope observations of M32: The color-magnitude diagram. Astron. J. 112, 1975-1987 (1996).
27. Johnston, K. V., Spergel, D. N. & Hernquist, L. The disruption of the Sagittarius dwarf galaxy. Astrophys. J. 451, 598-606 (1995).
28. Ibata, R. A. & Lewis, G. F. Galactic indigestion: Numerical simulations of the Milky Way’s closest neighbor. Astrophys. J. 500, 575-590 (1998).
29. Schwarzschild, M. Mass distribution and mass-luminosity ratio in galaxies. Astron. J. 59, 273-284 (1954).
30. Ibata, R. A. & Razoumov, A. O. Archer of the Galactic disk? The effect on the outer HI disk of the Milky Way of collisional encounters with the Sagittarius dwarf galaxy. Astron. Astrophys. 336, 130-136 (1998).

Acknowledgements. This paper is based on observations made with the Isaac Newton Telescope operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias.

Figure 1. Surface density of red-giant branch (RGB) stars over the outer southeastern halo of M31. The surface-density maps produced by our survey are shown as ’tiles’ overlying an optical image centred on M31 (see below). The over-density of stars revealed in the present study is seen as a stream extending out of M31 close to its minor axis. The density of extended background galaxies (approximately 25,000 per degree2 to $`V<24`$ mag) becomes comparable to the stellar density at a projected distance of 20-25 kpc from the centre of M31 along the minor axis. A fraction of these galaxies are misclassified by the reduction algorithm, thereby contaminating the stellar sample. There is a further contamination from foreground Galactic disk dwarfs, which number $`<10`$% of the M31 halo population over the survey region. However, the contribution of these contaminating sources to the number density maps is easily removed. Indeed, the stream can be detected at signal-to-noise ratio $`>5`$ over tiny individual regions of area 0.01 degree2. The gaps to the northeast and southwest of the map correspond to fields taken in poor seeing conditions where the image quality criteria of the survey were not attained. The surface density maps are superimposed upon a photographic Palomar Sky Survey image of M31; clearly visible are the two companion galaxies M32 and NGC 205. The major and minor axes of M31 are displayed, and an ellipse denoting the size of a flattened ellipsoid (aspect ratio 3:5) of semi-major axis length 55 kpc has been superimposed to show the spatial extent of the survey. The location of the Andromeda stream is labelled. The overall shape of the M31 halo revealed by our observations appears rather boxy, with indications of other possible substructures-but these results require further study to confirm their reality.

Figure 2. The relation between V-i colour and i-band magnitude. Data are shown for a stream field (a), and an off-stream field (b). a and b both cover 0.3 degree2 of sky, and were taken in similar observing conditions (0.97 arcsec in i and 0.93 arcsec in V for the stream field, and 1.08 arcsec for both i and V for the off-stream field). a and b are also at a similar distance ($`2^{}`$) from the centre of M31, with the off-stream field almost exactly on the minor axis of the galaxy. The red-giant branch (RGB) stars of M31 are seen at magnitudes fainter than $`i=20.5`$. To guide the interpretation of these diagrams, we have superimposed on each panel the RGB tracks (shifted to the distance of M31 and corrected for extinction) of four well-studied galactic globular clusters of different metallicity (\[Fe/H\]); these are, from left to right, NGC6397 $`([Fe/H]=1.91)`$, NGC1851 $`([Fe/H]=1.29)`$, 47 Tuc $`([Fe/H]=0.71)`$ and NGC6553 $`([Fe/H]=0.2)`$. The difference between the stellar populations in a and b is displayed in c in the form of a binned colour-magnitude diagram. With much of the contamination removed in this way, we see that the stream contains a broad range of stellar populations, from metallicities similar to that of NGC1851 through to solar abundance, with a mean metallicity slightly higher than that of 47 Tuc. It is particularly difficult to constrain the distance to the stream owing to this metallicity spread; to within an uncertainty of possibly as much as 0.5 mag the stream RGB appears at the same position as the M31 RGB, which suggests a distance of $`d800\pm 200`$ kpc. (The photometric uncertainties are approximately 0.2 mag at the faint limit of the data in the diagrams; 0.1 mag uncertainties occur at $`V=23.7`$ mag, $`i=22.5`$ mag).