text:
# Low Temperature Metallic State of Ultrathin Films of Bismuth

## Abstract

Measurements of resistance vs. temperature have been carried out on a sequence of quench-condensed ultrathin films of amorphous bismuth (a-Bi). The resistance below about $`0.1K`$ was found to be temperature independent in a range of films with thicknesses spanning the superconductor-to-insulator transition that would be inferred by analyzing data obtained above $`0.14K`$. Film magnetoresistance was temperature dependent in the same temperature range over which the resistance was temperature independent. This implies that the low temperature metallic regime is intrinsic, and not a consequence of failure to cool the electrons.

Research on ultrathin superconducting films has been focused on superconductor-to-insulator (SI) transitions , which are believed to be quantum phase transitions (QPTs). In a QPT, tuning an external parameter of the Hamiltonian alters the ground state. The external parameters in previous work have been disorder, varied by changing the thickness of quench-condensed films, and perpendicular and parallel magnetic fields .

A complication in studying QPTs is that they occur at zero temperature, whereas data is acquired down to the lowest accessible temperatures, which are naturally nonzero. Finite-size scaling is used to infer the existence of the QPT, and to extract critical exponents. A successful finite-size scaling analysis is evidence of a QPT, but is not a proof. If other physics intervenes at temperatures lower than those accessed in the measurements, one may draw incomplete, if not erroneous, conclusions . For this reason it is essential to strive to extend measurements to lower temperatures in search of the “true” ground state.

Metallic behavior at low temperatures was recently reported in studies of the field-driven transition of homogeneous MoGe  films that would have been considered to be superconducting in the $`T0`$ limit, based on a finite-size scaling analysis. Earlier, quench-condensed granular films exhibited metallic behavior in zero magnetic field, beginning at temperatures$`(2K)`$, just below those at which the resistance initially began to drop towards zero . In previous work on ultrathin films of a-Bi grown on an a-Ge layer, there were no indications of a metallic regime down to temperatures of order $`0.14K`$ .

In this letter we report measurements of $`R(T,H)`$ of a-Bi films over a range of thicknesses spanning the putative SI transition. We find metallic behavior in zero magnetic field, as evidenced by a nonzero temperature-independent resistance at the lowest measured temperatures, which were of order $`0.050K`$ . Considering only data taken at temperatures above $`0.14K`$, this same set of films would appear to be undergoing a thickness-tuned SI transition. A metallic regime was also present when the films were subjected to magnetic fields applied in the plane. Such fields mainly affect the amplitude of the order parameter, as they do not introduce vortices as perpendicular fields do.

Amorphous Bi films were deposited in situ at liquid helium temperatures onto single-crystal SrTiO3(100) substrates, pre-coated (in situ) with a $`6`$Å thick film of a-Ge . This was done in a chamber held at a pressure of $`10^{10}`$ Torr. To prevent annealing, substrate temperatures were held below $`12K`$ during growth, and below $`18K`$ during other processing and handling. Film thicknesses were increased in increments as small as $`0.05`$ Å, as measured using a calibrated quartz crystal monitor. Films processed in this fashion are believed to be homogeneously disordered . However, a previous experiment showed that a-Ge plays an active role in electrical transport. Critical features of the present experiments are the ability to change the nominal thickness of a film in tiny increments and to grow films that are homogeneous in thickness to one part in 10$`^4.`$ This was accomplished by using a large source-to-substrate distance (60cm) and employing Knudsen cells as vapor sources, and maintaining ultrahigh vacuum conditions during growth. All of the effects reported here occurred over a nominal thickness range of order 0.2Å out of approximately 9.0Å, and would not have been seen without such stringent control.

Resistance measurements were carried out in a Kelvinox 400 dilution refrigerator, using either DC or AC four-probe techniques. Electrical leads into the cryostat were filtered at room temperature using $`\pi `$-section filters with a cutoff frequency of about 500 Hz. Extra filtering was employed when a DC current source was used, bringing the cutoff frequency down to about 10 Hz. Power dissipation in the measurement process was kept below 1 pW. Only magnetic fields parallel to the film plane could be applied, with values extending up to $`12.5T`$. In the process of insuring adequate cooling of the films, rotation of the sample platform was restricted to a few degrees. Parallel alignment was adjusted very carefully at room temperature. In high fields, even a slight misalignment will result in a substantial perpendicular field component, which could dominate the high-field physics. For example, an error of $`1^0`$ at a field of $`12T`$ would result in a perpendicular component of about $`0.2T`$. This is greater than the observed critical tuning field in previous studies of the SI transition in a-Bi films. Given the fact that one full turn of the rotator control corresponded to about $`2^0`$, the actual error of the alignment in this work is estimated be of order $`0.1^0`$.

The evolution of $`R(T)`$ of eleven films with thicknesses ranging from 8.5 Å to 9.3 Å is shown in Fig 1. Thinner films and thicker films, grown in other runs (not shown) were insulating and glass-like in their responses, or fully superconducting, respectively. Focusing on temperatures above $`0.14K`$, insulating and superconducting films appear to be clearly separated into two groups. Although we do not show it here, resistance data of this set of films could be collapsed in the usual manner employing a finite size scaling form, $`R(\delta ,T)=R_cF(\delta /T^{\frac{1}{\nu z}})`$, with a critical exponent product $`\nu z`$ of $`1.1\pm 0.1`$ in agreement with previous work . The coherence length and dynamical critical exponents are $`\nu `$ and $`z`$ respectively, and $`\delta `$ is the control parameter given by $`dd_c`$, where $`d`$ is film thickness and $`d_c`$ is a critical thickness. This scaling breaks down when temperature-independent, metallic, resistance data below $`0.14K`$ are included in the analysis. It should be noted that this metallic regime was found for both insulator- and superconductor-like films, as might be inferred from the data above $`0.14K`$.

In interpreting a temperature independent $`R(T)`$, one must consider the possibility that as the measured temperature decreases, the electrons do not cool. This is an important issue for disordered ultrathin films, which, as well as being antennae for electromagnetic radiation, have tiny heat capacities. Apart from heating by the measuring current, which is unlikely because I-V characteristics were always linear, electromagnetic noise is likely the major source of heating. Noise can be external or internal to the cryostat. The latter would be Johnson noise in the electrical leads that are in equilibrium with black-body radiation. These leads are at temperatures higher than that of the film over much of their length. To be certain that the electrons have cooled would require measurement of some well-understood property of the film and its use as a thermometer, or the provision of a separate thermometer that is in certain thermal contact with the film that has a response to possible heating effects from the environment identical to that of the film. Since neither of these possibilities was available, we rely on indirect arguments. These will emerge in the course of the discussion of $`R(T,H)`$, which follows.

In Fig.2, $`R(H)`$ at various temperatures is plotted for the $`9.3\AA `$ thick film. At the highest temperatures, in the lowest fields, $`dR/dH`$ is slightly positive. With increasing field, a regime in which $`dR/dH`$ becomes negative is entered. This is followed by a minimum in $`R(H)`$, and finally an upturn at higher fields. In sufficiently high fields superconducting fluctuations are completely quenched. The regime in which $`dR/dH<0`$ is possibly associated with the suggestion of Kivelson and Spivak, that a local order parameter density can fluctuate from point to point in sign as well as in magnitude in disordered systems near the SI transition. This can bring about a negative magnetoresistance in a perpendicular magnetic field. As stated previously the small error in alignment can result in a small component of magnetic field perpendicular to the plane.

The systematics of the magnetoresistance provide evidence that the metallic regime is intrinsic. The negative magnetoresistance grows as T decreases as shown in Fig.2, but becomes abruptly positive when $`R(T)`$ becomes temperature independent at low temperatures. It continues to increase with further decrease of $`T`$. This can be seen in the inset of Fig. 2 which emphasizes the low field regime. The magnetoresistance in the metallic regime is a function of temperature, both in its initial change with field and for all values of field. If the saturation of the resistance were a consequence of electrons not cooling then these effects would not be expected.

The upturn of $`R(H)`$ in high parallel fields is hard to explain quantitatively as it is a consequence of the combined action of the large parallel component of the field and the very much smaller, but unknown perpendicular component. A linear dependence of $`R(H)`$ on field would be expected if there were flux flow resistance due to a perpendicular field component . It should be noted that at high fields there is a significant contribution to the resistance which is quadratic with field. This could be associated with the weakening of the order parameter amplitude fluctuations by the parallel component of field.

Further interesting effects are seen if the temperature dependence of the resistance is plotted for fields above the magnetoresistance minimum. The quenching of superconducting fluctuations in this regime is shown in Fig. 3 for the $`9.19\AA `$ thick film. Its low temperature behavior is always metallic. As the field is increased, the domain in temperature over which $`dR/dT=0`$ broadens significantly and reaches a maximum width of about $`100mK`$.. This maximum width can be seen much more dramatically in Fig. 4 which shows data from the $`9.09\AA `$ thick film, in which the temperature of the onset of the metallic regime is plotted as a function of its resistance. For this film, the resistance of the widest metallic regime was $`12,900\mathrm{\Omega }`$, which is remarkably close to twice the value of the quantum resistance for pairs. This resistance is actually equal to the zero-field resistance of the film in the low temperature limit. This particular film is the first of those in the sequence of films shown in Fig. 1 that exhibits superconducting fluctuations as evidenced by a downturn in $`R(H)`$ at higher temperatures. All of the films exhibiting superconducting fluctuations display a similar response to an applied parallel magnetic field. It is important to note that the systematics of this variation of resistance with field are different from what was reported by Mason and Kapitulnik..  This may be a consequence of the fact that the field in this work was parallel rather than perpendicular to the film plane.

One may compare the present results with those obtained in studies of quench-condensed a- Ga films which behaved in a manner very similar to the films of the present work, but became metallic at temperatures as high as $`2K`$ . In that work, the temperatures at which $`dR/dT`$ first fell to zero were dependent on the measuring current. In the subsequent work of Ref., on similar films, it was proposed that the nonlinearities and metallic behaviors in films that appeared to become insulating at high temperatures were a consequence of charge motion, whereas for those that appeared to become superconducting, they were due to vortex motion. This was arrived at because the systematics of the data suggested that insulating and superconducting states at $`T=0`$ would only be found in the limit of zero measuring current. Although the I-V characteristics were nonlinear, the systematics of the data could also be used to eliminate the possibility of heating. It is conceivable that similar nonlinear features of the I-V characteristics for a-Bi films exist, but their observation would require much lower current levels than were employed in the measurements.

Metallic regimes at low temperatures in ultrathin films that exhibit local superconductivity have been the subject of several theoretical works, which emphasize the bosonic nature of these systems as well as, in some instances, the role of dissipation . The data we have presented cannot validate any of these theories, some of which make very specific predictions, such as glass-like behavior. On the other hand, the data can be used to make a case for the metallic regime being intrinsic and not a consequence of failure to cool the electrons.

The values of either thickness or magnetic field that result in maximally wide metallic regimes could be evidence of the metallic quantum critical point of the SI transitions proposed in the context of phase-only Bose Hubbard models . We conjecture that when the tuning parameter of the quantum phase transition is close to its critical value, evidence of quantum critical fluctuations, in this case, a metallic regime, extends to high temperatures, resulting in behavior such as that exhibited in Fig. 4. In this scenario metallic behavior at other fields (or thicknesses) is seen at nonzero temperatures because the free motion of vortices and charges that respond to the measuring current masks the superconducting and insulating ground states, as discussed above. The actual value of the resistance of the separatrix for the thinnest film exhibiting superconducting fluctuations is very close to $`h/2e^2`$ which is double the quantum resistance for electron pairs. It turns out that this is the value of the critical resistance obtained from finite-size scaling of Monte Carlo data in simulations of the phase-only Bose Hubbard model. This scenario implies that there is a quantum critical point separating insulating and superconducting phases, with the metallic regime that obscures it resulting from a nonzero measuring current. Alternatively, there might not be a quantum critical point at all. In this instance, the striking magnetic field dependence of the temperature onset of the metallic regime such as shown in Fig. 4 would be just a feature of the phase diagram of the putative Bose metal regime.

In summary, the extension of measurements of $`R(T)`$ of ultrathin, quench-condensed films of a-Bi down to temperatures of order $`0.05K`$ have revealed metallic behavior over a range of film thicknesses near the SI transition. This suggests that the apparent superconducting and insulating ground states inferred from the analysis of data obtained at temperatures above $`0.14K`$ and discussed in many theories, may not be realized. Evidence was presented supporting the results being intrinsic, and not a consequence of the electrons in the film not being cooled. Lower temperatures and lower measuring current densities may be needed to find the “true” ground state, at least for films in the transition regime between insulating and superconducting behavior.

This work was supported by the National Science Foundation Condensed Matter Physics Program under grant Grant DMR-0138209.

….


format: text

text:
# Minijet Initial Conditions For Non-Equilibrium Parton Evolution at RHIC and LHC

## .1 Boost Invariant Initial Distribution Function From Minijets

First let us assume that a boost invariant description bjorken  is appropriate. Then the gluon distribution function depends only on the boost invariant parameters $`\tau `$, $`p_T`$ and $`\xi =\eta y`$. We will parametrize the $`\eta y`$ correlation as Gaussian:

$$f(p_T,\eta y,\tau _0)=f(p_T)e^{\frac{(\eta y)^2}{\sigma ^2(p_T)}},$$
(10)
where $`f(p_T)`$ will be determined from eq. (9) by using eq. (3) which is obtained by using pQCD. Using the above equation in eq. (9) we get:

$$\frac{dN^{jet}}{dp_T}=g_C2\pi ^2R_A^2\tau _0p_T^2f(p_T)𝑑y_{\mathrm{}}^{\mathrm{}}𝑑\eta \mathrm{cosh}(\eta y)e^{\frac{(\eta y)^2}{\sigma ^2(p_T)}}$$
(11)
which gives:

$$\frac{dN^{jet}}{dp_T}=g_C2\pi ^2\sqrt{\pi }R_A^2\tau _0p_T^2f(p_T)\sigma (p_T)e^{\sigma ^2(p_T)/4}𝑑y$$
(12)
where the maximum allowed value of the momentum rapidity at RHIC and LHC depends on $`p_T`$ and $`\sqrt{s}`$. Integrating over $`y`$ we get from the above equation:

$$f(p_T)=\frac{\frac{dN^{jet}}{dp_T}}{g_C2\pi ^2\sqrt{\pi }R_A^2\tau _0p_T^22\mathrm{ln}(\sqrt{s}/2p_T+\sqrt{s/4p_T^21})}\frac{e^{\sigma ^2(p_T)/4}}{\sigma (p_T)},$$
(13)
which gives a boost invariant initial distribution function:

$$f(p_T,\xi ,\tau _0)=\frac{\frac{dN^{jet}}{dp_T}}{g_C2\pi ^2\sqrt{\pi }R_A^2\tau _0p_T^22\mathrm{ln}(\sqrt{s}/2p_T+\sqrt{s/4p_T^21})}\frac{e^{\frac{\xi ^2}{\sigma ^2(p_T)}}}{\sigma (p_T)e^{\sigma ^2(p_T)/4}}.$$
(14)
Using this boost invariant initial phase-space distribution function of the gluon the initial minijet number density is given by:

$`n(\tau _0)=g_C{\displaystyle 𝑑\mathrm{\Gamma }p^\mu u_\mu f(p_T,\eta y,\tau _0)}=g_C{\displaystyle d^2p_Tp_T𝑑\xi \mathrm{cosh}\xi f(p_T,\xi ,\tau )}`$
$`=2\pi {\displaystyle 𝑑p_T\frac{\frac{dN^{jet}}{dp_T}}{2\pi ^2R_A^2\tau _02\mathrm{ln}(\sqrt{s}/2p_T+\sqrt{s/4p_T^21})}},`$ (15)
where $`d\mathrm{\Gamma }=\frac{d^3p}{p^0}`$. Similarly the initial energy density is given by:

$`ϵ(\tau _0)=g_C{\displaystyle 𝑑\mathrm{\Gamma }(p^\mu u_\mu )^2f(p_T,\eta y,\tau _0)}=g_C{\displaystyle d^2p_Tp_T^2𝑑\xi \mathrm{cosh}^2\xi f(p_T,\xi ,\tau )}`$
$`=2\pi {\displaystyle 𝑑p_Tp_T\frac{\frac{dN^{jet}}{dp_T}}{4\pi ^2R_A^2\tau _02\mathrm{ln}(\sqrt{s}/2p_T+\sqrt{s/4p_T^21})}(e^{\sigma ^2(p_T)/4}+e^{3\sigma ^2(p_T)/4})}.`$ (16)
## .2 Boost Non-Invariant Initial Distribution Function From Minijets

In the above we used a boost invariant distribution function $`f(\tau ,\xi ,p_T)`$ as the initial phase-space distribution function. However this is not consistent with the actual rapidity distribution function obtained from pQCD. To relax the boost invariant assumption, we will introduce a function $`f(p_T,y)`$ instead of $`f(p_T)`$ in eq. (10). The form of the boost non-invariant distribution we take is:

$$f(p_T,\eta ,y,\tau _0)=f(p_T,y)e^{\frac{(\eta y)^2}{\sigma ^2(p_T)}},$$
(17)
where $`f(p_T,y)`$ is obtained from $`\frac{dN^{jet}}{dydp_T}`$. Note that $`f(p_T)`$ in eq. (10) is obtained from $`\frac{dN^{jet}}{dp_Tdy}`$ after integrating over $`y`$. However, in the above form (eq. (17)), the $`y`$ dependence is more general which for small $`\sigma ^2`$ is close to that obtained by using pQCD. Using the above form in eq. (9) we get:

$$\frac{dN^{jet}}{dydp_T}=g_C2\pi ^2R_A^2\tau _0p_T^2f(p_T,y)_{\mathrm{}}^{\mathrm{}}𝑑\eta \mathrm{cosh}(\eta y)e^{\frac{(\eta y)^2}{\sigma ^2(p_T)}}$$
(18)
which gives:

$$\frac{dN^{jet}}{dydp_T}=g_C2\pi ^2\sqrt{\pi }R_A^2\tau _0p_T^2f(p_T,y)\sigma (p_T)e^{\sigma ^2(p_T)/4}.$$
(19)
From the above equation we get:

$$f(p_T,y)=\frac{\frac{dN^{jet}}{dydp_T}}{g_C2\pi ^2\sqrt{\pi }R_A^2\tau _0p_T^2\sigma (p_T)e^{\sigma ^2(p_T)/4}}$$
(20)
which gives a boost non-invariant initial phase-space gluon distribution function:

$$f(p_T,\eta ,y,\tau _0)=\frac{\frac{dN^{jet}}{dydp_T}}{g_C2\pi ^2\sqrt{\pi }R_A^2\tau _0p_T^2}\frac{e^{\frac{(\eta y)^2}{\sigma ^2(p_T)}}}{\sigma (p_T)e^{\sigma ^2(p_T)/4}}.$$
(21)
Using the above boost non-invariant gluon distribution function the initial minijet number density is given by:

$`n(\tau _0,\eta )=g_C{\displaystyle 𝑑\mathrm{\Gamma }p^\mu u_\mu f(p_T,\eta ,y,\tau _0)}=g_C{\displaystyle d^2p_Tp_T𝑑y\mathrm{cosh}(\eta y)f(p_T,\eta ,y,\tau _0)}`$
$`=2\pi {\displaystyle 𝑑p_T𝑑y\frac{\frac{dN^{jet}}{dydp_T}}{2\pi ^2\sqrt{\pi }R_A^2\tau _0}\frac{e^{\frac{(\eta y)^2}{\sigma ^2(p_T)}}}{\sigma (p_T)e^{\sigma ^2(p_T)/4}}\mathrm{cosh}(\eta y)},`$ (22)
and the initial minijet energy density given by:

$`ϵ(\tau _0,\eta )=g_C{\displaystyle 𝑑\mathrm{\Gamma }(p^\mu u_\mu )^2f(p_T,\eta ,y,\tau _0)}=g_C{\displaystyle d^2p_Tp_T^2𝑑y\mathrm{cosh}^2(\eta y)f(p_T,\eta ,y,\tau _0)}`$
$`=2\pi {\displaystyle 𝑑p_Tp_T𝑑y\frac{\frac{dN^{jet}}{dydp_T}}{2\pi ^2\sqrt{\pi }R_A^2\tau _0}\frac{e^{\frac{(\eta y)^2}{\sigma ^2(p_T)}}}{\sigma (p_T)e^{\sigma ^2(p_T)/4}}\mathrm{cosh}^2(\eta y)}.`$ (23)
## .3 Initial Energy Density and Number Density at RHIC and LHC

Using the above expressions for the minijet initial phase-space distribution functions we can now predict the initial energy density and number density of the minijet plasma at RHIC and LHC. As temperature can not be defined in a non-equilibrium situation we define an effective temperature by:

$$T_{eff}=[\frac{15}{8\pi ^2}ϵ]^{1/4},$$
(24)
with the non-equilibrium energy density obtained above.

Before we calculate the initial conditions we would like to find out the values of the unknown parameter $`\sigma ^2`$ appearing in the above equation. We will fix this unknown parameter by equating the first $`p_T`$ moment of the distribution function with the pQCD predicted $`E_T`$ distribution as given by eq. (4). We have for the boost invariant case:

$`{\displaystyle 𝑑p_T\frac{dE_T^{jet}}{dp_T}}=2g_C\pi ^2R_A^2\tau _0{\displaystyle 𝑑p_Tp_T^3f(p_T)}={\displaystyle 𝑑p_T\frac{p_T\frac{dN^{jet}}{dp_T}}{\sqrt{\pi }2\mathrm{ln}(\sqrt{s}/2p_T+\sqrt{s/4p_T^21})\sigma (p_T)e^{\sigma ^2(p_T)/4}}},`$ (25)
and for the boost non-invariant case:

$`{\displaystyle 𝑑p_T𝑑y\frac{dE_T^{jet}}{dydp_T}}=2g_C\pi ^2R_A^2\tau _0{\displaystyle 𝑑p_T𝑑yp_T^3f(p_T,y)}={\displaystyle 𝑑p_T𝑑y\frac{p_T\frac{dN^{jet}}{dydp_T}}{\sqrt{\pi }\sigma (p_T)e^{\sigma ^2(p_T)/4}}}.`$ (26)
We determine the value of $`\sigma ^2`$ by using the pQCD values of $`\frac{dE_T^{jet}}{dydp_T}`$ and $`\frac{dE_T^{jet}}{dp_T}`$ from eq. (4) in the left hand side of the above equations. Assuming $`\sigma ^2`$ independent of $`p_T`$ we find from the above equations: $`\sigma ^2`$ = 0.0034(0.0015) at RHIC(LHC) for boost invariant case and $`\sigma ^2`$ = 0.28(0.28) at RHIC(LHC) for boost non-invariant case. Note that for boost non-invariant case $`\sigma ^2`$ is same at RHIC and LHC as the normalization factor $`\sqrt{\pi }\sigma (p_T)e^{\sigma ^2(p_T)/4}`$ in eq. (26) is energy independent whereas the corresponding normalization factor $`\sqrt{\pi }2\mathrm{ln}(\sqrt{s}/2p_T+\sqrt{s/4p_T^21})\sigma (p_T)e^{\sigma ^2(p_T)/4}`$ in eq. (25) depends on energy. Using the above values for the boost invariant case in eq. (16) we find $`ϵ_0`$=31 GeV/fm3 with $`T_{eff}`$= 460 MeV at RHIC and $`ϵ_0`$=287 GeV/fm3 with $`T_{eff}`$= 800 MeV at LHC. The number density for the boost invariant case obtained from eq. (15) is $`n_0`$=20/fm3 at RHIC and 85/fm3 at LHC. Using the values of the $`\sigma ^2`$ for boost non-invariant case we find from eq. (23): $`ϵ_0`$=51 GeV/fm3 with $`T_{eff}`$=521 MeV at RHIC and $`ϵ_0`$=517 GeV/fm3 with $`T_{eff}`$=932 MeV at LHC for $`\eta `$=0. The number density for the boost non-invariant case obtained from eq. (22) is $`n_0`$=30/fm3 at RHIC and 138/fm3 at LHC for $`\eta `$=0. It can be noted that the energy density for the boost invariant case is less than that of the boost non-invariant case. This is because we have assumed a boost invariance in the entire rapidity range. However, the actual pQCD minijet rapidity distribution is not flat in the entire rapidity range as can be seen from Fig. 2. Therefore, in order to maintain the boost invariance in the entire rapidity range the density has to be smaller. However, in boost non-invariant case the rapidity distribution of the minijet is minimaly altered and the pQCD form of $`\frac{dN}{dp_Tdy}`$ is used throughout the calculation.

Let us mention briefly other approaches for obtaining initial condition at RHIC and LHC. Commonly used way is to obtain energy density from eq. (5) by using pQCD estimate of the minijet $`E_T`$ distribution hijing ; eskk ; gwa ; bhal . This formula assumes that the initial partons are uniformly distributed in the coordinate space so that one divides the volume to obtain energy density. In nayak  a boost invariant Boltzmann form $`f(x,p,t_0)=e^{p_T\mathrm{cosh}(\eta y)/T}`$ was used for the initial distribution function with $`\eta y`$ (boost invariant quantity) correlation taken into account. In the parton cascade model gegm  momentum space structure function along with coordinate distribution of the nucleons inside the nucleus was used directly in the transport equation which is different from the minijet calculation by using pQCD. In the approach of this paper we have obtained a phase-space initial non-equilibrium distribution function $`f(x,p)`$ for the boost non-invariant case with minimal altering the $`p_T`$ and $`y`$ distribution of the minijet partons obtained by using pQCD calculations.

To conclude, we have obtained a physically reasonable initial phase-space distribution function of the partons formed at the very early stage of the heavy-ion collisions at RHIC and LHC which can be used to study equilibration of the quark-gluon plasma bhal ; nayak ; cw . To obtain such phase-space distribution function we have considered several plausible parametrizations of $`f_i(x,p)`$ and have fixed the parameters by matching $`f(x,p)p^\mu 𝑑\sigma _\mu `$ to the invariant momentum space semi-hard parton distributions obtained using QCD (pQCD) as well as fitting low order moments of the distribution. Once the parameters of $`f(x,p)`$ are found, we then have determined the initial number density, energy density and the corresponding (effective) temperature of the minijet plasma at RHIC and LHC energies. For a boost non-invariant minijet phase-space distribution function we obtain $``$ 30(140) /fm3 as the initial number density, $``$ 50(520) GeV/fm3 as the initial energy density and $``$ 520(930) MeV as the corresponding initial effective temperature at RHIC(LHC). Our minijet initial conditions are obtained above $`p_T`$=1(2) GeV at RHIC(LHC). As partons below these momentum cut-off values might be described by the formation of a classical chromofield, the $`f(x,p,t_0)`$ obtained in this paper is the distribution function for the matter part. Initial conditions for lower transverse momentum has to be supplied in form of the strength of the chromofield or in the form of field energy density in the non-abelian transport equation to study production and evolution of the quark-gluon plasma at RHIC and LHC.

###### Acknowledgements.

We thank Ulrich Heinz for useful discussions.


text:
# Untitled Document

An Electrical Network Model of Plant Intelligence

Bikas K. Chakrabarti(1) and Omjyoti Dutta(2)

(1)Saha Institute of Nuclear Physics, 1/AF Bidhan Nagar, Kolkata-700064.

(2)Electrical Engineering Dept. (4th Yr.), Jadavpur University, Kolkata-700032.

Abstract

A simple electrical network model, having logical gate capacities, is proposed here for computations in plant cells. It is compared and contrasted with the animal brain network structure and functions.

————————————————

Talk at the Condensed Matter Days-2002, held in Bhagalpur University, Bhagalpur, during August 29-31, 2002 (to be published in Ind. J. Phys.).

Considerable investigations and efforts have been made in understanding how the animal and human brains compute and recognise various spatial and temporal patterns . The essential model consists of a network of large number (about 1012 in case of human, 106 in case of birds) of two state electrical devices called neurons which are capable of just summing over the various (milli-volt order) input electrical pulses for a short synaptic period (of milli-second order) collected by the (106 or so) dendrites of each neuron, and comparing this sum with a threshold. The synaptic interactions among the neurons develop during the “learning" process, and can be both excitory or inhibitory, rendering the network randomly frustrated. The computational capabilities emerge out of the collctive dynamics of the network, which is nonlinear (due to the threshold behaviour of each neuron). For symmetric interactions, one can define an energy function (or free energy at finite noise or “temperature" level) for the network and the local free energy minima corresponds to the various local attractor patterns or memory states of the network (Hopfield ). For long-range interactions, the statistical physics of such a network is analytically tractable to a large extent (Amit et al ). The processing of informations in such network models and their detailed analysis are now established (see e.g., Nishimori ). These demonostrated capabilities of such networks are of course very limited in their emerging computational abilities  and far short of anything like consciousness, where some aspects of quantum mechanics (entanglement in the molecules in microtubules of a single neuron) are speculated to be involved .

Are the plants around us intelligent? Do they also deserve our attention in this context of modelling for information processing and computation? Plants have remarkable adaptibility in changed environments. They survive in every landscape of this earth, representing almost 99$`\%`$ of its biomass. Such marvellous adaptive behaviour must be interpreted to be intelligent; although naive definitation of intelligence seem to involve movement of the animal (either bodily or part of it) and plants can not move (bodily) . How such intelligent behaviour of mindless plants, having no brain, compare with those of animals ? Plants do not have neuronal cells either.

Almost eighty years back, Bose detected electrical signalling between plant cells in coordinating its responses to the environment . Although the chemical diffusion of (uncharged) molecules is a dominent source of signalling between the plant cells, it is a very slow mode. It is now established  that some signals are trasmitted within the plants at much smaller time scale (with signal velocity about 30-400 mm/sec, depending on the plant and its environment). Such fast transmissions are due to electrical pulses, generated by ionic motions within the plant cells. Although not the dominent mode, except in some very sensetive plants like Desmodium or Mimosa , the electrical mode (due to migration of Ca+, K+, etc ions) generally present in the cells of all the plants . However, these electrically excitable plant cells do not have many dendrites, like for the neurons, nor are they connected by random excitory/inhibitory (frustrating) interactions.

In absence of the highly connected (frustrated) network of neuron-like units, as in the animal brains, the plants might be utilising the non-linear current (I)-voltage (V) characterestics of their cell membranes for the logical operations (gates). In fact, the plant vacuolar membrane current-voltage characterestics  is now established to be equivalent to that of a Zener diode, as indicated in Fig. 1.

Fig. 1. (a) Typical I-V curves of the plant vacuolar membrane fast-activating channel (from ). The current being mainly due to Ca+ ions and the (reversible) effect of divalant Putrescine (C4H14N$`{}_{}{}^{2^+}{}_{2}{}^{}`$) are shown at two different concentrations. (b) The equivalent Zener diode-like behaviour of the membrane, where the Putrescine concentrations ($`c`$) modulate the changes in the threshold voltage VT.

One can utilise such a threshold behaviour of the plant cell membranes to develop or model gates for perfoming simple logical operations. In Fig. 2, such a model network containing four such threshold units; one in the output and the other three in the input. Each of these threshold units is modelled as a binary unit, having two states: 0 and 1. The inter-unit connection strength is denoted here by the matrix $`W`$. The output $`O`$ of the network considerd can then be expressed as

$$O=\theta (I\varphi ),$$
$`(1)`$
where $`\theta `$ is the step function ($`\theta (x)=1`$ for $`x0`$; and 0 otherwise),

$$I=W_1I_1+W_2I_2+W_3I_3,$$
$`(1a)`$
and $`\varphi `$ is the threshold strength (determined by the threshold voltage $`V_T`$) for the output unit. $`I_1`$, $`I_2,`$ $`I_3`$ are the inputs to the three input units and $`W_1`$, $`W_2`$, $`W_3`$ are their connectivity strengths with the output unit, as indicated in Fig. 2.

Fig. 2. A simple network containing four threshold units (three in the input and one in the output) for performing logical operations like AND, OR, NAND, etc.

For different combinations of $`I_1`$, $`I_2`$ and $`I_3`$, the outputs for different logic gates are given in the Table 1. These can be easily achieved using the combinations of inter-cell connections and the output cell thresholds as:

$$W_1=W_2=W_3=1,\varphi =3,$$
$`(2a)`$
for the AND gate;

$$W_1=W_2=W_3=1,\varphi =1,$$
$`(2b)`$
for the OR gate; and

$$W_1=W_2=W_3=1,\varphi =2$$
$`(2c)`$
for the NAND gate.

Table I: The input-output (truth) table for the logic gates

|  | Inputs |  |  | Output O |  |
| --- | --- | --- | --- | --- | --- |
| $`I_1`$ | $`I_2`$ | $`I_3`$ | AND | OR | NAND |
| 0 | 0 | 0 | 0 | 0 | 1 |
| 0 | 0 | 1 | 0 | 1 | 1 |
| 0 | 1 | 0 | 0 | 1 | 1 |
| 0 | 1 | 1 | 0 | 1 | 1 |
| 1 | 0 | 0 | 0 | 1 | 1 |
| 1 | 0 | 1 | 0 | 1 | 1 |
| 1 | 1 | 0 | 0 | 1 | 1 |
| 1 | 1 | 1 | 1 | 1 | 0 |

These gate capabilities of simple networks of the plant cell membranes, using their nonlinear characterestics and cosequent threshold behaviour (with adjustable thresholds through changed concentrations of, for example, the putrescine and the interaction matrix $`W`$) would allow (cf. ) simple computations in the electrical channels of the plants. It may be noted that such networks here are much more local and tiny in structure, compared to the massively connected and parallelly working network of animal brains. Also, the network matrix $`W`$ elements are either all positive (excitory) or all negative (inhibitory) in eqns. (2). As such, they do not involve any frustration as in the animal brains and have got consequently several limitations in their computational capabilities; for example, they lack the distributed parallel computational capacity, associate memory, etc.

We are grateful to Indrani Bose, Arnab Chatterjee and Dibyendu Sengupta for many useful comments and suggestions.

References:

 J. J. Hopfield, “Neural Networks and Physical Systems with Emergent Computational Abilities”, Proc. Nat. Acad. Sc. (USA) 79 2554 (1982); D. J. Amit, H. Gutfreund and Sompolinsky, “Statistical Mechanics of Neural Networks Near Saturation”, Ann. Phys. 173 30 (1987).

 For reviews see e.g., J. Hertz, A. Krough and R. G. Palmer, Introduction to the Theory of Neural Computation, Addison-Wesley, Reading, MA (1991); B. K. Chakrabarti and P. K. Dasgupta, “Modelling Neural Networks”, Physica 186 33 (1992); H. Nishimori, Statistical Physics of Spin Glasses and Information Processing, Oxford Univ. Press, Oxford (2001).

 R. Penrose, The Emperor’s New Mind, Oxford Univ. Press, Oxford (1999); M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information, Cambridge Univ. Press, Cambridge (2000).

 A. Trewavas, “Mindless Mastery”, Nature 415 841 (2002).

 B. K. Chakrabarti, “Mathematics, Brain Modelling and Indian Concept of Mind”, Proc. Asiatic Soc., Kolkata (in press); also at arXiv:cond-mat/0205094 (2002).

 J. C. Bose, The Nervous Mechanism of Plants, Longmans, Green & Co, London (1923).

 V. A. Shepherd, “Bioelectricity and Rhythms of Sensitive Plants: The Biophysical Research of J. C. Bose”, Current Sc. 77 189 (1999).

 T. Genoud and J.-P. Mètraux, “Crosstalk in Plant Signalling”, Trends in Plant Sc., 4 503 (1999).

 L. I. Brüggemann, I. I. Pottosin and G. Schonknecht, “Cytoplasmic Polyamines Block Fast-Activating Vacuolar Cation Channel”, The Plant J. 16 101 (1998).


text:
# Dynamical model of financial markets: fluctuating ‘temperature’ causes intermittent behavior of price changes

## 1 Introduction

Financial returns are known to be non-gaussian and exhibit fat-tailed distribution fa65 -hkk01 . The fat tail relates to intermittency — an unexpected high probability of large price changes, which is of utmost importance for risk analysis. The recent development of high-frequency data bases makes it possible to study the intermittent market dynamics on time scales of less than a day ms00 -hkk01 .

Using foreign exchange (FX) intraday data, Müller et al. mdd97  showed that there is a net flow of information from long to short timescales, i.e., the behavior of long-term traders influences the behavior of short-term traders. Motivated by this hierarchical structure, Ghashghaie et al. gbp96  have discussed analogies between the market dynamics and hydrodynamic turbulence fr95 ; cgh90 , and claimed that the information cascade in time hierarchy exists in a FX market, which corresponds to the energy cascade in space hierarchy in a three-dimensional turbulent flow. These studies have stimulated further investigations on similarities and differences in statistical properties of the fluctuations in the economic data and turbulence ms96 -bgt00 . Differences have also emerged. Mantegra and Stanley ms96 ; ms97  and Arneodo et al. abc96  pointed out that the time evolution, or equivalently the power spectrum is different for the price difference (nearly white spectrum) and the velocity difference ($`f^{1/3}`$ spectrum, i.e., $`5/3`$-law for the spectrum of the velocity). Moreover, from a parallel analysis of the price change data with time delay and the velocity difference data with time delay (equivalent to the velocity difference data with spatial separation under the Taylor hypothesis fr95 ), it was shown that the time evolution of the second moment and the shape of the probability density function (PDF), i.e., the deviation from gaussian PDF are different in these two stochastic processes ms97 . On the other hand, non-gaussian character in fully developed turbulence fr95  has been linked with the nonextensive statistical physics ts88 -bls01 .

As dynamical foundation of nonextensive statistics, Beck recently proposed a new model describing hydrodynamic turbulence be01 ; be01a . The velocity difference $`\mathrm{\Delta }v`$ of two points in a turbulent flow with the spatial separation $`\mathrm{\Delta }r`$ is described by Brownian motion (an overdamped Langevin equation ri89 ) in a power-law potential. Assuming a $`\chi ^2`$-distribution for the inverse temperature, he obtained a Tsallis distribution ts88  for $`\mathrm{\Delta }v`$. However, if we take into account the almost uncorrelated behavior of the price change ms96 -ms97 , the picture by means of the Brownian motion seems to be more appropriate for the market data rather than turbulence. Moreover, the description by the Langevin equation is able to relate the PDF of the price change to that of the volatility, a quantity known as a measure of the risk in the market. Thus we applied the model to FX market dynamics by employing the correspondence by Ghashghaie et al. gbp96 .

## 2 Model

We substitute the FX price difference $`Z(t)y(t+\mathrm{\Delta }t)y(t)`$ with the time delay $`\mathrm{\Delta }t`$ for the velocity difference $`\mathrm{\Delta }v`$ with the spatial separation $`\mathrm{\Delta }r`$. Beck’s model for turbulence then reads

$$\frac{dZ}{dt}=\gamma F(Z)+R(t),$$
(1)
where $`\gamma >0`$ is a constant, and $`R(t)`$ is gaussian white noise corresponding to the temperature $`kT`$, satisfying $`R(t)R(t^{})=2\gamma kT\delta (tt^{})`$. The ‘force’ $`F=U(Z)/Z`$ is assumed to be obtained by a power-law potential $`U(Z)=C|Z|^{2\alpha }`$ with an exponent $`2\alpha `$, where $`0<\alpha 1`$ and $`C`$ is a positive constant. That is, the system is subject to a restoring force proportional to the power of the price difference, $`|Z|^{2\alpha }`$ besides the random force $`R(t)`$. Especially when $`\alpha =1`$, the restoring force is linear to $`Z`$. Under a constant temperature $`kT`$, Eq. (1) leads to a stationary (i.e., thermal equilibrium) distribution of $`Z`$ as

$$P_{\mathrm{\Delta }t}(Z|\beta )=\frac{e^{\beta U(Z)}}{e^{\beta U(Z)}𝑑Z}=\frac{\alpha }{\mathrm{\Gamma }(\frac{1}{2\alpha })}(C\beta )^{1/2\alpha }e^{\beta C|Z|^{2\alpha }},$$
(2)
where $`\beta 1/kT`$ is the inverse temperature be01 . The ‘local’ variance of $`Z`$, which is defined for a fixed value of $`\beta `$, is obtained from the conditional probability in Eq. (2) as

$$Z^2_\beta =Z^2P_{\mathrm{\Delta }t}(Z|\beta )𝑑Z=\frac{\mathrm{\Gamma }(\frac{3}{2\alpha })}{\mathrm{\Gamma }(\frac{1}{2\alpha })}(C\beta )^{1/\alpha }(kT)^{1/\alpha }.$$
(3)
We define the volatility by the square root of the local variance of $`Z`$ (see Eq. (9)). When $`\alpha =1`$ and $`C=1/2`$, $`Z^2_\beta `$ coincides with the temperature and the conditional PDF reduces to gaussian, while for $`\alpha 1`$, $`\{Z^2_\beta \}^\alpha `$ is proportional to $`kT`$.

Let us assume that the ‘temperature’ of the FX market is not constant and fluctuates in larger time scales, and $`\beta `$ is, just as in Beck’s model for turbulence, $`\chi ^2`$-distributed with degree $`n`$:

$$f_{\mathrm{\Delta }t}(\beta )\frac{1}{\mathrm{\Gamma }(\frac{n}{2})}\left(\frac{n}{2\beta _0}\right)^{n/2}\beta ^{n/21}\mathrm{exp}\left(\frac{n\beta }{2\beta _0}\right),n>2,$$
(4)
where $`\mathrm{\Gamma }`$ is the Gamma function, $`\beta _0`$ is the average of the fluctuating $`\beta `$ and $`n`$ relates to the relative variance of $`\beta `$:

$$\beta =\beta _0,\frac{\beta ^2\beta ^2}{\beta ^2}=\frac{2}{n}.$$
(5)
Equation (4) implies that the local variance $`Z^2_\beta v`$ fluctuates with the distribution $`p(v)v^{(\alpha +1)}f_{\mathrm{\Delta }t}(v^\alpha )`$. The conditional probability in Eq. (2) together with Eq. (4) yields a Tsallis-type distribution ts88 ; be01  for the ultimate PDF of $`Z`$:

$`P_{\mathrm{\Delta }t}(Z)`$ $`=`$ $`{\displaystyle P_{\mathrm{\Delta }t}(Z|\beta )f_{\mathrm{\Delta }t}(\beta )𝑑\beta }={\displaystyle \frac{1}{Z_q}}{\displaystyle \frac{1}{\left\{1+\frac{C\beta _02\alpha (q1)|Z|^{2\alpha }}{2\alpha (q1)}\right\}^{1/(q1)}}},`$ (6)
$`{\displaystyle \frac{1}{Z_q}}`$ $`=`$ $`\alpha \left\{{\displaystyle \frac{C\beta _02\alpha (q1)}{2\alpha (q1)}}\right\}^{1/2\alpha }{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{1}{q1}\right)}{\mathrm{\Gamma }(\frac{1}{2\alpha })\mathrm{\Gamma }(\frac{1}{q1}\frac{1}{2\alpha })}},`$ (7)
where Tsallis’ nonextensivity parameter $`q`$ is defined by

$$q1+\frac{2\alpha }{\alpha n+1},$$
(8)
which satisfies $`1<q<5/3`$ because of $`0<\alpha 1`$ and $`n>2`$. Since $`q>1`$, the distribution of $`Z`$ exhibits power-law tails for large $`Z`$: $`P_{\mathrm{\Delta }t}(Z)Z^{2\alpha /(q1)}=Z^{(\alpha n+1)}`$. Hence, the $`m`$th moment $`Z^m=Z^mP_{\mathrm{\Delta }t}(Z)𝑑Z`$ converges only for $`m<\alpha n`$. In the limit of $`q1`$, $`P_{\mathrm{\Delta }t}(Z)`$ in Eq. (6) reduces to the canonical distribution of extensive statistical mechanics: $`lim_{q1}P_{\mathrm{\Delta }t}(Z)=P_{\mathrm{\Delta }t}(Z|\beta _0)`$.

## 3 Results and discussions

We have applied the present model to the same FX market data set as used in Ref. gbp96  (provided by Olsen and Associates oa93  which consists of 1 472 241 bid-ask quotes for US dollar-German mark exchange rates during the period October 92 - September 93). The volatility is often estimated by the standard deviation of the price change in an appropriate time window ms00 . Employing this definition, we express the volatility in terms of the local standard deviation of the price change $`Z(t_i)y(t_i+\mathrm{\Delta }t)y(t_i)`$ as

$$V\sqrt{\frac{1}{N}\underset{j=1}{\overset{N}{}}\left\{Z(t_j)\right\}^2}.$$
(9)
Here the window size has been chosen as $`N=35`$. Since $`V^{2\alpha }`$ corresponds to $`\beta `$ (see Eq. (3)) which is $`\chi ^2`$-distributed, $`n`$ can be explicitly obtained from the relative variance of $`\beta `$ using Eq. (5). Thus, there is only one adjustable parameter among $`(\alpha ,q,n)`$, because we have another relation, Eq. (8). In other words, the PDF, $`P_{\mathrm{\Delta }t}(Z)`$ of the price change and the PDF, $`f_{\mathrm{\Delta }t}(\beta )`$ of the inverse power $`V^{2\alpha }\beta `$ of the volatility are determined simultaneously once the value of $`\alpha `$ has been specified.

The PDF’s with time delay $`\mathrm{\Delta }t`$ varying from five minutes up to approximately six hours are displayed in Fig.1 together with theoretical curves obtained from Eq. (6). As the time scale $`\mathrm{\Delta }t`$ increases, $`n`$ increases, while $`\alpha `$ and $`q`$ decrease. The nonextensivity parameter $`q`$ tends to the extensive limit: $`q1`$ as $`\mathrm{\Delta }t`$ increases. Using the same parameter values $`(\alpha ,q,n)`$, the PDF’s of $`\beta `$ are compared in Fig. 2.

The average ‘temperature’ in the market $`1/\beta _0Z^2^\alpha (\mathrm{\Delta }t)^{\alpha \zeta _2}`$ increases with $`\mathrm{\Delta }t`$ since $`\alpha \zeta _2`$ is positive. (We obtained the scaling exponent gbp96 ; fr95  $`\zeta 0.85`$, which is larger than 2/3 obtained for turbulence.) In contrast, the fluctuation of the temperature increases with decreasing $`\mathrm{\Delta }t`$ because the variance of the inverse temperature is proportional to $`2/n`$. The smaller values of $`n`$ imply the stronger intermittency which occurs in small time scales. The intermittent character in the price change can be seen as a fat tail of $`P_{\mathrm{\Delta }t}(Z)`$ in Fig. 1. Also in Fig. 2, the peak of $`f_{\mathrm{\Delta }t}(\beta )`$ shifts to smaller $`\beta `$ as $`\mathrm{\Delta }t`$ decreases (from d to a in Fig. 2), which means relatively high temperatures are realized more frequently in short time scales.

It should be noted that the PDF in Eq. (6) with $`\alpha =1`$ reduces to Student’s $`t`$-distribution ms00 , which has often been used to characterize the fat tails hkk01 . When $`\alpha =1`$, there is no adjustable parameter because $`n`$ is decided from Eq. (5). The FX market system is then subject to a restoring force linear to the price change and the volatility is proportional to the temperature. We have found that the PDFs of $`Z`$ and $`\beta `$ reproduce, although very roughly, the Olsen and Associates’ data points even if $`\alpha `$ is fixed at $`\alpha =1`$. However, adjusting the parameter $`\alpha `$ improved the line shape of $`P_{\mathrm{\Delta }t}(Z)`$ in a range close to $`Z=0`$. The data points in Fig. 1 exhibit a cusp at $`Z=0`$ for large time scales $`\mathrm{\Delta }t`$, which implies a singularity of the second derivative of the PDF at $`Z=0`$. The larger reduction of $`\alpha `$ from unity leads to the stronger singularity. (Note that the factor $`|Z|^{2\alpha 2}`$ arises from $`d^2P_{\mathrm{\Delta }t}(Z)/dZ^2`$.) Thus the better fitting for large $`\mathrm{\Delta }t`$ was obtained from a reduced value of $`\alpha `$. The trend of decrease in $`\alpha `$ with increasing $`\mathrm{\Delta }t`$ was observed for the turbulent flow be01  as well. However, the deviation from $`\alpha =1`$ is much smaller than the present case and the cusp is invisible.

Ghashghaie et al. gbp96  have used a model for turbulence by Castaing et al. cgh90 , in which a log-normal distribution has been assumed for the local standard deviation of the price change. The present model reduces to the model by Ghashghaie et al. if the stochastic process given in Eq. (1) is assumed with $`\alpha =1`$ (the local variance of $`Z`$ is then proportional to $`kT`$) and the $`\chi ^2`$-distribution for $`\beta `$ (the inverse of the local variance of $`Z`$) is replaced by the log-normal distribution for $`\sqrt{kT}`$ (the local standard deviation of $`Z`$). Although no analytic expression like Eq. (6) for $`P_{\mathrm{\Delta }t}(Z)`$ is obtained, a similar qualitative explanation can be applied to their model: The volatility (or equivalently, the square root of the temperature) fluctuates slowly with a log-normal distribution, and the smaller time scale corresponds to the larger variance of the logarithm of the volatility. (The variance is denoted by $`\lambda ^2`$ in Ref. 4.) However, the power-law behavior of the tail of the volatility distribution lgc99  can be better described by the $`\chi ^2`$-distribution for the inverse of the variance.

We have proposed the stochastic process described by Eq. (1) for FX market dynamics in small time scales. In fact, Eq. (1) is the simplest stochastic process which can realize the thermal equilibrium distribution, Eq. (2) in the power-law potential. More realistic and more complicated processes that assure convergence to Eq. (2) at local temperatures might be possible. Mantegna and Stanley have proposed a different stochastic model of the price change, which is described by a truncated Levy flight (TLF) ms95 ; ms00 . The model, yielding approximately a stable distribution, well reproduces the self-similar property of PDF at different time scales: $`\mathrm{\Delta }t=`$ 1 min to 1 000 min. However, the parameters $`\alpha `$ and $`\gamma `$ characterizing the stable distribution fluctuate for larger (monthly) time scales ms00 , where $`\gamma `$ gives a measure of the volatility. In other words, the ultimate distribution (let us denote it $`P_{\mathrm{\Delta }t}^{\mathrm{TLF}}(Z)`$) should be obtained, like Eq. (6), from the weighted average over these parameters. A difference between $`P_{\mathrm{\Delta }t}^{\mathrm{TLF}}(Z)`$ and $`P_{\mathrm{\Delta }t}(Z)`$ is that $`P_{\mathrm{\Delta }t}^{\mathrm{TLF}}(Z)`$ has no cusp at $`Z=0`$: $`P_{\mathrm{\Delta }t}^{\mathrm{TLF}}(Z)1AZ^2`$ for $`|Z|1`$, whereas $`P_{\mathrm{\Delta }t}(Z)1A|Z|^{2\alpha }`$, and the present FX date set indeed exhibits a cusp as seen in Fig. 1.

Finally, a fundamental question has been left open: How to derive theoretically the increasing trend of fat-tailed character with decreasing time scale, i.e., $`\mathrm{\Delta }t`$-dependence of the parameters ($`\alpha ,q,n`$). An attempt to derive the non-gaussian fat-tailed character in small time scales was made by Friedlich et al. recently fpr00 . They derived a multiplicative Langevin equation from a Fokker-Planck equation and showed that the equation becomes more multiplicative and hence fat-tailed as $`\mathrm{\Delta }t`$ decreases. Clarifying the relation between their multiplicative Langevin equation and Eq. (1) (the latter being rather simple although including the fluctuating temperature) should be important as well as interesting.

Acknowledgments

We would like to thank Professor T. Watanabe for valuable comments. The FX data set was provided by Olsen and Associates.


text:
# Thermal components in the early X-ray afterglow of GRBs

## 1 Introduction

Almost all the $`\gamma `$-ray bursts (GRBs) detected so far are associated to a transient X-ray afterglow (Lazzati et al. 2002a). This afterglow is supposed to be due to the early deceleration of the fireball by the interstellar medium (Meszaros & Rees 1997) and its radiation produced by non-thermal synchrotron. Such an interpretation is corroborated by spectral (van Paradijs et al. 2002) and polarimetric (Covino et al. 1999) observations.

More recently, several X-ray afterglows observed by Newton–XMM had a spectrum that can be better fit by an optically thin thermal bremsstrahlung model (Watson et al. 2002) rather than with an absorbed power-law (even though the latter model cannot be unambiguously ruled out on purely statistical grounds). In particular, several high ionization emission lines were detected in the early afterglow of GRB 011211 (Reeves et al. 2002, hereafter R02) which can be accounted for, together with the observed continuum, by a moderately enriched thermal plasma.

The possibility of a thermal origin of the X-ray lines detected in several afterglows was first discussed in Lazzati et al. (1999; see also Vietri et al. 1999 and Kumar & Narayan 2002), who pointed out how a 10 times solar plasma may produce a Fe $`K_\alpha `$ line with the observed large luminosity and equivalent width. In this paper we derive more rigorously the conditions for the environment and the heating mechanism that must be satisfied in order to observe line and continuum emission with the prescribed luminosities, equivalent widths (EWs) and variability time scales.

## 2 Continuum and line emission in thermal plasma

To compute the continuum emission from a thermal optically thin plasma we adopt the treatment and basic equations of Rybicki & Lightman (1979). To compute the luminosities of $`K_\alpha `$ lines from a thermal plasma, instead, we adopt the widely used code MEKAL (Mewe et al. 1985; Liedahl et al. 1995) as implemented in XSPEC (Arnaud 1996). We concentrate in particular the elements Mg, Si, S, Ar, Ca and Fe, for which emission lines have been detected in the afterglow of GRBs. In Fig. 1 we show the line production efficiency for an optically thin thermal plasma as a function of the plasma temperature and for several values of metallicity. We plot $`\eta _{\mathrm{line}}`$, i.e. the ratio of the $`K_\alpha `$ emission line for the six elements above, irrespective of their ionization state, over the total luminosity of the plasma. Gray shading highlights regions in which the equivalent width (EW) of the lines is less than 100 eV, a robust lower limit to any emission feature detected in the afterglow so far.

When using the MEKAL code to evaluate line emission luminosities, one has to remember that the code does not include any radiation transfer, since it assumes that the medium is optically thin to radiation. The actual optical depth of a cloud of plasma at temperature $`T`$ depends on the temperature and on the frequency of the radiation considered. When X-ray continuum radiation is considered, the optically thin approximation can be used up to column densities $`N_H\stackrel{<}{}\sigma _T^11.5\times 10^{24}`$ cm-2. If, however, line emission is concerned, it must be taken into account that intermediate-high $`Z`$ elements retain some electrons which may cause the plasma to be optically thick due to photoionization. In Fig. 2 we show the optical depth of a solar metallicity plasma with $`N_H=1.5\times 10^{24}`$ cm-2 as a function of frequency for a range temperatures between $`10^6`$ and $`10^8`$ K. The opacity of a cold gas is also shown for comparison. In the region of the considered emission lines, the plasma can be optically thick up to temperatures of several keV. This will limit the maximum line luminosity and EW: increasing the column density of the plasma will have no effect on the line luminosity since line photons will be able to escape freely only out of the optically thin surface layer of the medium. Instead of properly introducing radiation transfer in the optically thin MEKAL code, which is a formidable task, in the following we will assume that, after the plasma becomes optically thick to radiation at the line frequency, the line luminosity does not increase if the column density of the gas is increased (see Lazzati et al. 2002b for discussions on a similar assumptions in reflection models).

## 3 Thermal models for GRB afterglows

Consider a cloud of plasma of electron density $`n`$ at a temperature $`T`$ which covers a solid angle $`\mathrm{\Omega }=2\pi (1\mathrm{cos}\theta )`$ along the line of sight to the observer at a distance $`R`$ from the burst explosion site (see Fig. 3). The luminosity that is inferred from infinity, taking into account the light travel time effects, is given by111Here and in the following we indicate a quantity $`Q`$ as $`Q=10^xQ_x`$. We also adopt CGS units.:

$$L1.7\times 10^{23}T_8^{1/2}EI\frac{t_{\mathrm{heat}}+t_{\mathrm{cool}}}{R/c(1\mathrm{cos}\theta )+t_{\mathrm{heat}}+t_{\mathrm{cool}}}$$
(1)
where $`EI`$ is the emission integral ($`EI=_Vn_en_i𝑑V`$), $`t_{\mathrm{cool}}1.2\times 10^{15}T_8^{1/2}n^1`$ s is the plasma cooling time and $`t_{\mathrm{heat}}`$ is the heating time, i.e. the time during which heat (or energy) is supplied to the emitting plasma. It is assumed that during this heating time the temperature is held constant, i.e. the plasma is in equilibrium. In the following we will consider a uniform density plasma and approximate $`n_en_in`$. In this case $`EI=n^2V`$, where V is the volume of the emitting cloud.

In many astrophysically relevant situations (think e.g. to the IGM in galaxy clusters) $`R/ct_{\mathrm{cool}}`$ and both the heating and cooling time scales are irrelevant and Eq. 1 simplifies to the usual free-free equation:

$$L=1.7\times 10^{23}T_8^{1/2}EI$$
(2)
In GRBs this may not be the case. In order to reproduce afterglow observations in the non-cooling regime, the plasma cloud must satisfy simultaneously three conditions. First, it must produce the observed luminosity in the non-cooling regime

$$L=1.7\times 10^{23}T_8^{1/2}n^2\frac{4\pi }{3}R^3\eta _R\frac{\mathrm{\Omega }}{4\pi }10^{46}L_{46}\frac{\mathrm{erg}}{\mathrm{s}}$$
(3)
where the numeric value is typical for the early X-ray afterglow (see e.g. GRB 011211, R02) and $`\eta _R`$ is the volume filling factor of a possible shell-like or clumpy cloud. Second, it must fulfill the non-cooling condition:

$$R(1\mathrm{cos}\theta )ct_{\mathrm{cool}}=3.6\times 10^{25}T_8^{1/2}n^1$$
(4)
finally, it must require a total amount of energy smaller than the total energy of a GRB:

$$E=8.7\times 10^8T_8nR^3\eta _R\frac{\mathrm{\Omega }}{4\pi }10^{52}E_{52}\mathrm{erg}$$
(5)
Conditions 3 and 4 yield the density constraint:

$$n3.7\times 10^8T_8^2\eta _RL_{46}^1\frac{\mathrm{\Omega }}{4\pi }$$
(6)
while conditions 3 and 5 give:

$$n1.2\times 10^9T_8^{1/2}L_{46}E_{52}^1$$
(7)
The above conditions cannot clearly be satisfied simultaneously if all the parameters are taken equal to their fiducial values. Since $`\eta _R`$ and $`\mathrm{\Omega }/4\pi `$ are both numbers smaller than unity, a change in the geometry does not help, making the constraint of Eq. 6 more stringent. Also the temperature cannot be changed significantly, especially as long as X-ray lines must be taken into account.

In conclusion, a thermal component that contributes to the early afterglow of a typical GRB cannot be emitted in the non-cooling regime. In fact, in order to reach a luminosity large enough without involving a too large thermal energy, the plasma density must be so large to make the cooling time very short. As an example, R02 derived a cooling time of $`2`$ s for the afterglow of GRB 011211. We cannot therefore adopt the simple Eq. 2 but we must use the more complete Eq. 1.

We therefore consider in the following two limiting cases: $`R/c(1\mathrm{cos}\theta )t_{\mathrm{cool}}t_{\mathrm{heat}}`$ (flash heating) and $`R/c(1\mathrm{cos}\theta )t_{\mathrm{heat}}t_{\mathrm{cool}}`$ (steady heating). In principle also the case $`t_{\mathrm{heat}}R/c(1\mathrm{cos}\theta )t_{\mathrm{cool}}`$ may be interesting. It is however difficult to identify a heating source that can be active for a time comparable to the line emission time scale. Even if the emitting plasma would be heated by afterglow photons, the heating time could be at most a few per cent of the line emission time scale (Lazzati et al. 2002b).

### 3.1 Flash heating

Consider a portion of a spherical shell of plasma of radius $`R`$, temperature $`T`$ and uniform density $`n`$ covering a solid angle $`\mathrm{\Omega }=1\mathrm{cos}\theta `$ (see, again, Fig. 3). Its emission integral can be written as $`EI=n^2R^3\eta _R\mathrm{\Omega }`$ where $`\eta _R\mathrm{\Delta }R/R`$ is the ratio of the thickness of the shell to its radius222Note that if a clumpy medium is assumed instead of a spherical shell, $`\eta _R`$ is equal to the volume filling factor of the configuration that maximizes the emitted luminosity.. Considering the values $`\eta _{\mathrm{line}}`$ defined above (and shown in Fig. 1), the cloud will produce a line with luminosity333The continuum luminosity can be easily derived by adopting $`\eta _{\mathrm{line}}=1`$. (see Eq. 1):

$$L_{\mathrm{line}}=1.9\times 10^{27}\eta _{\mathrm{line}}RT_8\tau _T$$
(8)
where we have used the Thomson opacity $`\tau _T=nR\eta _R\sigma _T`$. Eq. 8 coupled with the efficiencies shown in Fig. 1 implies that when we observe a line with luminosity $`L_{\mathrm{line}}\stackrel{>}{}5\times 10^{44}`$ erg s-1 (Ghisellini et al. 2002) a thermal model in the flash heating regime requires a radius larger than $`R\stackrel{>}{}2\times 10^{17}`$ cm. Analogous results are obtained by imposing a ten times larger continuum luminosity. Such a large radius implies that the line luminosity should remain constant for a time $`t=R(1\mathrm{cos}\theta )/c`$, to be compared with the observed variability time scales $`t_{\mathrm{var}}2\times 10^4`$ s observed in e.g. GRB 011211 (R02). This constraint can be satisfied if the heated plasma covers a very small solid angle, of opening angle $`\theta 4.5^{}`$, at the limit but consistent with the smallest opening angles inferred for GRB jets from afterglow modelling (Frail et al. 2001; Panaitescu & Kumar 2002). In order to satisfy the cooling condition, however, the plasma must be dense ($`n\stackrel{>}{}10^{11}`$ cm-3) and therefore confined in a very thin shell (or small clumps) with $`\eta _R\stackrel{<}{}5\times 10^5`$. In addition, if we assume that the plasma is due to a spherical shell surrounding the GRB explosion site which is heated only in the small polar cap by afterglow photons, the shell mass would turn out to be $`M=4\pi R^2m_p\tau _T/\sigma _T600M_{}`$ for the case explored above. Both the geometry factor $`\eta _R`$ and the total mass required seem to be rather extreme, even though the whole picture is, in this case, completely self consistent.

### 3.2 Steady heating

Let us now consider a plasma in which a source of heating is active for a time scale $`t_{\mathrm{heat}}>1.2\times 10^{15}T_8/n=t_{\mathrm{cool}}`$. In case of radiative heating, this time scale may be for example the GRB duration or even longer, if the afterglow photons can contribute to the heating. In steady heating, the line luminosity (continuum luminosity for $`\eta _{\mathrm{line}}=1`$) can be written as:

$$L_{\mathrm{line}}=7.2\times 10^{36}\tau _T^2T_8^{1/2}\eta _{\mathrm{line}}\frac{t_{\mathrm{heat}}}{\eta _R}c_h\frac{t_{\mathrm{heat}}}{\eta _R}$$
(9)
where we have defined the parameter $`c_h`$ in order to emphasize the ratio of the heating time scale over the geometry parameter. In Fig. 4 we show the highest possible values of the coefficient $`c_h`$ as a function of temperature for five different values of metallicity (as in Fig. 1). The lines are plotted only when the equivalent width of the line is predicted to be larger than 100 eV, corresponding to the unshaded areas in Fig. 1.

To compare these results with afterglow data, consider a Fe line with $`L=10^{45}`$ erg s-1 (GRB 991216; Piro et al. 2000) or the S line in GRB 011211 (R02) with luminosity $`L=4\times 10^{44}`$ erg s-1. In both cases a ratio $`t_{\mathrm{heat}}/\eta _R\stackrel{>}{}10^9`$ is required, by comparison with the appropriate panel in Fig. 4. A typical heating time can be considered to be the fireball transit time, which is equivalent to the total duration of the prompt GRB emission. In both cases such duration is of the order of $`100`$ s, leaving us with an extreme requirement of anisotropy $`\eta _R\stackrel{<}{}10^7`$. Alternatively, the heating may be provided radiatively by the absorption of GRB and early afterglow photons. In that case, the duration of the heat supply may be of several per cent of the line emission time scale (Lazzati et al. 2002b). In that case the constraint would be relaxed to the (yet still challenging) value $`\eta _R\stackrel{<}{}10^6`$. The most convenient solution is represented by an (unknown) heating mechanism acting for a time comparable to the line emission time scale itself. In this case, the constrain on the geometric factor $`\eta _R`$ would be similar to the one derived above for the flash heating case. Such a heating mechanism is however presently obscure and will have to face the problem of stability discussed below (§4).

In a steady heating scenario, therefore, thermal lines and continua can dominate the early afterglow emission, only in case of extreme geometric conditions, in which the emission is produced either in a sheet like shell or an extremely clumpy medium. Interestingly, analogous extreme conditions were independently inferred for the environment of GRB 000210 (Piro et al. 2002) in order to account for the lack of ionization features in the soft X-ray afterglow spectrum. Piro et al. (2002, and references therein) argue also that such conditions may be realized in giant molecular clouds. Stability considerations (see § 4), however, show that more extreme conditions are required in this case. Differently from what derived in the flash heating condition, in this case the radius of the shell does not have to be very large, and therefore the total mass required is not huge.

## 4 Stability

In the above section we have derived some geometrical constraints to the emitting plasma in order to reproduce the observed features and continua. Since the results envisage particular geometric conditions, we here analyze their stability. First, since the heating energy is supposed to come from a central point (the GRB progenitor), one has to make sure that the emitting material is not accelerated to high velocities.

Consider a shell of mass $`M`$ absorbing energy from a relativistic outflow. If the material is radiative, it acquires a bulk velocity $`v=E/(Mc)`$. A bulk velocity $`v0.1c`$ was measured in GRB 011211 (R02) and GRB 991216 (Piro et al. 2000). Requiring that our optically thin shell is accelerated to a comparable or smaller speed implies a radius larger than:

$$R\frac{L\sigma _T}{4\pi m_pc^2\eta _{\mathrm{line}}\tau _Tv}=3\times 10^{13}\frac{L_{\mathrm{line},45}}{\eta _{\mathrm{line}}\tau _Tv_9}\mathrm{cm}$$
(10)
This is not a compelling limit, given the radii discussed above.

Provided that the emitting medium is not accelerated to relativistic speeds by the energy input, we also want that the thin emitting shell (or blobs) do not expand in a time scale smaller than the emission one (which can be either the heating or cooling time scale). The shell (or blobs) was in fact in equilibrium with the ambient medium when it was cold. Now that its temperature is increased it will tend to expand under the effect of the increased internal pressure. If it expands at the speed of sound $`c_s`$, its density will be sizably modified in a time scale $`t_{\mathrm{exp}}=R\eta _R/c_s`$. It is therefore required that the expansion time is longer than the emission time scale.

In the case of steady heating we obtain:

$$R>c_s\frac{t_{\mathrm{heat}}}{\eta _R}=10^8T_8^{1/2}\frac{t_{\mathrm{heat}}}{\eta _R}\stackrel{>}{}10^{17}T_8^{1/2}$$
(11)
Also in the steady heating case, therefore, the radius has to be large in order to allow for the production of thermal components in the early GRB X-ray afterglows. For such large radii, however, the density required to fulfill the steady heating conditions is $`n\stackrel{>}{}10^{15}/t_{\mathrm{heat}}\stackrel{>}{}10^{12}`$, several orders of magnitudes larger than what inferred by Piro et al. (2002). It seems therefore that the steady heating case requires more extreme conditions than the flash heating one (the same distance from the explosion site and limits on the total mass involved, but larger densities and smaller filling factors). Applying the same stability condition to a flash heating case, we are left with the more relaxed constraint

$$\tau _T>0.1$$
(12)
which is therefore not difficult to fulfill. These stability considerations suggest therefore that the only viable way to produce a sizable thermal afterglow component is by heating a small portion of a massive $`\tau _T1`$ (either very thin or clumpy) shell of material located at a relatively large distance from the burst explosion site (Fig. 3). It should however be emphasized that these stability considerations cannot be applied if the heating is provided hydrodynamically, in such a way that the source of heating is providing also the confining pressure. In this case also a steady heating scenario, with a much smaller shell, may be viable.

## 5 Clumpiness of GRB ambient medium

The above discussion shows that the emission lines detected in several GRB afterglows can have a thermal origin only if the GRB surroundings are extremely clumpy, with density and geometric contrasts of order of a hundred thousand or more. We now discuss the reliability of these conditions in different GRB progenitor scenario.

In the hypernova scenario, the GRB sets on simultaneously to a supernova explosion, and the ambient medium is the result of the interaction of the pre-SN star with its surroundings. In particular, the nearby ambient will be dominated by the late stages of the mass ejection history of the star. These stages are known to be unsteady and clumps or shell-like structures can be envisaged. If the mass ejection is caused by radiative effects, however, extreme structure cannot be produced, since the stellar luminosity varies on time scales comparable to the Kelvin time-scale, which is of the order of hundreds of years, yielding a thick shell. In the wind environment of SN1998bw, for example, Li & Chevalier (1999) found that inhomogeneities up to a factor of a few were present. A more appealing scenario, in this perspective, is the supranova model by Vietri and Stella (1998), where the GRB is supposed to explode several weeks to years after a supernova. The supernova explosion, having a much smaller time scale, can generate a more extreme geometry. The clumpiness of SN ejecta has been investigated by numerical simulations. It is found that Rayleigh-Taylor instabilities can produce high density clumps with angular scales $`\stackrel{>}{}1^{}`$ (Böttcher et al. 2002). Such structures have $`\eta _Ra/R10^2`$ where $`a`$ is the clump radius, much larger than the value required to produce sizable line emission from collisional excitation. It is however possible, in the supranova scenario, that the remnant is illuminated by a super-Eddington relativistic wind (Vietri & Stella 1998, Konigl & Granot 2002) in the time span between the SN and GRB explosions. The interaction of this wind with the ejected SN shell may increase the inhomogeneities originally present in the shell (Guetta & Granot 2003; Lazzati & Rees in preparation).

## 6 Discussion

In at least three XMM-Newton spectra of X-ray GRB afterglows, the presence of a thermal component has been claimed (R02; Watson et al. 2002). We have studied the geometrical conditions that the thermal material must satisfy in order to contribute significantly to the early X-ray afterglow. We first concluded that the time scale during which the emission is observed must be set by the light crossing time of the emitting region. This imply that standard free-free equations that relate the emission integral to the luminosity cannot be used and a more general equation (Eq. 1) has to be adopted. The implications of this are of great importance. In fact, the use of the standard formalism (Eq. 2) to compute the EI led R02 to underestimate the emission integral for GRB 011211 by a factor $`10^4`$, and conclude that thermal emission from a shell with $`R10^{15}`$ cm can be a self consistent solution. We showed that, with the correct treatment of light-crossing effects, this is not the case. We then studied the possibility of (i) flash heating and (ii) steady heating. In both cases, we conclude that the thermal material must be extremely clumped in order to contribute significantly to a typical X-ray afterglow. In the steady heating case, the thermal material is continuously heated for a time longer than the cooling one, reaching a stationary state (albeit for a time smaller than the light-crossing one). We showed that in this case the material has to be so extremely clumped that the clumps are dissolved by the increased temperature in a very short time scale, so that an additional source of confinement must be envisaged. In the case of flash heating, a self consistent solution can be found, requiring a less extreme (even if yet compelling) clumping. In this case, the thermal material must be located far from the source and therefore, in order to preserve the fast time variability, only a small portion of it has to be heated. Assuming that the same density is spread all over the GRB explosion site, we derive that the total mass involved is large ($`600M_{}`$), larger than even a type-II supernova remnant. A moderate asymmetry of the remnant could however ameliorate this requirement. It is interesting to note that, even though such clumpy and dense environments seem extreme, they were already inferred independently to account for the lack of ionization features in GRB 000210 (Piro et al. 2002) and to explain the haigh energy emission in GRB 940217 (Katz 1994).

To date there is not yet a compelling evidence that there are indeed thermal components in the early X-ray afterglow of GRBs. It is however intriguing then in three out of four afterglows observed by XMM-Newton a thermal model yields to a better fit than an absorbed power-law (R02, Watson et al. 2002a). In the case of GRB 011211 (R02) is particularly difficult to explain the absence of an iron or nickel line, while reflection models seem to give a better explanation to the observed line ratios (Lazzati et al. 2002b).

Nonetheless, should a thermal component be confirmed in the X-ray afterglow of at least a sub-fraction of GRBs, it would be difficult to avoid the conclusion that (i) the GRB explosion site is surrounded by a massive shell of matter, similar to a 1–2 years old supernova remnant (as expected in some versions of the Supranova model, Vietri & Stella 1998); (ii) that this remnant is extremely clumped, suggesting that some precursor activity has been taking place. Whether such a clumpiness can help making the prompt emission in an external shock scenario (Dermer et al. 2000) heavily depends on the presence of clumping at smaller scales, a test that cannot be done with thermal emission from the clumps, since would contribute in any case a negligible flux to the early afterglow emission.

###### Acknowledgements.

I am grateful to G. Ghisellini, P. Kumar, M. J. Rees and F. Tavecchio for useful comments and discussions. I acknowledge financial support from the PPARC.


text:
# Isolating Clusters with Wolf-Rayet Stars in I Zw 181

## 1. INTRODUCTION

I Zw 18 (Mrk 116) represents an important local fiducial in our understanding of galaxy formation and evolution. It is a blue compact galaxy with the smallest known abundance of heavy elements in its ionized gas ($``$2% of solar; V$`\stackrel{´}{\mathrm{i}}`$lchez & Isglesias-P$`\stackrel{´}{\mathrm{a}}`$ramo 1998 and references therein). For the past three decades, the debate has focused on the age of its stellar populations. Early analyses of Wide Field Planetary Camera 2 (WFPC2) images suggested that star formation began several tens of Myr ago (Dufour et al. 1996; Hunter & Thronson 1995), while later analyses of the same WFPC2 images (Aloisi, Tosi, & Greggio et al. 1999) and new images from the Near Infrared Camera and Multiobject Spectrometer ($`\ddot{\mathrm{O}}`$stlin 2000) implied the presence of an older population (0.5–1 Gyr) from a small number of faint red stars.

Two independent sets of ground-based observations have recently discovered Wolf-Rayet (WR) emission within the brighter (NW) component of I Zw 18: Izotov et al. (1997) found evidence for 17 WN and 5 WC stars, while Legrand et al. (1997) found evidence for 1–2 WC stars. While neither group had the spatial information to tie this emission to stellar clusters, the existence of WR stars in I Zw 18, especially of subtype WC, is surprising, because canonical (non-rotating and non-binary) evolutionary theory predicts few WN stars and no WC stars to form at the low metallicity of I Zw 18 (Cervi$`\stackrel{~}{\mathrm{n}}`$o & Mas-Hesse 1994; Meynet 1995). de Mello et al. (1998) subsequently searched for these WR stars in narrow-band WFPC2 images, and found faint peaks in the He II $`\lambda `$4686 emission they attributed to either nebular emission or 5–9 WN stars.

Using the Space Telescope Imaging Spectrograph (STIS) on board the Hubble Space Telescope (HST), we recently obtained a two-dimensional far-UV spectroscopic map of I Zw 18, along with far-UV and near-UV images, in

1Based on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. These observations are associated with proposal 9054.

order to further understand its star-formation history. These data reveal two clusters containing WC stars. While neither of these clusters is associated with the WR detections of Legrand et al. (1997) and de Mello et al. (1998), one of the clusters might be responsible for the WC emission seen by Izotov et al. (1997). Our STIS data also provide a high-resolution H I map, based on the Ly-$`\alpha `$ absorption; we show that the H I column depth ($`N_{HI}`$) is sharply peaked within the SE component.

## 2. OBSERVATIONS

We obtained far-UV spectra of I Zw 18 from 2002 Feb 9 to 2002 Feb 12, using the G140L grating and the $`52\mathrm{}\times 0.5\mathrm{}`$ slit on STIS (for a full description of the instrument, see Woodgate et al. 1998 and Kimble et al. 1998). The field-of-view of the STIS UV detectors is $`25\mathrm{}\times 25\mathrm{}`$. We placed the slit at seven adjacent positions in I Zw 18, parallel to the axis of the main body (position angle = 145o), resulting in spatially resolved spectra of its bright star-forming regions (Figure 1). The exposure time at four positions was 6058 s (the central 3 positions and the SW extreme), while at the remaining three positions it was 5376 s. At each slit position, half the exposures were offset from the others by $`0.5\mathrm{}`$ (20 pixels) along the slit, in order to allow masking of bad, hot, and shadowed pixels in the coaddition, and to smooth over variations in the detector response. The spectra were processed through the standard calibration pipeline, except that we also subtracted an appropriately scaled profile of the dark rate, from a sum of 500 ks of dark exposures. The subtraction of this dark “glow” (which peaks in the upper left-hand quadrant and can be 20 times higher than the nominal $`7\times 10^6`$ cts s-1 pix-1 dark rate) is important when studying faint extended sources. Finally, we subtracted the airglow spectrum using a region in each spectral image free from source flux.

On 2002 Feb 25, we also obtained far-UV and near-UV images of I Zw 18 using the F25SRF2 and F25QTZ filters, respectively, with bandpasses centered at 1457 Å and 2365 Å, and exposure times of 5331 s and 5786 s. We again used $`0.5\mathrm{}`$ offsets to allow masking of pixels and smoothing of detector response. All of the frames for each camera were coadded, rescaled, and shifted with the DRIZZLE package in IRAF, in order to match the G140L spectral images described above.

In addition to our STIS data, we also obtained archival WFPC2 images at longer wavelengths, in the F555W ($`V`$), F439W ($`B`$), F336W, and F469N (He II $`\lambda 4686`$) filters, for comparison with earlier work. The WFPC2 data were cleaned of cosmic rays, coadded, and drizzled to register with the STIS data. Figure 1 shows a false color image produced from our STIS data in the blue and green channels, and the WFPC2/F555W data in the red channel.

Fig. 1– Far-UV spectra and broad-band images of I Zw 18. Panel a: The spectrum with the broadest Ly-$`\alpha `$ absorption and thus the highest $`N_{HI}`$, with (black) and without (blue) airglow subtraction. Panels b and e: Spectra of regions with previously reported WN detections (de Mello et al. 1998); we find no evidence for WN stars. Panels c and d: Spectra of clusters with WC stars; note the C IV and He II emission lines. Panel f: Composite HST image of I Zw 18 in 3 bandpasses (blue = STIS/FUV/F25SRF2, green = STIS/NUV/F25QTZ, red = WFPC2/F555W). The seven adjacent positions of the STIS slit are marked (white lines). Many of the old stars reported by Aloisi et al. (1999) appear as faint red stars here, but some are UV-bright, indicating a hotter temperature or a hot companion. Yellow circles mark two regions where we find possible narrow He II emission but no associated C IV emission. The de Mello et al. (1998) WN detections are encircled in red.

## 3. ANALYSIS

### 3.1. WC Stars

The CALSTIS pipeline software produces two-dimensional rectified spectral images that are useful for examining the spectral energy distributions (SEDs) of extended objects. These images are linear in wavelength along the x-axis, and linear in spatial extent along the y-axis. Upon examination of these images, we saw a compact source with obvious WR/WC emission

features at C IV $`\lambda \lambda 1548,1551`$ and He II $`\lambda 1640`$ (shown in Figure 1$`c`$). We then performed a systematic search for other such objects in the data.

To find the locations of individual spectra in these images, we collapsed the spectra along the dispersion axis and convolved with a Gaussian normalized to unit area and zero sum (similar to the finding algorithm of the DAOPHOT package; Stetson 1987). The local maxima in these convolutions gave the y-positions of compact sources along the slit. Returning to the original spectral images, we extracted spectra at each of these positions, using row sums 5 pixels high. We also extracted spectra from the two regions with significant WR/WN detections in previous analyses of WFPC2 data (de Mello et al. 1998; those not marked with “?” in their Table 3), although these regions are not associated with compact sources in our spectral images, or with sources in the UV and optical HST images (compare their Figures 2 and 3$`a`$). Each spectrum was then shifted to account for the position of the object within the slit, using a centroid on the corresponding region in the far-UV image.

After extracting 205 spectra at the local maxima in our spectral images, we performed a systematic search for broad He II $`\lambda 1640`$ emission, using as a guide the spectrum of the source with obvious WR/WC features discussed above. We fit a linear continuum to each spectrum, using spectral regions free of strong emission and absorption features (1400–1500 Å, 1590–1620 Å, and 1670–1710 Å), summed the net flux exceeding this continuum in the region 1635–1665 Å (the heliocentric velocity of I Zw 18 is 751 km s-1; Thuan et al. 1999), and returned a list of $`>3\sigma `$ detections of excess flux near 1640 Å. The automated search detected the obvious source we had already found by eye, plus one additional source of WR emission; both spectra also show broad C IV emission, indicating the presence of WC stars (Figure 1, panels $`c`$ and $`d`$). Note that the emission lines are much broader than a resolution element in our spectra ($`1`$ Å for point sources and $`12`$ Å for sources filling the slit). We also found two spectra with narrow emission features near 1640 Å that could be noise or nebular emission; the spectra do not have C IV emission, and their associated clusters are marked in Figure 1$`f`$ (yellow circles).

Our WC detections are located $`1.3\mathrm{}`$ SE and $`0.5\mathrm{}`$ W of the center of the NW starforming region. Legrand et al. (1997) report WC emission $`1.5\mathrm{}`$ SW of the central cluster, so our detections are distinct. Izotov et al. (1997) only restrict their detections of WR stars to the NW starforming region. None of the detections of WN stars by de Mello et al. (1998) correspond to the positions of our WC stars; in Figure 1$`b`$ and 1$`e`$, we show spectra from the regions of their putative detections.

To estimate the number of WR/WC stars needed to produce the observed He II emission, we used the tabulations of Schaerer & Vacca (1998). In the middle of the range of WC stars, WC5 stars produce an average of He II luminosity of $`5\times 10^{36}`$ erg s-1. Assuming a distance of 12.6 Mpc ($`\ddot{\mathrm{O}}`$stlin 2000) and a foreground extinction of $`E(BV)=0.03`$ mag (Schlegel, Finkbeiner, & Davis 1998) but no reddening from the patchy internal dust (Cannon et al. 2002), we can convert our He II fluxes to luminosity. For the spectra shown in Figure 1$`c`$ and 1$`d`$, we measure $`L_{1640}=(2.4\pm 0.5)\times 10^{37}`$ erg s-1 and $`L_{1640}=(3.2\pm 1.0)\times 10^{37}`$ erg s-1, respectively, corresponding to 5$`\pm 1`$ and 6$`\pm 2`$ WR/WC stars. Using the 1500 Å continuum flux and a 5 Myr isochrone (Bertelli et al. 1994), we estimate that these clusters contain 26 and 119 stars, respectively, with with $`M>`$$`M_{}`$, suggesting that the ratio WR/(WR+O) in these clusters is respectively 0.2 and 0.05. The stronger WR/WC detection in Figure 1$`c`$ is due to the higher fraction of WR/WC stars in that cluster; note the prominent C IV and Si IV absorption in Figure 1$`d`$, where the hot underlying population of O stars is more dominant. Our measured WR/(WR+O) ratios fall within the range observed in other galaxies, but are unusually high given the metallicity of I Zw 18 (see Maeder & Meynet 1994). Of course, the ratio is far lower for the galaxy as a whole, given the many UV-bright clusters showing no WR emission.

The ratio of C IV to He II emission in Figure 1$`c`$ is $``$2.3, which is typical for the WC subtype (Niedzielski & Rochowicz 1994; Gr$`\ddot{\mathrm{a}}`$fener, Koesterke, & Hamann 2002). However, this ratio is only $`0.7`$ in Figure 1$`d`$; the smaller ratio is partly due to the C IV absorption from the hot stellar population and interstellar medium, but it could also indicate that one or two of the WR stars in this cluster are of type WN instead of WC.

de Mello et al. (1998) claimed detections of 1.5 and 3 WNL stars, respectively, where we are showing the spectra in Figure 1$`b`$ and 1$`e`$. Assuming $`L_{1640}=1.6\times 10^{36}`$ erg s-1 (Schaerer & Vacca 1996), individual WNL stars would be difficult to detect in our data, but we do not see any signs of these stars. Where de Mello et al. (1998) found 1.5 WNL and 3 WNL stars, we formally find $`0.5\pm 3`$ WNL stars and $`+0.5\pm 1.6`$ WNL stars, respectively. Note that there are no point sources in our UV images (or in the extant WFPC2 images) at the locations of the de Mello et al. detections, and de Mello et al. note that their detections are consistent with both nebular or WR emission. Unfortunately, their F469N image has a hot pixel exactly where we have our strongest detection of WR/WC emission, which probably precluded their detection of these sources.

### 3.2. Old Clusters or Large Interstellar Hydrogen Column?

Compared to the spectra in the NW half of I Zw 18, our spectra in the SE half generally show broader Ly-$`\alpha `$ $`\lambda `$1216 absorption, but it is especially broad in one particular spectrum (see Figure 1$`a`$). Without other data, the absorption could arise from two situations: either a high column density of interstellar H I, or a predominately older ($`500`$ Myr) and cooler population of stars. We investigated this latter possibility by examining the SED with a longer wavelength baseline than that available in our far-UV spectra. Note that if the spectrum in Figure 1$`a`$ came from a 500 Myr population, it would be have to be $``$60 times more massive than the cluster shown Figure 1$`d`$.

Using the HST images in the far-UV, near-UV, $`B`$, and $`V`$ bands, we measured the flux in the region of our spectral extraction (an area $`0.5\mathrm{}\times 0.12\mathrm{}`$). This flux comes from a string of compact sources aligned with the long axis of the extraction box, so we estimated the aperture correction required by using TinyTim models of the HST PSF in the different bandpasses. We do not expect perfect agreement with the G140L spectrum (due to wavelength-dependent aperture correction uncertainties in the images and spectrum), but they agree at the $`10`$% level (see Figure 2). Furthermore, it is difficult to quantify exactly the dominant effective temperature of the stellar population, given the reddening uncertainties and the fact that we are looking in the Rayleigh-Jeans tail. However, the SED of this region is clearly dominated by a population of hot stars. We show this by integrating the synthetic spectra of Kurucz (1993) over the isochrones of Bertelli et al. (1994), for two different ages (Figure 2). A population of about 500 Myr (or a single star of $`\mathrm{T}_{\mathrm{eff}}12,000`$ K) is needed to explain the width of Ly-$`\alpha `$ absorption without a significant contribution from interstellar H I. However, such a population would produce far more flux at longer wavelengths than is actually seen. Instead, a much younger population (e.g., 16 Myr) reddened by $`E(BV)=0.06`$ mag (SMC extinction law; Witt & Gordon 2000) can reproduce the spectral energy distribution from the far-UV to the optical, but this requires that the Ly-$`\alpha `$ absorption come from a high H I column: $`N_{HI}2\times 10^{22}`$ atoms cm-2. This result is independent of the dominant $`\mathrm{T}_{\mathrm{eff}}`$ in the underlying population, for $`\mathrm{T}_{\mathrm{eff}}15,000`$ K. The spectra show much narrower Ly-$`\alpha `$ absorption as one moves away from this position; e.g., looking along a sightline to a bright cluster 1$`\mathrm{}`$ to the SE, the H I column must be at least ten times smaller (with the exact amount depending upon the assumed age/temperature of the underlying population).

van Zee et al. (1998) produced high spatial resolution (5$`\mathrm{}`$) maps of the H I distribution in I Zw 18, using the Very Large Array (VLA), and found that the column density peaks in the SE half of the galaxy at $`3\times 10^{21}`$ atoms cm-2. Because each of the two starbursting regions (NW and SE) subtends $`3\mathrm{}`$ on the sky, a sharp peak in the H I distribution would be much smaller than a VLA resolution element, and thus the column density could peak at a significantly larger value than the peak in the van Zee et al. (1998) map.

### 3.3. Summary

We have presented high spatial resolution images and spectra of the starburst galaxy I Zw 18. These spectra show strong WR/WC emission features originating in two compact clusters. From the He II $`\lambda `$1640 emission, we estimate that 6 WR/WC stars reside in a very UV-bright cluster situated well within the NW starbursting region, and that 5 WR/WC stars lie in a somewhat fainter cluster just outside of this region. The fact that any WR/WC stars exist at all in I Zw 18 is surprising, because canonical theory for the evolution of single, massive, non-rotating stars predicts few WR stars and no WC stars to form at the low metallicity of I Zw 18 (Meynet 1995; Cervi$`\stackrel{~}{\mathrm{n}}`$o & Mas-Hesse 1994).

The existence of WR/WC stars at the low metallicity of I Zw 18 demonstrates the importance of sophisticated models of massive star evolution. Maeder & Meynet (2001) find that the inclusion of rotation in their models decreases the minimum mass for a single star to become a WR star, allowing the formation of WR stars at lower metallicity. Furthermore, there is both theoretical and observational evidence for the importance of close binary evolution in the formation of WR stars, especially as one moves to lower metallicity (see Bartzakos, Moffat, & Niemela 2001); the role of binaries in WR formation has a long history (see Paczy$`\stackrel{´}{\mathrm{n}}`$ski 1967), even though much of the WR literature focuses on single-star evolution.

Our data also show a sharp peak in the H I column density of I Zw 18, situated within the SE half of the galaxy. This peak falls within the peak of previous high spatial resolution maps of the H I taken with the VLA (van Zee et al. 1998). However, our measurement of 2$`\times 10^{22}`$ atoms cm-2 exceeds the earlier estimates by an order of magnitude, because of the higher resolution in our HST data.

Support for proposal 9054 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. The authors are grateful to C. Leitherer for advice and useful discussions. A. Aloisi kindly provided her catalog of WFPC2 photometry.
Fig. 2– Our spectrum showing the broadest Ly-$`\alpha `$ absorption (see Figure 1$`a`$), here over a longer wavelength baseline, combining the G140L data (black curve) and HST images (circles). The flux decline at longer wavelengths and the small Balmer jump show clearly that this region is dominated by a population of hot stars. To demonstrate, we show the integrated flux from a lightly reddened 16 Myr population (light gray curve) and the flux from an unreddened 500 Myr population (dark gray curve). An H I column of $`2\times 10^{22}`$ atoms cm-2 is required (not shown for clarity) to produce the width of the Ly-$`\alpha `$ absorption if a young (hot) population of stars is assumed.


text:
# The Infrared Spectrum of R Doradus

## 1. Observations and analysis

Our goal is to discern whether a moderately varying M giant, such as R Doradus, could be modelled satisfactory with a hydrostatic model photosphere. We observed the $`2.63.7\mu \text{m}`$ region with the ISO SWS 06 at a resolution of $`R=20002500`$ and the full $`2.445\mu \text{m}`$ region with the SWS 01 at $`R=30`$. The low-resolution SWS01 observation is used to ensure that the different SWS06 observations in this region were properly aligned.

We have calculated hydrostatic model photospheres with the marcs code in spherical geometry. Synthetic spectra were calculated using these model photospheres. The synthetic spectra, also calculated in spherical geometry, include lines from H2O (Partridge & Schwenke 1997), CO (Goorwitch 1994), SiO (Langhoff & Bauschlicher 1993), CH (Jørgensen et al. 1996), CN (Jørgensen & Larsson 1990 and Plez 1998, priv. com.), OH (Goldman et al. 1998), C2 (Querci et al. 1971 and Jørgensen 2001, priv. com.), and CO2 (Hitemp; Rothman et al. 1992).

We are able to fit the observed range well with a hydrostatic model photosphere, see Figs. 1 and 2. Our best fit model have the fundamental parameters of $`T_{eff}=3000`$ K, $`\mathrm{log}(g)=0`$, and solar metallicity. We find that the wavelength region investigated is sensitive to the effective temperature used in the modeling of the photosphere. The agreement between the synthetic spectrum and the ISO observations is encouraging, especially in the wavelength region of $`2.83.7\mu \text{m}`$, suggesting that a hydrostatic model photosphere is adequate for the calculation of synthetic spectra in the near infrared for this moderately varying red giant star. However, an additional absorption component is needed at $`2.62.8\mu \text{m}`$. Thus, we also find that the model does not account for this extra absorption component at the beginning of the interval. For a discussion of this discrepancy see Ryde & Eriksson (2002). The spectral signatures are dominated by water vapour in the stellar photosphere, but several photospheric OH, CO, and SiO features are also present.

Figs. $`\text{3}\text{9}`$ show which molecules are important contributors to the emergent spectrum (up to $`5\mu \text{m}`$) of R Doradus. The synthetic spectra shown there are all calculated using our R Dor model atmosphere. However, the spectra are synthesized with only one molecule at a time.

## 2. Concluding remarks

R Dor is a moderately varying M Giant. We have modelled its 3 micron ISO spectrum (resolving power is ca 2000) successfully with a hydrostatic model photosphere. Dynamical models are, however, certainly needed for the Mira stars, which show larger variations with time. Line opacities contributing at $`2.63.7`$ microns are H2O, CO, OH, and SiO. The region is found to be very temperature sensitive, and it is found that an extra absorption is needed at 2.6 - 2.8 microns.

### Acknowledgments.

We should like to thank Drs. B. Gustafsson and D. L. Lambert for inspiration and enlightening discussions. This work was supported by the P.E. Lindahl Foundation Fund of the Royal Swedish Academy of Sciences and the Swedish Foundation for International Cooperation in Research and Higher Education.

## References

Goldman et al. 1998, JQSRT 59, 453

Goorwitch 1994, ApJS 95, 535

Jørgensen & Larsson 1990, A&A238, 424

Jørgensen et al. 1996, A&A 315, 204

Langhoff & Bauschlicher 1993, Chem. Phys. L. 211, 305

Partridge & Schwenke 1997, J. Chem. Phys. 106, 4618

Querci et al. 1971, A&A 15, 256

Rothman et al. 1992, JQSRT 48, 469

Ryde & Eriksson 2002, A&A 386, 874