text:
# Evolution of photoionization and star formation in starbursts and H ii galaxies

## 1 Introduction

The “starburst” phenomenon calls to mind a class of objects dominated by the radiation of massive stars (Gallego et al. 1995) embedded in a dusty H ii region. Modelling starbursts consistently is the key to interpret the properties of the actively star-forming galaxies detected at a redshift $`z>2`$ (Giavalisco et al. 1996; Madau et al. 1996; Steidel et al. 1996; Lowenthal et al. 1997). To this purpose, the extensive study of local samples is a requisite step.

Local starbursts span a wide range of types, from individual H ii regions in spiral arms to blue compact H ii galaxies and huge kpc-scale nuclear starbursts (Coziol et al. 1998). Types differ from one another in their main characteristics (line ratios, equivalent widths, colors), and presumably have different stellar populations and physical conditions. So, is it possible to find correlations between the basic properties: age, initial mass function, metallicity, relative distributions of stars, gas and dust? Among these, which one dominates spectral properties? In which range do parameters vary? Answering these questions requires a consistent model of starbursts coupling the evolution of stars, dust and gas.

The first studies dealing with emission lines of H ii regions (Shields 1974, 1978; Stasińska 1978, 1980; McCall et al. 1985) relied on single-star photoionization models. Average physical properties of large starburst samples were deduced from these pioneering works. As an example, a relation between the metallicity $`Z`$ and the ionization parameter $`U`$ was derived by Dopita & Evans (1986) by fitting the emission line ratios from large samples of extragalactic H ii regions.

More realistic stellar populations of star clusters were computed by McGaugh (1991) with the Salpeter (1955) initial mass function (IMF), to calibrate metallicity dependant line ratios. Further improvements were the implementation of theoretical tracks of massive stars, allowing to follow the evolution in time of a starburst (e.g. García-Vargas & Díaz 1994). Thanks to the computation of evolutionary tracks for various metal abundances, the influence of the metallicity could also be studied (Cid-Fernandes et al. 1992; Cerviño & Mas-Hesse 1994; Olofsson 1997; García-Vargas et al. 1995a,b; Stasińska & Leitherer 1996). Evans (1991) emphasized the importance of using up-to-date stellar atmosphere models to calibrate nebular diagnostics. The most recent models take into account modern, updated stellar physics (Leitherer & Heckman 1995; Leitherer et al. 1999), in particular the impact on stellar spectra of line blanketing and of departures from local thermodynamic equilibrium (Stasińska & Schaerer 1997).

Till now, only a few large datasets covering a wide range of physical conditions have been analyzed with state-of-the-art models, coupling evolving stellar populations and photoionization. Stasińska & Leitherer (1996), for example, restricted their analysis to metal-poor ($`Z<Z_{}/4`$) objects. As a matter of fact, statistical properties of starbursts are not definitely established. A metal-dependent IMF was early proposed as an explanation for the \[O iii\]λ5007/H$`\beta `$ decrease at high $`Z`$ (Terlevich 1985; Shields & Tinsley 1976). Although a standard IMF seems in agreement with the bulk of observations (Leitherer 1998), the IMF slope (see e.g. Eisenhauer et al. 1998; Greggio et al. 1998) and mass cut-offs (Goldader et al. 1997) are still under debate.

Our aim is hereafter to study the evolution of the spectral properties of aging starbursts with the help of models taking into account the effects of metallicity, geometry, dust and the excitation level of the gas. To avoid normalization problems, we selected relative properties independent of distance (colors, equivalent widths and line ratios). We focused on emission line ratios at close wavelengths – thus reducing reddening effects – to study the ionizing spectrum and the excitation of the gas. Equivalent widths and colors are used preferentially to study age effects, as well as the impact of the spatial distribution of stars, gas and dust.

The observational sample is presented in Sect. 2. The coupling of the codes PÉGASE, to model the evolution of star formation and stellar emission, and CLOUDY, for a consequent photoionization of the gas by massive stars in a given geometry, is described in Sect. 3. The relation between the metallicity $`Z`$ and the ionization parameter $`U`$ is analyzed in Sect. 4, while equivalent widths and colors are considered in Sect. 5. The contribution of an underlying population is analyzed in Sect. 6. Discussion and conclusion are respectively in Sect. 7 and 8.

## 2 A selected sample of starbursts and H ii galaxies

Our data set (Table 1) is a selection of 754 objects from observational samples of galaxies having properties characteristic of starbursts: intense emission lines (French 1980; Veilleux & Osterbrock 1987; Terlevich et al. 1991), strong far-infrared emission (Leech et al. 1989; Veilleux et al. 1995), or excess of blue or ultraviolet emission (Storchi-Bergmann et al. 1995); in addition, Contini et al. (1998) require galaxies to have a bar. Objects dominated by an active galactic nucleus (AGN) are excluded from our sample, though some of them may be faintly contaminated by nuclear activity. Effects of aperture and differences in data reduction processes are considered in Sect. 2.2.

### 2.1 Spectral classification of the original sample

Objects listed in the source papers are classified in two families. The first one includes metal-poor H ii galaxies (H iiGs) and extra-nuclear H ii regions, both likely devoid of active nucleus. The second one gathers galactic central regions where the ionized-gas emission could be partly due to shocks (Kim et al. 1998) or AGNs (see Heckman 1991 and references therein). To select only starburst nucleus galaxies (SBNGs) in this last class, we keep the classification of the authors, which is based on the pioneering works of Baldwin et al. (1981, hereafter BPT) and Veilleux & Osterbrock (1987, hereafter VO).

The Terlevich et al. (1991) sample is mainly composed of H iiGs with some possible SBNGs. The authors used two of the BPT criteria based on the \[O iii\]λ5007/H$`\beta `$ and \[O ii\]λ3727/\[O iii\]λ5007 ratios. We also analyzed their sample with the VO method and obtained the same classification.

Contini et al. (1998) adopted the VO criteria. Their sample includes SBNGs and giant H ii regions located in IRAS barred spiral galaxies, selected from Mazzarella & Balzano (1986) and listed in the Lyon-Meudon Extragalactic Database (LEDA).

Two samples contain only SBNGs located inside IRAS galaxies. Veilleux et al. (1995) followed both BPT and VO criteria, while Leech et al. (1989) used only some of the criteria of BPT. Actually, an important part of their sample lies far beyond the limit reported by VO between H ii-like regions and AGNs in the \[O iii\]/H$`\beta `$ vs. \[O i\]/H$`\alpha `$ and \[O iii\]/H$`\beta `$ vs. \[N ii\]/H$`\alpha `$ diagrams. These objects are clearly misclassified and have been excluded from our selection. Storchi-Bergmann et al. (1995) reported observations of emission line galaxies of various types. Following Coziol et al. (1998), all compact galaxies are hereafter classified as H ii galaxies, whether dwarf or not. Finally, we include in our selection the observations of five starbursts from the sample of Balzano (1983) reported by VO, and the data of French (1980) on 14 H iiGs.

### 2.2 Homogeneity of reddening, aperture and absorption corrections

Most line ratios considered here are computed from close lines; so, reddening effects should be small and will be neglected in the following.

In most samples, the contamination, through large apertures, of the starburst light by the host galaxy population is likely to increase the continuum emission and to dilute the equivalent width of emission lines. To avoid this problem, we will restrict the analysis of colors and equivalent widths to the sample of Contini et al. (1998). These data were acquired with a long-slit spectrograph, but the spectra presented by the authors correspond to H$`\alpha `$ emitting regions exclusively. Hence, the contamination of the starburst continuum by the environment of the host galaxy should be very weak.

The host galaxy can also modify the apparent line ratios involving Balmer lines through the presence of stellar absorption lines. Emission line fluxes are corrected for this effect as follows:

$$F_{\mathrm{line}}^{\mathrm{corr}}=F_{\mathrm{line}}^{\mathrm{obs}}\left(1+\frac{W_{\mathrm{line}}^{\mathrm{abs}}}{W_{\mathrm{line}}^{\mathrm{obs}}}\right),$$
(1)
where $`F_{\mathrm{line}}^{\mathrm{corr}}`$ and $`F_{\mathrm{line}}^{\mathrm{obs}}`$ are respectively the absorption-corrected and the apparent emission line fluxes, and $`W_{\mathrm{line}}^{\mathrm{abs}}`$ and $`W_{\mathrm{line}}^{\mathrm{obs}}`$ are the absorption and measured emission equivalent widths. Not all the samples were corrected (Table 1). The correction should be negligible for the H$`\alpha `$ line because $`W_{\mathrm{H}\alpha }^{\mathrm{abs}}`$ $``$ $`W_{\mathrm{H}\alpha }^{\mathrm{obs}}`$. It is stronger for H$`\beta `$ (up to 100 % in Leech et al. 1989) and the ratio \[O iii\]/H$`\beta `$ may be affected. This ratio is systematically higher in the samples of French (1980) and Terlevich et al. (1991), mainly composed of HIIGs, than in the other ones, which include mainly SBNGs. Is this trend real, or is it due to the fact that these two data sets are not corrected for underlying stellar lines? Note that Storchi-Bergmann et al. (1995) already noticed that \[O iii\]/H$`\beta `$ is higher in H ii galaxies than in SBNGs. Moreover, the correction for absorption lines should be very small in H iiGs, since these objects have simultaneously large $`W_{\mathrm{H}\beta }^{\mathrm{obs}}`$ (67 Å in Terlevich et al. 1991) and low $`W_{\mathrm{H}\beta }^{\mathrm{abs}}`$. Hence, we conclude that the trend we observe is real and not due to the discrepancies between the reduction procedures.

## 3 The coupling of spectral evolution and photoionization

The evolutionary synthesis code we use, PÉGASE, takes into account metallicity and dust effects. PÉGASE computes the stellar spectral energy distributions (SEDs) and the metallicities of starbursts and galaxies of the Hubble sequence at any stage of evolution, within the metallicity range $`Z=10^4`$ to $`10^1`$. Typical parameters of PÉGASE are the star formation rate (SFR) and the initial mass function. The new version111This code is available at *http://www.iap.fr/users/fioc/PEGASE.html* or by anonymous ftp at *ftp.iap.fr* in */pub/from$`\mathrm{\_}`$users/pegase*. used hereafter, PÉGASE.2, is based on the evolutionary tracks of Girardi et al. (1996), Fagotto et al. (1994a,b,c) and Bressan et al. (1993). The AGB to post-AGB phases are computed following the prescriptions of Groenewegen & de Jong (1993) models. The synthetic stellar spectral library is from Kurucz (1992), modified by Lejeune et al. (1997) to fit observed colors. For details, see Fioc & Rocca-Volmerange (1999; 2000, in preparation)

The photoionization code CLOUDY (version 90.04, Ferland 1996) predicts the spectra of low- to high-density astrophysical plasmas in the most extreme astrophysical sites, such as starbursts and quasar environments. It takes into account recent changes in atomic databases and new numerical methods. In addition to its physical performances, the code was chosen for its easy Web access, its documentation and its friendly on-line help. CLOUDY proposes options to explore the influence of geometrical factors: type of geometry – spherical or plane parallel –, covering factor $`\mathrm{\Omega }/(4\pi )`$, filling factor $`f`$, and the distance $`R`$ from the source to the plasma.

The coupling of PÉGASE.2 with CLOUDY allows to compute coherently the stellar and nebular emission, and to fit simultaneously the continuum and lines from stellar populations embedded in a gas cloud. The metallicity $`Z(t)`$, and the stellar SED – in particular the energy distribution of the ionizing photons – provided by PÉGASE are used as inputs by CLOUDY. Hereafter, we call $`Q_{\mathrm{H}^0}`$ the emission rate of ionizing photons. This quantity is used to normalize the SEDs computed by PÉGASE.

The results of the current version of the code have been compared to those obtained by similar “coupled models”. They are presented in Table 2 for two sets of predictions: the model 16 from García-Vargas et al. (1995b), and the $`Z_{}`$-model computed by Stasińska & Leitherer (1996) for a $`10^6M_{}`$-cluster. In both cases, we used inputs matching the reference model (IMF, age, geometrical parameters and chemical abundances). For the first comparison test, the agreement is very good (Table 2, top). Our results also agree relatively well with those of Stasińska & Leitherer (1996), given that the evolutionary tracks, the spectral libraries and the photoionization codes are different in the two models.

### 3.1 Star formation parameters

An instantaneous star formation episode (“pure starburst”) is assumed for starbursts – a scenario well suited to study emission line properties. In Sect. 6, we will also consider an underlying stellar population formed either during a previous burst or continuously. We tentatively tested several IMFs of the form d$`N`$/d$`mm^\alpha `$, $`m[M_{\mathrm{low}},M_{\mathrm{up}}]`$: $`\alpha =2.35`$, $`m[0.1,120]M_{}`$ (Salpeter 1955); $`\alpha =3.35`$, $`m[0.1,120]M_{}`$; and $`\alpha =2.35`$, $`m[3,30]M_{}`$. The last one is truncated, as suggested by various infrared studies (Rieke et al. 1993; Lançon & Rocca-Volmerange 1996). The initial metallicity is left as a free parameter.

### 3.2 Physical properties of the gas

The emission line spectrum of an ionized nebula depends on the combination of the ionizing spectrum, the chemical composition of the gas and the so-called ionization parameter $`U`$.

The elemental abundances of the gas have been made consistent with the metallicity of the ionizing stars. The solar abundances used here are given in Table 3. All the abundances were scaled linearly according to the metallicity, except for nitrogen and helium. For the likely secondary element N, we adopted the law proposed by Coziol et al. (1999), multiplied by a factor 1.5 to account for the excess of nitrogen observed in starbursts. For He, we followed a prescription proposed by Dopita & Kewley (private communication; see Table 3 for details).

The ionization parameter $`U`$ is defined as:

$$U=\frac{Q_{\mathrm{H}^0}}{4\pi R_\mathrm{s}^2n_\mathrm{H}c},$$
(2)
where $`n_\mathrm{H}`$ is the hydrogen density, $`c`$ the speed of light, and $`R_\mathrm{s}`$ the Strömgren radius (i.e. the radius of the ionization front). Taking into account the inner radius $`R`$ (e.g. Evans & Dopita, 1985), $`R_\mathrm{s}`$ is roughly equal to:

$$R_\mathrm{s}=\left(\frac{3Q_{\mathrm{H}^0}}{4\pi \alpha _\mathrm{B}n_\mathrm{H}^2f}+R^3\right)^{1/3},$$
(3)
where $`\alpha _\mathrm{B}`$ is the case B temperature-dependent recombination coefficient and $`f`$ is the filling factor. This formula highlights the impact of each parameter on $`U`$, but is only approximate, since it assumes a constant temperature through all the nebula. The effective Strömgren radius is obtained only at the end of a CLOUDY calculation.

The spherical geometry (Stasińska & Leitherer 1996; González-Delgado et al. 1999) and plane-parallel geometry (García-Vargas & Díaz 1994; García-Vargas et al. 1997) are extreme situations for which the right-hand side of eq. (3) is respectively dominated by the first and the second term. The adopted geometry was spherical. The hydrogen density $`n_\mathrm{H}`$ was 100 cm-3 in most models, but we also computed a few cases for $`n_\mathrm{H}=10`$ and $`n_\mathrm{H}=1000`$ cm-3, to assess the specific impact of density on emission line ratios. We set the number of ionizing photons to $`Q_{\mathrm{H}^0}=10^{52}`$ photons s-1, corresponding to a maximum luminosity of $``$ $`10^{40}`$ erg s-1 in the H$`\alpha `$ emission line – a value characteristic of intense star formation sites (Shields 1990). For emission lines, the exact value of $`Q_{\mathrm{H}^0}`$ is not important as long as the density does not exceed the critical limit of collisional de-excitation, below which models with the same $`U`$, ionizing spectrum and metallicity predict the same line ratios. To obtain reasonable values for $`U`$, we varied the volume filling factor between $`10^5`$ and 1 and set $`R`$ to 4 pc.

We assume that the gas is distributed in a spherical shell covering a solid angle $`\mathrm{\Omega }`$ around the star cluster. The covering factor $`\mathrm{\Omega }/(4\pi )`$ was allowed to vary in the interval \[0.1, 1\]. The continuum luminosity from the starburst $`F_\lambda `$ (erg s-1 Å-1) at any wavelength $`\lambda `$ is:

$$F_\lambda =F_\lambda ^{\mathrm{trans}}+F_\lambda ^{\mathrm{neb}}+F_\lambda ^{\mathrm{stell}}.$$
(4)
The transmitted stellar continuum $`F_\lambda ^{\mathrm{trans}}`$ and the nebular continuum $`F_\lambda ^{\mathrm{neb}}`$ both scale with $`\mathrm{\Omega }/(4\pi )`$. $`F_\lambda ^{\mathrm{stell}}`$, scaling as $`(1\mathrm{\Omega }/(4\pi ))`$, refers to the part of the stellar spectrum not crossing the gas (and thus not contributing to the ionization process).

### 3.3 Extinction of the spectrum

Assuming a dust screen at the outer edge of the H ii region, the continuum flux is reddened as follows:

$$F_\lambda =(F_\lambda ^{\mathrm{trans}}+F_\lambda ^{\mathrm{neb}})e^{\tau _\lambda ^{\mathrm{neb}}}+F_\lambda ^{\mathrm{stell}}e^{\tau _\lambda ^{\mathrm{stell}}}.$$
(5)
Emission line fluxes are also multiplied by $`e^{\tau _\lambda ^{\mathrm{neb}}}`$. The optical depth $`\tau _V^{\mathrm{neb}}`$ was a free parameter, as well as the $`\tau _V^{\mathrm{stell}}/\tau _V^{\mathrm{neb}}`$ ratio ($`1`$). This option takes into account a possible dust concentration outside the nebula rather than on the line-of-sight not covered by the gas. The variation of the optical depth with wavelength is derived from the metallicity-dependent extinction laws of the Magellanic Clouds and our Galaxy.

## 4 Results on emission line ratios for instantaneous starburst models

We have explored the impact of the IMF, the age of the starburst, the metallicity and the ionization parameter on the following optical line ratios: \[O iii\]λ4363/\[O iii\]λ5007, \[O i\]λ6300/H$`\alpha `$, \[S ii\]λλ6716,6731/H$`\alpha `$, \[N ii\]λ6584/H$`\alpha `$ and \[O iii\]λ5007/H$`\beta `$. The value of $`f`$ was varied between $`10^5`$ and 1. The corresponding $`\mathrm{log}(U)`$ (as defined in Sect. 3.2) belong to $`[5.5,1.5]`$. Note that similar results could also be derived by varying the inner radius $`R`$ (in a plane-parallel geometry), the density $`n_\mathrm{H}`$ or a combination of both.

### 4.1 Metallicity effects

Emission line ratios are sensitive to the metallicity through two competing effects. On the one hand, the opacity increases with $`Z`$, lowering the stellar effective temperatures and thus the hardness of the ionizing spectrum. A trend towards lower excitation levels is therefore expected in high-$`Z`$ starbursts, even if there is no dependence of the IMF or the geometry on the metallicity. This effect is reinforced at high $`Z`$ by the efficient cooling by metals and thus the low electronic temperature $`T_\mathrm{e}`$.

On the other hand, the energy in metal emission lines tends naturally to increase with higher elemental abundances. For oxygen lines, this process dominates at $`Z1/5Z_{}`$, while the opposite happens at higher metallicity (e.g. McGaugh 1991). As a consequence, two models with metallicities respectively below and above this limit can predict comparable emission lines. For this reason, we present hereafter our predictions for $`Z`$ in the range $`[1/50,3/2]Z_{}`$ if no degeneracy happens in the considered diagram, but in $`[1/5,3/2]Z_{}`$ otherwise.

### 4.2 Effects of the geometrical filling factor $`f`$

For a given density $`n_\mathrm{H}`$, varying the filling factor $`f`$ changes the optical depth per unit length and the ionization structure: the smaller $`f`$, the smaller the mean ionization level, and the lower the luminosities of high excitation lines such as \[O iii\]λ5007. Low-excitation lines behave inversely.

Models of pure instantaneous starbursts, with metallicities in the range $`[1/5,3/2]Z_{}`$, are compared to the data in Fig. 1 for six values of the filling factor $`f`$: $`10^5`$, $`10^4`$, $`10^3`$, $`10^2`$, $`10^1`$, and 1. The slope of the IMF is -1.35, for a mass range of 0.1-120 M. The effects of dust are not taken into account in the models plotted on Fig. 1. The data are largely covered by the models when at least two parameters vary – here chemical and geometrical. Interestingly, different types of objects seem to occupy different zones in the $`Z`$-$`U`$ plane: H ii galaxies are in better agreement with low-$`Z`$ (1/50 $`Z_{}`$ to 3/4 $`Z_{}`$), high-$`U`$ models ($`f[10^3,1]`$); SBNGs are compatible with low values of $`f`$ ($`10^5`$ to $`10^2`$), but encompass a wide range of metallicities.

We also report on Fig. 2 the comparison of predictions with observations in following diagrams: \[O ii\]λ3727/H$`\beta `$ vs. \[O iii\]λ5007/H$`\beta `$ and \[O ii\]λ3727/\[O iii\]λ5007 vs. \[S iii\]λ9069/\[S ii\]λλ6716,6731. Note that the \[O ii\]/H$`\beta `$, \[O ii\]/\[O iii\] and \[S iii\]/\[S ii\] ratios are sensitive to reddening and, for this reason, we do not include them in the discussion. We can only notice that, given the uncertain impact of dust on these ratios, our results are compatible with the observations.

### 4.3 The apparent relation between $`Z`$ and $`U`$: degeneracy problems

A mean $`U`$-$`Z`$ relation fitting the data could easily be derived from Fig. 1, as in Dopita & Evans (1986). However, this relation would suffer from too many uncertainties: first, it would not take into account the high $`Z`$-low $`Z`$ degeneracy; second, age effects are very similar to those induced by the metallicity. This is illustrated in Fig. 4, where \[O i\]/H$`\alpha `$ and \[S ii\]/H$`\alpha `$ vs. \[O iii\]λ5007/H$`\beta `$ are plotted for $`Z_{}/5`$ at $`t=0`$, 4, 5 and 6 Myr. Variations of the IMF give rise to similar uncertainties, as shown on Fig. 6: the lower number of massive stars for a steeper or truncated IMF causes a strong decrease of both \[O iii\]/H$`\beta `$ and \[O i\]/H$`\alpha `$, just as if the metallicity was higher.

Due to this age-metallicity-IMF degeneracy, a meaningful $`Z`$-$`U`$ analytical relation cannot be established on the basis of line ratio diagrams. Moreover, Fig. 1 clearly shows that the dispersion of $`U`$, rather than the mean value of this parameter, is linked to the metallicity: while starbursts exist at any $`U`$ at low $`Z`$, only low-$`U`$ starbursts are observed at high metallicity. We consider a possible cause of this trend in the next paragraph.

### 4.4 The upper envelope of the data sequences

In the \[O i\]/H$`\alpha `$ vs. \[O iii\]/H$`\beta `$ and \[S ii\]/H$`\alpha `$ vs. \[O iii\]/H$`\beta `$ diagrams, many data points lie above predictions, even when the misclassified objects of Leech et al. (1989) are excluded from the sample. To check whether this discrepancy can be solved by varying the density, we have performed two series of calculations at $`Z=Z_{}/5`$ with $`n_\mathrm{H}=10`$ and 1000 cm-3. We chose this metallicity because the corresponding sequence defines the upper envelope of our model grid in these diagrams. The results are shown in Table 4. The main one is that model sequences with different densities largely overlap: clearly, the outliers cannot be explained by a variation of $`n_\mathrm{H}`$.

Many authors have proposed that a faint contribution from shocks or AGNs could explain the intensity of the \[O i\] line in starbursts (see e.g. Stasińska & Leitherer 1996). We consider here an alternative possibility where there are, inside the starburst, both high- and low-excitation zones, due for example to a distribution of gas and stars more complex than the classical geometry adopted in our calculations. To test this hypothesis, we computed a series of sequences obtained by varying the relative weight of high- and low-excitation models ($`f=1`$ and $`f=10^5`$, respectively) for five metallicities. The results (Fig. 3) are in much better agreement with the data than pure $`U`$ sequences, as they explain both the trend and the dispersion of the observations. Another advantage of this explanation is that no other source of ionization, such as AGNs or shocks, is required.

If this hypothesis is correct, the evolution of the apparent $`U`$ with $`Z`$ in fact corresponds to the presence in low-metallicity objects of a high-$`U`$ component, which is absent in high-$`Z`$ starbursts. We come back to this point in the discussion.

### 4.5 Constraints on age and IMF

Fig. 5 shows the evolution, up to 6 Myr, of the \[O iii\]λ5007/H$`\beta `$ ratio for the models plotted on Fig. 1 and 2. The constraints on age shown by this plot are very strong, though noticeably different for H ii galaxies and SBNGs: H ii galaxies, of high-excitation level, need the high energy photons emitted by massive stars, so that only young ($``$3 Myr) models are compatible with observations (see Sect. 7); SBNGs, requiring a lower excitation level, span a wider range of ages. In any case, the rapid drop of the \[O iii\]/H$`\beta `$ ratio after 5 Myr requires the presence of stars younger than 6 Myr.

Line ratios depend on the hardness of the spectrum, and so on the IMF. The \[O i\] line, emitted from the partly ionized zone, requires high-energy photons and is thus particularly sensitive to the high-mass end of the IMF. Fig. 6 presents the sensitivity of \[O i\]λ6300/H$`\alpha `$ and \[O iii\]/H$`\beta `$ to the initial mass function. A Salpeter IMF is favored by our results, in good agreement with recent studies (e.g. Stasińska and Leitherer 1996; García-Vargas et al. 1995). An interesting point is that the steeper and truncated IMFs, while clearly excluded by the data for high-$`U`$ H iiGs, are acceptable for metal-rich starbursts.

## 5 Results on equivalent widths and colors: covering factor and extinction

In the following, we aim to keep our previous fits of emission line ratios and, simultaneously, to reproduce $`W(\mathrm{H}\alpha )`$ and $`BV`$ and $`VR`$ colors from the sample of Contini et al. (1998). The models plotted on Fig. 1 and 2 predict H$`\alpha `$ equivalent widths above 500 Å, while the observed $`W(\mathrm{H}\alpha )`$ seldom exceed 300 Å. Predicted equivalent widths in excess to the observations are a well-known problem (Bresolin et al. 1999 and references therein). The hypothesis of a “density bounded” nebula (Leitherer et al. 1996) is ruled out by the observed intensities of the low-ionization lines; the reduction of the thickness of the partly-ionized zone would lead to lower \[S ii\]/H$`\alpha `$ and \[O i\]/H$`\alpha `$ than in Fig. 1, and the models would not fit any more these line ratios.

A low covering factor $`\mathrm{\Omega }/(4\pi )`$ reduces the equivalent widths by increasing the stellar emission relative to the nebular one, without changing line ratios. A similar explanation of the faint $`W(\mathrm{H}\beta )`$ of blue compact and irregular galaxies was suggested by Mas-Hesse & Kunth (1999).

Fig. 7 shows the influence of dust and age in the $`W(\mathrm{H}\alpha )`$ vs. $`BV`$ and $`W(\mathrm{H}\alpha )`$ vs. $`VR`$ planes, for $`\mathrm{\Omega }/(4\pi )=0.1`$. For comparison, we also show $`t=0`$ and $`t=6`$ Myr sequences for $`\mathrm{\Omega }/(4\pi )=1`$. We plot only $`Z=Z_{}`$ models, since the metallicity has only weak effects on $`BV`$, $`VR`$ and $`W(\mathrm{H}\alpha )`$. Similarly, for such a low value of $`\mathrm{\Omega }/(4\pi )`$, $`f`$ hardly affects the colors and the only models on Fig. 7 are for $`f=1`$. Ages between 0 and 6 Myr, $`\tau _V^{\mathrm{stell}}/\tau _V^{\mathrm{gas}}=1`$ and $`\tau _V^{\mathrm{gas}}`$ in the range  are suggested by the $`W(\mathrm{H}\alpha )`$ vs. $`BV`$ diagram. However, $`\tau _V4`$ is excluded by the $`VR`$ data, and ages $`6`$ Myr are incompatible with line ratios (see also Sect. 7.2).

A ratio $`\tau _\lambda ^{\mathrm{stell}}/\tau _\lambda ^{\mathrm{gas}}0.5`$ (following Calzetti et al. 1997; see also Fanelli et al. 1988; Keel 1993; Mas-Hesse & Kunth 1999) can simultaneously fit $`BV`$ and $`VR`$, but requires very high extinctions, up to $`\tau _V^{\mathrm{neb}}=8`$ (Fig. 8), while Contini et al. (1998) derived a maximum value of $``$ 4 from the Balmer decrement. The $`BV`$ of the reddest objects are considered as anomalous (Contini, private communication) and are not reproduced by any model.

## 6 The contribution of an evolved population

The distribution of the data in the $`W(\mathrm{H}\alpha )`$ vs. $`BV`$ or $`VR`$ diagrams is very difficult to explain with pure instantaneous starburst models. The main problem is the very low covering factor required. A value of 0.1 seems in contradiction with previous estimates of the escaping fraction of ionizing ohotons in HII regions (see Sect. 7.2). An alternative hypothesis, able to explain the equivalent widths and colors of the Contini et al. (1998) sample with the much more reasonable value $`\mathrm{\Omega }/(4\pi )=0.5`$, is the presence of evolved stars inside the starburst. Such a population, created during past formation episodes, can account for the observed weakness of equivalent widths by increasing the continuum (e.g. McCall et al. 1985; Díaz et al. 1991). It may also explain the red colors of the sample without requiring huge amounts of dust.

We considered two possible underlying populations. The first one consists of an instantaneous burst which occured 100 Myr ago. This old burst is assumed to have been ten times more intense than the current one. The predictions of these models are plotted on Fig. 9. As in Fig. 8, $`\tau _\lambda ^{\mathrm{stell}}/\tau _\lambda ^{\mathrm{neb}}`$ is set to 0.5. However, these models face two problems: first, it is still necessary to extend the age of the young burst to $`t=6`$ Myr, unless the covering factor is much lower than the value adopted on Fig. 9 (0.5); second, the sum of the evolved and the young population does not reproduce the colors. To reconcile this type of models with the data, it would be necessary to assume that the old burst was systematically $`100`$ times more intense than the current one.

The second type of underlying population is assumed to have formed continuously. We computed a series of scenarios matching the observations of standard Hubble sequence galaxies at $`z=0`$, and added a starburst. For the sake of consistency, the initial metallicity of the burst is the current metallicity $`Z(t)`$ of the host galaxy and is provided by PÉGASE.

The input parameters are chosen so that the starburst metallicities belong to the range $`[1/5Z_{},3/2Z_{}]`$ derived from line ratios. Each scenario is defined by three parameters: the age $`t_{\mathrm{gal}}`$ of the host galaxy when the starburst occurs; the infall timescale $`t_\mathrm{c}`$ used to parameterize the accretion rate of the galaxy; the gas-to-star conversion efficiency, $`\nu `$, relating the star formation rate $`\zeta (t)`$ to the gas density $`\sigma (t)`$ by the Schmidt law $`\zeta (t)=\nu \sigma (t)`$.

Five models (A, B, C, D, E) were computed for each scenario, corresponding to five different contributions of the evolved population to the total emission: A is the pure host galaxy, E is the pure starburst and B, C and D are intermediate cases. Note that the star formation does not stop after the burst, but returns to the quiescent regime. The fraction of ionizing photons due to the underlying population immediately after the burst is given in Table 5, as well as the values of $`t_{\mathrm{gal}}`$, $`t_\mathrm{c}`$, and $`\nu `$ for the considered scenario.

Fig. 10 presents the impact of the parameters on colors and equivalent widths. We compare the results of model D at $`Z=Z_{}/2`$ to the sample of Contini et al. (1998). We also show the results of model C for $`Z=1/2Z_{}`$, and the predictions of model D for $`3/2Z_{}`$. The filling factors are $`10^2`$ at $`1/2Z_{}`$ and $`5\times 10^4`$ at $`3/2Z_{}`$. As in Fig. 8, $`\tau _\lambda ^{\mathrm{stell}}`$/$`\tau _\lambda ^{\mathrm{neb}}`$ is set to 0.5. However, the covering factor (0.5) is significantly higher, and the $`\tau _V^{\mathrm{neb}}`$ range (0 to 4) is much more reasonable.

Another important point is that, contrary to the pure starburst case, ages of 6 Myr and more are compatible with emission line ratios. This is due to the presence of massive stars, formed after the burst in the underlying population. The presence of such stars ensures that the line ratios are virtually unchanged compared to their values during the burst. Contrary to the case of a burst that occured 100 Myr ago, these scenarios are able to explain simultaneously the colors, the equivalent widths and the line ratios in our sample. They do not require very high levels of extinction, and, maybe more important, do not imply a covering factor lower than 0.5. For all these reasons, the hypothesis of an underlying population formed continuously is our preferred one.

## 7 Discussion

### 7.1 A physical link between the ionization parameter and the metallicity?

Dopita & Evans (1986) proposed an analytical relation between the metallicity and the ionization parameter to describe the observational distribution of line ratios in their sample of extragalactic H ii regions. In Sect. 4.3, we saw that such a relation suffers from many degeneracies due to age, IMF and metallicity. Moreover, it seems that the quantities which are really related are the dispersion of $`U`$ and the metallicity.

What could be the physical link between $`Z`$ and $`U`$? According to one hypothesis, this is a link through geometrical effects. Following Castor et al. (1975) and assuming their equations are valid also for a star cluster, the inner radius $`R`$ of a gaseous shell surrounding a massive star or cluster is:

$$R\left(\frac{\dot{E}t^3}{\rho }\right)^{1/5},$$
(6)
where $`\dot{E}`$ is the energy deposition rate of a star (or cluster) of age $`t`$ in the interstellar medium (ISM) of density $`\rho `$. Leitherer et al. (1992) found that $`\dot{E}`$ is constant for an instantaneous burst until $``$ 6 Myr and scales linearly with $`Z`$, regardless of the age. Hence, high-$`Z`$ H ii regions expand faster.

We therefore propose the following scenario to explain the distribution of starbursts in the line ratio diagrams. For some reason, a site of intense star formation appears in a galaxy with a metal-poor interstellar medium. This starburst can be composed of many individual H ii regions surrounding star clusters. As these sub-components evolve with age, the H ii regions expand because of stellar winds. At the same time, the enrichment of the ISM by massive-star ejecta begins. If star formation goes on, very young components will coexist inside the starburst with more evolved components. The former will dominate the nebular emission during the first few Myr, and the latter after.

Two additional effects can explain the lack of high-$`Z`$ and high-$`U`$ objects: first, at high metallicity, individual H ii regions expand faster, leading to a lower ionization parameter; second, according to this scenario, the bulk of nebular emission in high-$`Z`$ starbursts is produced by evolved, low-excitation H ii regions. Checking the validity of this hypothesis requires to model also the expansion of nebulae. This is the subject of a forthcoming paper.

### 7.2 The equivalent-width problem

The problems encountered by pure starburst scenarios to account for the weakness of the H$`\alpha `$ equivalent widths from the sample of Contini et al. (1998) are several. The main one is that a very low covering factor is required. From Fig. 7 and 8, we deduce that a value $`\mathrm{\Omega }/(4\pi )0.1`$ is necessary to reproduce both $`W(\mathrm{H}\alpha )`$ and the line ratios in the case of a pure starburst. This means that 90 % of the Lyman continuum photons do not contribute to ionization, but are directly absorbed by dust or escape from the starburst. This is in contradiction with many previous studies using H$`\alpha `$ or H$`\beta `$ luminosities to estimate the mass of stellar clusters inside starbursts. In many cases, the agreement between such estimations and those based on ultraviolet continuum luminosities is excellent (González-Delgado et al. 1999). Some discrepancies appear frequently (Leitherer et al. 1996; Vacca et al. 1995), but they are however weak enough to exclude an absorption fraction of the ionizing photons by the gas as low as 10 %. Moreover, observations of diffuse H$`\alpha `$ emission in nearby galaxies lead to an estimation of 15–50 % for the fraction of ionizing photons escaping from individual H ii regions (Kennicutt 1998 and references therein).

Such a systematic low value of the covering factor is therefore clearly ruled out by the observations. Note also that pure starbursts models, even with $`\mathrm{\Omega }/(4\pi )=0.1`$, cannot reproduce simultaneously the emission line ratios and the equivalent widths (see Sect. 5). Hence, the presence of an evolved stellar population coexisting with a starburst, as analysed in Sect. 6, seems mandatory to explain equivalent widths and colors, even for the data of Contini et al. (1998), who strictly limitated their apertures to H$`\alpha `$-emitting areas. A similar conclusion for galaxies observed through 5-arcsec. apertures was drawn by Lançon & Rocca-Volmerange (1996) from the near-infrared spectral synthesis of starbursts.

### 7.3 Nature of the underlying population

Including underlying populations in the models allows to fit the data in a satisfactory way in the $`W(\mathrm{H}\alpha )`$ vs. $`BV`$ and $`VR`$ diagrams with $`\mathrm{\Omega }/(4\pi )=0.5`$, which means that $`50\%`$ of the ionizing photons escape from the gaseous nebula. This value is close to the 60 % fraction suggested by Mas-Hesse & Kunth (1999). The nature of this population remains however elusive.

We have tested an old (100 Myr old) instantaneous burst as a possible underlying population. The agreement with the data, at first sight, is relatively good (Fig. 9). In this case, however, the age of the youngest burst has to extend up to 8 Myr to explain the faint end of the W(H$`\alpha `$) observational distribution. At this age, our models are incompatible with line ratio data. Note that there is no Wolf-Rayet(W-R) star SED in our spectral library; instead, the spectra of the hottest stars available in the library are used during W-R stages. The implementation of W-R SEDs in the library could in principle extend the range of ages compatible with line ratios up to 6 Myr or more. As an example, González-Delgado et al. (1998) obtained an age between 6 and 9 Myr for IRAS 0833+6517. However, this conclusion is based on models with an IMF truncated at $`M_{\mathrm{up}}=30M_{}`$. Such a low value is clearly excluded for the bulk of our sample (Fig. 6). Moreover, the role of W-R stars in starburst nebular emission is still under debate. According to the most recent studies, the bulk of nebular emission is due to stellar populations younger than 3 Myr (Bresolin et al. 1999).

We have assumed that the old burst was only ten times more intense than the current one. Changing this value to one hundred, for example, would move the predictions on Fig. 9 toward the bottom-right and might explain the bulk of the data with ages lower than, or equal to, 4 Myr. In this case, however, such a difference between past and present star formation remains to be explained.

In this context, the models plotted on Fig. 10 appear as the most attractive ones. Scenarios with an underlying population providing from 2.5 up to 20 % of the number of ionizing photons are in acceptable agreement with the data and explain naturally all the observables (emission lines, colors and equivalent widths). The range of $`\tau _V`$ is consistent with the Balmer decrement measurements of Contini et al. (1998). Moreover, the age and metallicity of the underlying population are typically those of normal spirals, suggesting that starbursts are normal events in galaxy evolution. These scenarios draw a coherent picture of the starburst phenomenon, and for this reason, remain the most plausible ones.

## 8 Conclusion

We have analyzed the emission line and continuum properties of a large sample of starbursts and put new constraints on their stellar populations, metallicity, photoionization state and dust content. Our tool was the spectral evolution code PÉGASE coupled to the photoionization code CLOUDY.

From the comparison of a dataset of emission lines with pure instantaneous starburst scenarios, we have been able to constrain the evolution of the ionization parameter $`U`$ with the metallicity $`Z`$. At high abundance, the high-excitation component clearly detected at low $`Z`$ is absent or very weak. A preliminary analysis based on energy deposition rate considerations is unsufficient to explain this trend, and a more refined analysis of the dependency of the H ii region expansion with $`Z`$ is needed.

The second important result concerns underlying populations inside starbursts; they are detected even through small-aperture observations. This contribution is needed to reproduce simultaneously all the observables of our sample (line ratios, colors and equivalent widths). The underlying populations are simulated through the whole computation of the past star formation history of the galaxy hosting the starburst. In the scenarios favored by our analysis, the burst itself is supposed to be a brief, intense episode, preceded and followed by a quiescent regime. After the burst, this quiescent formation rate is sufficient to keep the line ratios compatible with the data.

Future prospects should consider the problem of metal depletion on dust grains, since depletion changes the composition of nebular gas. However, the question of the presence of dust *inside* H ii regions will be solved only with more constraining data, e.g. those of the ISO satellite. Finally, the large grids computed with the help of a coupled photoionization and stellar evolutionary synthesis model will be also used to study other samples with better statistics.

###### Acknowledgements.

We are pleased to thank Gary Ferland for his multiple consultation and help for the interface with CLOUDY. M. F. acknowledges support from the National Research Council through the Resident Research Associateship Program.


format: text

text:
# Closed-Time Path Integral Formalism and Medium Effects of Non-Equilibrium QCD Matter

## I introduction

Much effort has been given to detect a new form of matter known as quark-gluon plasma (QGP) for decades. The relativistic heavy-ion collider experiments at RHIC (Au-Au collisions at $`\sqrt{s}`$ = 200 GeV) and LHC (Pb-Pb collisions at $`\sqrt{s}`$ = 5.5 TeV) will provide the best opportunity to study such a state of matter in the laboratory. Perturbative QCD estimates that the energy density of produced jets and minijets might be larger than 50 and 1000 GeV/$`fm^3`$ at RHIC and LHC  which is much larger than the required energy density to produce quark-gluon plasma. As QGP lives for a very short time (several fermi) in a small volume ($``$ 100 fm3) a direct detection of this phase is not possible. Hence various indirect signatures are proposed for its detection. The prominent among them are: 1) $`J/\mathrm{\Psi }`$ suppression , 2) strangeness enhancement , dilepton and direct photon production . However, many uncertainties exist which make it difficult to claim the existence of the quark-gluon plasma. The main uncertainty lies in the lack of an accurate determination of the space-time evolution of the quarks and gluons produced just after the collision of two nuclei at RHIC and LHC. While the quark-gluon plasma during the equilibrium stage is described by Bjorken’s hydrodynamic evolution equations, it is much more difficult to determine the space-time evolution of partons in the pre-equilibrium stage. An accurate study of the pre-equilibrium stage will determine the equilibration time and initial conditions for hydrodynamic evolutions in the equilibrium stage. This study is also very important for the accurate determination of various signatures of the quark-gluon plasma.

To describe the space-time evolution of quark-gluon plasma at RHIC and LHC one needs to know how the partons are formed in these high energy nuclear collisions. The hard parton (jets and minijets) production can be calculated by using pQCD. The pre-equilibrium evolution of these hard partons can be studied by solving relativistic transport equations with secondary collisions among these partons taken into account . However, soft parton production can not be computed within pQCD formalism. There are coherent effects for the soft partons and they may be described by formation of a classical chromofield . For simplicity, we will only consider the evolution of minijet plasma which can be studied by solving relativistic transport equation:

$$p^\mu _\mu f(x,p)=C(x,p),$$
(1)
with secondary collision among the partons taken into account. In the above equation

$`C(x,p)={\displaystyle \frac{d^3p_2}{(2\pi )^3p_2^0}\frac{d^3p_3}{(2\pi )^3p_3^0}\frac{d^3p_4}{(2\pi )^3p_4^0}|M(pp_2p_3p_4)|^2\delta ^4(p+p_2p_3p_4)}`$ (2)
$`[f(x,p_3)f(x,p_4)(1+f(x,p))(1+f(x,p_2))f(x,p)f(x,p_2)(1+f(x,p_3))(1+f(x,p_4))]`$ (3)
is the collision term for a partonic scattering process $`pp_2p_3p_4`$. In the above expression $`p`$, $`p_2`$, and $`p_3`$, $`p_4`$ are the four momentum of the partons before and after the collision. Different partonic scattering processes which can be considered are like $`gggg`$, $`qqqq`$, $`q\overline{q}q\overline{q}`$ and $`gqgq`$ etc.. As the gluons are the dominant part of the minijet production we only consider $`gggg`$ process in the following. The squared matrix element for this process is given by:

$$|M(\widehat{s},\widehat{u},\widehat{t})|^2=\frac{9\pi ^2\alpha _s^2}{8}[3\frac{\widehat{u}\widehat{t}}{\widehat{s}^2}\frac{\widehat{u}\widehat{s}}{\widehat{t}^2}\frac{\widehat{s}\widehat{t}}{\widehat{u}^2}],$$
(4)
where $`\widehat{s}=(p+p_2)^2=(p_3+p_4)^2`$, $`\widehat{t}=(pp_3)^2=(p_2p_4)^2`$, $`\widehat{u}=(pp_4)^2=(p_2p_3)^2`$ are the Mandelstam variables. For massless gluon they are related by

$$\widehat{t}=\frac{\widehat{s}}{2}[1cos\theta ],\widehat{u}=\frac{\widehat{s}}{2}[1+cos\theta ],$$
(5)
where $`\theta `$ is the center of mass scattering angle which goes from $`0\frac{\pi }{2}`$ for identical partons in the final state. When one puts $`|M(\widehat{s},\widehat{u},\widehat{t})|^2`$ from Eq. (4) in Eq. (3) one encounters divergence in the collision term at small angle ($`\theta `$ 0) or small momentum transfer $`\widehat{t}`$ 0. This divergence is inevitable as long as we use the free propagator in vaccuum to evaulate the Feynman diagrams. However, this infrared divergence can be removed when one uses medium modified propagators instead of the vacuum propagator to evaluate the collision term in medium . The medium modified propagtaors have been obtained in the thermal field theory for the case of an equilibrium plasma. However, these finite temperature calculations are valid only in equilibrium where there is a static temperature such as in the case of a heat bath and the system is isotropic. For realistic situations in the high energy heavy-ion collisions at RHIC and LHC the partons formed at the initial time are in non-equilibrium and finite temperature QCD is not applicable at this stage. One has to use closed-time formalism to compute various quantities in non-equilibrium . In this paper, using closed-time path integral formalism, we derive the medium modified resummed gluon propagator in non-equilibrium to the one loop order in the covariant gauge which is necessary to obtain finite collision term to study equilibration of expanding minijet plasma at RHIC and LHC.

The paper is organized as follows. In section II we describe the formulation of closed-time path integral of SU(3) guage theory. In section III we derive the resummed gluon propagator to one loop level. We summarize and conclude our main results in section IV.

## II Closed-Time Path Integral Formalism in Gauge Theory

We consider SU(3) pure gauge theory which is QCD without quarks. In high energy heavy-ion collisions at RHIC and LHC most of the parton formed are gluons. Hence we concentrate on gluons only. However, extending close-time path integral formalism to quarks is straightforward. As the two nuclei travel almost at a speed of light at RHIC and LHC the system is dynamically evolving and many quantities have to be formulated in Boost invariant way . For this expanding system of partons we work in the covariant gauge in this paper. The QCD action without quark is given by:

$$S=d^4x[\frac{1}{4}F^{a\mu \nu }F_{\mu \nu }^a\frac{1}{2\xi }(G^a)^2_{FP}],$$
(6)
where the gluon field tensor:

$$F^{a\mu \nu }=^\mu A^{a\nu }^\nu A^{a\mu }+gf^{abc}A^{b\mu }A^{c\nu },$$
(7)
the ghost Lagrangian density:

$$_{FP}=[_\mu \overline{C}]D^\mu [A]C$$
(8)
and the gauge fixing term given by

$$G^a=_\mu A^{a\mu }.$$
(9)
In Eq. (6) $`\xi `$ is the gauge fixing parameter in covariant gauge. In the closed-time path integral formalism the gauge field $`A`$, the ghost field $`C`$ and the corresponding sources $`j,\chi `$ are defined in both the time branches. To make the formulas simpler we denote the fields and the corresponding sources by following common notations:

$$Q=(A,C,\overline{C}),J=(j,\overline{\chi },\chi ).$$
(10)
Denoting $`Q^+`$, $`Q^{}`$ and $`J^+`$,$`J^{}`$ the fields and the sources on the upper and lower branch of the time path, the in-in generating functional becomes:

$$Z[J_+,J_{},\rho ]=DQ^r<Q^+,t_0|\rho |Q^{},t_0>e^{i(S[Q^r]+J_rQ^r)},$$
(11)
where $`S[Q^r]=S[Q^+]S^{}[Q^{}]`$ with $`r=+,`$. More explicitly, the above equation can be written as:

$`Z[J_+,J_{},\rho ]={\displaystyle [dQ^+]d[Q^{}]<Q^+,t_0|\rho |Q^{},t_0>e^{i[S_0[Q]+S_{int}[𝐐]+Tr𝐉𝐐]}},`$ (12)
where

$`S_0[Q]={\displaystyle \underset{r,s=+,}{}}{\displaystyle \frac{d^4p}{(2\pi )^4}[\frac{1}{2}A_\mu ^r(p)[G_{\mu \nu }^1(p)]^{rs}A_\nu ^s(p)+\overline{C}^r(p)[S^1(p)]^{rs}C^s(p)]},`$ (13)
and

$$Tr𝐉𝐐=\frac{d^4p}{(2\pi )^4}[j_\mu ^+(p)A_\mu ^+(p)+\overline{C}^+(p)\chi ^+(p)+\overline{\chi }^+(p)C^+(p)+(+)].$$
(14)
In the Eq. (13) $`G_{\mu \nu }(p)`$ and $`S(p)`$ are gluon and ghost free propagators and in Eq. (12) $`S_{int}[𝐐]=S_{int}[Q^+]S_{int}^{}[Q^{}]`$.

In covariant gauge the further complications arise because of the presence of the unphysical ghost fields. In equilibrium one can define the ghost distribution function (BE) and hence can deal with ghost fields in the medium even in covariant gauge. However, in non-equilibrium there is no easy procedure to obtain a ghost distribution function in the QCD medium at RHIC and LHC. Quark and gluon distribution functions in non-equilibrium situations at RHIC and LHC can be obtained from minijets by using pQCD or via other methods . For this reason we work in the Landshoff and Rebhan scheme of frozen ghost formalism , where the gauge theory in the covariant gauge is obtained by restricting the space of initial state $`|Q^r,t_0>`$ to the physical one. This means gluons with spatially transverse polarization will contribute to the trace. Hence the non-local Kernel $`K`$ appearing in the path integral of the generating function  will couple only to the transverse component of the gauge field $`A`$.

Now we consider a cylindrically symmetric expanding system of partons in 1$``$1 dimension. For this purpose we introduce the flow velocity of the medium

$$u^\mu =(\mathrm{cosh}\eta ,0,0,\mathrm{sinh}\eta ),$$
(15)
where $`\eta =\frac{1}{2}\mathrm{ln}\frac{t+z}{tz}`$ is the space-time rapidity and $`u_\mu u^\mu =1`$. We define the four symmetric tensors :

$`T_{\mu \nu }(p)`$ $`=`$ $`g_{\mu \nu }{\displaystyle \frac{(up)(u_\mu p_\nu +u_\nu p_\mu )p_\mu p_\nu p^2u_\mu u_\nu }{(up)^2p^2}},`$ (16)
$`L_{\mu \nu }(p)`$ $`=`$ $`{\displaystyle \frac{p^2}{(up)^2p^2}}\left(u_\mu {\displaystyle \frac{(up)p_\mu }{p^2}}\right)\left(u_\nu {\displaystyle \frac{(up)p_\nu }{p^2}}\right),`$ (17)
$`C_{\mu \nu }(p)`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2[(up)^2p^2]}}}\left[\left(u_\mu {\displaystyle \frac{(up)p_\mu }{p^2}}\right)p_\nu +\left(u_\nu {\displaystyle \frac{(up)p_\nu }{p^2}}\right)p_\mu \right]\mathrm{and}`$ (18)
$`D_{\mu \nu }(p)`$ $`=`$ $`{\displaystyle \frac{p_\mu p_\nu }{p^2}},`$ (19)
which are required for a cylindrically symmetric system in 1$``$1 expanding plasma at RHIC and LHC in the very early stage.

Here $`T^{\mu \nu }`$ is transverse with respect to the flow-velocity but $`L^{\mu \nu }`$ and $`D^{\mu \nu }`$ are mixtures of space-like and time-like components. These tensors satisfy the following transversality properties with respect to $`p^\mu `$:

$$p_\mu T^{\mu \nu }(p)=p_\mu L^{\mu \nu }(p)=0,p_\mu p_\nu C^{\mu \nu }(p)=0.$$
(20)
In addition to this, the above tensors satisfy the following properties:

$$\begin{array}{c}TL=TC=TD=0,T+L+D=1,\\ TT=T,LL=L,CC=\frac{1}{2}(L+D),\\ TrCL=TrCD=0.\end{array}$$
(21)
Any symmetric tensor $`S^{\mu \nu }`$ can be written in terms of the above four tensors:

$$S^{\mu \nu }=aT^{\mu \nu }+bL^{\mu \nu }+cC^{\mu \nu }+dD^{\mu \nu }$$
(22)
with

$$a=\frac{1}{2}TrTS,b=TrLS,c=TrCS,d=TrDS.$$
(23)
In terms of this tensor basis the gluon propagator in the covariant gauge is given by:

$`G_{\mu \nu }(p)_{ij}`$ $`=`$ $`iT_{\mu \nu }(p)\left(\left[G(p)\right]_{ij}^{vac}+\left[G(p)\right]_{ij}^{med}\right)i(L_{\mu \nu }(p)+\xi D_{\mu \nu }(p))\left[G(p)\right]_{ij}^{vac}`$ (24)
$`=`$ $`i\left(g_{\mu \nu }+(\xi 1)D_{\mu \nu }(p)\right)\left[G(p)\right]_{ij}^{vac}iT_{\mu \nu }(p)\left[G(p)\right]_{ij}^{med},`$ (25)
where $`i,j=+,`$. The forms of $`G(p)_{ij}`$ are

$$\left[G(p)\right]_{ij}^{vacuum}=\left(\begin{array}{cc}\frac{1}{p^2+iϵ}& 0\\ 0& \frac{1}{p^2iϵ}\end{array}\right)2\pi i\delta (p^2)\left(\begin{array}{cc}0& \theta (p_0)\\ \theta (p_0)& 0\end{array}\right),$$
(26)
and

$$\left[G(p)\right]_{ij}^{medium}=2\pi i\delta (p^2)\stackrel{~}{f}(\stackrel{}{p})\left(\begin{array}{cc}1& 1\\ 1& 1\end{array}\right),\stackrel{~}{f}(\stackrel{}{p})=f(\stackrel{}{p})\theta (p_0)+f(\stackrel{}{p})\theta (p_0).$$
(27)
Here $`f(\stackrel{}{p})`$ is the distribution function. Note that in this Landshoff-Rehban scheme the transverse component of the gauge propagator proportional to $`T^{\mu \nu }`$ contains the medium effect which is gauge-parameter independent. As the initial density matrix contains only sum over transverse polarizations of gluons the ghost propagator is simply given by:

$$S(p)_{ij}=i\left[G(p)\right]_{ij}^{vac},$$
(28)
because ghost fields do not couple to the kernel $`K`$ in the generating function.

## III Resummed Gluon Propagator at One Loop Level in Non-Equilibrium

In the last section we derived the gauge field propagator in the medium for the free part of the action $`S_0`$. In this section we will consider the interaction term of the action $`S_{int}`$. For such a situation one will have to consider the full propagator instead of the free propagator in the medium. This implies one has to solve the full Schwinger-Dyson equation which is practically impossible. Usually one has to truncate the series at one or two loop level of the self-energy. Such a truncation demands that the QCD coupling constant should not be too large. We consider here the evolution of the gluon-minijet plasma at RHIC and LHC where the average transverse momentum of the partons are found to be large at the earlry stage . Therefore the coupling constant corresponding to such average transverse momentum is found to be small. In this situation we truncate the Schwinger-Dyson equation at one loop level of the self energy and consider the corresponding resummed two-point Green’s function $`\stackrel{~}{G}`$.

The resummed two-point Green’s function $`\stackrel{~}{G}`$ of gluon can be decomposed as:

$$\stackrel{~}{G}_{\mu \nu }(p)=iT_{\mu \nu }(p)\stackrel{~}{G}^T(p)iL_{\mu \nu }(p)\stackrel{~}{G}^L(p)i\xi D_{\mu \nu }(p)\stackrel{~}{G}^D(p),$$
(29)
where $`\stackrel{~}{G}^T`$, $`\stackrel{~}{G}^L`$, $`\stackrel{~}{G}^D`$ correspond to $`T`$, $`L`$ and $`D`$ components respectively. The last part $`\stackrel{~}{G}^D(p)`$ is identical to the vacuum part  and hence we do not consider it any more. It can be mentioned that there are separate Dyson-Schwinger equation for different components which does not couple with each other.

Dyson-Schwinger equations for different components can be written as

$$\left[\stackrel{~}{G}^{T,L}(p)\right]_{ij}=\left[G^{T,L}(p)\right]_{ij}+\underset{l,k}{}\left[G^{T,L}(p)\right]_{il}\left[\mathrm{\Pi }^{T,L}(p)\right]_{lk}\left[\stackrel{~}{G}^{T,L}(p)\right]_{kj}.$$
(30)
Please remember that $`i,j,k,l=+,`$ and suppression of Lorentz and color indices in the above equation is understood.

Instead of using the matrix form of this equation, we prefer to use the retarded, advanced and symmetric Green’s functions. The retarded and advanced resummed Green’s function are found to be:

$$\stackrel{~}{G}_{R,A}^{T,L}(p)=G_{R,A}^{T,L}(p)+G_{R,A}^{T,L}(p)\mathrm{\Pi }_{R,A}^{T,L}(p)\stackrel{~}{G}_{R,A}^{T,L}(p).$$
(31)
The straightforward solution of the above equation gives:

$$\stackrel{~}{G}_{R,A}^{T,L}(p)=\frac{G_{R,A}^{T,L}(p)}{1G_{R,A}^{T,L}(p)\mathrm{\Pi }_{R,A}^{T,L}(p)}=\frac{1}{p^2\mathrm{\Pi }_{R,A}^{T,L}(p)\pm isgn(p_0)ϵ},$$
(32)
where the self-energy contains the medium effects. Similar but more complicated equation is obtained for the resummed symmetric Green’s function:

$$\begin{array}{c}\stackrel{~}{G}_S^{T,L}(p)=G_S^{T,L}(p)+G_R^{T,L}(p)\mathrm{\Pi }_R^{T,L}(p)\stackrel{~}{G}_S^{T,L}(p)\\ +G_S^{T,L}(p)\mathrm{\Pi }_A^{T,L}(p)\stackrel{~}{G}_A^{T,L}(p)+G_R^{T,L}(p)\mathrm{\Pi }_S^{T,L}(p)\stackrel{~}{G}_A^{T,L}(p).\end{array}$$
(33)
After some algebra it can be shown that

$$\begin{array}{c}\stackrel{~}{G}_S^{T,L}(p)=\left[1+2f(\stackrel{}{p})\right]sgn(p_0)[\stackrel{~}{G}_R^{T,L}(p)\stackrel{~}{G}_A^{T,L}(p)]\\ +\left(\mathrm{\Pi }_S^{T,L}(p)(1+2f(\stackrel{}{p}))sgn(p_0)[\mathrm{\Pi }_R^{T,L}(p)\mathrm{\Pi }_A^{T,L}(p)]\right)\times \stackrel{~}{G}_R^{T,L}(p)\times \stackrel{~}{G}_A^{T,L}(p).\end{array}$$
(34)
Let us evaluate the various components of the self energy. For short hand notation we denote Gl for gluon loop, Ta for tapole diagram (see Fig. 1) and FP represents ghost loop diagram. Explicitly the expressions for various self energies are given by:

$`\left[\mathrm{\Pi }_{Gl;\mu \nu }^{ab,}(p)\right]_{kl}`$ $`=`$ $`(k\times l){\displaystyle \frac{g^2}{2i}}\delta _{ab}N_c{\displaystyle \frac{d^4q}{(2\pi )^4}\left[G^{\alpha \gamma }(q)\right]_{kl}\left[G^{\beta \delta }(pq)\right]_{kl}}`$ (35)
$``$ $`\left[g_{\mu \alpha }(pq)_\beta +g_{\alpha \beta }(2qp)_\mu +g_{\beta \mu }(2pq)_\alpha \right]`$ (36)
$``$ $`\left[g_{\nu \gamma }(p+q)_\delta +g_{\gamma \delta }(p2q)_\nu +g_{\nu \delta }(q2p)_\gamma \right].`$ (37)
$`\left[\mathrm{\Pi }_{Ta;\mu \nu }^{ab}(p)\right]_{kl}`$ $`=`$ $`g^2N_c\delta _{ab}\delta _{kl}{\displaystyle \frac{d^4q}{(2\pi )^4}(g_{\mu \lambda }g_{\nu \sigma }g_{\sigma \lambda }g_{\mu \nu })G^{\sigma \lambda }(q)_{kl}}.`$ (38)
$`\left[\mathrm{\Pi }_{FP;\mu \nu }^{ab}(p)\right]_{kl}`$ $`=`$ $`(k\times l)ig^2\delta _{ab}N_c{\displaystyle \frac{d^4q}{(2\pi )^4}\left[G^{vac}(q)\right]_{kl}\left[G^{vac}(pq)\right]_{kl}q_\mu (qp)_\nu }.`$ (39)
It has to be remembered that $`k,l(=(+,))`$ are not contracted in the right hand side of the above equations. So it is obvious that $`\mathrm{\Pi }_{Ta;\mu \nu }^{ab}(p)_{kl}=0`$, when $`kl`$. In the following calculation we neglect the vacuum part and concentrate on the medium part of the self energy. Note that the distribution function is contained only in the medium part of the self energy. The divergence of the vacuum part of the self energy is absorbed in the redefinition of the bare quantities which is well known. Like the Green’s function the self energy can be written as :

$$\mathrm{\Pi }_{\mu \nu }(p)=T_{\mu \nu }(p)\mathrm{\Pi }^T(p)+L_{\mu \nu }(p)\mathrm{\Pi }^L(p).$$
(40)
Using Eq. (19) $`\mathrm{\Pi }`$’s are decomposed as:

$`\mathrm{\Pi }^{a;T}(p)`$ $`=`$ $`\mathrm{\Pi }_{\mu \nu }^a(p){\displaystyle \frac{T^{\mu \nu }(p)}{2}},a=R,A,S,\mathrm{and}`$ (41)
$`\mathrm{\Pi }^{a;L}(p)`$ $`=`$ $`\mathrm{\Pi }_{\mu \nu }^a(p)L^{\mu \nu }(p),a=R,A,S.`$ (42)
Total self energy of the medium part is the sum of gluon loop and tadpole contributions (see Fig. 1) because the ghost loop contribution is present in the vacuum sector in the frozen ghost formalism. The real and imaginary part of the total self energy can be written as:

$$Re\mathrm{\Pi }^{T,L}(p)_a=Re\mathrm{\Pi }^{T,L}(p)_{Gl;a}+\mathrm{\Pi }^{T,L}(p)_{Ta;a},a=R,A,S,$$

$$Im\mathrm{\Pi }^{T,L}(p)_a=Im\mathrm{\Pi }^{T,L}(p)_{Gl;a},a=R,A,S.$$

One can check that the retarded and advanced self energy are related by:

$`Re\mathrm{\Pi }_{Gl;R}^{T,L}(p)=Re\mathrm{\Pi }_{Gl;A}^{T,L}(p);Im\mathrm{\Pi }_{Gl;R}^{T,L}(p)=Im\mathrm{\Pi }_{Gl;A}^{T,L}(p),`$ (43)
$`Re\mathrm{\Pi }_{Gl;R}^{T,L}(p)=Re\mathrm{\Pi }_{Gl;A}^{T,L}(p);Im\mathrm{\Pi }_{Gl;R}^{T,L}(p)=Im\mathrm{\Pi }_{Gl;A}^{T,L}(p)\mathrm{and}`$ (44)
$`\mathrm{\Pi }_{Ta;R}^{T,L}(p)=\mathrm{\Pi }_{Ta;A}^{T,L}(p).`$ (45)
Furthermore we have $`\mathrm{\Pi }_{Ta;S}^{T,L}(p)`$=0. Therefore we only need calculate $`Re\mathrm{\Pi }_{Gl;R}^{T,L}`$, $`Im\mathrm{\Pi }_{Gl;R}^{T,L}`$, $`\mathrm{\Pi }_{Gl;S}^{T,L}`$ and $`\mathrm{\Pi }_{Ta;R}`$ which we compute in this paper for any non-equilibrium gluon distribution function $`f(p)`$. Simplifying Eq. (39) and using Eq. (42) the explicit expression of retarded self energy of the gluon loop diagram is found to be:

$`\mathrm{\Pi }_{Gl;R}^{T,L}(p)`$ $`=`$ $`{\displaystyle \frac{g^2}{2}}\delta _{ab}N_c{\displaystyle }{\displaystyle \frac{d^4q}{(2\pi )^3}}[{\displaystyle \frac{\stackrel{~}{f}(\stackrel{}{q})\delta (q^2)H^{T,L}(q,p)}{(pq)^2+isgn(p_0q_0)ϵ}}`$ (46)
$`+`$ $`{\displaystyle \frac{\stackrel{~}{f}(\stackrel{}{p}\stackrel{}{q})\delta ((pq)^2)H^{T,L}(pq,p)}{q^2+isgn(q_0)ϵ}}],`$ (47)
where

$`H^T(q,p)`$ $`=`$ $`8{\displaystyle \frac{(uq)(qp)(up)}{(up)^2p^2}}4{\displaystyle \frac{(qp)^2}{(up)^2p^2}}4{\displaystyle \frac{p^2(qu)^2}{(up)^2p^2}}`$ (48)
$``$ $`\left[(p+q)^2\right]\left[1{\displaystyle \frac{(qp)(up)}{(uq)((up)^2p^2)}}+{\displaystyle \frac{(qp)^2}{2(uq)^2((up)^2p^2)}}+{\displaystyle \frac{p^2}{2((up)^2p^2)}}\right]`$ (49)
$``$ $`4p^2+8{\displaystyle \frac{(qp)(up)}{(uq)}}4{\displaystyle \frac{(qp)^2}{(uq)^2}}`$ (50)
$`+`$ $`{\displaystyle \frac{(\xi 1)(p^4((qp)^2+2(qu)(pu)(qp)+(qu)^2(p^22(pu)^2)))}{2(qu)^2((qp)p^2)(p^2(pu)^2)}}`$ (51)
and

$`H^L(q,p)`$ $`=`$ $`{\displaystyle \frac{8p^2}{(up)^2p^2}}[((uq){\displaystyle \frac{(up)(qp)}{p^2}})^2]`$ (52)
$``$ $`\left[(p+q)^2\right]\left[{\displaystyle \frac{2(qp)(pu)^3}{(qu)p^2((pu)^2p^2)}}{\displaystyle \frac{(qp)^2(pu)^2}{(qu)^2p^2((pu)^2p^2)}}{\displaystyle \frac{(pu)^2}{((pu)^2p^2)}}\right]`$ (53)
$``$ $`4p^2+8{\displaystyle \frac{(qp)(up)}{(uq)}}4{\displaystyle \frac{(qp)^2}{(uq)^2}}`$ (54)
$`+`$ $`(\xi 1){\displaystyle \frac{(pu)^2((qp)^22(qu)(pu)(qp)+(qu)^2p^2)}{(qu)^2(p^2(pu)^2)}}.`$ (55)
Using the $`\delta `$ function we get:

$`\mathrm{\Pi }_{Gl;R}^{T,L}(p)`$ $`=`$ $`{\displaystyle \frac{g^2}{2}}\delta _{ab}N_c{\displaystyle }{\displaystyle \frac{d^3q}{(2\pi )^3}}{\displaystyle \frac{1}{2|\stackrel{}{q}|}}[{\displaystyle \frac{f(\stackrel{}{q})H^{T,L}(q,p)|_{q_0=|\stackrel{}{q}|}}{(p_0|\stackrel{}{q}|)^2|\stackrel{}{q}\stackrel{}{p}|^2+isgn(p_0|\stackrel{}{q}|)ϵ}}`$ (56)
$`+`$ $`{\displaystyle \frac{f(\stackrel{}{q})H^{T,L}(q,p)|_{q_0=|\stackrel{}{q}|}}{(p_0+|\stackrel{}{q}|)^2|\stackrel{}{q}\stackrel{}{p}|^2+isgn(p_0+|\stackrel{}{q}|)ϵ}}]`$ (57)
$`+`$ $`(qpq).`$ (58)
Similar but more complicated calculation for symmetric self energy yields:

$`\mathrm{\Pi }_{Gl;S}^{(T,L)}(p)`$ $`=`$ $`{\displaystyle \frac{g^2}{2i}}\delta _{ab}N_c{\displaystyle \frac{d^3q}{(2\pi )^2}\delta \left((|\stackrel{}{q}|p_0)^2|\stackrel{}{q}\stackrel{}{p}|^2\right)\frac{K^{(T,L)}(q,p)|_{q_0=|\stackrel{}{q}|}}{2|\stackrel{}{q}|}f(\stackrel{}{q})\stackrel{~}{f}(\stackrel{}{p}\stackrel{}{q})}`$ (59)
$`+`$ $`\delta \left((|\stackrel{}{q}|+p_0)^2|\stackrel{}{q}\stackrel{}{p}|^2\right){\displaystyle \frac{K^{(T,L)}(q,p)|_{q_0=|\stackrel{}{q}|}}{2|\stackrel{}{q}|}}f(\stackrel{}{q})\stackrel{~}{f}(\stackrel{}{p}\stackrel{}{q}),`$ (60)
where $`K^{T,L}(q,p)`$ are given by:

$`K^T(q,p)={\displaystyle \frac{p^2(p^2+4(qu)((qu)(pu)))}{4(qu)^2(qupu)^2(p^2(pu)^2)}}[8(qu)^416(pu)(qu)^3`$ (61)
$``$ $`12(p^22(pu)^2)(qu)^2+4(3p^2(pu)4(pu)^3)(qu)+8(pu)^4+p^48p^2(pu)^2]`$ (62)
and

$`K^L(q,p)`$ $`=`$ $`{\displaystyle \frac{p^2[(pu)2(qu)]^2}{4(qu)^2(qupu)^2(p^2(pu)^2)}}[8(qu)^416(pu)(qu)^3`$ (63)
$`+`$ $`4(2(pu)^2+p^2)(qu)^24p^2(pu)(qu)+p^4].`$ (64)
The expressions for the transverse and longitudinal part of the self energy coming from the Tadpole part (Eq. (39) and (42)) are found to be:

$`\mathrm{\Pi }_{Ta;R}^T(p)=g^2\delta _{ab}N_c{\displaystyle }{\displaystyle \frac{d^3q}{(2\pi )^3}}{\displaystyle \frac{f(\stackrel{}{q})}{2|\stackrel{}{q}|}}[1+{\displaystyle \frac{(up)(qp)}{2(uq)((up)^2p^2)}}{\displaystyle \frac{(qp)^2}{2(uq)^2((up)^2p^2)}}`$ (65)
$``$ $`{\displaystyle \frac{p^2}{2((up)^2p^2)}}]|_{q_0=|\stackrel{}{q}|}+{\displaystyle \frac{f(\stackrel{}{q})}{2|\stackrel{}{q}|}}[1+{\displaystyle \frac{(up)(qp)}{2(uq)((up)^2p^2)}}{\displaystyle \frac{(qp)^2}{2(uq)^2((up)^2p^2)}}`$ (66)
$``$ $`{\displaystyle \frac{p^2}{2((up)^2p^2)}}]|_{q_0=|\stackrel{}{q}|}`$ (67)
and

$`\mathrm{\Pi }_{Ta;R}^L(p)=g^2\delta _{ab}N_c{\displaystyle \frac{d^3q}{(2\pi )^3}\frac{f(\stackrel{}{q})}{2|\stackrel{}{q}|}[32\frac{(up)(qp)}{p^2(uq)}+\frac{p^2}{(up)^2p^2}(1\frac{(up)(qp)}{p^2(uq)})^2]}|_{q_0=|\stackrel{}{q}|}`$ (68)
$`+`$ $`{\displaystyle \frac{f(\stackrel{}{q})}{2|\stackrel{}{q}|}}[32{\displaystyle \frac{(up)(qp)}{p^2(uq)}}+{\displaystyle \frac{p^2}{(up)^2p^2}}(1{\displaystyle \frac{(up)(qp)}{p^2(uq)}})^2]|_{q_0=|\stackrel{}{q}|}`$ (69)
respectively. To summarize Eqs. (58), (60), (67) and (69) contain all the expressions for different parts of the self energy in non-equilibrium which will be used in the medium modified propagator to obtain a finite collision term. In this paper we make hard momentum loop approximation in non-equilibrium which is simillar to hard thermal loop approximation in equilibrium . In the hard momentum loop approximation the self energies are found to be gauge parameter independent (see Appendix).

To obtain the collision term for the process $`gggg`$ in the medium we have to use the resummed Feynman propagator $`[\stackrel{~}{G}(p)_{++}]`$ which can be obtained from the resummed advanced, retarded and symmetric propagators via the relation:

$$\left[\stackrel{~}{G}(p)\right]_{++}=\frac{1}{2}\left[\stackrel{~}{G}_S(p)+\stackrel{~}{G}_A(p)+\stackrel{~}{G}_R(p)\right],$$
(70)
where

$$\stackrel{~}{G}_{R,A}(p)=\frac{1}{p^2\mathrm{\Pi }_{R,A}(p)\pm isgn(p_0)ϵ},$$
(71)
and

$`\stackrel{~}{G}_S(p)`$ $`=`$ $`\left[1+2f(\stackrel{}{p})\right]sgn(p_0)[\stackrel{~}{G}_R(p)\stackrel{~}{G}_A(p)]`$ (72)
$`+`$ $`\left(\mathrm{\Pi }_S(p)(1+2f(\stackrel{}{p}))sgn(p_0)[\mathrm{\Pi }_R(p)\mathrm{\Pi }_A(p)]\right)\times \stackrel{~}{G}_R(p)\times \stackrel{~}{G}_A(p).`$ (73)
Different parts of the self energy appearing in the above equations are given in Eqs. (58), (60), (67) and (69). Finally, using the relation between various parts of the self energy (see Eq. (45)) we find from the above equations:

$$\stackrel{~}{G}_{++}(p)=\frac{p^2Re\mathrm{\Pi }_R(p)+\frac{1}{2}\mathrm{\Pi }_S(p)}{(p^2Re\mathrm{\Pi }_R(p))^2+(Im\mathrm{\Pi }_R(p))^2}$$
(74)
which is the required expression for the medium modified resummed Feynman gluon propagator in non-equilibrium at one loop level of the self energy.

Our main purpose is to remove infrared divergence appearing in the small angle partonic scatterings which plays an important role in the production and equilibration of minijet plasma. For this purpose it is necessary to study the infrared behaviour of the self energies obtained in this paper for non-equilibrium situtations. We note that in the static limit ($`p_0=0,|\stackrel{}{p}|0`$) one obtains $`\mathrm{\Pi }^L(p_0=0,|\stackrel{}{p}|0)=m_D^2`$ (the Debye screening mass) and $`\mathrm{\Pi }^T(p_0=0,|\stackrel{}{p}|0)=m_g^2`$ (the magnetic screening mass). It is widely believed that while Debye screening mass is non-zero the magnetic screening mass is zero at one loop level of the self energy. For this reason one expects that while the electric field is screened, the magnetic field is not screened and one still has infrared divergence in the magnetic sector at one loop level. However, this is true for a system where the gluon distribution function is isotropic (in momentum space) or in equilibrium. This is not true for non-isotropic gluon distribution function which is the case in the early stage of the heavy-ion collisions at RHIC and LHC. The static limit result of the transverse part of the self energy derived in this paper is not zero for non-isotropic (in momentum space) gluon distribution function. The typical values of the magnetic screening masses at one loop level is found to be 257 MeV at RHIC and 330 MeV at LHC by using non-equilibrium gluon-minijet distribution function  at the initial time. Hence the medium modified gluon propagator derived in this paper is safe from infrared divergences both in electric and magnetic sector in non-equilibrium situations at RHIC and LHC. Therefore these propagators can be used to obtain finite collision terms to study equilibration of minijet plasma at RHIC and LHC. Note that the above values are obtained by using the minijet distribution function at the initial time. The time evolution of these screening masses have to be determined by solving relativistic transport equations with Bjorken’s boost invariance picture taken into account. As these involves extensive numerical work (see ) we will report the self-consistent space-time evolution study elsewhere.

## IV Conclusions

In this paper we have applied the closed-time path integral formalism to study medium effects of the minijet plasma in non-equilibrium. In particular we have derived the medium modified resummed gluon propagator to the one loop order of self energy which is necessary to obtain a finite collision term to study equilibration of parton plasma at RHIC and LHC. These medium modified propagators have been studied in more detail in finite temperature QCD formulations. However, at RHIC and LHC the parton momentum distributions at early stage are anisotropic and finite temperature QCD formulations can not be applied to these non-equilibrium situations. This is because a parton inside the nucleus (which travels almost at the speed of light at RHIC and LHC) carries mostly longitudinal momentum before an ultra relativistic nuclear collisions. After jets and minijets are formed and suffer secondary collisions the isotropy between longitudinal and transverse momentum may be achieved. According to Bjorken’s proposal  many quantities are expected to be expressed in terms of boost invariant parameters. For this reason and to have a covariant formulation we have worked in the covariant gauge which is a suitable gauge for expanding plasma. We give the result of the resummed gluon propagator upto one loop level of self energy in non-equilibrium. Furthermore it is shown that these propagators are infrared divergence free, both in electric and magnetic sector.

The medium modified resummed gluon propagator we derived in this paper will be used to obtain finite collision term for the $`gggg`$ scattering process in non-equilibrium to solve the relativistic kinetic equation (Eq. (1)) to study equilibration of minijet plasma at RHIC and LHC. In this way one does not have to put ad-hoc values for the momentum transfer cut-off which crucially changes all the properties and hence determination of all the signatures of the quark-gluon plasma. In future we hope to use our results of the resummed gluon propagators in the collision term to study production and equilibration of the minijet plasma at RHIC and LHC. As the solution of the relativistic transport equation involves extensive numerical work (see ) we will report it elsewhere.

###### Acknowledgements.

We thank Fred Cooper for useful discussions. We thank Joerg Ruppert for his help in using Feynman package. C-W. K. and G. C. N. acknowledge the financial support from Alexander von Humboldt Foundation. Appendix A: The explicit forms of self energies in Hard Loop Momentum Approximation:
We use hard momentum loop (HML) approximation in non-equilibrium which is equivalent to hard thermal loop (HTL) approximation in equilibrium. In the hard momentum loop approximation the loop momentum $`q`$ is harder than external momentum $`p`$ . In HML approximation we find that the whole expressions are independent of gauge-fix parameters $`\xi `$. In the hard momentum loop approximation we find:    $`Re\mathrm{\Pi }_{Gl;R}^T(p)`$ $``$ $`{\displaystyle \frac{g^2}{2}}\delta _{ab}N_c{\displaystyle \frac{1}{(up)^2p^2}}{\displaystyle }{\displaystyle \frac{d^3q}{(2\pi )^3}}{\displaystyle \frac{1}{2|\stackrel{}{q}|}}[{\displaystyle \frac{f(\stackrel{}{q})}{(|\stackrel{}{q}|p_0)^2|\stackrel{}{q}\stackrel{}{p}|^2}}`$ (75)
$``$ $`[8(uq)(qp)(up)4(qp)^24p^2(qu)^2`$ (76)
$``$ $`2(pq)(up)^2+p^2(pq)+2{\displaystyle \frac{(pq)^2(up)}{uq}}2{\displaystyle \frac{(qp)^3}{(uq)^2}}]|_{q_0=|\stackrel{}{q}|}`$ (77)
$`+`$ $`{\displaystyle \frac{f(\stackrel{}{q})}{(|\stackrel{}{q}|+p_0)^2|\stackrel{}{q}\stackrel{}{p}|^2}}[8(uq)(qp)(up)4(qp)^24p^2(qu)^2`$ (78)
$``$ $`2(pq)(up)^2+p^2(pq)+2{\displaystyle \frac{(pq)^2(up)}{uq}}2{\displaystyle \frac{(qp)^3}{(uq)^2}}]|_{q_0=|\stackrel{}{q}|}]`$ (79)
$`+`$ $`(qpq),`$ (80)
and

$`Im\mathrm{\Pi }_{Gl;R}^T(p)`$ $``$ $`{\displaystyle \frac{g^2}{2}}\delta _{ab}N_c{\displaystyle \frac{1}{(up)^2p^2}}{\displaystyle \frac{d^3q}{(2\pi )^2}sgn(p_0|\stackrel{}{q}|)\delta \left((p_0|\stackrel{}{q}|)^2|\stackrel{}{q}\stackrel{}{p}|^2\right)\frac{f(\stackrel{}{q})}{2|\stackrel{}{q}|}}`$ (82)
$`[8(uq)(qp)(up)4(qp)^24p^2(qu)^2`$
$``$ $`2(pq)(up)^2+p^2(pq)+2{\displaystyle \frac{(pq)^2(up)}{uq}}2{\displaystyle \frac{(qp)^3}{(uq)^2}}]|_{q_0=|\stackrel{}{q}|}`$ (83)
$`+`$ $`sgn(p_0+|\stackrel{}{q}|)\delta \left((|\stackrel{}{q}|+p_0)^2|\stackrel{}{q}\stackrel{}{p}|^2\right){\displaystyle \frac{f(\stackrel{}{q})}{2|\stackrel{}{q}|}}`$ (85)
$`[8(uq)(qp)(up)4(qp)^24p^2(qu)^2`$
$``$ $`2(pq)(up)^2+p^2(pq)+2{\displaystyle \frac{(pq)^2(up)}{uq}}2{\displaystyle \frac{(qp)^3}{(uq)^2}}]|_{q_0=|\stackrel{}{q}|}`$ (86)
$`+`$ $`(qpq),`$ (87)
Similarly for longitudinal components we obtain:

$`Re\mathrm{\Pi }_{Gl;R}^L(p)`$ $``$ $`{\displaystyle \frac{g^2}{2}}\delta _{ab}N_c{\displaystyle \frac{1}{(up)^2p^2}}[{\displaystyle }{\displaystyle \frac{d^3q}{(2\pi )^3}}{\displaystyle \frac{f(\stackrel{}{q})}{2|\stackrel{}{q}|}}{\displaystyle \frac{1}{(p_0|\stackrel{}{q}|)^2|\stackrel{}{q}\stackrel{}{p}|^2}}[8p^2((uq){\displaystyle \frac{(up)(qp)}{p^2}})^2`$ (88)
$``$ $`4{\displaystyle \frac{(pq)^2(up)^3}{(uq)p^2}}+2{\displaystyle \frac{(qp)^3(up)^2}{(uq)^2p^2}}+2(up)^2(qp)]_{q_0=|\stackrel{}{q}|}`$ (89)
$`+`$ $`{\displaystyle }{\displaystyle \frac{d^3q}{(2\pi )^3}}{\displaystyle \frac{f(\stackrel{}{q})}{2|\stackrel{}{q}|}}{\displaystyle \frac{1}{(p_0+|\stackrel{}{q}|)^2|\stackrel{}{q}\stackrel{}{p}|^2}}[8p^2((uq){\displaystyle \frac{(up)(qp)}{p^2}})^2`$ (90)
$``$ $`4{\displaystyle \frac{(pq)^2(up)^3}{(uq)p^2}}+2{\displaystyle \frac{(qp)^3(up)^2}{(uq)^2p^2}}+2(up)^2(qp)]_{q_0=|\stackrel{}{q}|}]`$ (91)
$`+`$ $`(qpq),`$ (92)
and

$`Im\mathrm{\Pi }_{Gl;R}^L(p)`$ $``$ $`{\displaystyle \frac{g^2}{2}}\delta _{ab}N_c{\displaystyle \frac{8p^2}{(up)^2p^2}}{\displaystyle \frac{d^3q}{(2\pi )^2}\frac{f(\stackrel{}{q})}{2|\stackrel{}{q}|}}`$ (94)
$`sgn(p_0|\stackrel{}{q}|)\delta ((|\stackrel{}{q}|p_0)^2|\stackrel{}{q}\stackrel{}{p}|^2)[8p^2((uq){\displaystyle \frac{(up)(qp)}{p^2}})^2)`$
$``$ $`4{\displaystyle \frac{(pq)^2(up)^3}{(uq)p^2}}+2{\displaystyle \frac{(qp)^3(up)^2}{(uq)^2p^2}}+2(up)^2(qp)]_{q_0=|\stackrel{}{q}|}]`$ (95)
$`+`$ $`{\displaystyle \frac{f(\stackrel{}{q})}{2|\stackrel{}{q}|}}sgn(p_0+|\stackrel{}{q}|)\delta ((|\stackrel{}{q}|+p_0)^2|\stackrel{}{q}\stackrel{}{p}|^2)[8p^2((uq){\displaystyle \frac{(up)(qp)}{p^2}})^2)`$ (96)
$``$ $`4{\displaystyle \frac{(pq)^2(up)^3}{(uq)p^2}}+2{\displaystyle \frac{(qp)^3(up)^2}{(uq)^2p^2}}+2(up)^2(qp)]_{q_0=|\stackrel{}{q}|}]`$ (97)
$`+`$ $`(qpq)`$ (98)
For the symmetric part of the self energy we find:

$`\mathrm{\Pi }_{Gl;S}^{(T)}(p)`$ $`=`$ $`{\displaystyle \frac{g^2}{2i}}\delta _{ab}N_c{\displaystyle \frac{4p^2}{(up)^2p^2}}{\displaystyle \frac{d^3q}{(2\pi )^2}\delta \left((|\stackrel{}{q}|p_0)^2|\stackrel{}{q}\stackrel{}{p}|^2\right)|\stackrel{}{q}|(u_0\stackrel{}{u}\widehat{q})^2f(\stackrel{}{q})\stackrel{~}{f}(\stackrel{}{p}\stackrel{}{q})}`$ (99)
$`+`$ $`\delta \left((|\stackrel{}{q}|+p_0)^2|\stackrel{}{q}\stackrel{}{p}|^2\right)|\stackrel{}{q}|(u_0\stackrel{}{u}\widehat{q})^2f(\stackrel{}{q})\stackrel{~}{f}(\stackrel{}{p}\stackrel{}{q}).`$ (100)
It is interesting to see that in the Hard Loop Momentum Approximation we have $`\mathrm{\Pi }_{Gl;S}^{(T)}(p)\mathrm{\Pi }_{Gl;S}^{(L)}(p)`$.


text:
# 1 Introduction

## 1 Introduction

The dynamics of spatially extended systems can be measured by sets of multi-detector arrays. Most spatio-temporal analysis methods fitting multi-dimensional dynamical models  consider data over the full time range. In , a method was described for partitioning spatio-temporal signals into time segments, in which the signal can be modeled by deterministic ordinary differential equations near fixed points. Each dynamical system is determined by a nonlinear spatio-temporal analysis . The earlier proposed algorithm in  works with an arbitrary number of clusters $`k`$. Since results are dependant on $`k`$, an objective criterion for the number of clusters is necessary. The present report introduces a different segmentation algorithm and aims to derive an indirect criterion for the number of temporal segments.
In the present report, we use the K-Means algorithm  for segmenting data, which addresses each data point to a single cluster. Since K-Means works with an arbitrary number of clusters and this number is crucial to clustering results, we derive a probability function representing the degree of membership of a data point at time $`t`$ to a cluster. It addresses data to clusters or transition parts between clusters and hence determines the number of necessary clusters. Applications to simulated non-stationary data illustrates the probability measure.

## 2 Fixed Point Clustering (FPC)

In the following, a signal trajectory is assumed as compound of a sequence of segments governed by saddle point dynamics. Under the hypothesis, that these segments comprise the main functionality of the underlying system, we aim to extract them from the signal. Trajectories approach saddle points along their stable manifolds whereas they leave the vicinity of the fixed points along the unstable manifolds. The signal points accumulate close to the fixed points if the signal is sampled at a constant rate. This accumulation may also be regarded as a point cluster in data space. Subsequently, stable manifolds in multi-dimensional signals lead to point clusters and their detection can be treated as a recognition problem in data space .

### The clustering algorithm

A $`N`$-dimensional spatio-temporal signal can be described by a data vector $`𝐪(t)\mathrm{}^N`$, where the component $`q_j(t_i)`$ represents a data point at time $`i`$ and detection channel $`j`$. The clustering algorithm aims at cluster centers {$`𝐤_k`$}, whose mean Euclidean distance to a set of data points $`𝐪(t_i)`$ is minimal. The presented implementation follows Moody et al. and is sketched in Fig. 1.
Cluster centers $`𝐤^0`$ are initialized at random locations in the data and their Euclidean distances to each data point are calculated. K-Means defines memberships of data points to a cluster by the smallest Euclidean distance to its center. Thus, data are segmented into $`k`$ clusters and new cluster centers $`𝐤^1`$ are calculated as means of clustered data points. Distances between data points and centers $`𝐤^n`$ are re-estimated until a convergence condition is fulfilled. This criterion can be set either as a upper Euclidean distance limit between sequential cluster centers $`𝐤^n,𝐤^{n+1}`$ or as number of iterations. We choose to limit the number of iterations to $`25`$.

### Simulated spatio-temporal data and results

Now, a low-dimensional simulated signal $`𝐀(t)`$ is introduced describing amplitudes of multi-dimensional spatial patterns $`𝐯_i`$ by

$`𝐪(t)={\displaystyle \underset{i}{}}A_i(t)𝐯_i.`$
This superposition describes a spatio-temporal signal $`𝐪(t)`$.
The dataset $`𝐀(t)`$ is generated by

$`\dot{A}_1=ϵA_1A_1[A_1^2+(2+b)A_2^2+(2b)A_3^2]+\mathrm{\Gamma }(t)`$
$`\dot{A}_2=ϵA_2A_2[A_2^2+(2+b)A_3^2+(2b)A_1^1]+\mathrm{\Gamma }(t)`$ (1)
$`\dot{A}_3=ϵA_3A_3[A_3^2+(2+b)A_1^2+(2b)A_2^2]+\mathrm{\Gamma }(t).`$
Parameters are set to $`ϵ=1`$, $`b=2`$ and $`\mathrm{\Gamma }(t)[0.05\mathrm{}0.05]`$ represents additive noise following a uniform deviate. Equations 1 describe the convection onset of a Rayleigh-Benard-experiment in the presence of rotation.
A 3-dimensional trajectory $`𝐀(t)`$ is calculated by $`2200`$ integration steps with the initial condition $`𝐀(t=0)=(0.03,0.2,0.8)`$, see Fig. 2. The trajectory passes the saddle points $`𝐀_3^0=(0,0,1),𝐀_1^0=(1,0,0)`$ and $`𝐀_2^0=(0,1,0)`$ in this sequence, and then returns to $`𝐀_3^0`$.
The K-Means algorithm is applied on the simulated data for different number of clusters $`k=2,..,7`$. In Fig. 3, the Euclidean distances from each data point to the determined cluster centers are plotted in temporal sequence for each $`k`$. When a trajectory approaches or moves from a cluster center, its Euclidean distance to the center decreases resp. increases. These changes can be observed in Fig. 3. For fixed number of clusters, each data point is considered to be member of a cluster, whose center is closest to the data point.
Comparing obtained clustering results for different $`k`$, clustered time windows \[0;$``$350\], \[$``$350;$``$1050\], \[$``$1160;$``$1610\] and \[$``$1740;2200\] are recognized, which borders remain similar for different $`k`$.

## 3 The cluster criterion

Although there might be only a limited number of clusters $`k_d<k`$ in the data, K-Means determines $`k`$ clusters also including void clusters. In Fig. 3, small clustered time windows are visible, whose occurences and temporal widths strongly depend on $`k`$. They are considered as invalid clusters. Conversely, a first qualitative criterion for valid clusters may be formulated as:

* cluster widths and locations in time remain independent of $`k`$ and
* the Euclidean distances of clustered data points to centers is obviously smaller then the Euclidean distances of points to the next nearest cluster center and
* the width of the clustered time window is not too small.

Although these criteria are rather heuristic than formal, they proved to be useful in practice . Now, we try to evolve them quantitatively. The first item can be formulated as a sum over all clustering results: valid contributions are additive if they occur for all $`k`$, others vanish in the sum as small contributions. Thus the contribution of a valid cluster to the sum should be large, not reliable clusters should contribute with small values. A good quantity for these contributions is the area between the curves of the signal-nearest cluster-distance and signal-next cluster-distance. This definition allows the analytical formulation of the second item and is outlined in Fig. 4. Each data point $`t_i`$ obtains an index corresponding to the cluster $`j`$ it is member of. The index is equal the relativ area $`\frac{A_j^{(k)}(t_i)}{T_jA_j^{(k)}}`$, where $`T`$ denotes the number of data points. By summing up the indices over $`K`$ cluster realizations for every data point, a degree of membership $`M(t)`$ for every data point is obtained:

$`A_r^{(k)}(t_i)`$ $`=`$ $`{\displaystyle \frac{A_j^{(k)}(t_i)}{_jA_j^{(k)}}}`$
$`M(t_i)`$ $`=`$ $`{\displaystyle \frac{_{k=2}^{KT}A_r^{(k)}(t_i)}{_{i=1}^T_{k=2}^{KT}A_r^{(k)}(t_i)}}.`$
$`M(t)`$ represents a probability, that a data point at time $`t`$ belongs to a cluster. The application to simulated data with $`K=30`$ leads to results of $`M(t)`$ shown in Fig. 5. Clustered time windows can be recognized as regions of high values of $`M(t)`$. Four plateaus of $`M(t)`$ are recognized at \[0;250\], \[340;1010\], \[1180;1560\] and \[1750;2200\] with borders at drastic value changes. Regions between these plateaus are considered as non-functional transitions parts. Comparing time windows in the original signal(Fig. 2) and detected clustered time windows in Fig. 3 and Fig. 5, good accordance of time windows near fixed points and cluster results are recognized.

## 4 Conclusion

The present brief report extends the fixed point clustering method by introducing a probability function $`M(t)`$. High values of $`M(t)`$ indicate clustered data points. Fixed point clustering relates temporal dynamics near fixed points showing attractive and repelling properties with clusters in dataspace. By the presented extension, regions in data space near such fixed points can be determined independant of the number of clusters. Applications to spatio-temporal signals in hydrodynamics, metereology or brain science  are possible.