text:
# Source of Atomic Hydrogen in the Atmosphere of HD 209458b

## 1 INTRODUCTION

Since the discovery of the first extrasolar planet, 51 Peg b, in 1995 (Mayor & Queloz, 1995), a total of 102 planets have so far been discovered (e.g., Butler et al. 2003; Udry et al. 2002, and references therein) and analyzed statistically in order to characterize the formation environment (Fischer et al., 2002; Santos et al., 2003). The formation of gas giants is thought to be complete in 10 Myr, before the disappearance of the gaseous stellar accretion disk, at distances $`>5`$ AU from the parent star. They are then pulled to their present positions by tidal interaction between the gas disk and planet (e.g., Pollack et al. 1996; Ward 1997).

An edge-on planet provides a unique opportunity to investigate the planetary atmosphere. HD 209458b is such a planet, providing the first extrasolar planetary detection using the light curve obtained during a planetary transit of its parent star (Charbonneau et al., 2000; Henry et al., 2000). The orbital parameters were accurately determined by Charbonneau et al. (2000), Henry et al. (2000), and Mazeh et al. (2000). Strong absorption lines are required to make an atmospheric detection and Seager & Sasselov (2000) theoretically characterized the most prominent absorption features, viz., Na I and K I doublet resonance and He I $`2^3S2^3P`$ triplets. Charbonneau et al. (2002) detected the Na I doublet at 589.3 nm in HD 209458b with $`4\sigma `$ confidence level. Following this, Vidal-Madjar et al. (2003) made the first observation of the extended upper atmosphere of HD 209458b with a $`4\sigma `$ detection of the H I atomic hydrogen absorption of the stellar Lyman-$`\alpha `$ line. They reported an absorption of $`15\pm 4`$% and claimed this should be taking place beyond the Roche limit, thus implying hydrodynamic escape of hydrogen atoms from HD 209458b’s atmosphere.

The temperature and UV flux of close-in planets are high. This motivates us to study the chemistry that may be important in this ”hot jupiter”. In this paper, we consider a simple hydrocarbon/oxygen chemistry model to determine a source of atomic hydrogen in the atmosphere of HD 209458b and represents the first effort to investigate the UV enhanced chemical processes in ”hot jupiters”.

## 2 MODEL

Our model is based on the four parent molecules H2, CO, H2O, and CH4 and is a derivative of the Caltech/JPL KINETICS model for the Jovian atmosphere. HD 209458b is orbiting at a distance of 0.05 AU. As HD 209458 is a G0 solar-type dwarf star, it is justified to use the solar spectrum. We expect the atmosphere to have a temperature, $`1000`$ K, and UV flux ($`<1800`$ Å), $`2\times 10^{15}`$ photons cm-2 s-1. By comparison the UV flux at Jupiter is a factor of $`10^4`$ lower.

Hydrocarbons The hydrocarbon photochemical scheme used here is a simplified version of the Jovian atmospheric model described in Strobel (1973) and Gladstone et al. (1996). The photodissociation of CH4 and the subsequent reactions of the species with hydrogen produce all the other hydrocarbons present in a Jovian-type atmosphere. For lower temperatures and weaker stellar irradiation, the main source of H is from H2 and CH4 photodissociation and the main sink via C2H2, which acts as a catalyst in recombining H. HD 209458b receives much greater stellar irradiation and is therefore much hotter than Jupiter. In this case, the formation of H is greatly enhanced by photolysis of H2O and reactions between O and OH radicals and H2. The sink for H is more complex (see $`\mathrm{\S }`$ 3.3).

Oxygen O is similar in abundance to C and represents a cosmic abundance of these species. The atmospheric H2O abundance will be controlled mainly by the comparative richness of these two species. The amount of H2O and CO governs the amount of atomic oxygen present and the related reactions are of great interest. The oxygen related reactions are taken from Moses et al. (2000).

Model atmosphere Three models are investigated in this paper. Our standard reference model using solar abundances is shown in Fig. 1 (Model A). We also consider 2 other cases, Model B and Model C, in which H2O and CO abundances are 10 times lower, respectively. We have taken the 1 bar level to be ”0” km, and all heights are referenced from this level. The temperature-pressure profile and chemical abundances are based on Seager et al. (2000). The temperature decreases from the bottom to the top of the atmosphere. The abundances of CO and H2O are, $`3.6\times 10^4`$ and $`4.5\times 10^4`$, respectively. These values are similar to solar abundances. The CH4 abundance is determined by the thermodynamic equilibrium chemistry in the deep atmosphere. We adopt the value $`3.9\times 10^8`$, which is the lower end of the model by Seager & Sasselov (2000). The temperature-pressure profiles are not certain, because global circulation and high temperature condensation (Seager & Sasselov, 2000; Sudarsky et al., 2000) are not included in generating the model atmosphere. Nevertheless, the present standard reference model is accurate enough for a first order understanding of the chemistry and characterizes a source of H in the upper atmosphere. The eddy diffusion is proportional to n (n is number density) where $`\alpha `$ is taken to be $`0.60.7`$.

## 3 RESULTS

The principle results for our three models will be described in three sections. Section 3.1 deals with OH and O radicals, followed by CO2 and CH4 in section 3.2. The important question of hydrogen production in relation to H2O is addressed in section 3.3.

### 3.1 OH and O radicals

Fig. 2 shows the OH and O radicals in our models. OH and O are the most important radicals as they drive most of the chemical reactions. O is the most important element in facilitating the formation of OH radicals in the water-poor atmosphere. These species are produced photochemically, for example,

$$\text{CO}+\text{h}\nu \text{C}+\text{O},$$

$$\text{H}\text{2}\text{O}+\text{h}\nu \text{H}+\text{OH}.$$

CO photolysis is an important source of O and H2O is the main source of OH. Though O(1D) is not as abundant as O ($`10^5`$ \[O\]), it can react with H2 to produce a similar amount of OH. O is increasing with altitude as a consequence of H2O and CO photolysis. OH is increasing with altitude until it starts decreasing at $`10`$ nbar. The decline of OH above 10 nbar is due to OH photodissociation. We see the mixing ratio of OH radicals is not sensitive to the abundance of CO and H2O. With an order of magnitude change in H2O (Model B) or CO (Model C), OH is changed only by a factor of $`3`$. However, O is sensitive to both CO and H2O concentrations and preferentially forms OH. From Fig. 2 we see that OH is not sensitive to H2O abundance. The abundance of H2O depends on the comparative richness of cosmic C and O. Under high stellar UV irradiation, a fraction of CO will be photodissociated. The resulting O will react with H2 to form OH, which eventually forms H2O by reacting with H2. Therefore, the abundance of OH radicals is determined by the amount of O in the system. It does not matter whether O is in the CO or H2O reservoir.

### 3.2 CO2 and CH4

CO2 is formed via the reaction of CO and OH,

$$\text{OH}+\text{CO}\text{CO}\text{2}+\text{H}$$

Fig. 3 shows the vertical profiles of CO2 for our three models. The CO2 mixing ratio is enhanced in the upper atmosphere. At pressures of $`10`$ nbar for Model A, the CO2 mixing ratio is only about 2 orders of magnitude less than its progenitor, CO. The CO2 abundance in the model is rather insensitive to the abundance of H2O. An order of magnitude decrease in H2O results in only a factor of $`3`$ decrease in CO2 abundance (cf. Models A and B in Fig. 3). However, CO2 abundance varies approximately linearly with the abundance of CO (cf. Models A and C in Fig. 3).

The formation of CH4 is initiated by the downward flux of C atoms produced in the photolysis of CO in the upper atmosphere. This obtains the following sequence of reactions$`:`$

$$\text{CO}+\text{h}\nu \text{C}+\text{O},$$

$$\text{C}+\text{H}\text{2}+\text{M}\text{3}\text{CH}\text{2}+\text{M},$$

$$2\text{3}\text{CH}\text{2}\text{C}\text{2}\text{H}\text{2}+2\text{H},$$

$$\text{C}\text{2}\text{H}\text{2}+\text{H}+\text{M}\text{C}\text{2}\text{H}\text{3}+\text{M},$$

$$\text{C}\text{2}\text{H}\text{3}+\text{H}\text{2}\text{C}\text{2}\text{H}\text{4}+\text{H},$$

$$\text{C}\text{2}\text{H}\text{4}+\text{H}+\text{M}\text{C}\text{2}\text{H}\text{5}+\text{M},$$

$$\text{C}\text{2}\text{H}\text{5}+\text{H}2\text{CH}\text{3},$$

$$\text{CH}\text{3}+\text{H}+\text{M}\text{CH}\text{4}+\text{M}.$$

Fig. 3 shows the vertical profiles of CH4 in our models. We see that the CH4 mixing ratio is increasing by a factor of 5-100 from the bottom to the 0.1 mbar level (cf. Fig. 3). Above this level, CH4 rapidly decreases due to photodissociation. Some of the C is eventually converted to CO2, whose mixing ratio increases while that of CH4 decreases.

The CH4 mixing ratio is increased by a factor $`2`$ when we lower the H2O abundance by an order of magnitude (cf. Models A and B in Fig. 3). We suggest this increase is due to less UV shielding by water above. The CH4 mixing ratio is decreased by an order of magnitude when we lower the CO abundance by an order of magnitude. The reason is that CO photolysis is the source of the C in CH4.

### 3.3 H and H2O

Fig. 3 shows the mixing ratios for H and H2O. The most striking features are the production of H (all three models) and the production of H2O (Model B). In our 1-D model, the production rate of atomic hydrogen is not sensitive to the exact abundances of CO and H2O. With an order of magnitude change in the abundance of either CO or H2O, the atomic hydrogen changes by only a small factor, $`12`$. This implies the production of H in three models is limited by the availability of UV photons. The H production is also not sensitive to the abundance of CH4. CH4 abundance has been increased to be as high as CO, and the H mixing ratio is only changed by a small factor. A more comprehensive discussion will be given in a separate paper (Liang et al. 2003, in preparation). The mixing ratio of H exceeds 10% at the top of atmosphere. At the top of the atmosphere ($`<1`$ nbar), H2 will be photolyzed and is a source of H. This atomic hydrogen will fuel the hydrodynamic loss process, as observed by Vidal-Madjar et al. (2003) and is discussed in more detailed in the companion paper by Parkinson et al. (2003, in preparation).

In a water-poor atmosphere (e.g., Model B), CO will be driving the photochemical reactions to form H2O$`:`$

$$\text{CO}+\text{h}\nu \text{C}+\text{O},$$

$$\text{O}+\text{H}\text{2}\text{OH}+\text{H},$$

$$\text{OH}+\text{H}\text{2}\text{H}\text{2}\text{O}+\text{H}.$$

H2O and H are the net products. The produced H2O will be recycled and is an important source of OH radicals and H atoms. Near the top of atmosphere, a large fraction of H2O is destroyed due to the high UV bombardment. However, CO is more stable.

Fig. 4 shows the production rate of H and the photolysis rate of H2O. The rates are not sensitive to the abundance of CO, but are sensitive to the abundance of H2O. By comparing the H2O photolysis rate with H production rate, it is evident that production of H is mainly driven by H2O photolysis and the reaction of OH with H2. Below $`1`$ mbar, HCO plays a role in the removal of H via

$$\text{CO}+\text{H}+\text{M}\text{HCO}+\text{M},$$

$$\text{HCO}+\text{H}\text{CO}+\text{H}\text{2}.$$

In the upper atmosphere, H atom recombination and reactions with CH and 3CH2 will drive the loss of H via

$$2\text{H}+\text{M}\text{H}\text{2}+\text{M},$$

and

$$\text{H}+\text{CH}\text{C}+\text{H}\text{2},$$

$$\text{H}+\text{3}\text{CH}\text{2}\text{CH}+\text{H}\text{2},$$

$$\text{C}+\text{H}\text{2}+\text{M}\text{3}\text{CH}\text{2}+\text{M}.$$

## 4 CONCLUSION

We have considered a series of possible chemical reactions using various models for a ”hot jupiter”. We have shown the mechanism for producing the atomic hydrogen. The production of H is not sensitive to the abundances of CO, H2O, and CH4. Lowering H2O or CO an order of magnitude changes the concentration of H by only a factor of $`2`$. However, the production rate of H is sensitive to the temperature profile. A 30% change in the temperature will result in $`50\%`$ change in the H concentration.

Our calculations show that the H mixing ratio at $`1`$ mbar is $`10^3`$ and exceeds 10% in the top of the atmosphere. Being less gravitationally bound, the atomic hydrogen formed at the top of atmosphere can escape hydrodynamically as putatively suggested by observations of Vidal-Madjar et al. (2003). Since these close-in gas-rich giant planets are probably tidally locked, circulation may be important in transporting heat because the temperature gradient can be as high as 1000 K over the entire planet, which in turn implies wind speeds of a few km s-1 (Showman & Guillot, 2002; Cho et al., 2003). Therefore, it is of interest to simulate the differences in chemical processes between the day and night sides and the global transport of heat and mass as well.

We thank Dr. R. L. Shia and Dr. M. Gerstell for helpful comments and discussions. The support of NASA Grant NAG5-6263 to the California Institute of Technology is gratefully acknowledged. This work represents one phase of research carried out at the Jet Propulsion Laboratory, California Institute of Technology, under contract to the National Aeronautics and Space Administrations.


format: text

text:
# Low Temperature Thermal Conductivity of High Purity YBa2Cu3O6.99 in the Vortex State

## Abstract

The low temperature thermal conductivity of a $`d`$-wave superconductor is considered for arbitrary strength of the impurity potential. A random distribution of vortices is included within the usual semiclassical approach. The vortex-quasiparticle scattering is accounted for phenomenologically, in addition to impurity scattering. Application of the theory to recent data on high purity YBa2Cu3O6.99 shows excellent agreement with experiment, both in zero field and in the presence of an applied magnetic field. We find that the strength of the impurity potential required to describe the data deviates qualitatively from either the unitary or Born limits.

Recent measurements of the in-plane thermal conductivity $`(\kappa )`$ of a high purity YBa2Cu3O6.99 single crystal hill  have revealed disagreement with the standard $`d`$-wave theory in several respects. The very low temperature $`(T)`$ variation of $`\kappa `$ fits neither the unitary nor the Born limit of impurity scattering theory. In addition, when a magnetic field $`(H)`$ is applied up to $`13`$ Tesla, the $`T=0`$ data disagree with the predictions of the usual semiclassical theory of the Doppler shift which accounts for the supercurrents circulating around the vortex cores. $`\kappa `$ initially increases very rapidly reaching almost twice its universal value $`\kappa _{00}`$ as $`H`$ increases to a fraction of a Tesla. Beyond this value, $`\kappa `$ remains almost unchanged up to $`13`$ Tesla while the semiclassical theory of Kübert and Hirschfeld kh  predicted it should continue to rise as $`\sqrt{H}`$.

Quasiparticles can be created not only thermally but also by the magnetic field. This additional creation of quasiparticles gives rise to a rapid increase of $`\kappa `$. However, as the density of vortices is increased the quasiparticles should also experience additional vortex-quasiparticle scattering due to Andreev reflection, and this should have the corresponding effect of decreasing $`\kappa `$. It is well-known that vortex scattering has important effects on the thermal conductivity of a type II superconductor cleary . However, it was not included in the work of Ref.kh , which applies only to the dirty limit, where it is assumed that impurity scattering is dominant. To treat the pure limit where the vortex scattering plays a significant role, we follow a semi-phenomenological approach. We use Matthiessen’s rule to add the vortex scattering to the usual impurity term which is left unchanged.

A random distribution of vortices is assumed and a configuration average is done as usual. After the averaging the vortex scattering is dependent only on the single vortex cross section. This quantity has been evaluated by Yu et al. yu , who gave a simple formula which is proportional to the magnetic energy $`E_H`$. Here the proportionality constant will be treated as a phenomenological parameter. Franz franz  has also given a dimensional argument for this form. However, his approach does not give the correct behavior of $`\kappa `$ at low temperatures because in his calculation $`\kappa \kappa _{00}`$ as $`T0`$ regardless of $`H`$, which is inconsistent with experiments. While Vekhter and Houghton vekhter  did include Andreev scattering in their work, they start, however, with an Abrikosov lattice which introduces extra complications associated with Bloch’s theorem. Their model does not apply to a random distribution and it also predicts a $`\sqrt{H}`$ behavior at low temperature. Moreover, their theory is not well-justified for $`HH_{c2}`$. Arguing that the Dirac nodes remain unaffected even in the vortex state (except for renormalized Fermi and gap velocities), Franz and Vafek fv  concluded that for $`H_{c1}HH_{c2}`$ the thermal conductivity is independent of $`H`$. This may happen only in a perfect vortex lattice; however, it is debatable if vortices form a lattice in a rigorous sense even in a very pure sample.

In the work of Hill et al. hill  for $`H=0`$, $`\kappa /T`$ rises out of its universal limit very fast as a function of temperature, and reaches 4 times this universal value within $`400mK`$. This behavior cannot be understood within either of the two limits extensively studied in the past graf : Born and unitary limits, corresponding to the impurity potential $`U_i0`$ and $`U_i\mathrm{}`$, respectively note1 . In this paper we consider explicitly the intermediate case. A dimensionless parameter $`c`$ is introduced, as usual, with $`c^2=1/(U_iN(0))`$, where $`N(0)`$ is the normal state density of states at the Fermi level. Within a $`t`$-matrix approximation, the effective scattering in the superconducting state can be characterized by $`c`$ and $`\mathrm{\Gamma }=n_i/(\pi N(0))`$, where $`n_i`$ is the impurity density. We find that the data of Hill et al. hill  can be fit naturally with a value of $`c0.014`$, which does not belong to either of the two limiting cases. The rapid rise in $`\kappa /T`$ with $`T`$ also reflects the high purity of the sample (i.e small value of $`\mathrm{\Gamma }`$).

An important feature of the impurity scattering is that the scattering rate $`\gamma (\omega )`$ in the superconducting state is strongly dependent on frequency and that our fit to data for YBa2Cu3O6.99 with $`T_c90K`$ gives $`\gamma (0)=0.29K`$, which is very small in comparison to, for example, the value of the magnetic energy at $`H_{c1}`$. This value is of order $`10^2\mathrm{\Delta }_0`$ while $`\gamma (0)0.0015\mathrm{\Delta }_0`$, where $`\mathrm{\Delta }_0`$ is the amplitude of a $`d`$-wave gap and is $`2.14T_c`$ for a $`d`$-wave superconductor. Thus when we include a magnetic field $`H`$ into the theory we will need to work in the regime $`E_H>\gamma (0)`$. The limit $`E_H<\gamma (0)`$ is unrealistic for such pure samples. Nevertheless, provided $`E_H/\mathrm{\Delta }_0`$ and $`T/\mathrm{\Delta }_0`$ remain small, we can apply the nodal approximation durst  to reduce complications associated with the calculations.

The superconducting $`(2\times 2)`$ matrix Green’s function in the Nambu notation takes the form $`\widehat{G}(𝐤,i\omega _n)=(\stackrel{~}{\omega }\widehat{\tau }_0+\mathrm{\Delta }_𝐤\widehat{\tau }_1+\xi _𝐤\widehat{\tau }_3)/(\stackrel{~}{\omega }^2\xi _𝐤^2\mathrm{\Delta }_𝐤^2),`$ where $`\widehat{\tau }`$’s are the Pauli matrices in spin space, $`\mathrm{\Delta }_𝐤`$ a $`d`$-wave order parameter, and $`\xi _𝐤`$ is the electronic energy dispersion in the normal state. The renormalized Matsubara frequencies $`\stackrel{~}{\omega }(i\omega _n)=i\omega _n\mathrm{\Sigma }(\stackrel{~}{\omega })`$ include the self-energy $`\mathrm{\Sigma }(\stackrel{~}{\omega })`$, which depends on what interactions are taken into account. Here we consider only the scalar component of the self-energy assuming that other components can be neglected or absorbed into $`\xi _𝐤`$ or $`\mathrm{\Delta }_𝐤`$durst . In our consideration for a vortex state, it has two parts: the impurity part $`\mathrm{\Sigma }_i`$ and the vortex part $`\mathrm{\Sigma }_v`$. For $`\mathrm{\Sigma }_i`$ we use the self-consistent $`t`$-matrix approximation while for $`\mathrm{\Sigma }_v`$ we take a more phenomenological approach due to the lack of a detailed microscopic footing.

In the absence of a magnetic field, $`\stackrel{~}{\omega }(\omega +i\delta )=\omega \mathrm{\Sigma }_{i,ret}(\stackrel{~}{\omega })`$. Introducing the spectral functions for the diagonal $`\left[G(𝐤,i\omega )\right]`$ and off-diagonal $`\left[F(𝐤,i\omega )\right]`$ components of $`\widehat{G}(𝐤,i\omega )`$, $`A(𝐤,\omega )=2\text{Im}\left[G(𝐤,\omega +i\delta )\right]`$ and $`B(𝐤,\omega )=2\text{Im}\left[F(𝐤,\omega +i\delta )\right]`$, respectively, we obtain the thermal conductivity at $`T`$

$`{\displaystyle \frac{\kappa (T)}{T}}`$ $`=`$ $`{\displaystyle \frac{1}{T^2}}{\displaystyle \underset{𝐤}{}}(v_f^2+v_g^2)`$ (1)
$`\times `$ $`{\displaystyle \frac{d\omega }{2\pi }\left[A^2(𝐤,\omega )B^2(𝐤,\omega )\right]\omega ^2\left(\frac{f}{\omega }\right)},`$
where $`v_{f(g)}`$ is the Fermi (gap) velocity and $`f(\omega )`$ is the Fermi function. As we mentioned earlier, application of the nodal approximation allows us to carry out the summation over $`𝐤`$ and yields

$$\frac{\kappa (T)}{T}=\frac{1}{\pi ^2}\frac{1}{T^2}\left(\frac{v_f}{v_g}+\frac{v_g}{v_f}\right)𝑑\omega 𝒜(\omega )\omega ^2\left(\frac{f}{\omega }\right),$$
(2)
where $`𝒜(\omega )=1+\left[h(\omega )+\frac{1}{h(\omega )}\right]\mathrm{arctan}[h(\omega )].`$ Here $`h(\omega )`$ is defined as $`h(\omega )=(\omega \text{Re}\left[\mathrm{\Sigma }_{i,ret}\right])/\gamma (\omega )`$ with $`\gamma (\omega )=\text{Im}\left[\mathrm{\Sigma }_{i,ret}\right]`$. In order to calculate the thermal conductivity in Eq. (2), we need to determine the self-energy $`\mathrm{\Sigma }_{i,ret}`$ self-consistently. In the $`t`$-matrix approximation $`\mathrm{\Sigma }_{i,ret}`$ can be written as $`\mathrm{\Sigma }_{i,ret}(\stackrel{~}{\omega })=\mathrm{\Gamma }G_0(\stackrel{~}{\omega })/\left[c^2G_0^2(\stackrel{~}{\omega })\right]`$, where $`G_0=[2\pi N(0)]^1_𝐤\text{Tr}[\widehat{\tau }_0\widehat{G}_{ret}]`$. As we described, $`c`$ here controls the strength of the impurity scattering. Applying the nodal approximation, we obtainewald

$$G_0(\stackrel{~}{\omega })\frac{2}{\pi }\left[\frac{\stackrel{~}{\omega }}{\mathrm{\Delta }_0}\mathrm{ln}\frac{\stackrel{~}{\omega }}{4\mathrm{\Delta }_0}i\frac{\pi }{2}\frac{\stackrel{~}{\omega }}{\mathrm{\Delta }_0}\right].$$
(3)
Substituting Eq. (3) into $`\mathrm{\Sigma }_{i,ret}`$, we compute $`\text{Re}[\mathrm{\Sigma }_{i,ret}]`$ and $`\text{Im}[\mathrm{\Sigma }_{i,ret}]`$ self-consistently. In the presence of a magnetic field, the results for $`\mathrm{\Sigma }_{i,ret}`$ remain unchanged.

Since $`𝒜(\omega )2`$ as $`T0`$, it is easy to show that $`\kappa (T)/T`$ approaches its universal value, $`(2/3)\left[v_f/v_g+v_g/v_f\right]\kappa _{00}/T`$. In Fig. 1 we show results for $`\kappa (T)/T`$ normalized by the universal value for different values of $`c`$. The figure is motivated by the data of Hill et al. on a very pure sample of YBa2Cu3O6.99. Their experimental fit appears as solid circles corresponding to $`(0.16+3.0T^2)/0.16`$ hill . We should note that this data requires subtraction of a phonon contribution; this can be done in a variety of sensible ways, which leads to an uncertainty in the electronic part of several percent taillefer . The solid curve shown, a theoretical calculation with $`c=0.014`$, is thus in excellent agreement with the data. While there are two parameters in the impurity theory, $`\mathrm{\Gamma }`$ and $`c`$, a constraint was applied to all fits: namely, $`\gamma (0)`$ is set to be $`0.29K`$. This constraint comes from consideration of finite as well as zero magnetic field data as will be explained later. As seen in Fig. 1 the unitary limit does not fit the data at finite $`T`$. On the other hand we can achieve reasonable agreement with the data for $`c=0.014`$ and $`\mathrm{\Gamma }/\mathrm{\Delta }_05.05\times 10^5`$ over the entire temperature range shown in Fig. 1. This corresponds to a very clean sample. If $`\mathrm{\Gamma }`$ were to be interpreted as the normal state impurity scattering rate, the corresponding mean free path would be anomalously large. We take this as an indication that for this sample the underlying normal state, on which the superconductivity is built, may not be an ordinary Fermi liquid.

We note that the theory does not give an exact straight line while the experimental data has been parameterized with a linear fit hill . However, it should be remembered that to obtain the electronic thermal conductivity from experiment it was necessary to subtract a considerable phonon contribution. Uncertainty remains about this subtraction as we will elaborate on later and so the small disagreement that remains between theory and experiment is not a serious limitation on the fit obtained. In the inset of Fig. 1 we show $`\gamma (\omega )`$ obtained for $`c=0.014`$ (solid curve) compared with the unitary limit (dotted curve). The curves are completely different. They both start with $`\gamma (0)/\mathrm{\Delta }_0=0.0015`$ by design but for finite $`\omega `$ the solid curve rises while the dotted curve drops monotonically. These differences lead to the large difference in the $`T`$ dependence of $`\kappa (T)/T`$ shown in Fig. 1.

In the presence of a magnetic field. Next we wish to include the effect on the thermal conductivity of an external magnetic field $`H`$ applied along the $`c`$-axis. To this end we need to know the self-energy of a quasiparticle due to vortices $`\mathrm{\Sigma }_v`$ as well as for the impurities. Formally speaking, the semiclassical approximation takes into account an approximate real part of the self-energy: $`\text{Re}\left[\mathrm{\Sigma }_v\right]𝐯_s(𝐫)𝐤`$ where $`𝐯_s(𝐫)`$ is the superfluid velocity at a position $`𝐫`$ measured from the center of a vortex core. Here $`𝐫`$ is restricted to the vortex unit cell and an average over $`𝐫`$ is to be carried out. For the imaginary part of $`\mathrm{\Sigma }_v`$, we consider the Andreev reflection as a possible mechanism. When a quasiparticle moves from $`𝐫`$ to $`𝐫^{}`$, the Andreev reflection takes place effectively if $`\sqrt{\xi _𝐤^2+\mathrm{\Delta }_𝐤^2}+𝐯_s(𝐫)𝐤=|\mathrm{\Delta }_𝐤|+𝐯_s(𝐫^{})𝐤`$. It has been shown yu  that for nodal quasiparticles the vortex scattering rate $`\gamma _v`$ is determined mainly by the intervortex distance $`R`$; namely, $`\gamma _vv_f/R`$. This implies that within the nodal approximation a vortex plays the role of a defect or a scattering center. Since the intervortex distance $`R\sqrt{\mathrm{\Phi }_0/H}`$, where $`\mathrm{\Phi }_0`$ is a flux quantum, and the magnetic energy $`E_H`$ is given by $`E_H=a\sqrt{\pi /2}v_f\sqrt{H/\mathrm{\Phi }_0}`$, where $`a`$ is a parameter of order unitykh , we can write $`\gamma _v=bE_H`$, which defines a vortex scattering parameter $`b`$. Instead of evaluating the microscopic quantities which determine the size of $`b`$ from first principles, we propose to treat it as a phenomenological parameter which we can use to fit the zero temperature data of the thermal conductivity as a function of $`H`$.

An averaged value of the thermal conductivity over the vortex distribution can be evaluated as usual: $`\kappa (H,T)=𝑑ϵ𝒫(ϵ)\kappa (ϵ,T)`$, where $`𝒫(ϵ)=A^1𝑑𝐫\delta \left[ϵ𝐯_s(𝐫)𝐤\right]`$ is the vortex distribution within the area of a vortex unit cell $`A`$. The appropriate average (series or parallel) depends on the direction of the field as discussed by Kübert and Hirschfeldkh . We will not deal directly with these issues here. We made numerical calculations assuming series and parallel arrangements for the thermal conductivity but have found no qualitative differences. Now the thermal conductivity in a vortex state becomes

$`{\displaystyle \frac{\kappa (H,T)}{T}}`$ $`=`$ $`{\displaystyle \frac{1}{\pi ^2}}{\displaystyle \frac{1}{T^2}}\left({\displaystyle \frac{v_f}{v_g}}+{\displaystyle \frac{v_g}{v_f}}\right){\displaystyle 𝑑ϵ𝒫(ϵ)}`$ (4)
$`\times `$ $`{\displaystyle 𝑑\omega 𝒜(\omega ,ϵ)\omega ^2\left(\frac{f}{\omega }\right)},`$
where $`𝒜(\omega ,ϵ)=1+\left[h(\omega ,ϵ)+\frac{1}{h(\omega ,ϵ)}\right]\mathrm{arctan}[h(\omega ,ϵ)].`$ Note the similarity between Eqs. (2) and (4). The effects of the Doppler shift and the Andreev scattering appears in $`h(\omega ,ϵ)`$; $`h(\omega ,ϵ)=\left(\omega \text{Re}\left[\mathrm{\Sigma }_{i,ret}\right]ϵ\right)/\gamma _t(\omega )`$, where $`\gamma _t(\omega )=\gamma (\omega )+bE_H`$.

Taking the limit of $`T=0`$ in Eq. (4) leads to a simplified expression as follows: $`\kappa (H,0)=(\kappa _{00}/2)𝑑ϵ𝒫(ϵ)𝒜(0,ϵ)`$ with $`𝒜(0,ϵ)=1+\left[ϵ/\gamma _t(0)+\gamma _t(0)/ϵ\right]\mathrm{arctan}\left[ϵ/\gamma _t(0)\right]`$. To proceed further some vortex distribution function $`𝒫(ϵ)`$ needs to be specified. We consider the Gaussian distribution; $`𝒫(ϵ)=\mathrm{exp}\left[ϵ^2/E_H^2\right]/(\sqrt{\pi }E_H)`$, which is believed to be favorable for Andreev scattering distribution . Now we obtain the thermal conductivity as $`T0`$:

$`\kappa (0,H)`$ $`=`$ $`{\displaystyle \frac{\kappa _{00}}{2}}{\displaystyle 𝑑x\frac{e^{x^2}}{\sqrt{\pi }}}`$ (5)
$`\times `$ $`\left\{1+\left[h(0,x)+{\displaystyle \frac{1}{h(0,x)}}\right]\mathrm{arctan}\left[h(0,x)\right]\right\}`$
where $`h(0,x)=x/\left[\gamma (0)/E_H+b\right]`$. For $`E_H\gamma (0)`$, $`h(0,x)x/b`$ and, thus, $`\kappa (0,H)/\kappa _{00}g(b)`$, where $`g(b)=(1/2\sqrt{\pi })𝑑xe^{x^2}\left[1+(x/b+b/x)\mathrm{arctan}(x/b)\right]`$. Consequently, the magnetic field saturated value of $`\kappa (0,H)`$ is $`g(b)\kappa _{00}`$ and not $`\kappa _{00}`$ as found in Ref.franz .

The function $`g(b)`$ is shown in the inset of Fig. 2. Note that as $`b0`$, $`g(b)`$ goes to infinity as $`\sqrt{\pi }/(4b)`$ and no saturation occurs in this case. This also corresponds to the dirty limit case considered in Ref.kh  where it is assumed that $`\gamma (0)`$ is much bigger than $`bE_H`$. In the main frame of Fig. 2 we show our results for the normalized $`\kappa (H,T)/T`$ vs $`H`$ at $`T=0`$ for the two samples considered by Hill et al. , a very pure sample denoted by BZO with $`T_c90K`$ (triangle) and an impure sample denoted by YSZ (circle). In both cases the agreement with the data is excellent (compare with Fig. 3 of Ref. hill ).

For the BZO sample with $`\gamma (0)0.29K`$ the value of $`b`$ obtained in the fit is $`b=0.3`$ assuming $`v_f10^7cm/s`$ and $`a=1/2`$ so that $`E_H/\mathrm{\Delta }_00.036\sqrt{H}`$T-1/2chiao . Using $`H_{c1}30`$mT of YBa2Cu3Ox with $`T_c91K`$ handbook , we have $`E_H/\mathrm{\Delta }_00.006\mathrm{\Delta }_0`$, which is greater than $`\gamma (0)`$. For the dirty sample (YSZ) with $`T_c93K`$, the estimated value of $`\gamma (0)`$ is about $`13K`$ hill . Using these parameters with $`v_f2\times 10^7cm/s`$, we found that no vortex scattering parameter needs to be introduced to fit the data; namely, $`b=0`$. In this case the magnetic energy at a few Tesla is still of the order of $`\gamma (0)`$ and so the data is not sensitive to a small value of $`b`$.

In Fig. 3 we show results for the normalized $`\kappa (H,T)/T`$ in units of the universal value as a function of $`(T/T_c)^2`$. The solid curve reproduces our $`H=0`$ result for $`c=0.014`$. The different symbols represent different values of $`H`$; namely, $`H=1`$T (open circle), $`4`$T (open square), and $`9`$T (open diamond). Note that within $`10\%`$ or less the saturated value of $`\kappa (H,T)/T`$ due to $`H`$ is temperature independent. This justifies the process used in Hill et al. to subtract out the phonon contribution from their raw data which includes both an electronic and a phonon part. We also emphasize that the crossing with the $`H=0`$ curve occurs around $`(T/T_c)^2=5\times 10^5`$ or $`T0.2K`$, which agrees well with experiment. We point out that to get this agreement we insisted on $`\gamma (0)0.29K`$ and the crossing temperature $`T^{}0.2K`$ is close to $`\gamma (0)`$ so that it is the crossing that determines this parameter in our approach. For $`b=0`$ the crossing cannot occur. This implies that both the crossing itself and the saturated value of $`\kappa (0,T)`$ show the significance of the vortex scattering. In the inset we show a few more results for small values of $`H(H_{c1})`$. This inset shows that for a weak field the thermal conductivity retains significant temperature dependence although considerably reduced from the $`H=0`$ case. As one can see in this inset or the main frame of Fig. 2, at $`T=0`$ the thermal conductivity depends strongly on $`H`$ for $`H_{c1}H<1`$T. This can be seen in the experiments.

In summary we have given a theory of the low temperature thermal conductivity of a $`d`$-wave superconductor using a self-consistent treatment of impurity scattering intermediate between the unitary and Born limits. We have also included the effects of a magnetic field in the range $`H_{c1}<HH_{c2}`$ within the semiclassical approximation which takes into account the Doppler shift due to the circulating supercurrent around the vortex cores. Assuming a random distribution of vortices, the Andreev scattering from the vortices is considered and this leads to saturation in the thermal conductivity at $`T=0`$ for a sufficiently large value of $`H`$. We find good qualitative and quantitative agreement in all respects with recent data on a very pure sample of YBa2Cu3O6.99. We also verify that the phonon contribution can be subtracted out using the thermal conductivity for a sufficiently high magnetic field.

###### Acknowledgements.

We wish to thank R. Hill and L. Taillefer for helpful correspondence and for sending us a preliminary draft of Ref.hill  before publication. This work was supported in part by the Natural Sciences and Engineering Research Council of Canada (NSERC), by ICORE (Alberta), and by the Canadian Institute for Advanced Research (CIAR).


text:
# 1 Introduction

## 1 Introduction

Data mining is the exploration and analysis, by automatic or semiautomatic means, of large quantities of data in order to discover meaningful patterns and rules . Its goal is to find *patterns* or *relationships* in data using techniques to synthesize models which are abstract representations of reality. There are two main types of data mining models : *descriptive* and *predictive*. *Descriptive* models represent patterns in data and are generally used to build meaningful groups or clusters. *Predictive* models are used to forecast explicit values, on the basis of samples with known results.

At present, astrophysics is a discipline in which the exponential growth and heterogeneity of data require the use of data mining techniques. Due to the advances in telescopes, detectors and computer technologies, astrophysics has become a domain extremely rich of scientific data. The primary source of astronomical data are the systematic sky surveys over a wide photon energy range (from $`10^7`$ eV to $`10^{13}`$ eV) . Large archives and digital sky surveys with dimensions of $`10^{12}`$ bytes currently exist, while in the near future they will reach sizes of the order of $`10^{15}`$ bytes. Numerical simulations are also producing comparable volumes of information.

Therefore, the use of data mining techniques is necessary to maximize the information extraction from such a growing quantity of data. Data mining has reached a certain degree of maturity in data exploration for decision support mostly in domains like marketing, sales and customer care. In the astrophysics domain the role of data mining is to help researchers building or verifying new physical models based on observational data.

A first issue in applying these techniques is the heterogeneity of astronomical data, due in part to their high dimensionality including both spatial and temporal components, due in part to the multiplicity of instruments and projects. Another issue is that currently astronomical data are organized in a traditional operational system, in which the emphasis is on data normalization; data mining techniques, instead, require putting more emphasis on efficiency and simplicity of data access, even if this entails some redundancy. Gathering data from multiple astronomical datasets to perform multi-wavelength analysis necessitates both using an informational system, or data warehouse, as a model for data management, and the definition of a common set of metadata to guarantee the interoperability between different archives and a simpler data exploration.

Data mining techniques are rather general and can be employed in different application domains in which there is an intensive use of data. Such techniques include :

1. Clustering techniques, such as Expectation Maximization (EM) with mixture models or Self-Organizing Maps (SOM), to find regions of interest, produce descriptive summaries and build density estimates of large datasets, or methods like Support Vector Clustering (SVC) to single-out rare objects or anomalous events.
2. Classification techniques, such as decision trees, nearest-neighbor classifiers, neural networks and statistical learning methods like Support Vector Machines (SVM), to categorize objects or clusters of objects of interest. The classification result is further analyzed to verify whether physically meaningful objects or groups of objects have been identified, and if these objects are present in some catalog or they are new.
3. Techniques to improve clustering and classification algorithms, such as genetic algorithms and Principal Component Analysis (PCA). These methods allow to find the best features and reduce the dimensionality of the domain space.
4. Software agents for automatic or semi-automatic information search and analysis.
5. Data visualization and presentation techniques (exploratory data mining). These techniques allow to present multidimensional information in a way easily understandable by a human user.

## 2 Data management in astrophysical databases

At present, several multi-wavelength projects are underway, for example SDSS, GALEX, GSC-2, POSS2, ROSAT, FIRST and DENIS. In the next years new spatial missions will be launched; two of them, AGILE and GLAST, will observe gamma-rays on a wide energy range. Besides, ground based telescopes sensitive to gamma-rays are being tested or starting taking data (MAGIC, HESS, VERITAS, CANGAROO III).

Therefore, astrophysicists will need a uniform interface to access all these data . Data gathered by all missions are heterogeneous as they are mission oriented and dependent on the particular platform or instrument (including hardware components information, quality flags decided inside the mission, derived measures for particular analysis). Several scientific research fields require to perform the analysis on multiple energy spectra and consequently to get the data from different missions. Typically, astrophysicists want to retrieve multi-spectral data for specific objects, classes of objects (i.e. AGN, HII region) or selected regions of the sky. However, metadata in mission archives are not designed to answer these queries .

### 2.1 The data warehousing approach

Most of the online resources available to the astrophysicists community are simple data archives. Typically, users can perform queries based on observational parameters (detector, type of the observation, coordinates, astronomical object, exposure time, etc.) to obtain images which are then processed by standard analysis tools. Many astronomical catalogs can be accessed online, even if it is still difficult to correlate objects in different archives or access multiple catalogs simultaneously. Some advances, in this direction, have been accomplished by projects like Vizier, Aladin and SkyView .

To identify objects and parameters which allow to answer directly to particular scientific issues it is necessary to build a *scientific archive* containing the results of data analysis - scientific measurements - rather then the data itself . The users of this archive should be able to perform queries based on scientific parameters (magnitude, redshift, spectral indexes, morphological type of galaxies, etc.), easily discover the object types contained into the archive and the available properties for each type, and define the set of objects which they are interested in by constraining the values of their scientific properties along with the desired level of detail.

Data mining applied to large astrophysical databases can involve the execution of complex queries and multiple scans of large quantity of data. Therefore, it is opportune to put more emphasis on data access efficiency rather than on data normalization.

All the aforesaid requirements can be satisfied organizing data in a data warehouse. A data warehouse can be defined as a *subject-oriented*, *integrated*, *time varying* and *non-volatile* data collection . Subject-oriented means that in a data warehouse data are collected and organized with the aim of a particular analysis. The second property is surely the most important one; in fact, a data warehouse has to integrate with the multiplicity of standards used by the sources it gathers data from (i.e. multiple astronomical catalogs). This data integration process can involve conversion of types, formats or units and the addition of derived types (i.e. several statistical measures). A data warehouse is time varying because its time horizon usually oscillates between 5 and 10 years and along this period of time data collected are a series of snapshots taken at fixed times . It is non-volatile because data updating, and the resulting loss of information, doesn’t take place within it.

In a data warehouse, data are arranged in a structure that can be easily explored and queried, with fewer tables and keys than the equivalent relational model. You start from a relational model, but some restrictions are introduced by using *facts*, *dimensions*, *hierarchies* and *measures* in a characteristic star structure called *star schema* . The central table is called “fact” table and it is the highest dimensional table of the scheme. It can represent a particular phenomenon that we want to study. This table is surrounded by a number of tables, called “dimensions”, which represent entities related to the phenomenon to be studied and connected to the central table, forming the ends of the star. Within the dimensions, attributes are arranged in hierarchies, determining the “drill-down” and “roll-up” operations available on each dimension: the result is a tree that the user can visit from the root to the leaves, refining his query (drill-down) or generalizing it (roll-up).

### 2.2 The metadata role

Accessing data into a set of continuously evolving catalogs states the problem of accessing and understanding parameters available in each catalog. A typical problem is understanding whether a catalog contains some specific data type, what is the reliability of these data, if they are written in a standard format, if they are taken from other publications or catalogs, how the associated data file can be processed. All these details describe data - they are metadata - and traditionally are presented in the introduction of printed catalogs or detailed in several publications analyzing the catalog data.

Metadata play an important role: a researcher has to obtain information about the environment in which data have been gathered in order to understand the respondence to the project requirements; for instance: date and/or data acquisition method, internal or external error estimates, aim of data. Computing systems have to access metadata to merge or compare data from different sources. For instance, it is necessary that units are expressed unambiguously to allow comparisons between data with different units.

The astrophysicists community, in addition to using the FITS (Flexible Image Transport System) exchange format, is currently considering alternatives like XML. Some attempts to define a common standard are XSIL (eXtensible Scientific Interchange Language), XDF (eXtensible Data Format) and VOTable .

## 3 Spatial and multidimensional data structures

Spatial DataBase Systems (SDBS) are designed to handle spatial data and the associated non-spatial information. Spatial data are characterized by a complex structure (a spatial object can be a single point or a set of polygons arbitrarily distributed). They are usually dynamic (requiring robust data structures for frequent insertions, deletions and updates), tend to be large (sky maps can reach sizes of Terabytes) requiring the integration of the secondary storage. There is no standard spatial algebra, that is the set of spatial operators depends on the specific application. Another important property is that there is no total ordering on spatial objects preserving spatial proximity . This characteristic makes difficult to use traditional indexing methods, like B-trees or linear hashing.

*Spatial data mining* analyzes the relationships between the attributes of a spatial object stored into the database and the attributes of the neighboring ones. Typical queries required by this kind of analysis are: *point queries*, to find all objects overlapping the query point; *range queries*, to find all objects having at least one common point with a query window; *nearest neighbor queries*, to find all objects that have a minimum distance from the query object.

These principles can be generalized to multidimensional data. Access methods to multidimensional databases can be classified in *Point Access Methods* (PAM) and *Spatial Access Methods* (SAM). PAM are designed to perform searches on point databases. They usually arrange data in buckets, each one corresponding to a disk page. The buckets are indexed by either flat or hierarchical data structures: flat structures are used in multidimensional hashing methods like the grid file and EXCELL; hierarchical structures are used in hierarchical access methods like quadtree, KD-tree and KD-B-tree. SAM manage objects with spatial properties like area and shape. Access methods in SAM are often extensions of PAM ones to handle objects with a spatial extent. Such methods include R-tree, R-tree and Multi-layer grid file.

### 3.1 Quadtree

The term *quadtree* is used to describe a class of hierarchical data structures based on the principle of recursive decomposition of the space. They can be distinguished by the following elements :

* the type of data they are used to represent;
* the principle guiding the decomposition process;
* the resolution (variable or not).

Until recent experiments, astronomical observations were restricted to a selected region of the sky, making a planar projection of the observed region adequate. However, the next generation experiments, like GLAST and AGILE, will provide a detailed observation of the whole sky and thus they will require the handling of data with a spherical distribution. In the astrophysical field, two methods for indexing the sky based on quadtrees have been designed: the Hierarchical Triangular Mesh (HTM) and the Hierarchical Equal Area isoLatitude Pixelization (HEALPix).

HTM  maps triangular regions of the sphere to unique identifiers. The technique for subdividing the sphere in spherical triangles is a recursive process. At each level of the recursion, the area of the resulting triangles is roughly the same (see figure 1). In areas with a larger data density, the recursion process can be applied with a greater level of detail than in areas with lower density. The starting point is a spherical octahedron which identifies 8 spherical triangles of equal size.

HEALPix  is a curvilinear subdivision of the sphere in quadrilaterals (pixel) of equal area (but variable shape). The contour of a pixel is defined by the equation $`cos\theta =a+b\times \varphi `$ on the equator and by $`cos\theta =a+b/\varphi ^2`$ on the polar regions. This structure makes more efficient the execution of operations typically performed on the sky maps including: convolution with local and global kernels, Fourier analysis with spherical harmonics, nearest-neighbor searches.

### 3.2 KD-tree and its variants

The KD-tree  is a binary tree that stores points of a $`k`$-dimensional space. In each internal node, the KD-tree divides the $`k`$-dimensional space into two parts with a $`(k1)`$-dimensional hyperplane. The direction of the hyperplane, that is the dimension on which the division in performed, alternates between the $`k`$ possibilities from one tree level to the following one. The subdivision process is recursive and terminates when the size of a node (its longer side) or the number of points contained into it is below a certain threshold. Given $`N`$ data points, the average cost of an insertion operation is $`O(\mathrm{log}_2N)`$. The tree structure and the resulting hierarchical division of the space depends on the *splitting rule*.

A drawback of KD-trees is that they have to be completely contained into the main memory. With large datasets this is not feasible. KD-B-trees  and hB-trees  combine properties of KD-trees and B-trees to overcome this problem.

### 3.3 R-tree and its variants

The R-trees  are hierarchical data structures meant to efficiently index multidimensional objects with a spatial extent. They are used to store not the real objects but their minimum bounding box (MBB). Each node of the R-tree corresponds to a disk page. Similar to B-trees, the R-trees are balanced and they guarantee an efficient memory usage. Due to the overlapping between the MBBs of sibling nodes, in an R-tree a range query can require more than one search path to be traversed.

To solve the overlapping problem, the R<sup>+</sup>-tree access method introduced in  uses a clipping operation to avoid the intersection between intervals at the same tree level. Objects intersecting more than one MBB at a specific level are clipped and copied in several pages. This way, a single search path is traversed for an exact match query. However, insertion operations are more complex.

In the R-trees, search performances depend on the insertion algorithms. In  an improved version of the R-tree, called R-tree, has been proposed. This version uses a new insertion policy which significantly improves search performances. The main target of this policy is to minimize the overlapping between MBBs of sibling nodes to reduce the number of search path to be traversed during a query operation.

## 4 Clustering algorithms on large datasets

Clustering algorithms have to locate regions of interest in which to perform more detailed analysis and point out correlations between objects. An important issue, in large datasets, is the efficiency and scalability of the clustering algorithms with respect to the dataset size.

Many scalable algorithms have been proposed in the last ten years, including: BIRCH , CURE , CLIQUE .

In particular, BIRCH is a hierarchical clustering algorithm. The main idea behind the algorithm is to compress data into small subclusters and then to perform a standard partitional clustering on the subclusters. Each subcluster is represented by a *clustering feature* which is a triplet summarizing information about the group of data objects, that is the number of points contained into the cluster and the linear sum and the square sum of the data points. This algorithm has a linear cost with respect to the number of data points.

CURE is an hierarchical agglomerative algorithm. Instead of using a single centroid or object, it selects a fixed number of well-scattered objects to represent each cluster. The distance between two clusters is defined as the distance between the closest pair of representatives points and at each step of the algorithm, the two closest clusters are merged. The algorithm terminates when the desired number of clusters is obtained. To reduce the computational cost of the algorithm, these steps are performed on a data sample (using suitable sampling techniques). Its computational cost is not worse than the BIRCH one.

CLIQUE has been designed to locate clusters in subspaces of high dimensional data. This is useful because generally, in high dimensional spaces, data are scattered. CLIQUE partitions the space into a grid of disjoint rectangular units of equal size. The algorithm is made up of three phases: first, it finds subspaces containing clusters of dense units, than identifies the clusters, and finally generates a minimum description for each cluster. Also this algorithm scales linearly with the database size.

## 5 Supervised learning and classification

Classification algorithms are required in order to identify objects belonging to known classes. In case of scientific (and in particular astrophysics) data, care has to be taken on the interpretability of the classification results. For this reason, one of the most popular methods to classify scientific data (in addition to neural networks) is the algorithm based on decision trees . In fact, with this method, the learning algorithm produces a binary tree which performs the classification by means of value ranges on the data attributes.

Recently, Support Vector Machines (SVMs) are an active research domain within the field of machine learning.

### 5.1 SVM for classification and novelty detection

SVM and the related kernel methods are becoming popular for data mining tasks like classification, regression and novelty detection. This approach is systematic, reproducible, and properly grounded by statistical learning theory .

In its simplest form, an SVM is able to perform a binary classification finding the “best” separating hyperplane between two linearly separable classes. There are infinite hyperplanes properly separating the data. So, the SVM finds this plane maximizing the distance, or *margin*, between the *support* planes for each class (see  for theoretical foundations). A plane supports a class if all points in that class are on one side of that plane. This problem is formulated as a quadratic programming problem (QP) and can be solved by effective robust algorithms. If the data is not linearly separable, *slack* variables are introduced into the QP problem to accept outliers. Finally, a further non-linearity is introduced using kernel functions (satisfying the Mercer’s condition ) to map data to a higher dimensional space.

In many real problems, the task is not classifying but novelties or anomalies detecting. In astrophysics, possible applications are the research of anomalous events or new astronomical sources. An approach is modeling the *support* of a distribution (rather than estimating the density function of the data). A method to solve this problem is represented by the Support Vector Clustering (SVC) algorithm , in which data are mapped to a higher dimensional space by means of a Gaussian kernel function. In the new space, the algorithm finds the minimum sphere enclosing the data. The Mapping of the sphere to the original input space generates a set of contours enclosing the data and corresponding to the support of the distribution.

## 6 Conclusions

In this work we have studied some data management and mining issues related to astrophysical data, aiming at a complete data mining framework. In particular, we have justified the need for a data warehousing approach to handle astrophysical data and we have focused on multidimensional access methods to efficiently index spatial and multidimensional data. A second issue concerns clustering techniques on large datasets, and we have discussed about some scalable algorithms with linear computational complexity. Finally, we have focused on classification algorithms introducing an increasingly popular method named Support Vector Machine, whose applications include the tasks of classification, regression and novelty detection.


text:
# Neuromorphometric characterization with shape functionals

## I introduction

Many natural phenomena are defined or influenced by the geometrical properties of the involved elements, and vice-versa. Examples of such a shape-function relationship include the chemical properties of proteins, the aerodynamic efficiency of wings, and the oxygen exchanges through elaborated bronchic structures. The close relationship between geometry and function provides strong motivation for the geometrical analysis of natural objects. Yet, the state-of-the-art of geometrical characterization, an area sometimes called morphometry or morphology, remains in a relatively incipient stage where several competing, and often divergent, approaches coexist. While powerful methods have been used in physics to express relevant geometrical properties, often including differential measurements such as curvature, the situation remains particularly challenging in biology.

The relationship between neuronal shape and function has attracted increasing attention due to its far reaching implications for basic neuroscience and for medical applications. Neuron morphology has the special characteristic, that it evolves during the developmental stage of the cell, being influenced by its molecular environment and the history of synaptic activity Whitford et al. (2002). The mature neuronal shape, together with its membrane electrical properties, determine the electric conductance of the cell Koch and Segev (2000) and account for part of its electrophysiological characteristics such as firing patterns and computational abilities van Ooyen et al. (2002); Fukuda et al. (1984); Wassle et al. (1981); Krapp et al. (1998). Software packages, such as NEURON, are available for modelling neuronal activity with basis on cable theory Rall (1977) which can be useful for analysis of real or virtual neurons Coelho and da F. Costa (1997); Ascoli (2002); G.Ascoli (2002). At the same time, neuronal shape can vary for different tissues, depending on a number of extracellular factors Kandel and Schwartz (1985); Whitford et al. (2002). Ultimately, it determines the patterns of connectivity and, consequently, the overall network computational abilities. On the application side, tasks such as automated morphological characterization of neuronal shape and the diagnosis of abnormalities, deserve further investigation.

While the use of morphological tools was severely constrained until recently by the cost of the relatively sophisticated systems needed to process images and geometry, the continuing advances in computer software and hardware have paved the way for an ever widening range of possible applications. Consequently, more effective morphological concepts and methods have been developed and reported in the literature, including the use of differential geometry concepts such as multi-scale curvature and bending energy Cesar and da F. Costa (1997); da F. Costa and Velte (1999); Cesar and da F. Costa (1998), methods from mathematical morphology such as skeletonization da F. Costa (2000); Falcao et al. (2002), as well as the recently reported framework known as Integral-Geometry Morphological Image Analysis (MIA) Michelsen and de Raedt (2001, 2000); Mecke (1998). The latter approach involves the use of additive shape functionals, i.e. mappings that take shapes to single scalar values, in terms of a parameter usually related to the spatial scale or time. As far as neuroscience is concerned, the contour of a neuronal cell has been shown to possess a fractal structure Morigiwa et al. (1989) and its multiscale fractal dimension been used to characterize different morphological classes of neuronal cells da F Costa et al. (2002).

Primarily motivated by the possibility of applying the methodology proposed in Michelsen and de Raedt (2001, 2000) as a novel and potentially useful tool for addressing the problem of neuronal shape characterization and classification, the present work also provides an assessment of those measures considering a database of real biological data, namely camera lucida images of cat ganglion neuronal cells. In order to provide a comparative reference, the multi-scale fractal dimension da F Costa et al. (2002), itself a shape functional, is also considered as a measure for shape characterization.

Integral geometry provides an adequate mathematical framework for morphological image analysis, having a core of useful theorems and formulae that in some cases leads to analytical results for averages of image functionals Mecke (1998) while also being quite efficient to implement computationally Michelsen and de Raedt (2000). Here, the class of functionals involved are restricted to additive, motion invariant and continuous, called Minkowski functionals. These functionals are related to usual geometric quantities, for instance, in the Euclidean plane, to area, perimeter and connectivity or the Euler number, which expresses the number of holes in a connected pattern such as the image of a neuronal cell.

In order to describe geometrically one object, a set of measures (functionals) is taken and the behaviour of these measures is monitored, as some control parameter is varied. In this work we compute additive functionals in the plane as the contour of a neuronal cell image is inflated by a parallel set dilation of radius $`r`$, the control parameter. The non-additive multi scale fractal dimension is derived from one of those computed additive functionals, at each radius of dilation, giving important da F Costa et al. (2001) complementary information.

This article starts by presenting the adopted methodology, the considered shape functionals, and the statistical procedure for cell identification. The results, which are presented subsequently, clearly indicate that the new measures are efficient for distinguishing morphologically the two functional classes alpha and beta as well as revealing a strong morphological coherence in one of the classes.

## II Methodology

### II.1 Additive shape functionals

The morphological characterization in Euclidean plane by means of shape functionals explores simple properties of convex sets. For these basic geometric objects, such as triangle or an ellipse, we may evaluate a change in area while the object undergoes a morphological dilation with the knowledge of its initial geometry. For example the change in area of a convex body $`K`$, after a parallel set dilation using a 2D ball of radius $`r`$, can be expressed as

$$A(K_r)=A(K)+U(K)r+\pi r^2,$$
(1)
where $`A(K)`$ and $`U(K)`$ stand for the initial area and perimeter of the object $`K`$ and $`r`$ is the dilation parameter. The process of taking parallel sets generalises naturally to higher dimensions, while the change in hyper volume preserves the general form (1) and is given by the Steiner formula

$$v^d(K_r)=\underset{\nu =0}{\overset{d}{}}\left(\genfrac{}{}{0pt}{}{d}{\nu }\right)W_\nu ^{(d)}(K)r^\nu ,$$
(2)
where the coefficients $`W_\nu ^{(d)}`$ are referred to as quermassintegrals or Minkowski functionals, Michelsen and de Raedt (2001). These functionals, as a generalisation of known geometric quantities, are additive, motion invariant and continuous. Moreover, a Theorem by Hadwiger Michelsen and de Raedt (2001); Mecke (1998) states that these functionals form a complete set of measures, with the above properties, on the set of convex bodies

$$\varphi (K)=\underset{j=0}{\overset{d}{}}c_jW_j^{(d)}(K).$$
(3)
Notwithstanding, the change of any, additive, motion invariant and continuous functional can be expressed, using a generalised Steiner formula Mecke (1998); Santaló (1976), in terms of the initial geometric information

$$\varphi (K_r)=\underset{j=0}{\overset{d}{}}\underset{k=0}{\overset{dj}{}}c_j\left(\genfrac{}{}{0pt}{}{dj}{k}\right)W_{k+j}^{(d)}(K)r^k.$$
(4)
The notion of connectivity number or Euler characteristic $`\chi `$ is central in establishing the aforementioned properties of Minkowski functionals. The usual definition of the connectivity from algebraic topology in two dimensions is the difference between the number of connected $`n_c`$ components and the number of holes $`n_h`$,

$$\chi (K)=n_cn_h,$$
(5)
while in three dimensions distinction should be made between two kind of holes namely cavities $`n_{hc}`$ and handles(tunnels) $`n_{hh}`$

$$\chi (K)=n_cn_{hh}+n_{hc}.$$
(6)
Integral Geometry provides an equivalent definition for connectivity number of a convex set $`K`$ which is given by

$$\chi (K)=\{\begin{array}{cc}1\hfill & K\mathrm{}\hfill \\ 0\hfill & K=\mathrm{}.\hfill \end{array}$$
(7)
Of great importance is its property of additivity

$`\chi (A)=`$ $`\chi (_{i=1}^lK_i)={\displaystyle \underset{i}{}}\chi (K_i){\displaystyle \underset{i<j}{}}\chi (K_iK_j)+\mathrm{}`$
$`\mathrm{}+(1)^{l+1}\chi (K_1\mathrm{}K_l).`$ (8)
Additivity and motion invariance is inherited by the Minkowski functionals as they are related to the connectivity number by the formulae

$`W_\nu ^{(d)}(A)={\displaystyle _𝒢}\chi (AE_\nu )𝑑\mu E_\nu \nu =0,\mathrm{},d1`$
$`W_d^{(d)}(A)=\omega _d\chi (A),\omega _d=\pi ^{d/2}/\mathrm{\Gamma }(1+d/2).`$ (9)
In the above expression $`E_\nu `$ stands for an $`\nu `$-dimensional plane in $`\mathrm{R}^d`$. The integral is to be taken for all positions, induced by isometries $`𝒢`$, of $`E_\nu `$ weighted by $`d\mu (E_\nu )`$, the kinematical density which is in turn related to the Haar measure on the group of motions $`𝒢`$, see Michelsen and de Raedt (2001); Mecke (1998); Santaló (1976).

To sum up, the Minkowski functionals $`W_\nu ^{(d)}(A)`$, as a generalisation of the usual procedure for volume determination, counts the number of possible intersection of a $`\nu `$-dimensional plane with the domain $`A`$.

If one is to take advantage of the above additivity property all intersection in (II.1) must be taken into account. When working on a lattice, there is a more expedient route, exploring the discrete nature of the images and the additivity of the Minkowski functionals, which consists of a decomposition of the 2D body $`A`$ into a disjoint collection of interior bodies, open edges and vertices. Following the usual nomenclature we denote the interior of a set $`A`$ by $`\stackrel{˘}{A}=A/A`$. For an open interior of a $`n`$-dimensional body embedded in a $`d`$-dimensional Euclidean space there is the following expression for the Minkowski functionals, Michelsen and de Raedt (2001),

$$W_\nu ^{(d)}(\stackrel{˘}{A})=(1)^{d+n+\nu }W_\nu ^{(d)}(A),\nu =0,\mathrm{},d.$$
(10)
We may then apply additivity and the lack of connectivity of open sets on the lattice to determine the functionals for the body as a whole

$$W_\nu ^{(d)}(𝒫)=\underset{m}{}W_\nu ^{(d)}(\stackrel{˘}{N}_m)n_m(𝒫),\nu =0,\mathrm{},d.$$
(11)
Where $`n_m(𝒫)`$ stands for the number of building elements of each type $`m`$ occurring in the pattern $`𝒫`$. For a two dimensional space which is our interest for the present neuron images, we display in Table 1 the value of Minkowski functionals for the building elements on a square lattice of pixels and their direct relation to familiar geometric quantities on the plane. Using the information presented in Table 1 and equation (11) we have

$$A(𝒫)=n_2,U(𝒫)=4n_2+2n_1,\chi (𝒫)=n_2n_1+n_0.$$
(12)
So the procedure of calculating Minkowski functionals of a pattern $`𝒫`$ has been reduced to the proper counting of the number of elementary bodies of each type that compose a pixel (squares, edges and vertices) involved in the make up of $`𝒫`$.

In section III we describe typical results for the evaluation of the above presented additive functionals using an actual neuron image. The procedure involves the implementation of, first an algorithm for the proper parallel set dilation throughout all permitted radius on the square lattice and, second, of an algorithm for the calculation of Minkowski functionals by counting disjoint building elements based on the formulae (12). An efficient routine for undertaking the latter is described in detail in Michelsen and de Raedt (2001, 2000).

### II.2 Multi scale fractal dimension

As an example of a related non-additive functional, we add to the previous measures the multi scale fractal, an approach which has been applied successfully to neuromorphometry da F Costa et al. (2001). The notion of multi-scale fractal dimension refers to a quantitative characterization of complexity and the degree of self-similarity at distinct spatial scales, see da F Costa et al. (2001, 2002). Intuitively, the fractal dimension indicates how much the curve extends itself throughout the space. As a consequence, more intricate curves will cover the surround space more effectively and will display a higher fractal dimension. This quantity is calculated via the derivative of the logarithm of the changing interior area as the neuron cell image undergoes a dilation. As such it is immediately derivable from the first additive functional (Area) in a pixel based approach as opposed to a curvature approach, see da F Costa et al. (2002).

### II.3 Implementation

We have conducted the evaluation of the functionals described above on a 800MHz ordinary personal computer running Linux. Both the algorithm for exact dilations on the square lattice, according to da F. Costa and Manoel (2001) and the pixel based algorithm for the estimation of the Minkowski functionals Michelsen and de Raedt (2001) was implemented in SCILAB-2.6. It took approximately 40s to calculate the functionals for each dilation radius.

## III Results

We start by describing a measure that is perhaps the simplest one but which holds important information not only on its own, but through relation with the multi scalar fractal dimension described in section II.2. In Figure 3 we show the typical monotonically increasing curve of the interior area of a neuron cell as its contour is inflated by a parallel set procedure. To capture the gross structure of this measure we calculate the area under the curve, $`\mathrm{𝗌𝗎𝗆}^a`$, and in a size independent manner, the radius $`𝖱_{\frac{1}{2}}^a`$ at which the area below this curve reaches half of its value at the end of dilation. The fine structure is given by its standard deviation, $`\mathrm{𝗌𝗍𝖽}^a`$.

In Figure 4 we show a typical curve of the perimeter of the evolving frontier of the neuron cell as the contour is inflated by a parallel set procedure. It is important to observe that this measure is affected by errors introduced by the discrete nature of the lattice of pixels and the low resolution of the original image, which becomes more acute in beta cells because of their reduced size. Nonetheless, this effect tends to be less important after the tenth radii or so. To capture the gross structure of this measure, we calculate the area under the curve $`\mathrm{𝗌𝗎𝗆}^p`$ and the size independent radius $`𝖱_{\frac{1}{2}}^p`$ at which the area below this curve reaches half of its overall value. The fine structure is given simply by the standard deviation of the data, $`\mathrm{𝗌𝗍𝖽}^p`$.

In contrast to the preceding measure the connectivity or the Euler number of the neuron shape as it undergoes the dilation processes is independent of the resolution of the image. It is a measure restricted to the topology of the shape counting essentially the number of holes at each radius of dilation no matter the holes are perfectly round or not. The alterations of connectivity are subtle from step to step, which is reflected by the complex distribution of cusps in Figure 5. The vast amplitude of scale for which there is an abrupt change of connectivity is a measure of complexity of this type of cell. Figure 6 shows the particular behaviour of the connectivity, for samples from the two classes, as the cell shape undergoes a dilation. We take the gross information of this measure by extracting the area under the interpolated curve, $`\mathrm{𝗌𝗎𝗆}_{\mathrm{𝗂𝗇𝗍𝖾𝗋𝗉}}`$, and a monotonicity index given by

$$i_s=\frac{s}{s+d+p},$$
(13)
where $`s,d`$ and $`p`$ counts respectively the number of times that the curve increase, decrease or reach a plateau. This index characterizes a perfect monotonically increasing curve when its value is 1 and reaching its minimum value for a curve of high variability. This measure is designed specially to explore the multi-scale nature of the neuronal complexity. The finer structure is captured solely by the standard deviation $`\mathrm{𝗌𝗍𝖽}_{\mathrm{𝖽𝗂𝖿𝖿}}`$ of the difference between the interpolated curve and the original data.

As a measure of complexity, the multi-scale fractal dimension has been experimentally found to be related to the connectivity of a shape. Although this relationship is not straightforward, there might be a correlation between these measures as commented bellow. For this image functional, we evaluate the maximum fractal dimension, the mean fractality and the standard deviation, respectively max, mean and $`\mathrm{𝗌𝗍𝖽}^f`$. A typical curve for the multi scale dimension is shown in Figure 7 for the same neuron appearing in Figure 1.

Among all the considered measurements, we found a good separation of alpha and beta type cells for the feature space defined by the perimeter half integral radius, $`𝖱_{\frac{1}{2}}^p`$, and the standard deviation for the fine structure of connectivity, $`\mathrm{𝗌𝗍𝖽}_{\mathrm{𝖽𝗂𝖿𝖿}}`$. Figure 8 shows the obtained alpha and beta clustering with both class exhibiting similar dispersion. Another good result for morphological characterization was obtained for the connectivity or the Euler number of the cell shape. A feature space involving the index of monotonicity and the integral of the interpolated connectivity curve is presented in Figure 9. More efficient than the former in separating the two classes, these measures produce a well localised clustering of beta cells characterising its geometrical intricacies. This result suggests that alpha class is indeed a more homogeneous category while the beta class may have a morphological subclass structure.

Table 2 shows the correlation coefficients for the twelve measurements considered in this work. Of special interest is the correlation between fractality and connectivity, measures that have experimentally been found to represent complementary but not redundant measures of complexity. Figure 10 shows a combination of two measures to form a feature space which in this case shows a poor separation of classes which is in accordance with the high anti-correlation of the involved measures presented in Table 2. Unusually high correlation appears between some measures notably as occurring between area and perimeter, suggesting a specific tendency that seems to be particular to the type of data (neurons) and not a general rule.

## IV Conclusions

The use of additive shape functionals has been recently considered for the characterization of the geometrical properties of several physical objects Michelsen and de Raedt (2001); Mecke (1998). The current article explored the use of a representative set of such functionals, namely the area, perimeter and connectivity, for the characterization of neural shapes represented in terms of a whole set of parallel expansions. The multi scale fractal dimension, a non-additive shape functional, was also considered as a standard for comparison.

All the adopted shape functionals consist of functions of the dilating radius of each parallel body. For the sake of efficience, the following compact subset of global features was selected: the area under the functionals, the value of radius where the area reaches its half-value, the standard deviation of the functionals, as well as a novel measurements expressing the monotonicity of the functionals and the decomposition of the connectivity functional in terms of a low and high variation signals.

All two-by-two combinations of these measures were investigated visually in order to identify the combinations of features leading to more pronounced separations between the two classes of considered neural cells, namely cat retinal ganglion cells of type alpha and beta. The coefficients of correlation of each pair of measures were also estimated and analysed, indicating decorrelation between several of the considered features. The obtained results confirmed a differentiated potential of each measurement for neural cell clustering, with the features derived from the connectivity functional accounting for the best separation between classes. However, further investigations and comparison to more established and recent methods (such as in Mizrahi et al. (2000)) based explicitly on dendritic morphology will be necessary to reveal the best realization of the novel methodology proposed in this article. The biological implications of such results are that the two type of cells differ in the distribution of holes (defined by the respective dendritic arborizations) for different spatial scales. The obtained clusters indicated that the alpha cells exhibit less uniform geometrical properties than the beta, suggesting the existence of morphological subclasses.

###### Acknowledgements.

The authors are grateful to FAPESP (processes 02/02504-01, 99/12765-2 and 96/05497-3) and to CNPQ (process 301422/92-3) for financial support.


text:
# Charmonium Production in 𝛾-A, p-A and A-A

## 1 Introduction

The NA50-collaboration  observed quite different $`J/\psi `$ and $`\psi ^{}`$-nucleon cross sections. Analogous results have been found at Fermilab : Both experiments found that the charmonium-nucleon cross section is larger in the target fragmention region than at midrapidity. This is in good agreement with models that assume that the $`\psi ^{}`$ is produced as colorless, point-like $`c\overline{c}`$ pair, which expands with time to its full size and that there exists a relationship between the spatial distribution of color in a hadron and the cross section of its interaction with a nucleon. Such a relation is proved in pQCD (perturbative QCD) . For the nonperturbative regime it is a well known experimental fact that spatially larger hadrons have larger interaction cross sections. For example, 3.5 mb for the inelastic $`J/\psi `$-nucleon cross section was found at SLAC  while the inelastic $`\pi `$-nucleon cross section at this energies is 20 mb. From charmonium models, e.g. in ref.  it is known that the different charmonium states ($`J/\psi `$, $`\chi `$ and $`\psi ^{}`$) have different spatial sizes. We concluded in ref.  that these cross sections are dominated by non-perturbative contributions.

The agreement of this scenario with the negative $`x_F`$ Fermilab data  with was first demonstrated in ref. . Another attempt to describe these data with such an expansion scenario was done in ref. . However, in ref.  20%-50% of the suppression of charmonium states is due to a color octet state added to the model to explain also the large $`x_F`$ regime. This state as defined in ref.  is in variance with QCD, because of the following reason. The eigen life time of this state adjusted to the Fermilab data is only 0.06 fm. A gluon emitted in such a short time has to have a momentum of $`1/(0.06\mathrm{fm})3.3`$ GeV relative to the $`c\overline{c}`$ pair . In ref.  this gluon was assumed to be massless.

In ref.  it was shown that the production of $`J/\psi `$’s in pA collisions can be described, if one takes into account the production and the subsequent decays of higher resonances ($`\chi `$,$`\psi ^{}`$) into $`J/\psi `$’s. This leads to a significant increase of the absorption of $`J/\psi `$’s as compared to the propagation of genuine $`J/\psi `$-states. In ref.  it was shown also that the production of $`J/\psi `$’s in AA collisions is additionally suppressed by the final state interaction of charmonium states with newly produced particles like $`\pi `$’s, $`\rho `$’s and so on.

In sect. 2 the used models are introduced and in sect. 3 they are compared to the data. A more detailed description can be found in ref. . We summarize in sect. 4.

## 2 Model Description

### 2.1 Semiclassical Glauber Approximation

The suppression factor $`S`$ for minimum bias $`pA`$ collisions can be evaluated within the semiclassical approximation (cf. ) as

$`S_A`$ $`=`$ $`{\displaystyle \frac{\sigma (pAX)}{A\sigma (pNX)}}`$ (1)
$`=`$ $`{\displaystyle \frac{1}{A}}{\displaystyle \mathrm{d}^2Bdz\rho (B,z)\mathrm{exp}\left(_z^{\mathrm{}}\sigma (XN)\rho (B,z^{})dz^{}\right)}.`$ (2)
Here $`\rho (B,z)`$ is the local nuclear ground state density (we used the standard parametrisation from ). $`\sigma (XN)`$ is the interaction cross section of the charmonium state $`X`$ with a nucleon and $`\sigma (pAX)`$ ($`\sigma (ppX)`$) is the production cross section of the state $`X`$ in a pA (pp) collision. We want to draw attention to the fact that this cross section changes with time due to the space-time evolution of color fluctuations. Therefore, it is necessary to keep $`\sigma `$ under the integral. The suppression factor $`S`$ of $`J/\psi `$’s produced in the nuclear medium is calculated as:

$$S=0.6(0.92S^{J/\psi }+0.08S^\psi ^{})+0.4S^\chi .$$
(3)
Here $`S^X`$ are the respective suppression factors of the different pure charmonium states $`X`$ in nuclear matter. Eq. (3) accounts for the decay of higher resonances after they left the target nucleus into $`J/\psi `$’s. The fractions of $`J/\psi ^{}s`$ that are produced in the decays of higher resonances in eq. (3) are taken from ref. . However, in ref.  it is assumed that the different charmonium states interact with nucleons with the same cross section, which is in disagreement with the data from the refs. .

In line with the above discussion we want to stress here that Eqs. (2) and (3) are applicable at CERN energies for central and negative rapidities, but have to be modified, if applied already at $`y_{c.m.}0`$ at RHIC or higher energies, because at higher energies charmonium states can be produced outside of the nucleus and the $`c\overline{c}`$ pairs propagate through the whole nucleus without forming a hadron. Data are often presented in the form $`\alpha `$ vs. A, where $`\alpha `$ is given by

$$\sigma _{pA}=\sigma _{pp}A^\alpha \text{ so that }S=A^{\alpha 1}.$$
(4)
In the semiclassical Glauber approximation, we take into account these color fluctuations in an effective way as described in ref. . We assume that charmonium states are produced at $`z`$ as small $`c\overline{c}`$ configurations predominantly through gluon-gluon-fusion, then they evolve – during the formation time $`t_f`$ – to their full size. Please note that there is up to now no theoretical or experimental proof for the assumption that charmonium states are produced in point-like configurations as predicted in pQCD. A way to test this experimentally was suggested recently in ref. . It’s based on the idea that the $`\psi ^{}`$ and the $`\psi ^{\prime \prime }`$ are mixed states of the $`2S`$ and the $`1D`$ charmonium states. The names of these states comes from a comparison of nonrelativistic charmonium models with the nonrelativistic wave functions of the positronium (actually even the name charmonium was given to pronounce the similarity with positronium). Namely

$`|\psi ^{}`$ $`=`$ $`\mathrm{cos}\theta |2S+\mathrm{sin}\theta |1D,`$ (5)
$`|\psi ^{\prime \prime }`$ $`=`$ $`\mathrm{cos}\theta |1D\mathrm{sin}\theta |2S.`$ (6)
Since only the $`S`$-wave contributes to the decay of $`\psi `$ states into $`e^+e^{}`$-pairs (at least in nonrelativistic charmonium models) the value of $`\theta =19\pm 2^o`$ can be determined from the data on the $`e^+e^{}`$ decay widths of $`\psi ^{}`$ and $`\psi ^{\prime \prime }`$. If the production of these charmonium states is pointlike, then only the $`S`$-wave is produced. These leads to an universal ratio of

$$\frac{\sigma (\psi ^{\prime \prime })}{\sigma (\psi ^{})}\frac{\mathrm{\Gamma }(\psi ^{\prime \prime }l^+l^{})}{\mathrm{\Gamma }(\psi ^{}l^+l^{})}0.1.$$
(7)
$`\sigma (\psi ^{\prime \prime })`$ resp. $`\sigma (\psi ^{})`$ are here the production cross sections of the $`\psi ^{\prime \prime }`$ and resp. $`\psi ^{}`$ in various processes. Predictions for different processes can be found in ref. .

If the formation length of the charmonium states, $`l_f`$, becomes larger than the average internucleon distance ($`l_f>r_{NN}1.8`$ fm), one has to take into account the evolution of the cross sections with the distance from the production point . Here we assume motivated by quantum diffusion that the cross section increases linearly with time. The formation length of the $`J/\psi `$ is given by the energy denominator $`l_f\frac{2p}{m_\psi ^{}^2m_{J/\psi }^2}`$, where $`p`$ is the momentum of the $`J/\psi `$ in the rest frame of the target. With $`p=30`$ GeV, the momentum of a $`J/\psi `$ produced at midrapidity at SPS energies ($`E_{lab}=200`$ AGeV), this yields $`l_f3`$ fm, i.e. a proper formation time of $`\tau _f=0.3`$ fm. As formation time of the $`\psi ^{}`$ in its rest system we use here the radius given by nonrelativistic charmonium models, e.g. see the refs. . This radius is $`r=0.45`$ fm for the $`\psi ^{}`$. A larger value of $`\tau _f`$ for the $`\psi ^{}`$ is supported also by the extraction of the formation time of the $`J/\psi `$ . Finally the formation time is $`t_f=\gamma \tau _f`$, where $`\gamma `$ is the the Lorentz-factor of the charmonium state relativ to the nuclear target. For higher gamma factors, i.e. at higher energies, the formation time becomes larger than the nuclear targets. In this regime a hadronic description during the formation time is questionable. A partonic model for this energy range was proposed in ref. .

### 2.2 Vector Dominance Model for $`\gamma `$A

The VDM (Vector Dominance Model) takes into account only the direct diffractive production of the $`J/\psi `$ and the $`\psi ^{}`$, while the GVDM (Generalized Vector Dominance Model) accounts also for the non-diagonal transitions ($`\psi ^{}+NJ/\psi +N`$ and $`J/\psi +N\psi ^{}+N`$). The later are needed, because in photoproduction the particles are produced as point like configurations and develop then to their average size. In a hadronic model like the GVDM this is taken into account in form of the interference due to the nondiagonal matrix elements. In the GVDM the photoproduction amplitudes $`f_{\gamma \psi }`$ and $`f_{\gamma \psi ^{}}`$ for the $`J/\psi `$ and the $`\psi ^{}`$ are given by

$`f_{\gamma \psi }`$ $`=`$ $`{\displaystyle \frac{e}{f_\psi }}f_{\psi \psi }+{\displaystyle \frac{e}{f_\psi ^{}}}f_{\psi ^{}\psi }`$ (8)
$`f_{\gamma \psi ^{}}`$ $`=`$ $`{\displaystyle \frac{e}{f_\psi }}f_{\psi \psi ^{}}+{\displaystyle \frac{e}{f_\psi ^{}}}f_{\psi ^{}\psi ^{}}`$ (9)
Here $`f_\psi `$ and $`f_\psi ^{}`$ are the $`J/\psi \gamma `$ and the $`\psi ^{}\gamma `$ coupling and $`f_{VV^{}}`$ are the amplitudes for the processes $`V+NV^{}+N`$, where $`V`$ and $`V^{}`$ are the $`J/\psi `$ and the $`\psi ^{}`$ respectively. In the VDM the non-diagonal amplitudes with $`VV^{}`$ are neglected. The importance of the nondiagonal transitions is evident,because the left hand side of eq. (9) is small. If it is neglected as a first approximation , then $`f_{\psi ^{}\psi }=\frac{f_\psi }{f_\psi ^{}}f_{\psi \psi }1.7f_{\psi \psi }`$. And due to the CPT-theorem $`f_{\psi ^{}\psi }=f_{\psi \psi ^{}}`$.

In ref.  the GVDM yields approximately $`8\pm 2`$ mb for the $`\psi ^{}`$-nucleon interaction cross section at SPS-energies as can be seen in Fig. 1. $`\omega `$ is the laboratory energy of the photon. The $`J/\psi `$-nucleon interaction cross section at SPS-energies approximately 3.5 - 4 mb is used as input into the analysis of ref. . The accuracy of such GVDM in predicting the $`\psi ^{}`$N cross sections is not clear. The above calculation demonstrates that implementing color transparency leads to significantly larger cross sections of the $`\psi ^{}`$ N interaction.

Figure 1: The energy dependence of the elementary charmonium-nucleon cross sections found in the GVDM. The filled areas show the variation of the cross sections due to the uncertainty of the experimental $`J/\psi `$N cross section.

## 3 Comparison with NA50 and NA51 Data

In Fig. 2 we show a comparison between calculations with different cross sections and different expansion times and the NA50 data  for pA collisions and the NA51 data  for pp and pD collisions for the cross section of $`\psi ^{}`$ N interaction vs. the mass of the target. The y-axis shows $`B_{\mu \mu }\sigma _\psi ^{}/A`$ where $`B_{\mu \mu }`$ is the branching ratio for the decay of the $`\psi ^{}`$ into dimuon pairs. and $`\sigma _\psi ^{}`$ is the production cross section. The ”5.1 mb, instant formation” curve in Fig. 2 is the fit of the NA50 collaboration to their data. Instant formation means that they assumed that the $`\psi ^{}`$ is produced with the full cross section and not as a point like particle as in the description of this paper. (Note the NA50 collaboration fitted $`B_{\mu \mu }\sigma _\psi ^{}/\sigma _{DY}`$, where $`DY`$ means Drell-Yan, we multiplied this fit with the $`DY`$ cross section in pp collisions measured by NA51).

The ”8 mb, $`t_f=0.45`$ fm” curve is the eye-ball fit of the model described in this paper. For the comparison with the data we need the production cross section of the $`\psi ^{}`$ in pp collisions as input. We used here the average of the pp and pD data of the NA51 collaboration. The value of 8 mb agrees well with the model parametrisations discussed in the sections 2. However, we compare also with the calculation with the parameters of ref. , i.e. $`\sigma (\psi ^{}N)=20`$ mb and $`t_f=0.6`$ fm. For this comparison we used the production cross section of the $`\psi ^{}`$ in pD collisions divided by two as input. This is also close to the value of the NA50 fit. One can see in Fig. 2 that the calculation with these parameters is also in good agreement with the data. A value for $`\sigma (\psi ^{}N)`$ of the size of 20 mb is favored by the nucleus-nucleus data as shown in ref.  and in Fig. 3. Plotted is $`B_{\mu \mu }\sigma (J/\psi )/\sigma (DY)`$ with the absorption cross sections from ref.  and the NA50 data  for PbPb collisions vs. the transverse energy $`E_t`$, a measure for the centrality of the collision. The calculation agrees well with the data. The calculation for the $`\psi ^{}`$ underestimates the data. However it is not understood, if this is due to the high value of $`\sigma (\psi ^{}N)=20`$ mb, or if nondiagonal transitions like in sect. 2.2 should be taken into

Figure 2: Plotted are the $`B_{\mu \mu }\sigma _\psi ^{}/A`$ values extracted from calculations with different absorption cross sections and different formation times, the NA51 data for pp and pD, and the NA50 data for pBe, pCu, pAg, and pW vs. the mass $`A`$ of the target.

account in AA collisions, too.

The value of 8 mb is smaller than the theoretical estimate 20 mb of ref. . This is because in ref.  a formation time of 0.6 fm was chosen for the $`\psi ^{}`$, while we used her 0.45 fm, the radius of the $`\psi ^{}`$ given by the charmonium models. The fact that the formation time is not known very well is another theoretical uncertainty. Further uncertainty comes from using diffusion model of expansion at the distances comparable to the scale of the soft interaction. Within the error bars the $`\psi ^{}`$-nucleon cross section extracted from these pA data and the prediction of the GVDM, discussed in sect. 2.2, are qualitatively similar. However, further data are needed to learn more about this cross section.

Figure 3: Plotted is $`B_{\mu \mu }\sigma (J/\psi )/\sigma (DY)`$ with the absorption cross sections from ref.  and the NA50 data for PbPb collisions vs. the transverse energy $`E_t`$

## 4 Conclusions

The new data of the NA50-collaboration  and the data of the E866-collaboration  prove that the $`\psi ^{}`$-nucleon cross section is much larger than the $`J/\psi `$-nucleon cross section. This is in agreement with the photoproduction data for these charmonium states as discussed in the framework of the GVDM in section 2.2. This confirms the QCD prediction that the strength of hadron-hadron interactions depends on the volume occupied by color.

Within the assumption that charmonium states are produced as point like white states, we demonstrated that the data  can be fitted with a $`\psi ^{}`$-nucleon cross section of $`\sigma (\psi ^{}N)8`$ mb. However, a much larger cross section of e.g. $`\sigma (\psi ^{}N)20`$ mb is not ruled out by the data. Due to the large experimental errors we conclude that the data and the QCD-motivated models agree, but further data with higher accuracy and covering larger rapidity range are needed.

This cross section will be measured soon in proton-nucleus collisions at HERA B at an energy of $`E_{lab}=920`$ GeV. The advantage of this experiment is that it covers a larger range of Feynman-$`x_F`$, especially in the negative $`x_F`$ region. In this region effects due to the formation time of the hadron will be less important and the genuine cross sections will be measured.

I thank Leonid Frankfurt, Mark Strikman, Horst Stöcker, Ramona Vogt, and Mihail Zhalov for discussions, Ramona Vogt for Fig. 3, and acknowledge support by the Minerva Foundation and the organizers of the SQM2003.

## References