text:
# Epidemic spreading in correlated complex networks

## Abstract

We study a dynamical model of epidemic spreading on complex networks in which there are explicit correlations among the node’s connectivities. For the case of Markovian complex networks, showing only correlations between pairs of nodes, we find an epidemic threshold inversely proportional to the largest eigenvalue of the connectivity matrix that gives the average number of links that from a node with connectivity $`k`$ go to nodes with connectivity $`k^{}`$. Numerical simulations on a correlated growing network model provide support for our conclusions.

Statistical physics has witnessed in recent years a renewed interest in graph theory due to the discovery that many natural and artificial systems can be described in terms of complex networks, in which the nodes represent typical units, and the links represent the interactions between pairs of units Albert and Barabási (2002); Dorogovtsev and Mendes (2002). The term complex network has been coined to refer to networks that typically exhibit two distinct properties: (i) A scale-free (SF) connectivity distribution. If we define the connectivity distribution, $`P(k)`$, as the probability that a node is connected to $`k`$ other nodes, then SF networks are characterized by a power-law behavior $`P(k)k^\gamma `$, where $`\gamma `$ is a characteristic exponent. This property renders classical models of random graphs Erdös and Rényi (1960), described by an exponentially bounded connectivity, inappropriate for the description of many real networks. (ii) The small-world property Watts and Strogatz (1998), that is defined by an average path length—average distance between any pair of nodes—increasing very slowly (usually logarithmically) with the network size $`N`$.

Apart from the empirical characterization of real complex networks and the development of new models, accounting for the observed properties, the interest in this field has been also driven by the discovery of the profound and subtle effects that the connectivity has on the behavior of dynamical systems defined on top of complex networks. These effects are particularly interesting in the study of disease transmission in SF networks, relevant for the understanding of the spreading of computer viruses Pastor-Satorras and Vespignani (2002) and sexually transmitted diseases Liljeros et al. (2001). Indeed, it was first noted in Ref. Pastor-Satorras and Vespignani (2001a) that in uncorrelated SF networks with a connectivity exponent $`\gamma 3`$, epidemic processes do not posses an epidemic threshold below which diseases cannot produce a macroscopic epidemic outbreak or the development of an endemic state. This feature, observed in several epidemic models Pastor-Satorras and Vespignani (2001a, b); May and Lloyd (2001); Moreno et al. (2002); Newman (2001), is deeply rooted in the presence of very large connectivity fluctuations in infinite SF networks.

The study of epidemic spreading in uncorrelated complex networks (that is, in graphs in which the connectivity of any node is independent of the connectivity of its neighbors) has been proved to be extremely successful, providing, for instance, the first satisfactory explanation of the long-standing problem of the generalized low prevalence of computer viruses without assuming any global tuning of microscopic parameters Pastor-Satorras and Vespignani (2001a); White (1998). Nevertheless, it represents a first approximation to real networks, that neglects the possibility of correlations in the connectivity of the nodes. The existence and importance of connectivity correlations has been recently pointed out in the literature. In fact, it has become clear that these correlations are critical in the understanding of the hierarchical structure of the Internet Pastor-Satorras et al. (2001); Vázquez et al. ; Goh et al. (2002). On the other hand, some growing network models have been proposed Callaway et al. (2001), in which correlations are spontaneously generated and have important effects in the percolation transition.

In this paper we shall present a study of epidemic spreading in complex random networks in which there are explicit correlations among the node’s connectivities. We will consider in particular the subset of undirected Markovian random networks, that are completely defined by their connectivity distribution $`P(k)`$ and the conditional probability $`P(k^{}|k)`$ that a node of connectivity $`k`$ is connected to a node of connectivity $`k^{}`$. These two functions can have any form (including SF behavior), and are assumed to be normalized ($`_kP(k)=_k^{}P(k^{}|k)=1`$) and restricted by the connectivity detailed balance condition

$$kP(k^{}|k)P(k)=k^{}P(k|k^{})P(k^{})kP(k,k^{}),$$
(1)
where the symmetric function $`(2\delta _{kk^{}})P(k,k^{})`$ is the joint probability that two nodes of connectivity $`k`$ and $`k^{}`$ are connected. The Markovian nature of this class of networks implies that all higher order correlations can be expressed as a function of $`P(k)`$ and $`P(k^{}|k)`$, allowing an exact treatment of epidemic models at the mean-field (MF) level. It is worth noticing, however, that a more detailed description, in terms of a Langevin equation (to be reported elsewhere Boguñá and Pastor-Satorras ), yields exactly the same results, confirming the accuracy of the MF description. In this framework, the topologically relevant magnitude is the connectivity matrix $`C_{kk^{}}=kP(k^{}|k)`$, that measures the average number of links that go from a node with connectivity $`k`$ to nodes with connectivity $`k^{}`$. We will show that the epidemic threshold is related to the largest eigenvalue of this matrix. Extensive numerical simulations on a random correlated network model Callaway et al. (2001) confirm the predictions of the present analysis. During the completion of this work we became aware of two recent preprints Newman ; Eguíluz and Klemm  in which it is also highlighted the general role of correlations in spreading and percolation in complex networks.

In order to study the effects of connectivity correlations in epidemic spreading, we will focus in the standard susceptible-infected-susceptible (SIS) model Diekmann and Heesterbeek (2000). All the results, however, can be easily extended to the more general susceptible-infected-removed-susceptible (SIRS) model Boguñá and Pastor-Satorras . In the SIS model each node in the network represents an individual, and each link represents a connection along which the infection can propagate. Susceptible (healthy) nodes become infected with probability $`\nu `$ if they are connected to one or more infected nodes. On the other hand, infected nodes recover spontaneously with probability $`\delta `$. The ratio of this two rates defines an effective spreading rate $`\lambda =\nu /\delta `$ (without lack of generality, we set $`\delta =1`$). For homogeneous networks, in which each node has more or less the same number of connections, $`kk`$, a general result states the existence of a finite epidemic threshold, separating an infected (endemic) phase, with a finite average density of infected individuals, from a healthy phase, in which the infection dies out exponentially fast. In terms of the average density of infected individuals $`\rho (t)`$ (the prevalence) we can describe the SIS model in homogeneous networks at a MF level by the following rate equation Pastor-Satorras and Vespignani (2001b)

$$_t\rho (t)=\rho (t)+\lambda k\rho (t)\left[1\rho (t)\right].$$
(2)
In this equation we have neglected higher order terms, since we are interested in the onset of the endemic state, close to the point $`\rho (t)0`$. Also, we have neglected correlations among nodes. That is, the probability of infection of a new node—the second term in Eq. (2)—is proportional to the infection rate $`\lambda `$, to the probability that a node is healthy, $`1\rho (t)`$, and to the probability that a link in a healthy node points to an infected node. This last quantity, assuming the homogeneous mixing hypothesis Anderson and May (1992), is approximated for homogeneous networks as $`k\rho (t)`$, i.e. proportional to the average number of connection and to the density of infected individuals, and independent of the connectivity. From Eq. (2) it can be proved the existence of an epidemic threshold $`\lambda _c=k^1`$ Diekmann and Heesterbeek (2000), such that $`\rho =0`$ if $`\lambda <\lambda _c`$, while $`\rho (\lambda \lambda _c)`$ if $`\lambda \lambda _c`$. In this context, it is easy to recognize that the SIS model is a generalization of the contact process model, widely studied as the paradigmatic example of an absorbing-state phase transition to a unique absorbing state Marro and Dickman (1999).

For general complex networks, in which large connectivity fluctuations and correlations might be allowed, we must relax the homogeneous hypothesis made in writing Eq. (2) and work instead with the relative density $`\rho _k(t)`$ of infected nodes with given connectivity $`k`$; i.e. the probability that a node with $`k`$ links is infected. Following Refs. Pastor-Satorras and Vespignani (2001a, b), the rate equation for $`\rho _k(t)`$ can be written as

$$\frac{d\rho _k(t)}{dt}=\rho _k(t)+\lambda k\left[1\rho _k(t)\right]\mathrm{\Theta }_k(t).$$
(3)
In this case, the creation term is proportional to the spreading rate $`\lambda `$, the density of healthy sites $`1\rho _k(t)`$, the connectivity $`k`$, and the variable $`\mathrm{\Theta }_k(t)`$, that stands for the probability that a link emanating from a node of connectivity $`k`$ points to an infected site. In the case of an uncorrelated random network, considered in Refs. Pastor-Satorras and Vespignani (2001a, b), the probability that a link points to a node of connectivity $`k^{}`$ is independent of the connectivity $`k`$ of the node from which the link is emanating. Therefore, in this case $`\mathrm{\Theta }_k=\mathrm{\Theta }^{\mathrm{nc}}`$ is independent of $`k`$ and can be written as

$$\mathrm{\Theta }^{\mathrm{nc}}=\frac{1}{k}\underset{k^{}}{}k^{}P(k^{})\rho _k^{}(t),$$
(4)
since the probability that a node is pointing to a node of connectivity $`k^{}`$ is proportional to $`k^{}P(k^{})`$. Substituting the expression (4) into Eq. (3), one can solve for the steady state solution and find the existence of an epidemic threshold $`\lambda _c`$, below which there are no solutions with a nonzero value of $`\mathrm{\Theta }^{\mathrm{nc}}`$. The expression of the epidemic threshold for uncorrelated random networks is

$$\lambda _c^{\mathrm{nc}}=\frac{k}{k^2}.$$
(5)
For infinite SF networks with $`\gamma 3`$, we have $`k^2=\mathrm{}`$, and correspondingly $`\lambda _c^{\mathrm{nc}}=0`$. Finally, from the solution of $`\rho _k`$, one can compute the total prevalence $`\rho `$ using the relation $`\rho =_kP(k)\rho _k`$.

For a general network in which the connectivities of the nodes are correlated, the above formalism is not correct, since we are not considering the effect of the connectivity $`k`$ into the expression for $`\mathrm{\Theta }_k`$. This effect can be taken into account, however, for Markovian networks, whose correlations are completely defined by the conditional probability $`P(k^{}|k)`$. In this case, it is easy to realize that the correct factor $`\mathrm{\Theta }_k`$ can be written as

$$\mathrm{\Theta }_k(t)=\underset{k^{}}{}P(k^{}|k)\rho _k^{}(t),$$
(6)
that is, the probability that a link in a node of connectivity $`k`$ is pointing to an infected node is proportional to the probability that any link points to a node with connectivity $`k^{}`$, times the probability that this node is infected, $`\rho _k^{}(t)`$, averaged over all the nodes connected to the original node. Eqs. (3) and (6) define together the MF equation describing the SIS model on Markovian complex networks,

$$\frac{d\rho _k(t)}{dt}=\rho _k(t)+\lambda k\left[1\rho _k(t)\right]\underset{k^{}}{}P(k^{}|k)\rho _k^{}(t).$$
(7)
The exact solution of Eq. (7) can be difficult to find, depending on the particular form of $`P(k^{}|k)`$. However, it is possible to extract the value of the epidemic threshold by analyzing the stability of the steady-state solutions. Of course, the healthy state $`\rho _k=0`$ is one solution. For small $`\rho _k`$, we can linearize Eq. (7), getting

$$\frac{d\rho _k(t)}{dt}\underset{k^{}}{}L_{kk^{}}\rho _k^{}(t).$$
(8)
In the previous equation we have defined the Jacobian matrix $`𝐋=\{L_{kk^{}}\}`$ by

$$L_{kk^{}}=\delta _{kk^{}}+\lambda kP(k^{}|k),$$
(9)
where $`\delta _{kk^{}}`$ is the Kronecker delta function. The solution $`\rho _k=0`$ will be unstable if there exists at least one positive eigenvalue of the Jacobian matrix $`𝐋`$. Let us consider the connectivity matrix $`𝐂`$, defined by $`C_{kk^{}}=kP(k^{}|k)`$. Using the symmetry condition Eq. (1), it is easy to check that if $`v_k`$ is an eigenvector of $`𝐂`$, with eigenvalue $`\mathrm{\Lambda }`$, then $`P(k)v_k`$ is an eigenvector of the transposed matrix $`𝐂^T`$ with the same eigenvalue. From here it follows immediately that all the eigenvalues of $`𝐂`$ are real. Let $`\mathrm{\Lambda }_m`$ be the largest eigenvalue of $`𝐂`$. Then, the origin will be unstable whenever $`1+\lambda \mathrm{\Lambda }_m>0`$, which defines an epidemic threshold

$$\lambda _c=\frac{1}{\mathrm{\Lambda }_m},$$
(10)
above which the solution $`\rho _k=0`$ is unstable, and another nonzero solution takes over as the actual steady-state—the endemic state.

It is instructive to see how this general formalism recovers previous results Pastor-Satorras and Vespignani (2001a, b), implicitly obtained for random uncorrelated networks. For any random network, in which there are no correlations among the connectivities of the nodes, we have that the connectivity matrix is given by $`C_{kk^{}}^{\mathrm{nc}}=kP(k^{}/k)kk^{}P(k^{})/k`$, since the probability that a link points to a node of connectivity $`k^{}`$ is proportional to $`k^{}P(k^{})`$. It is easy to check that the matrix $`\{C_{k^{}k}^{\mathrm{nc}}\}`$ has unique eigenvalue $`\mathrm{\Lambda }_m^{\mathrm{nc}}=k^2/k`$, corresponding to the eigenvector $`v_k^{\mathrm{nc}}=k`$, from where we recover the now established result Eq. (5).

In order to check the theoretical prediction Eq. (10), we have performed numerical simulations of the network model proposed by Callaway et al. Callaway et al. (2001). This model is very simply defined: each time step we add a new node, and, with probability $`\delta `$, two nodes are randomly selected and joined with a link. This recipe yields a random network with a exponential connectivity distribution, $`P(k)=(2\delta )^k/(1+2\delta )^{k+1}`$, which shows correlations among the connectivities of the nodes. These correlations can be analytically computed by means of the joint probability $`P(k,k^{})`$, which fulfills the recursion relation Callaway et al. (2001)

$`P(k,k^{})`$ $`=`$ $`{\displaystyle \frac{2\delta }{1+4\delta }}[P(k1,k^{})+P(k,k^{}1)]`$
$`+`$ $`{\displaystyle \frac{P(k1)P(k^{}1)}{1+4\delta }}.`$
Even though it is not possible to solve for $`P(k,k^{})`$ in a closed form, we can numerically estimate the largest eigenvalue of the connectivity matrix $`𝐂`$ by generating a finite matrix $`P(k,k^{})`$ and using the relation (1) to find $`P(k^{}|k)`$. The correctedness of this numerical estimate is ensured by the exponential decay of $`P(k,k^{})`$ for large values of $`k`$ and $`k^{}`$, as can be seen from the recursion relation. The estimate of the largest eigenvalue of the connectivity matrix for Callaway’s model is $`\mathrm{\Lambda }_m6.47656`$, which yields an epidemic threshold $`\lambda _c0.15`$. This value is to be compared with the prediction for an uncorrelated network, $`\lambda _c^{\mathrm{nc}}=k/k^2=0.20`$.

Figure 1 shows the results of numerical simulations of the SIS model on Callaway’s network, as well as on random networks with the same connectivity distribution, generated using the Molloy and Reed (MR) algorithm Molloy and Reed (1995, 1998). The MR algorithm generates a random network with a prescribed connectivity distribution and no correlations among nodes, and thus it is expected to yield an epidemic threshold given by $`\lambda _c^{\mathrm{nc}}`$. Simulations were performed for a fixed value $`\delta =1`$ in networks of size up to $`N=10^7`$, averaging over at least $`100`$ different starting configurations, performed on at least $`10`$ different realizations of the network. Figure 1 depicts the steady state prevalence as a function of the spreading rate $`\lambda `$. For the MR network, the function $`\rho (\lambda )`$ shows a clear linear behavior. The epidemic threshold estimated from a least-squares fitting is $`\lambda _c=0.21\pm 0.01`$, in excellent agreement with the prediction for an uncorrelated network. On the other hand, Callaway’s network exhibits a very different behavior, which might be compatible with the presence of a transition of infinite order Marro and Dickman (1999). In fact, Fig. 1 is reminiscent of the behavior found in Ref. Callaway et al. (2001) for the giant component of the network as a function of the parameter $`\delta `$. In that work, the size of the giant component was fitted to an stretched exponential with exponent $`1/2`$. Guided by this intuition, in Fig. 2 we perform a fit of the stationary prevalence in Callaway’s model to the form $`\rho (\lambda )\mathrm{exp}[\alpha (\lambda \lambda _c)^{1/2}]`$ Callaway et al. (2001). The fit yields a prefactor $`\alpha =1.52\pm 0.02`$, and an epidemic threshold $`\lambda _c=0.11\pm 0.02`$, smaller by a factor $`2`$ than the value corresponding to an uncorrelated random network, and in quite good agreement with the prediction from the largest eigenvalue of the connectivity matrix.

In summary, we have shown that, in the presence of correlations, the epidemic threshold in complex networks is determined by the connectivity matrix $`𝐂`$, and not by the connectivity distribution $`P(k)`$, as happens in uncorrelated networks. This fact implies that the previously predicted null epidemic threshold for SF networks with $`\gamma 3`$ might be shifted in correlated graphs, attaining a positive value depending on the nature of the correlations as given by the connectivity matrix. At this respect, it might be surprising that some complex networks, such as the Barabási-Albert (BA) graph Barabási and Albert (1999), are exactly described at the uncorrelated level given by Eq. (3) with $`\mathrm{\Theta }_k`$ independent of $`k`$ Pastor-Satorras and Vespignani (2001a, b). This fact must be taken as an evidence of the lack of correlations in the BA model; lack that, on the other hand, has been already checked numerically in Ref. Vázquez et al. . The formalism presented in this paper represents a refinement over previous works because it includes the effects of correlations between pairs of nodes and, in this sense, it is exact for Markovian networks. Real networks, such as the Internet, however, posses a more complex correlation structure. Our formalism will provide an improved approximation to epidemic dynamics in these cases, but it still remains the task of ascertaining the effects of higher order correlations. Further work is necessary in this direction.

###### Acknowledgements.

This work has been partially supported by the European Commission - Fet Open project COSIN IST-2001-33555. R.P.-S. acknowledges financial support from the Ministerio de Ciencia y Tecnología (Spain). We thank A. Vespignani for helpful comments and discussions.


format: text

text:
# EFFECTS OF A BURST OF FORMATION OF FIRST-GENERATION STARS ON THE EVOLUTION OF GALAXIES

## 1 INTRODUCTION

In the standard big bang cosmology, the amount of heavy elements synthesized in the universe until the decoupling of matter from primordial radiation is negligibly small; i.e., $`Z<10^{10}`$ where $`Z`$ is the mass fraction of heavy elements (Wagoner, Fowler, & Hoyle 1967). On the other hand, gas clouds in the present-day galaxies have an appreciable amount of heavy elements $`Z0.01`$; e.g., the solar heavy element abundance is $`Z_{}0.02`$. This means that chemical and thermal properties of gas clouds in primeval galaxies at high redshift are significantly different from those in the present-day galaxies.

Just after the recombination of the universe, the hydrogen gas could be in atomic phase at a temperature of $`10^4`$ K. Then the formation of primordial hydrogen molecules can proceed through the gas phase reaction (e.g., Peebles & Dicke 1968). Since the rotational and vibrational transitions in hydrogen molecules can cool the gas further, the first generation stars can be born in such gas clouds at redshift $`10`$ \- 30. Since the temperature is still much higher than a typical kinetic temperature of molecular gas clouds in the present-day galaxies (e.g., $`10`$ K), the initial mass function (IMF) in primeval galaxies may be different from that in the present-day galaxies; i.e., a Salpeter-like IMF (Salpeter 1955; Scalo 1986). Therefore, one of fundamental questions on the star formation in primeval galaxies is; “What kind of stars (i.e., Population III stars) were made in such extremely metal deficient gas clouds ?” \[Matsuda, Sato, & Takeda 1969; Yoneyama 1972; Yoshii & Saio 1986 (hereafter YS86); Silk 1977a; Tegmark et al. 1997; Nakamura & Umemura 1999, 2001a, 2001b; Bromm et al. 1999, 2001, 2002; Coppi et al. 2001\].

Applying the opacity-limited fragmentation theory \[the fragmentation of gas clouds due to the gravitational instability continues until the fragments become to be optically thick to the cooling radiation (Silk 1977b, 1977c)\], YS86 first derived an IMF for a hydrogen-helium gas cloud. They found that typical masses of the first-generation stars are from several $`M_{}`$ to 10 $`M_{}`$. Nakamura & Umemura (1999) obtained similar results based on their one-dimensional hydrodynamical calculations together with nonequilibrium processes for hydrogen molecule formation. Abia et al. (2001) pointed out that the very large C and N enhancements observed in the extremely metal-poor stars in our Galaxy favor a Population III IMF peaked at intermediate-mass stars.

Recently, on one hand, some groups investigated the collapse and fragmentation of primordial, metal-free gas by numerical simulations (Abel, Bryan, & Norman 2000, 2002; Bromm et al. 1999, 2001, 2002; Coppi et al. 2001). They showed that the first generation gaseous clumps have masses of $`1000M_{}`$ which corresponds to the Jeans mass of the gas with temperature $`T`$ a few 100 K and the density $`n10^310^4\mathrm{cm}^3`$. As Bromm et al. (2002) assumed, if the clump masses are indicative of the final stellar mass, Population III stars may be more massive than $`100M_{}`$. Nakamura & Umemura (2001a, 2001b) have developed their hydrodynamical simulations and found that the initial mass function of Population III stars may be bimodal with peaks of $`1`$$`10M_{}`$ and $`100M_{}`$. The mass depends on the initial density and the initial fraction of molecular hydrogen ($`x_{\mathrm{H}_2,0}`$) (Nakamura & Umemura 2001b). If $`x_{\mathrm{H}_2,0}`$ is larger than $`3\times 10^3`$ for low-density ($`n10^5`$ cm-3), the effective HD cooling makes the high mass peak small to $`10M_{}`$. However, if $`x_{\mathrm{H}_2,0}`$ is smaller than $`3\times 10^3`$, very massive stars with $`100M_{}`$ are formed preferentially.

In this way, the mass of first-generation stars is controversial. It seems also hard to examine observationally which mode is more feasible; i.e., Population III stars are either intermediate-mass star or very massive ones. If most Population III stars are as massive as $`100M_{}`$, they could explode quickly and then contribute to the chemical enrichment of either galaxies or intergalactic medium or both (e.g., Madau, Ferrara, & Rees 2001 and references therein). Since usual Population II stars could be formed in cold gas clouds polluted chemically by first Population III stars, we cannot examine observationally whether or not such very massive stars were really formed in the primordial gas clouds. Furthermore, according to Nakamura & Umemura (2001a), proto-galactic clouds with a higher density tend to form intermediate-mass stars. If such proto-galactic systems are present, their evolution can be governed by evolution of the intermediate-mass stars whose lifetime is $`10^8`$ yr, being much longer than those of very massive stars. We investigate what happens and follows in such bursts of formation of intermediate-mass stars and discuss effects on the further evolution of galaxies, adopting the IMF of YS86. In particular, intermediate-mass stars experience inevitably the planetary-nebula-nucleus (PNN) phase, providing a lot of high-energy photons (e.g., Vassiliadis & Wood 1994). Therefore, it is expected that the galaxies appear to be luminous emission-line galaxies during the PNN phase. We discuss how they look and compare our model results with observations of some interesting high-redshift galaxies.

## 2 EVOLUTION OF POPULATION-III STARBURSTS IN FORMING GALAXIES

### 2.1 Initial Mass Function for Population III Stars derived by YS86

In many previous galaxy evolution models, it has been often assumed that stars are formed with a Salpeter-like initial mass function (IMF) in a galaxy;

$$\varphi (m)=\beta m^\mu $$
(1)
where $`m`$ is the stellar mass in units of $`M_{}`$ and $`\beta `$ is a normalization constant determined by the relation

$$_{m_l}^{m_u}\varphi (m)𝑑m=1,$$
(2)
which leads to

$$\beta =\frac{(\mu 1)m_l^{\mu 1}}{1(m_l/m_u)^{\mu 1}}.$$
(3)
In this formulation, there are three free parameters; the power index ($`\mu `$), and the upper and lower mass limits of the IMF ($`m_u`$ and $`m_l`$). In order to describe the evolution of the solar neighborhood, one adopts conventionally $`\mu `$ = 1.35, $`m_l=0.1M_{}`$, and $`m_u=60M_{}`$ or $`100M_{}`$. However, it is highly uncertain whether or not this Salpeter-like IMF is applicable to the star formation history in primeval galaxies.

YS86 investigated the IMF for a hydrogen-helium gas cloud in detail based on the opacity-limited fragmentation theory. They found that typical masses of the first-generation stars are from several $`M_{}`$ to 10 $`M_{}`$, depending on the assumed mass-luminosity relation for protostars derived by Silk (1977c); assuming that the opacity sources within the protostar is due to grains, he obtained the protostar luminosity

$$L(m)=ϵ10^{2.47}(m/M_{})^\eta L_{},$$
(4)
where $`ϵ=(Z/Z_{})^1`$ and $`\eta =3`$. Since this equation cannot be applicable to the luminosity of stars with $`Z10^4Z_{}`$, YS86 treated $`ϵ`$ as well as $`\eta `$ as free parameters in the ranges $`10ϵ100`$ and $`1.5\eta 3`$. Then YS86 found that the mass at the peak of the IMF is $`4M_{}`$ for $`\eta =3`$ and $`10M_{}`$ for $`\eta =1.5`$. Another important conclusion in YS86 is that stars more massive than these peaks are formed less efficiently than those expected from the Salpeter IMF with $`\alpha =1.35`$. In addition, stars less massive than these peaks are also formed less efficiently than those expected from the same Salpeter IMF (see Figure 3 in YS86). Therefore, the IMFs derived by YS86 are significantly different from the Salpeter IMFs. We therefore adopt the IMFs of YS86 and investigate luminosity and spectral evolution of galaxies for various sets of the parameters.

For convenience, we present a parametric form of the YS86 IMF in the mass range between 1 $`M_{}`$ and 60 $`M_{}`$ for $`\eta =3`$ (YS86a),

$$\varphi (m)=\{\begin{array}{ccc}1.5944\times 10^2m^{2.0}\hfill & & 1m/M_{}<2.5\hfill \\ 4.1728\times 10^2m^{0.95}\hfill & & 2.5m/M_{}<3.7\hfill \\ 1.4462\times 10^1m^0\hfill & & 3.7m/M_{}<5\hfill \\ 1.0813m^{1.25}\hfill & & 5m/M_{}<7.5\hfill \\ 2.0081\times 10m^{2.7}\hfill & & 7.5m/M_{}<11\hfill \\ 2.2088\times 10^2m^{3.7}\hfill & & 11m/M_{}60\hfill \end{array}$$
(5)
and that for $`\eta =1.5`$ (YS86b),

$$\varphi (m)=\{\begin{array}{ccc}7.3036\times 10^4m^{3.5}\hfill & & 1m/M_{}<2.9\hfill \\ 4.4629\times 10^3m^{1.8}\hfill & & 2.9m/M_{}<5\hfill \\ 8.0865\times 10^2m^0\hfill & & 5m/M_{}<8\hfill \\ 1.2849m^{1.33}\hfill & & 8m/M_{}<12\hfill \\ 1.4671\times 10m^{2.31}\hfill & & 12m/M_{}60\hfill \end{array}$$
(6)
Here the above numerical coefficients are determined by the following normalization for each case,

$$_1^{60}\varphi (m)𝑑m=1.$$
(7)
In Figure 1, we show the shapes of IMF which we use in this paper. Since the massive end slope of YS86a IMF is steepest, the number of massive stars ($`>20M_{}`$) is smallest among the three, although the number of intermediate-mass star exceeds that of Salpeter IMF.

### 2.2 How Populous are Intermediate-Mass Stars ?

Prior going to investigating the luminosity and spectral evolution of galaxies using the YS86 IMFs given in equations (5) and (6), we estimate how populous intermediate-mass stars are in the YS86 IMFs with respect to those expected from usual Salpeter-like IMFs. Since it is expected that the most spectacular effect of intermediate-mass enhanced star formation is brought by a number of luminous PNNs with very high temperature $`10^6`$ K (e.g., Vassiliadis & Wood 1994). Since the mass of such high-temperature PNN progenitors lies in a range between $`4M_{}`$ and $`5M_{}`$, we estimate how many stars within this mass range are formed.

If we adopt a Salpeter IMF, the number of stars with a mass range $`m_1m_{}m_2`$ formed in the gas with 1 $`M_{}`$ is estimated as

$$N(m_1m_{}m_2)=_{m_1}^{m_2}\frac{\varphi (m)}{m}𝑑m.$$
(8)
Using $`\beta `$ given in equation (3), we re-write equation (8) as

$$N(m_1m_{}m_2)=\left(\frac{\beta }{\mu }\right)(m_1^\mu m_2^\mu )\mathrm{stars}M_{}^1.$$
(9)
Taking $`m_1=4M_{}`$ and $`m_2=5M_{}`$, we obtain the numbers for the following two cases; 1) $`\mu `$ = 1.35, $`m_l=0.1M_{}`$, and $`m_u=60M_{}`$, and 2) $`\mu `$ = 1.35, $`m_l=1M_{}`$, and $`m_u=60M_{}`$ (see Table 1).

If we adopt a YS86 IMF, we can estimate the number of stars with a mass range $`4M_{}m_{}5M_{}`$ by integrating equations (5) and (6) for $`\eta =1.5`$ and $`\eta =3`$, respectively. The results are shown in Table 1. We also show the number if all stars have a mass of 5 $`M_{}`$; i.e., the IMF is described as a $`\delta `$ function with $`m_{}=5M_{}`$.

As shown in Table 1, the number of intermediate-mass stars with $`m_{}=4M_{}`$$`5M_{}`$ for the YS86 IMF with $`\eta =1.5`$ is nearly the same as that for the Salpeter IMF with $`m_l=1M_{}`$. However, that for the YS86 IMF with $`\eta =3`$ is larger by a factor of 2.4 than that for the Salpeter IMF with $`m_l=1M_{}`$. For both cases, the ratio of the number of intermediate-mass stars with $`m_{}=4M_{}`$$`5M_{}`$ for YS86 IMFs to that for Salpeter IMF is within a factor of 3. If stars with mass of $`5M_{}`$ are only formed, the number is larger by a factor of 15 than that for the Salpeter IMF with $`m_l=1M_{}`$ (see section 3).

### 2.3 Galaxy Evolution Models

In this section, we show the evolution of bolometric luminosity, ionizing photon production rate, and emission line ratios of Population III star clusters with the YS86 IMFs. We adopt the $`10^7`$ yr burst model as the star formation history. We calculate the spectral energy distributions (SEDs) of star clusters by using the evolutionary population synthesis code PEGASE2 (Fioc & Rocca-Volmerange 1997, 2000). Although the code can be used to calculate the photometric evolution of galaxies being the effect of the chemical evolution into account, we calculate that of a single metallicity population with $`Z=0.0001`$; note that this metallicity is the poorest metallicity available in PEGASE2. The duration of the starburst is longer than the lifetime of the most massive star in our model ($`60M_{}`$), $`5\times 10^6`$ yr (Girardi et al. 1996). Since the lifetime is a half of the starburst duration, the metallicity of a half of stars formed in the starburst is $`Z<10^{10}`$.

Figure 2 shows the SED evolution for each model: Salpeter IMF (solid lines), YS86a IMF (dotted lines) and YS86b IMF (dashed lines). We show the SEDs at 0.01 Gyr, 0.1 Gyr and 0.5 Gyr. Note that the ordinate is the flux in units of erg s-1 $`M_{}^1`$, corresponding to the case in which the total gas mass is $`1M_{}`$. The SEDs are mainly determined by the upper main sequence stars and are not sensitive to the IMF. During the starburst, the ionizing photons are radiated from massive main sequence stars whose effective temperature is lower than about 60,000 K (Girardi et al. 1996). Therefore there is a large break at 228 Å; the photons with wavelength shorter than this are capable of ionizing the singly ionized helium. At an age of $`>10^8`$ yr, hottest stars in the star cluster are not main sequence stars but post-asymptotic giant branch (AGB) stars whose progenitors are stars with a degenerate carbon core, e.g., intermediate-mass stars. The effective temperature of such post-AGB stars (i.e., PNNs) is about 200,000 K (Blöcker 1995). Therefore the break of 228 Å of such stars becomes to be small and thus the SED shape of far-UV is very hard.

The global shape of SED depends slightly on the IMF (Figure 3). During a period of the starburst, between $`4\times 10^6`$ yr and $`1\times 10^7`$ yr, colors ($`UB`$, $`BV`$, and $`VK`$) become redder because of red supergiants. In this phase, the star cluster with the Salpeter IMF is reddest and that with the YS86a is bluest. This is because the number ratio of red supergiants to upper main sequence stars is highest for the Salpeter IMF ($`m^{1.35}`$) and lowest for the YS86a IMF ($`m^{3.7}`$). On the other hand, for ages longer than $`10^8`$ yr, the star cluster with the Salpeter IMF is bluest and that with the YS86b is reddest. The reason for this is that the number ratio of redgiants to upper main sequence stars is lowest for the Salpeter IMF ($`m^{1.35}`$) and highest for the YS86b IMF ($`m^{3.5}`$).

Figure 4 shows the evolution of bolometric luminosity of a star cluster of $`1M_{}`$ for various IMFs. The bolometric luminosity reaches its maximum at the end of starburst ($`10^7`$ yr). The star cluster with the YS86b IMF is brightest while that with the YS86a IMF is faintest among the models. After that, following the death of massive stars, the bolometric luminosity of stellar systems decreases with time. The decreasing rate of the bolometric luminosity also depends on the functional form of IMF. Since the YS86a IMF has a peak at $`m4M_{}`$ whose lifetime is about 0.15 Gyr, the decreasing rate of the bolometric luminosity is small till the time, and becomes large after that. Since the lifetime of stars with $`1M_{}`$ is about 5.6 Gyr, all the stars formed during the starburst have died till the time. As we show below, the dominant sources of the ionizing photons become the PNN at the age of $``$ 0.1 Gyr. The bolometric luminosity to mass ratios decrease from $`66L_{}/M_{}`$ to $`3.2L_{}/M_{}`$ during this period.

In the upper panel of Figure 5, we show the evolution of production rate of photons capable of ionizing neutral hydrogen, $`N(\mathrm{H}^0)`$, which is calculated by integrating the spectra below 912 Å, and that of photons capable of ionizing ionized helium, $`N(\mathrm{He}^+)`$, which is calculated by integrating the spectra below 228 Å. First, we mention how the production rate of the hydrogen ionizing photons evolves. As the bolometric luminosity increases, the ionizing photon production rate reaches its maximum at the end of starburst, because the lifetime of early B stars of $`10M_{}`$ is longer than $`10^7`$ yr. The evolution of $`N(\mathrm{H}^0)`$ depends on IMF. During the starburst, the ionizing photon production rate per unit mass is largest for the star cluster with the YS86b IMF but smallest for that with the YS86a IMF. In our calculations, we adopt the mass of most massive stars $`M_u=60M_{}`$ whose lifetime is about $`5\times 10^6`$ yr. When the age of star clusters reaches the lifetime of the most massive stars, the number of hydrogen ionizing photons begins to decrease with time, following the death of stars radiating hydrogen ionizing photons. At the age of $`10^8`$ yr, the intermediate-mass stars evolve to PNNs and they become dominant ionization sources for hydrogen. In the star cluster with the Salpeter IMF, the number of PNNs hardly changes and thus the production rate of ionizing photons is kept nearly constant. On the other hand, for the star cluster whose IMF is either YS86a or YS86b, the number of PNNs and the production rate of ionizing photons decrease with decreasing number of low mass stars. After the age of $``$ 1 Gyr, the ionizing photon production rate is highest for the star cluster with the Salpeter IMF and smallest for that with the YS86b IMF since it is proportional to the number of PNNs. Next, we comment on the evolution of the production rate of ionized helium ionizing photons. The lifetime of stars with $`m>10M_{}`$ is shorter than a few $`\times 10^7`$ yr. After the massive stars have died, the remaining main sequence stars cannot ionize the neutral helium. When the age of the star cluster become $`1\times 10^8`$ yr, PNNs come into existence. The effective temperature of PNNs is very high enough to radiate He+ ionizing photons.

The lower panel of Figure 5 shows the evolution of the ratio of $`N(\mathrm{He}^+)`$ to $`N(\mathrm{H}^0)`$ which indicates the hardness of UV spectra. During the starburst, the ratio is lower than $`10^3`$. However, after PNNs become the dominant source of UV radiation, the ratio exceeds $`10^3`$. The stellar mass of a progenitor of PNN decreases with time. Therefore, the effective temperature of PNN decreases and thus the ratio of $`N(\mathrm{He}^+)`$ to $`N(\mathrm{H}^0)`$ also decreases (Vassiliadis & Wood 1994).

We investigate ultraviolet emission-line properties of the gaseous nebula photoionized by PNN clusters during the course of galaxy evolution. We use the photoionization code CLOUDY94 (Ferland 1997), which solves the equations of statistical and thermal equilibrium and produces a self-consistent model of the run of temperature as a function of depth into the nebula. Here we assume that a uniform density, dust-free gas cloud with plane-parallel geometry is ionized by various continuum sources. The parameters for the calculations are (1) the hydrogen density of the cloud, $`n_\mathrm{H}`$, (2) the ionization parameter, $`U=Q(\mathrm{H}^0)(4\pi r^2n_\mathrm{H}c)^1`$, where $`Q(\mathrm{H}^0)`$ is the number of ionizing photons, $`r`$ is the distance from the ionizing source, and $`c`$ is the speed of light (Osterbrock 1989), (3) the SED of the ionizing radiation, and (4) the chemical composition. The SED of the ionizing radiation is derived as noted before using PEGASE2 (see Figure 2). We perform our photoionization calculations as a function of the ionization parameter between $`\mathrm{log}U=4`$ and $`1`$ with a logarithmic interval of 0.5. We adopt $`n_\mathrm{H}=10^2\mathrm{cm}^3`$. We also adopt a metallicity of $`Z=0.001=0.05Z_{}`$. The calculations were stopped when the temperature fell to 3000 K, below which little optical emission is expected.

For reference, we also calculate AGN photoionization models and shock-heating models in the following ways. a) AGN photoionization models are calculated by using CLOUDY94. We adopt that the input SED has a power-law form with $`\alpha =1.5`$ ($`f_\nu \nu ^\alpha `$). b) Shock-heating models are taken from results of Dopita & Sutherland (1996). Their models assume solar metallicities, shock velocities from $`150\mathrm{km}\mathrm{s}^1`$ to $`500\mathrm{km}\mathrm{s}^1`$, and the magnetic parameter $`0B/\sqrt{n}4\mu \mathrm{G}\mathrm{cm}^{3/2}`$. The PNN photoionization model results are shown together with these models in Figure 6; the UV emission-line diagnostic diagrams proposed by De Breuck et al. (2000).

Since high-energy photons supplied from the PNN clusters can cause high-ionization lines such as C iv, the UV emission-line properties of the PNN photoionization models are quite similar to those of the AGN photoionization. In particular, the model locus predicted for $`10^7`$ yr-burst stellar populations with age of 0.5 Gyr cannot be distinguished from that for the AGN models with a power-law continuum with $`\alpha =1.5`$. This similarity has been already utilized to show that some LINERs (low-ionization nuclear emission-line regions: Heckman 1980) can be explained by photoionization by PNN clusters (Taniguchi, Shioya, & Murayama 2000). On the other hand, the shock-heating models appear to show different properties from those of the PNN and AGN models.

Figure 7 shows the $`L_{\mathrm{HeII}}/L_{\mathrm{Bol}}`$, Ly$`\alpha `$/He ii, and N v/He ii ratios as a function of $`\mathrm{log}U`$ for PNN models with age of 0.1 Gyr (open circles) and 0.5 Gyr (filled circles). These ratios hardly depend on the IMF but depend on the age.

Recently, the evolutionary tracks of zero-metal stars have been published (e.g., Marigo et al. 2001). Based on these new models, Schaerer (2002) demonstrated that massive population III stars can radiate the He+ ionizing photons because of both the high effective temperature ($`T_{\mathrm{eff}}>80,000`$ K) of zero-metal stars and the effect of the non-LTE stellar atmosphere (see also Tumlinson & Shull 2000). The small $`N(\mathrm{He}^+)/N(\mathrm{H}^0)`$ ratio during the starburst shown in Fig.5 is attributed to the use of the population synthesis model with metallicity of $`Z=0.0001`$ instead that with $`Z=0`$. We note, however, H ii regions around massive population III stars cannot radiate metal lines in which we are interested since the gas around massive Population III stars is still primordial.

### 2.4 Comparison with high-z star forming galaxies

If the star formation in some forming galaxies takes place with YS86 IMF, such galaxies could experience a PNN-dominated phase $`10^8`$ yr after the initial starburst. Therefore it seems important to look for such galaxies at high redshift. As is shown in section 2.3, the SED of star clusters in PNN-dominated phase is so hard that the emission-line properties of the gaseous nebula photoionized by PNN clusters are indistinguishable from those photoionized by AGN. Considering the stellar evolution of intermediate mass stars, the PNN phase follows the post-AGB phase when stars eject a dusty gaseous envelope. It is therefore expected that the PNN cluster is surrounded by dusty gas (see the next section). Until the PNN-dominated phase begin, dust formed in the ejecta of Type II supernovae is also present in the interstellar medium (Todini & Ferrara 2001). Therefore, a part of ionizing photons could be absorbed by such dust grains. However, we neglect this effect in the following discussion for simplicity.

The above expected features of the PNN star cluster are very interesting because the PNN star cluster may be observed as a dusty, ultraluminous infrared galaxy with high-ionization emission lines at high redshift. It is known that some high-$`z`$ forming galaxy candidates share nearly the same properties; e.g., IRAS F 10214$`+`$4724 at $`z=2.3`$ (Rowan-Robinson et al. 1991) and SMM 02399$``$0136 at $`z=2.8`$ (Ivison et al. 1998). Their observed huge luminosities ($`L_{\mathrm{bol}}10^{12}L_{}`$) have been often considered to be powered by massive stars. However, high-ionization emission lines such as C iv $`\lambda `$1549 Å and N v $`\lambda `$1240 Å are also observed in them. Since these emission lines cannot be produced by photoionization by ordinary massive stars, it has been often considered that active galactic nuclei (AGNs) also reside in them (Ivison et al. 1998; Serjeant et al. 1998). In particular, the presence of hidden broad emission lines reinforces this interpretation for IRAS F10214$`+`$4724 (Goodrich et al. 1996). Spectropolarimetry of SMM J02399$``$0136 reveals that the polarization behavior is similar to BAL quasars (Vernet & Cimatti 2001). Even in the local universe, bursts of star formation (starbursts) are sometimes associated with galaxies hosting AGNs and thus the simultaneous presence of both starburst and AGNs may not be unusual (e.g., Taniguchi 1999; Mouri & Taniguchi 2001; Storchi-Bergmann et al. 2001). However, the evidence for AGNs in them has been obtained often by rest-frame ultraviolet and optical emission lines which are subject to misunderstanding in some cases (Taniguchi et al. 1999; Lutz, Veilluex, & Genzel 1999). Since the effective temperatures of PNNs exceed $`10^5`$ K, high energy photons are much more numerous than those from ordinary massive stars. Here we examine whether or not the PNN cluster is responsible for the formation of high-ionization emission lines.

In Figure 6, we plot the observed positions of IRAS F10214+4714 (Serjeant et al. 1998) and the high-$`z`$ powerful radio galaxies studied by De Breuck et al. (2000). Although the high-$`z`$ powerful radio galaxies have an AGN, at least some of them are considered to be very young galaxies (e.g., Eales et al. 1993). We also plot available data of planetary nebulae (Feibelman 2000). It is found that the emission line ratios of IRAS F10214+4724 as well as the high-$`z`$ powerful radio galaxies in these UV diagnostics can be explained either by the AGN models with $`2\mathrm{log}U1.5`$ or by the PNN model with an age of 0.5 Gyr. The similar behaviors between the PNN and AGN photoionization models unable us to distinguish the two models unfortunately.

It is noted that both the AGN and PNN models appear to underpredict the C ii\]/C iii ratio as shown in the right-lower panel in Figure 6. This discrepancy suggests that the observed UV emission-line ratios of IRAS F10214+4714 and the high-$`z`$ powerful radio galaxies can be explained either by a combination between the shock heating and the AGN photoionization or by a combination between the shock heating and the PNN photoionization.

Next, we study whether or not the star clusters dominated by PNNs can explain the luminosity of high-$`z`$ forming galaxy candidates. As shown in section 2.2 the bolometric luminosity to mass ratio changes from $`66L_{}/M_{}`$ to $`3.2L_{}/M_{}`$. Therefore, the mass of stellar system needed to reproduce the bolometric luminosity of IRAS F10214+4724 ($`10^{13}L_{}`$), is $`10^{11}10^{12}M_{}`$. This value is comparable to the stellar mass of massive galaxies. It is impossible to radiate the huge luminosity from only PNNs of first generation stars for this galaxy. There are two reasons why it is hard to explain the huge luminosity by a lot of PNNs. As shown in section 2, even for the YS86’s IMF the number of intermediate-mass stars (e.g., $`4M_{}<M_{}<5M_{}`$) is only a few time larger than that for the Salpeter IMF. The other reason is related to the difference of the lifetime of progenitor stars. The lifetime of stars with mass of $`4M_{}`$ is $`1.54\times 10^8`$ yr and that of $`5M_{}`$ is $`9.9\times 10^7`$ yr (Girardi et al. 1996). Therefore, the mass range of stars which evolve to PNNs at the same time is very small and as a result the number of those stars is very small. To solve this discrepancy, we consider an extreme model in the next section.

## 3 AN EXTREME MODEL

### 3.1 Scenario

Here, we consider an extreme case that the formation of intermediate-mass stars occurs during a period of $`\tau _{\mathrm{SF}}`$. Most luminous forming galaxy candidates at high redshift contain a lot of molecular gas (e.g., $`M_{\mathrm{gas}}10^{11}M_{}`$) and their inferred star formation rates (SFR) are of the order of $`10^3M_{}`$ yr-1 (Frayer et al. 1998, 1999). Therefore, we assume that a galaxy has a gas mass of $`10^{11}M_{}`$ and the initial star formation occurs at a rate of $`SFR\dot{M}_{\mathrm{SF}}=10^3M_{}`$ yr-1. We adopt a duration of the star formation $`\tau _{\mathrm{SF}}=10^7`$ yr. The gas mass converted to stars is estimated as $`M_{\mathrm{SF}}=\dot{M}_{\mathrm{SF}}\times \tau _{\mathrm{SF}}=10^{10}M_{}`$. Thus in this model, the star formation efficiency is $`\eta _{\mathrm{SF}}=`$ 10%. The star formation may take place in a central region (say, within a radius of $``$ 1 kpc) of galaxies. We assume for simplicity that all stars formed in this initial star formation have a mass of 5 $`M_{}`$ (i.e., a typical intermediate-mass star). Therefore, the SFR adopted above means that 200 stars with $`m_{}=5M_{}`$ are formed in a year; $`\dot{N}_{}=\dot{M}_{\mathrm{SF}}/m_{}=200`$ stars yr-1.

Here we consider the evolution of a cluster of such stars. The bolometric luminosity of each star in the main sequence phase is $`L_{}1\times 10^3L_{}`$ and the age of the star is $`\tau _{\mathrm{age}}1\times 10^8`$ years. The total number of stars formed in the duration is $`N_{}`$(total) = $`\dot{N}_{}\times \tau _{\mathrm{SF}}=2\times 10^9`$ stars. Therefore, the total luminosity of the star cluster increases during the star formation period ($`10^7`$ years) at a rate of $`\dot{N}_{}\times L_{}2\times 10^5L_{}`$ yr-1 and reach the maximum bolometric luminosity $`L_{\mathrm{bol}}(\mathrm{max})2\times 10^{12}L_{}`$ at $`t=10^7`$ years. This maximum luminosity keeps as long as the stars are alive; i.e., up to $`t=\tau _{\mathrm{age}}=10^8`$ years. Shortly after this age, the stars evolve from the main sequence to PNNs through AGB.

Since one PNN has a bolometric luminosity of $`L_{\mathrm{PNN}}10^4L_{}`$ (e.g., Vassiliadis & Wood 1994), the total PNN luminosity amounts to $`L_{\mathrm{PNN}}(\mathrm{total})=\dot{N}_{}\tau _{\mathrm{PNN}}L_{\mathrm{PNN}}=2\times 10^{11}L_{}`$ where $`\tau _{\mathrm{PNN}}10^5`$ yr is the lifetime of the PNN phase. Nevertheless this short lifetime, note that the PNN phase in the star cluster continues during nearly the same period as that of star formation; i.e., $`10^7`$ yr. In the PNN phase, the total luminosity of the star cluster is mostly dominated by unevolved main sequence stars while the ionizing photons are supplied mainly by PNNs.

Another important property of the PNN phase is that a lot of gas is ejected from the stars. The core mass of each star is estimated to be $`m_\mathrm{c}0.6M_{}`$. Therefore, each star ejects a metal-enriched gas of $`m_{\mathrm{gas},\mathrm{PN}}=m_{}m_\mathrm{c}5M_{}0.6M_{}4.4M_{}`$. The total ejected gas mass amounts to $`M_{\mathrm{gas},\mathrm{PN}}=m_{\mathrm{gas},\mathrm{PN}}N_{}(\mathrm{total})8.8\times 10^9M_{}`$. Since the gas-to-dust mass ratio in such gas may be $``$ 100, the star cluster can produce a lot of dust grains with a mass of $`M_{\mathrm{dust}}10^8M_{}`$.

In summary, the PNN phase of the star cluster has the following properties; 1) the bolometric luminosity is $`L_{\mathrm{bol}}10^{12}L_{}`$; 2) the ionizing luminosity is $`2\times 10^{11}L_{}`$; 3) the ionizing continuum is so hard that it creates highly-ionized regions around the star cluster; 4) the dust mass is $`M_{\mathrm{dust}}10^8M_{}`$; 5) the heavy element abundance is increased up to $`ZM_{\mathrm{dust}}/M_{\mathrm{gas}}0.001`$, and 6) this phase lasts for $`10^7`$ yr given that the duration of star formation is $`10^7`$ yr.

### 3.2 Emission-line properties of the PNN cluster

As shown in section 2.3, the PNN cluster is capable of the formation of the high-ionization lines such as C iv $`\lambda `$1549 and C iii\] $`\lambda `$1909 emission lines. For example, the observed C iv/He ii ratios are 5.7 and 1.9 for SMM 02399$``$0136 and IRAS F10214$`+`$4724, respectively (Ivison et al. 1998; Serjeant et al. 1998). However, our models cannot explain the observed higher ratios of N v$`\lambda `$1240/He ii; 8.5 for SMM 02399$``$0136 (Ivison et al. 1998) and 1.4 for IRAS F10214$`+`$4724 (Serjeant et al. 1998) (see Figure 7). In our model, we have assumed that each heavy element abundance is five percent of its solar value. If only the helium burning occurs in the core of PNN, carbon could be enriched. However, if the hydrogen burning also occurs, nitrogen could be enriched because nitrogen is formed as the end product of CNO cycle in the core. If such nitrogen enrichment occurs, the observed strong N v emission may be explained. It is worthwhile noting that the overabundance of nitrogen is also suggested for the broad emission-line gas in high-$`z`$ quasars (Hamann & Ferland 1993, 1999).

### 3.3 Further Evolution of the Star Cluster

We consider the final fate of the star cluster in our model given in section 3.1. The intermediate-mass stars formed in the initial starburst will evolve to white dwarf (WD) stars with a mass of $`m_{\mathrm{WD}}=m_\mathrm{c}0.6M_{}`$ through the phases of AGB and planetary nebula. Since the evolved WD stars have a luminosity of $`10^4L_{}`$ at their ages of $`10^{10}`$ yr, the star cluster formed in the starburst becomes finally very faint; i.e, $`L_{\mathrm{WD}}(\mathrm{total})=L_{\mathrm{WD}}N_{}(\mathrm{total})2\times 10^5L_{}`$. Since there is no such candidate in the local universe, it is expected that further star formation could have occurred in galaxies even if starbursts with only intermediate-mass stars occurred in the galaxies. After the death of PNN stars, dust grains with a mass of $`10^8M_{}`$ are left in a galaxy. Since the galaxy has the gas with a mass of $`10^{11}M_{}`$, the heavy element abundance increases at a level of $`Z0.001`$ after the PNN phase. Therefore, so-called Population II stars would form in the gas clouds polluted chemically by the intermediate-mass stars. In other word, galaxies with a PNN cluster may evolve to ordinary galaxies observed in the present-day universe.

### 3.4 Final Remarks

What we would like to propose in this section is not that some particular high-$`z`$ dusty galaxies are indeed the PNN galaxies but that candidates of the extreme PNN galaxies share nearly the same properties of the observed high-$`z`$ dusty galaxies. Recent deep surveys both at FIR and submillimeter wavelengths have discovered a significant number of hyperluminous dusty galaxies at high redshift (Rowan-Robinson et al. 1991; Kawara et al. 1998; Puget et al. 1999; Smail, Ivison, & Blain 1997; Hughes et al. 1998; Barger et al. 1998; Eales et al. 1999; Rowan-Robinson 2000). If some proto-galactic systems tend to form intermediate-mass stars preferentially as Population III stars, it is possible that candidates of the PNN galaxies may be found in some high-$`z`$ dusty galaxies. If present, such galaxies will be detected in future sensitive search for dusty galaxies made by Atacama Large Millimeter Array (ALMA; Kurz & Shaver 1999). Deep optical spectroscopy using new 30-m class telescopes will be also necessary to investigate detailed rest-frame UV spectra of such high-$`z`$ dusty galaxies.

We would like to thank Gary Ferland for making his code CLOUDY available to the public. We would like to thank Hideyuki Saio and Satoru Ikeuchi for useful discussion. We would also like to thank an anonymous referee for his/her useful comments and suggestions. This work was financially supported in part by Grant-in-Aids for the Scientific Research (Nos. 10044052, and 10304013) of the Japanese Ministry of Education, Culture, Sports, and Science. YS is JSPS fellow.


text:
# Entanglement between two macroscopic fields by coherent atom-mediated exchange of photons

## Abstract

Using two different criteria for continuous variable systems we demonstrated that pump and probe beams became quantum correlated in a situation of Electromagnetically Induced Transparency in a sample of 85Rb atoms. Our result combines two important features for practical implementations in the field of quantum information processing. Namely, we proved the existence of entanglement between two macroscopic light beams, and this entanglement is intrinsically associated to a strong coherence in an atomic medium.

In the last two decades, the philosofical and technological impacts of quantum correlations or entanglement in multipartite physical systems have been accentuated by the theoretical and experimental investigations leading to the development of the area of quantum information processing bennett . The existence of entanglement, which is an evidence of the non-local character of the quantum theory, has been confirmed in quantum optics experiments using dichotomics aspect  and continuous variable systems kimble1 . These experimental realizations made possible some implementations in the fields of quantum information kimble ; gisin  and computation haroche1 ; ions . Even so, the unavoidable interaction between the physical system used to process the quantum information and the environment leads to a loss of coherence, and consequently, to a loss of information that introduces important limitations to the practicality of the quantum information processing technology. Since it is impossible to have isolated systems, among the different approaches employed to reduce the influence of the environment we find the use of continuous variable systems, specifically, intense light beams as in the case of quantum teleportation based on the entanglement between twin beams issued from an OPO kimble . In that sense, there are propositions, for exemple, for quantum cryptography ralph , quantum computation lloyd , dense coding samuel , and quantum key distribution leuchs . Another approach is to use coherently-prepared atomic media and, in this case, the Electromagnetically Induced Transparency (EIT) is a good candidate as it has been demonstrated recently by the observation of very slow light pulse propagation paq  and light storage p . In this context, a natural question arises: can we produce entanglement between two intense ligth beams in a coherently-prepared atomic medium, combining in that way the two mentioned approaches? As we will see through out this paper, the answer to this question is yes and this result is very important due to the recent huge interest in the applications of coherently-prepared atomic media using EIT mastko .

In this paper, we show that three-level atoms can produce entanglement in two intense travelling light fields, in the EIT regime. We study this system from a theoretical point of view, and we show that for some particular parameters, the two fields present quantum correlations after interacting with the atoms, even if they are initially completely uncorrelated. We demonstrate the existence of entanglement between these two propagating fields according to two different criteria.

In our model, we consider three-level atoms in a closed $`\mathrm{\Lambda }`$ configuration (ground states $`|1`$ and $`|2`$, and excited state $`|0`$) interacting with two copropagating fields treated quantum-mechanically. In the Heisenberg picture, the electric field operator for the propagating mode $`j`$ (pumping laser $`j=1`$, probe laser $`j=2`$) is given by the expression

$$\widehat{\stackrel{}{E}}_j(t)=_{0\omega _{Lj}}\overline{ϵ}_je^{i\omega _{Lj}t}\widehat{A}_j(t)+h.c.,$$
(1)
where $`_{0\omega _{Lj}}`$ , $`\overline{ϵ}_j`$ and $`\omega _{Lj}`$ are the amplitude, the polarization direction and the angular frequency of mode (laser) $`j`$, respectively. $`\widehat{A}_j`$ ($`\widehat{A}_j^{}`$) is the annihilation (creation) operator and represents the slowly varying amplitude of the laser field. We take the hamiltonian

$$\widehat{H}_{Lj}=_{\mathrm{}}^+\mathrm{}𝑑\omega |G_j(\omega )|^2\mathrm{}\omega \widehat{a}_{j\omega }^{}\widehat{a}_{j\omega }$$
(2)
as the energy source of the interacting field where, $`\widehat{a}_{j\omega }`$ is the annihilation operator of the field inside the laser source cavity and, through the non-dimensional function $`G_j(\omega )`$, we take into account the influence of the external vacuum modes. This function is determined by the frequency-dependent reflectivity of the output mirror of the laser cavity and provides the laser linewidth $`\gamma _j`$, which we assume to be constant here, in accordance with the Markov approximation gardiner . Taking into account the theoretical and experimental studies about the laser sources, we took $`G_j(\omega )`$ as a lorentzian profile centered at the laser frequency $`\omega _{Lj}`$, allowing the lower integration limit in (2) be taken equal to $`\mathrm{}`$ instead of zero.

The coupling between the source and the propagating mode is given by the linear hamiltonian

$$\widehat{H}_{LjCj}=i\mathrm{}\sqrt{\frac{\gamma _j}{2\pi }}_{\mathrm{}}^+\mathrm{}𝑑\omega \widehat{a}_{j\omega }^{}\widehat{A}_j(t)+h.c.$$
(3)
We obtain the interaction hamiltonian for the two light beams and the atoms, within the usual dipole and rotating-wave approximations

$`\widehat{H}_{int}=\mathrm{}g_1\widehat{S}_1^+(t)\widehat{A}_1(t)+\mathrm{}g_2\widehat{S}_2^+(t)\widehat{A}_2(t)+h.c.,`$ (4)
where $`g_1`$ ($`g_2`$) is the atom – field 1 (field 2) coupling strength, and $`\widehat{S}_1^+`$ ($`\widehat{S}_2^+`$) the slowly varying envelope of the atomic polarization on the transition $`|1|0`$ ($`|2|0`$).

The dynamics of the system is determined by twelve coupled quantum Langevin equations derived from the Heisenberg equations of motion. Since we are dealing with macroscopic systems, the quantum fluctuations of the operators are studied by linearizing them around their steady-state values and the dynamics of these fluctuations is described by a matrix linear stochastic differential equation for the fluctuation operators tobe . Recently, intensity correlations between the pump and probe fields in EIT have been measured exper . These can be understood by inspection of the equations for the fluctuations of one field and for the corresponding atomic polarization:

$$\frac{d\delta \widehat{A}_1(t)}{dt}=\frac{\gamma _1}{2}\delta \widehat{A}_1(t)ig_1\delta \widehat{S}_1^{}(t)+\sqrt{\gamma _1}\delta \widehat{A}_{1in}(t),$$
(5)
$`{\displaystyle \frac{d\delta \widehat{S}_1^{}(t)}{dt}}=\left({\displaystyle \frac{\mathrm{\Gamma }_1+\mathrm{\Gamma }_2}{2}}i\delta _{L1}\right)\delta \widehat{S}_1^{}(t)`$
$`+ig_1w_1\delta \widehat{A}_1(t)+ig_1\alpha _1\delta \widehat{W}_1(t)`$
$`ig_2s_{12}^{}\delta \widehat{A}_2(t)ig_2\alpha _2\delta \widehat{S}_{12}^{}(t)+\widehat{F}_{S1}(t).`$ (6)
Here we define $`\widehat{A}_{1in}`$ the annihilation operator of the source or input field 1, $`\mathrm{\Gamma }_1`$ ($`\mathrm{\Gamma }_2`$) the spontaneous emission rate from $`|0|1`$ ($`|0|2`$), $`\delta _{L1}`$ the detuning between field 1 and the corresponding atomic transition, $`w_1`$ the steady-state inversion between states $`|0`$ and $`|1`$, $`\alpha _1`$ ($`\alpha _2`$) the steady-state amplitude of field 1 (field 2), $`\widehat{W}_1`$ the inversion (operator) between states $`|0`$ and $`|1`$, $`s_{12}^{}`$ the steady-state coherence between ground states $`|1`$ and $`|2`$, $`\widehat{S}_{12}^+`$ the coherence operator, $`\widehat{F}_{S1}`$ the Langevin fluctuation force. The notation $`\delta \widehat{S}_1^{}`$ means fluctuations of the corresponding operator. The fluctuations of the input $`\delta \widehat{A}_{1in}`$, the interacting $`\delta \widehat{A}_1`$ and the detected $`\delta \widehat{A}_{1out}`$ fields are related by the expression $`\delta \widehat{A}_{1out}=\delta \widehat{A}_{1in}\sqrt{\gamma _1}\delta \widehat{A}_1`$.

From Eq. (6), we notice that noise correlations between the fields are created, owing to the coherent effect in the atomic medium exper . In Fig. 1 we show the quadrature correlations of the fields in the frequency domain as a function of probe detuning for a resonant pump. This theoretical prediction (and the following too) corresponds to a system of N=$`10^8`$ atoms of 85Rb where the states of the $`\mathrm{\Lambda }`$ configuration are designated as follows: $`|0=|5P_{3/2},F^{^{}}=3`$, $`|1=|5S_{1/2},F=3`$ and $`|2=|5S_{1/2},F=2`$. The pump and probe lasers are taken linearly polarized with equal intensities (2.8 mW/cm2) and issued from two independent sources with quantum fluctuations corresponding to a coherent state. We took the analysis frequency $`\mathrm{\Omega }=\mathrm{\Gamma }/6`$, where $`\mathrm{\Gamma }=\mathrm{\Gamma }_1+\mathrm{\Gamma }_2=2\pi 6`$ MHz is the total decay rate of the rubidium excited state. The correlation, taking values in the interval \[-1;1\], is defined as the ratio between the fields covariance and the squared root of the product of the fields’ variances. Outside the EIT window the fields are completely uncorrelated and, in the EIT condition (zero probe detuning), there is the following correlation between the fluctuations of the fields quadratures: $`\delta \widehat{Y}_{1out,0}\delta \widehat{Y}_{2out,0}`$ and $`\delta \widehat{Y}_{1out,\pi /2}\delta \widehat{Y}_{2out,\pi /2}`$. The subscript 0 ($`\pi /2`$) stands for the field amplitude (phase) quadrature and the general quadrature fluctuation operator is, for the field 2, given by $`\delta \widehat{Y}_{2out,\varphi }(t)=\delta \widehat{A}_{2out}(t)e^{i\varphi }+\delta \widehat{A}_{2out}^{}(t)e^{i\varphi }`$.

The observed correlations can be interpreted from the propagation dynamics of the beams. When the beams have the same intensity, the role “pump” and “probe” is interchangeable and for a resonant coupling of both fields the atomic medium presents exactly the same absorptive and dispersive responses for these beams. This regime, that may be called electromagnetically mutual induced transparency (EMIT), is broken down when we introduce a detuning in one of the fields, creating in this way a phase difference between them leading to their decorrelation.

We use two criteria for continuous variable systems to distinguish between quantum and classical correlations. We begin our analysis with the criterion of the inferred variances, described theoretically in reid  and experimentally implemented in kimble1 ; kimble . Let us suppose that we are interested in the inferrence of the probe field amplitude ($`\varphi =0`$) and phase ($`\varphi =\pi /2`$) quadratures from measurements of the pump field quadratures. In this case, the inferred variances of the probe quadratures are defined by the equations

$`\mathrm{\Delta }_{inf}^2Y_{2,0}(t)\left(\widehat{Y}_{2out,0}^{^{}}(t)\eta _0\widehat{Y}_{1out,0}^{^{}}(t)\right)^2,`$ (7)
$`\mathrm{\Delta }_{inf}^2Y_{2,\pi /2}(t)\left(\widehat{Y}_{2out,\pi /2}^{^{}}(t)+\eta _{\pi /2}\widehat{Y}_{1out,\pi /2}^{^{}}(t)\right)^2,`$ (8)
where $`\widehat{Y}_{jout,0}^{^{}}(t)=\widehat{Y}_{jout,0}(t)\widehat{Y}_{jout,0}(t)`$ and $`\widehat{Y}_{jout,\pi /2}^{^{}}(t)=\widehat{Y}_{jout,\pi /2}(t)\widehat{Y}_{jout,\pi /2}(t)`$ with $`j=1,2`$. The parameters $`\eta _0`$ and $`\eta _{\pi /2}`$ take into account the non-perfect correlation between the fields and the non-ideal efficiency of the measurement procedure. The values of these parameters are taken in order to minimize the inferred variances (7) and (8), allowing the following criterion, in the frequency domain, for the entanglement of the pump and probe fields

$$[\mathrm{\Delta }_{inf}^2Y_{2,0}(\mathrm{\Omega })]_{min}[\mathrm{\Delta }_{inf}^2Y_{2,\pi /2}(\mathrm{\Omega })]_{min}<1.$$
(9)
That is to say, if the product of the inferred variances is less than 1, then the correlation between the fields has a quantum nature. In Fig. 2 we show the product of the inferred variances for the probe field. As expected, outside de EIT region, the inequality (9) is violated since the fields are uncorrelated (see Fig. 1). However, in the EIT condition, the pump and probe fields became quantum correlated.

The other criterion used to establish the nature of the fields correlation is the theorem of Duan et al. zoller , which we will abreviate as DGCZ. Taking $`a=1`$, introducing the equivalence between the fields quadrature operators and the operators defined in zoller  as $`\widehat{x}_1\widehat{Y}_{1out,0}`$, $`\widehat{p}_1\widehat{Y}_{1out,\pi /2}`$, $`\widehat{x}_2\widehat{Y}_{2out,0}`$ and $`\widehat{p}_2\widehat{Y}_{2out,\pi /2}`$, and using the commutation relation $`[\widehat{x}_j,\widehat{p}_j^{^{}}]=2i\delta _{jj^{^{}}}`$ derived from the definition of the quadrature operators, we find the following necessary condition to prove that the join state of the pump and probe fields is separable

$$(\mathrm{\Delta }\widehat{u})^2_\rho +(\mathrm{\Delta }\widehat{v})^2_\rho 4.$$
(10)
Since this last inequality provides a necessary condition for the separability of the join fields state, then its violation is a sufficient condition for the inseparability or entanglement between the fields. In Fig. 3 we plot the left hand side of the inequality (10) and, again, outside the EIT window the equality is satisfied since the fields are uncorrelated. In the EIT condition, the violation of (10) indicates that the pump and probe fields are entangled. We must point out that both criteria report about 40 % of entanglement of the fields. This amount of quantum correlation is limited by the decay rate of the coherence between the two ground states.

The predicted quantum correlation between the pump and probe fields in the EIT condition is associated to the existence of phase quadrature squeezing in both fields. This squeezing, produced in the coherent situation, depends on the intensities of the source fields and, as it can be observed from the bistability response of the detected fields, it is produced at the turning point of the bistability curve and is accompanied by an excess noise in the corresponding conjugate quadrature amplitude tobe .

These results may be somehow unexpected because it is believed that in the EIT condition the field fluctuations are not altered for field intensities higher or comparable to the saturation intensity of the atomic transition. As we showed theoretically and experimentally exper  this is not the case. In the EIT situation there is a coherent atom-mediated exchange of photons between the pump and probe fields that preserves their mean intensities and at the same time modifies their quantum noise properties creating a correlation between them. This modification of the field quantum fluctuations is a direct consequence of the strong coherence induced in the atomic medium by the two beams.

So far, the correlation properties of the pump and probe beams have not been extensively studied in the EIT experiments. Entanglement between two single-photon pulses (quantum fields) has been predicted before in a coherently-prepared medium by two classical beams lukin . In this paper, we show that the same intense beams used to prepare the transparent nonlinear medium, in particular circumstances, become entangled even when the investigated system is subjet to the influence of a reservoir and consequently the quantum correlations are predicted in a system that is not “pure”. Since we demonstrated the existence of entanglement between the light beams only, our result suggests that there exists stronger quantum correlations in the system atom – pump field – probe field. Another remarkable point is that the investigation of the correlations and quantum fluctuations of the light beams in the EIT situation provides a precise tool to determine the natural width of the EIT resonance, and this can be a very powerful method to attain the highest sensibility in detecting coherent effects in atomic media.

Given the possible technological applications of such intense entangled fields, it is our understanding that the statistical properties of these fields deserve further experimental investigations in the near future. Not only do they open the possibility to use such systems as a macroscopic resource for different quantum technologies, but they also help understanding the nature of entanglement and how it may arise from non-linear couplings in macroscopic media.

Finally, in Fig. 4 we present an experimental setup that can be employed to measure the entanglement between the pump and probe fields. Considering the light beams have orthogonal polarizations, they can be separated using polarizing cube beam splitters and then the variance of their quadratures and the correlation between them can be determined utilizing homodyne detectors. From practical considerations, the analysis frequency must be chosen as low as possible since in this case we have more sensibility to detect correlated photons.

M.F.S. would like to thank Prof. S. Salinas for his hospitality at the University of São Paulo. The authors acknowledge the financial support of the Brazilian agencies CAPES, FAPESP and CNPq.


text:
# Assessing System Agreement and Instance Difficulty in the Lexical Sample Tasks of Senseval-2

## 1 Introduction

This paper presents a post-mortem analysis of the English and Spanish lexical sample tasks of Senseval-2. Two closely related questions are considered. First, to what extent did the competing systems agree? Did systems tend to be redundant and have success with many of the same test instances, or were they complementary and able to disambiguate different portions of the instance space? Second, how much did the difficulty of the test instances vary? Are there test instances that proved unusually difficult to disambiguate relative to other instances?

We address the first question via a series of pairwise comparisons among the participating systems that measures their agreement via the kappa statistic. We also introduce a simple measure of the degree to which systems are complementary called optimal combination. We analyze the second question by rating the difficulty of test instances relative to the number of systems that were able to disambiguate them correctly.

Nearly all systems that received official scores in the Spanish and English lexical sample tasks of Senseval-2 are included in this study. There are 23 systems included from the English lexical sample task and eight from the Spanish. Table 1 lists the systems and shows the number of test instances that each disambiguated correctly, both by part of speech and in total.

## 2 Pairwise System Agreement

Assessing agreement among systems sheds light on whether their combined performance is potentially more accurate than that of any of the individual systems. If several systems are largely in agreement, then there is little benefit in combining them since they are redundant and will simply reinforce one another. However, if some systems disambiguate instances that others do not, then they are complementary and it may be possible to combine them to take advantage of the different strengths of each system to improve overall accuracy.

The kappa statistic \[Cohen, 1960\] is a measure of agreement between multiple systems (or judges) that is scaled by the agreement that would be expected just by chance. A value of 1.00 suggests complete agreement, while 0.00 indicates pure chance agreement. Negative values indicate agreement less than what would be expected by chance. \[Krippendorf, 1980\] points out that it is difficult to specify a particular value of kappa as being generally indicative of agreement. As such we simply use kappa as a tool for comparison and relative ranking. A detailed discussion on the use of kappa in natural language processing is presented in \[Carletta, 1996\].

To study agreement we have made a series of pairwise comparisons among the systems included in the English and Spanish lexical sample tasks. Each pairwise combination is represented in a 2 $`\times `$ 2 contingency table, where one cell represents the number of test instances that both systems disambiguate correctly, one cell represents the number of instances where both systems are incorrect, and there are two cells to represent the counts when only one system is correct. Agreement does not imply accuracy, since two systems may get a large number of the same instances incorrect and have a high rate of agreement.

Tables 2 and 3 show the system pairs in the English and Spanish lexical sample tasks that exhibit the highest level of agreement according to the kappa statistic. The values in the both–one–zero column indicate the percentage of instances where both systems are correct, where only one is correct, and where neither is correct. The top 15 pairs are shown for nouns and verbs, and the top 10 for adjectives. A complete list would include about 250 pairs for each part of speech for English and 24 such pairs for Spanish.

The utility of kappa agreement is confirmed in that system pairs known to be very similar have correspondingly high measures. In Table 2, duluth2 and duluth3 exhibit a high kappa value for all parts of speech. This is expected since duluth3 is an ensemble approach that includes duluth2 as one of its members. The same relationship exists between duluth7 and duluth8 in the Spanish lexical sample, and comparable behavior is seen in Table 3.

A more surprising case is the even higher level of agreement between the most common sense baseline and the lesk corpus baseline shown in Table 2. This is not necessarily expected, and suggests that lesk corpus may not be finding a significant number of matches between the Senseval contexts and the WordNet glosses (as the lesk algorithm would hope to do) but instead may be relying on a simple default in many cases.

In previous work \[Pedersen, 2001\] we propose a 50-25-25 rule that suggests that about half of the instances in a supervised word sense disambiguation evaluation will be fairly easy for most systems to resolve, another quarter will be harder but possible for at least some systems, and that the final quarter will be very difficult for any system to resolve. This same idea could also be expressed by stating that the kappa agreement between two word sense disambiguation systems will likely be around 0.50. In fact this is a common result in the full set of pairwise comparisons, particularly for overall results not broken down by part of speech. Tables 2 and 3 only list the largest kappa values, but even there kappa quickly reduces towards 0.50. These same tables show that it is rare for two systems to agree on more than 60% of the correctly disambiguated instances.

## 3 Optimal Combination

An optimal combination is the accuracy that could be attained by a hypothetical tool called an optimal combiner that accepts as input the sense assignments for a test instance as generated by several different systems. It is able to select the correct sense from these inputs, and will only be wrong when none of the sense assignments is the correct one. Thus, the percentage accuracy of an optimal combiner is equal to one minus the percentage of instances that no system can resolve correctly.

Of course this is only a tool for thought experiments and is not a practical algorithm. An optimal combiner can establish an upper bound on the accuracy that could reasonably be attained over a particular sample of test instances.

Tables 4 and 5 list the top system pairs ranked by optimal combination (1.00 - value in zero column) for the English and Spanish lexical samples. Kappa scores are also shown to illustrate the interaction between agreement and optimal combination. Optimal combination is maximized when the percentage of instances where both systems are wrong is minimized. Kappa agreement is maximized by minimizing the percentage of instances where one or the other system (but not both) is correct. Thus, the only way a system pair could have a high measure of kappa and a high measure of optimal combination is if they were very accurate systems that disambiguated many of the same test instances correctly.

System pairs with low measures of agreement are potentially quite interesting because they are the most likely to make complementary errors. For example, in Table 5 under nouns, the alicante system has a low level of agreement with all of the other systems. However, the measure of optimal combination is quite high, reaching 0.89 (1.00 - 0.11) for the pair of alicante and jhu. In fact, all seven of the other systems achieve their highest optimal combination value when paired with alicante.

This combination of circumstances suggests that the alicante system is fundamentally different than the other systems, and is able to disambiguate a certain set of instances where the other systems fail. In fact the alicante system is different in that it is the only Spanish lexical sample system that makes use of the structure of Euro-WordNet, the source of the sense inventory.

## 4 Instance Difficulty

The difficulty of disambiguating word senses can vary considerably. A word with multiple closely related senses is likely to be more difficult than one with a few starkly drawn differences. In supervised learning, a particular sense of a word can be difficult to disambiguate if there are a small number of training examples available.

Table 6 shows the distribution of the number of instances that are successfully disambiguated by a particular number of systems in both the English and Spanish lexical samples. The value under the # column shows the number of systems that are able to disambiguate the number of noun, verb, adjective and total instances shown in the row. The average number of training examples available for the correct answers associated with these instances is shown in parenthesis. For example, the first line shows that there were 59 noun instances that no system (of 23) could disambiguate, and that there were on average 16 training examples available for each of the correct senses for these 59 instances.

Two very clear trends emerge. First, there are a substantial number of instances that are not disambiguated correctly by any system (262 in English, 228 in Spanish) and there are a large number of instances that are disambiguated by just a handful of systems. In the English lexical sample, there are 1,277 test instances that are correctly disambiguated by five or fewer of the 23 systems. This is nearly 30% of the test data, and confirms that this was a very challenging set of test instances.

There is also a very clear correlation between the number of training examples available for a particular sense of a word and the number of systems that are able to disambiguate instances of that word correctly. For example, Table 6 shows that there were 174 English verb instances that no system disambiguated correctly. On average there were only 6 training examples for the correct senses of these instances. However, there were 28 instances that all 23 English systems were able to disambiguate. For these instances an average of 47 training examples were available for each correct sense.

This correlation between instance difficulty and number of training examples may suggest that future Senseval exercises provide a minimum number of training examples for each sense, or adjust the scoring to reflect the difficulty of disambiguating a sense with very few training examples.

Finally, we assess the difficulty associated with word types by calculating the average number of systems that were able to disambiguate the instances associated with that type. This information is provided for the English and Spanish lexical samples in Tables 7 and 8. Each word is shown with its part of speech, the number of test instances, and the average number of systems that were able to disambiguate each of the test instances.

The verb collaborate is the easiest according to this metric in the English lexical sample. It has 30 test instances that were disambiguated correctly by an average of 20.2 of the 23 systems. The verb find proves to be the most difficult, with 68 test instances disambiguated correctly by an average of 4.2 systems. A somewhat less extreme range of values is observed for the Spanish lexical sample in Table 8. The adjective claro had 66 test instances that were disambiguated correctly by an average of 7.6 of the 8 systems. The most difficult word was the verb conducir, which has 54 test instances that were disambiguated correctly by an average of 3.8 systems.

## 5 Conclusion

This paper presents an analysis of the results from the English and Spanish lexical sample tasks of Senseval-2. The analysis is based on the kappa statistic and a measure known as optimal combination. It also assesses the difficulty of the test instances in these lexical samples. We find that there are a significant number of test instances that were not disambiguated correctly by any system, and that there is some correlation between instance difficulty and the number of available training examples.

## 6 Acknowledgments

This work has been partially supported by a National Science Foundation Faculty Early CAREER Development award (#0092784).


text:
# 1 Introduction

## 1 Introduction

Non-perturbative renormalization techniques  have become a crucial ingredient in lattice computations of fundamental QCD parameters and hadronic matrix elements. They allow the limitations of lattice perturbation theory to be overcome in the determinations of the relevant renormalization constants by removing the uncertainties due to unknown higher order terms.

The RI/MOM technique proposed in Refs.  has been successfully applied to compute renormalization constants of composite fermion operators in many lattice regularizations: Wilson fermions , , Kogut–Susskind fermions , domain-wall fermions  and overlap fermions . The method imposes RI/MOM renormalization conditions on conveniently defined amputated Green functions computed non-perturbatively between off-shell quark states at large virtuality in a fixed gauge. In practical computations the renormalization scale $`\mu ^2`$, determined from the virtuality $`p^2`$ of the external states, must satisfy the condition $`\mathrm{\Lambda }_{\text{QCD}}\mu O(1/a)`$ to ensure the reliability of the results .

The RI/MOM renormalization constants depend in general on the external states and therefore on the gauge in which the RI/MOM renormalization conditions have been imposed. For gauge-invariant operators the coefficients needed to match in a given gauge-invariant scheme will cancel these dependences up to higher orders in the continuum perturbative expansion and up to discretization errors. A necessary assumption in this procedure is that, at large virtualities, the continuum limit of the renormalized gauge-dependent matrix elements, computed by imposing the non-perturbative Landau lattice gauge fixing, lead to the corresponding values computed in perturbation theory up to higher order . The existence of both continuum  and lattice Gribov copies  and the numerical noise that they can generate  (for a recent review see Ref. ) is an unsolved problem of the lattice non-perturbative gauge fixing. The real concern is, of course, the influence that these phenomena may have on the values of the renormalization constants, when computed using non-gauge-invariant quantities as in the RI/MOM scheme.

The aim of this paper is to study the systematics induced by the gauge-fixing procedure on the RI/MOM determinations of the renormalization constants. We have computed the quark field ($`Z_\psi `$), the axial-vector ($`Z_A`$) and scalar ($`Z_S`$) renormalization constants, by imposing the standard lattice Landau gauge and another realization of the Landau condition that, in the following, will be called Landau1 gauge. Whereas the two realizations impose the same gauge-fixing condition in the naive continuum limit, they are affected by different Gribov ambiguities. By comparing the different sets of results for the renormalization constants, we have found negligible differences within our statistical errors.

As will be discussed in section 2, a very interesting feature of the Landau1 gauge is that it can be generalized to impose a generic covariant gauge on the lattice, as proposed in Ref. . By exploiting this opportunity, we have extended the study of the gauge dependence to off-shell Green functions measured in a generic covariant gauge.

## 2 Gauge-Fixing Conditions

In this section we briefly describe the continuum analogous of the non-perturbative gauge-fixing conditions whose discretized forms have been used to obtain the numerical results presented below.

In the standard Landau gauge the expectation value of a gauge-dependent operator is given by

$$𝒪_F=\delta A_\mu 𝒪e^{S(A)}\mathrm{\Delta }_F[A]\delta (_\mu A_\mu ),$$
(1)
where the Faddeev–Popov factor is defined as

$$\mathrm{\Delta }_F[A]=det\left(_\mu D_\mu [A]\right)$$
(2)
and the gauge-fixing condition

$$_\mu A_\mu ^G=0\text{with periodic boundary conditions},$$
(3)
is enforced non-perturbatively by minimizing the Gribov functional

$$F_A[G]A^G^2=d^4x\text{Tr}\left[A_\mu ^GA_\mu ^G\right].$$
(4)
Since the functional in Eq. (4) can have multiple minima in the non-perturbative regime, the corresponding gauge-fixing condition is incomplete, a fact that is known as the Gribov copy problem .

In order to study the effects induced by Gribov’s ambiguities on gauge-dependent matrix elements, the Landau gauge can also be fixed by minimizing a different functional

$$H_A[G]d^4x\text{Tr}\left[(_\mu A_\mu ^G)^2\right],$$
(5)
which reaches its minima when

$$D_\nu _\nu _\mu A_\mu ^G=0,$$
(6)
where $`D_\nu `$ is the covariant derivative. Equation (6), supplemented with periodic boundary conditions, defines a new gauge which we call Landau1 in the following. In this case the expectation value of a gauge-dependent operator is given by

$$𝒪_H=\delta A_\mu 𝒪e^{S(A)}\mathrm{\Delta }_H[A]\delta (D_\nu _\nu _\mu A_\mu ).$$
(7)
where $`\mathrm{\Delta }_H[A]`$ is defined analogously as in Eq. (2) . By comparing Eqs. (1) and (7) it is straightforward to show that the two procedures lead to the same values of gauge-dependent matrix elements in the perturbative regime .

The solutions of Eq. (6) are all absolute minima of $`H_A[G]`$, at variance with the copies associated with the standard functional in Eq. (4). Therefore we expect a different distribution of Gribov’s copies for the two procedures, which could generate different effects in off-shell matrix elements. It is clear that Eq. (6) could also have “spurious” solutions, which correspond to zero modes of the operator $`D_\nu _\nu `$ and do not satisfy the gauge condition in Eq. (3). The numerical results of the exploratory study presented in Ref. , and, a posteriori, the new results presented below, indicate that such a possibility can be discarded.

The functional $`H_A[G]`$ has been proposed in the literature  because it can be easily generalized to the following form

$$H_{(A,\mathrm{\Lambda })}[G]d^4x\text{Tr}\left[(_\mu A_\mu ^G\mathrm{\Lambda })(_\nu A_\nu ^G\mathrm{\Lambda })\right],$$
(8)
where $`\mathrm{\Lambda }(x)`$ are matrices belonging to the Lie algebra of the group and are generated according to a Gaussian distribution with the variance $`\xi `$ fixed at a given value. It is the value of this parameter that determines the corresponding covariant gauge. The $`H_{(A,\mathrm{\Lambda })}[G]`$ minima satisfy the gauge condition

$$D_\nu _\nu (_\mu A_\mu ^G\mathrm{\Lambda })=0,\text{with periodic boundary conditions}.$$
(9)
The expectation value of a given operator is obtained by

$$𝒪_H_\xi =\delta \mathrm{\Lambda }e^{\frac{1}{\xi }{\scriptscriptstyle d^4x\text{Tr}[\mathrm{\Lambda }^2]}}\delta A_\mu 𝒪e^{S(A)}\mathrm{\Delta }_H[A]\delta (_\mu A_\mu \mathrm{\Lambda }).$$
(10)
For a complete description of the method, see Ref. . In the following we perform a further check of the gauge-dependence comparing the results obtained in the Landau gauge and in one of the lattice covariant gauges.

## 3 Numerical Results and Discussion

We have used a sample of $`80`$ SU(3) gauge configurations retrieved from the repository at the “Gauge Connection” (http://qcd.nersc.gov/), which were generated with the Wilson gluonic action at $`\beta =6.0`$ and $`V\times T=16^3\times 32`$ . By using the discretized version (see  for details) of the functionals in Eqs. (4), (5) and (8) we have rotated each configuration in the Landau gauge, in the Landau1 gauge and in the covariant gauge with the bare value of the gauge fixing parameter $`\xi =8`$ respectively.

The various gauges have been enforced with a quality factor $`\theta <10^6`$ (see  for details) with the Origin2000 computers at Boston University. The inversion of the Wilson–Dirac operator has been performed on an INFN APE100 machine for masses corresponding to the values of the $`k`$ parameter shown in Table 1.

Once the quark propagator $`S(x,0)=\psi (x)\overline{\psi }(0)`$ has been computed for each mass and gauge and Fourier-transformed, we have determined the inverse Euclidean bare-quark propagator $`S^1(pa)`$. We have computed the amputated Green functions of the local quark bilinears $`O_\mathrm{\Gamma }(x)=\overline{\psi }(x)\mathrm{\Gamma }\psi (x)`$ in the Fourier space 222In what follows we adopt the conventions of Ref. .

$$\mathrm{\Lambda }_\mathrm{\Gamma }(pa)=S^1(pa)G_\mathrm{\Gamma }(pa)S^1(pa),$$
(11)
where

$$G_\mathrm{\Gamma }(p)=d^4xd^4ye^{ip(xy)}\psi (x)O_\mathrm{\Gamma }(0)\overline{\psi }(y)$$
(12)
and $`\mathrm{\Gamma }`$ is a generic Dirac matrix corresponding to $`\mathrm{\Gamma }=\{\text{A},\text{S}\}=\{\gamma _\mu \gamma _5,\text{1I}\}`$.

In the RI/MOM scheme, the wave-function renormalization constant $`Z_\psi `$ is defined as

$$Z_\psi (\mu a)=\underset{m0}{lim}i\frac{1}{12}\text{Tr}\left[\frac{S^1(pa)}{\text{ / }p}\right]_{p^2=\mu ^2}.$$
(13)
In order to avoid numerical derivatives, it is convenient to use

$$Z_\psi ^{}(\mu a)=\underset{m0}{lim}i\frac{1}{12}\frac{\text{Tr}_{\mu =1,4}\gamma _\mu \mathrm{sin}(p_\mu a)S^1(pa)}{4_{\mu =1,4}\mathrm{sin}^2(p_\mu a)}|_{p^2=\mu ^2},$$
(14)
which differs from $`Z_\psi `$ by a finite term of order $`\alpha _s`$. The matching between the two renormalization constants can be computed in continuum perturbation theory. Since this correction is negligible at the level of the numerical precision of our calculation (see below), we will not take it into account in the following. Therefore throughout this paper we assume $`Z_\psi =Z_\psi ^{}`$.

The renormalization constant $`Z_\mathrm{\Gamma }^{\text{RI}}(\mu a,g_0)`$, which defines the renormalized operator $`\widehat{𝒪}_\mathrm{\Gamma }^{\text{RI}}=Z_\mathrm{\Gamma }^{\text{RI}}𝒪_\mathrm{\Gamma }`$, is fixed by imposing in the chiral limit the renormalization condition

$$Z_\mathrm{\Gamma }^{\text{RI}}(\mu a)Z_\psi ^1(\mu a)\text{Tr}\text{I}\text{P}_\mathrm{\Gamma }\mathrm{\Lambda }_\mathrm{\Gamma }(pa)|_{p^2=\mu ^2}=1,$$
(15)
where $`\text{I}\text{P}_\mathrm{\Gamma }`$ is a suitable projector on the tree-level amputated Green function .

Even though $`𝒪_A`$ and $`𝒪_S`$ are gauge-independent operators, their matrix elements between quark states acquire a gauge dependence. In Fig. 1 are shown the numerical results for $`Z_\psi `$, $`Z_A`$ and $`Z_S`$ calculated in the Landau and in the Landau1 gauge, as a function of the square lattice momenta

$$(ap_L^\mu )^2=4\mathrm{sin}^2\left(\frac{ap^\mu }{2}\right).$$
(16)
In the perturbative region, $`Z_\psi (a\mu )`$ has a flat dependence on the renormalization scale. This is expected since in the Landau gauge it has zero anomalous dimension at the leading order. Also $`Z_A(a\mu )`$ shows a flat behaviour due to the compatibility of the RI/MOM renormalization conditions with the axial Ward identities. The value of $`Z_S(a\mu )`$ increases with a logarithmic behaviour, as expected in perturbation theory (see also Fig. 3). Our data are in very good agreement with the results reported in the literature .

The results for all correlators, corresponding to Landau and Landau1 gauges, coincide within the statistical errors. As a consequence we can conclude that the data do not show any residual effect due to the presence of lattice Gribov’s copies in the statistical sample of configurations generated and for the lattice used. This represents one of the main results of this paper.

To show the amount of gauge dependence that can be found in our correlators and therefore our sensitiveness to the gauge, in the last part of this section we compare the results obtained in the Landau and in the covariant gauge with the bare gauge parameter fixed to $`\xi =8`$.

The gauge dependence of the primary quantities $`Z_\psi `$, $`\mathrm{\Lambda }_A`$ and $`\mathrm{\Lambda }_S`$ measured on the lattice is clearly visible in Fig. 2, where they are plotted for the Landau and the covariant gauge as a function of the momenta $`(ap)^2`$. All three quantities show a statistically significant gauge dependence while the renormalization constants $`Z_A`$, $`Z_S`$ shown in Fig. 3, obtained by computing the ratios as indicated in Eq. (15), do not exhibit any visible gauge dependence within the statistical precision reached in our simulation. For these quantities, the fluctuations of the simulation hide the weak gauge dependence that is expected in perturbation theory from the next-to-leading order terms in $`\alpha _s`$. These results indicate an upper limit to the numerical troubles that can be expected in the RI/MOM renormalization constants because of the lattice gauge fixing.

## 4 Conclusions

We have proposed a method to study the effects of Gribov’s copies on lattice gauge-dependent matrix elements. We have implemented it numerically in quenched QCD with Wilson fermions and we have computed the matrix elements needed to fix the quark field, the axial-vector and scalar renormalization constants in the RI/MOM non-perturbative renormalization technique. Our numerical results do not show any significant differences between data coming from the standard Landau and the Landau1 gauge. We do not find any relevant effect that can be traced back to the incompleteness of the gauge fixing for a number of gauge-field configurations and for lattice sizes of the order of those used in this calculation. The main indication of our paper is therefore that, at the level of the precision that we have reached, the systematics due to the numerical gauge-fixing machinery has negligible effect on the matrix elements needed for implementing the RI/MOM scheme.

## Acknowledgements

We warmly thank N. Cabibbo, M. L. Paciello, F. Rapuano, C. Rebbi, G.C. Rossi and M. Testa for many fruitful discussions on this subject. We thank the Center for Computational Science of Boston University, where part of the numerical computation has been done. We also gratefully acknowledge the use of the gauge configurations produced by the authors in Ref. . This work was partially supported by the M.U.R.S.T. and by the European Community’s Human Potential Programme under contract HPRN-CT-2000-00145, Hadrons/Lattice QCD.