text:
# Warm Molecular Layers in Protoplanetary Disks

## 1 Introduction

It is well established from millimeter and infrared observations that the birth of solar-mass stars is accompanied by the formation of a circumstellar disk (Beckwith & Sargent 1996, Natta et al. 2000). These disks are important both as reservoirs of material to be accreted onto growing stars and as sites of planetary formation. Because the gas and dust in the disk are the basic components from which future solar systems are built, studies of their chemistry are essential to investigate the link between interstellar and planetary matter. Moreover, the chemical abundances and molecular excitation depend on physical parameters in the disks such as temperature and density, and on processes such as radial and vertical mixing. Thus, studies of the chemistry in protoplanetary disks can help to constrain their physical structure.

Although CO millimeter lines are routinely used to trace the gas and Keplerian velocity field in disks around classical T Tauri stars (e.g., Kawabe et al. 1993, Koerner et al. 1993, Dutrey et al. 1994, 1996, Koerner & Sargent 1995, Saito et al. 1995, Guilloteau & Dutrey 1998, Thi et al. 2001), detections of other molecules are still rare. Dutrey et al. (1997) and Kastner et al. (1997) were the first to report observations of molecules such as HCO+, HCN, CN, HNC, H2CO, and C2H in the disks around GG Tau, DM Tau and TW Hya. Dutrey et al. found that the abundances of these species relative to hydrogen in the DM Tau and GG Tau disks are lower than those in molecular clouds by factors of 5-100 (Dutrey et al. 1994, Dutrey et al. 1997, Guilloteau et al. 1999). The abundance ratio of CN/HCN, on the other hand, is significantly higher than in molecular clouds in all three disks. More recently, Qi (2000), van Zadelhoff et al. (2001) and Thi et al. (in preparation) have reported observations of molecules other than CO in the disks around the T Tauri stars LkCa15 and TW Hya, and in those around the Herbig Ae stars MWC 480 and HD 163296, confirming the trends of high CN/HCN ratio and molecular depletion; the disk masses derived from CO (and its isotopes) are significantly smaller than those estimated from the dust continuum assuming a gas to dust ratio of 100.

Early models of the chemistry in disks considered mostly the one-dimensional radial structure in the cold midplane (e.g. Aikawa et al. 1997, Willacy et al. 1998). Subsequently, it has been recognized that the vertical stratification of the molecules is equally relevant. Aikawa & Herbst (1999a) investigated the two-dimensional chemical structure within the so-called Kyoto disk model, representative of the minimum mass solar nebula. This model has low midplane temperatures and is isothermal in the vertical direction. The molecular abundances were found to vary significantly with height $`Z`$ from the midplane. At large radii ($`R>`$100 AU), the temperature is so low that most species, except H2 and He, are frozen out onto the grains. This depletion is most effective in the midplane region ($`Z0`$) because of the higher density and hence shorter timescales for molecules to collide with and stick to the grains. In regions above and below the midplane, significant amounts of molecules can remain in the gas phase for longer periods because of the lower densities, and because of non-thermal desorption by cosmic rays and/or radiation (e.g., X-rays) from the interstellar radiation field and the central star. Aikawa & Herbst (1999a) suggested that the observed molecular line emission comes mostly from this region. There is a height distinction between stable molecules and radicals, however. In the surface region of a disk, radicals such as CN are very abundant because of photodissociation via ultraviolet radiation, whereas the abundances of the stable molecules such as HCN peak closer to the midplane. The molecular column densities obtained by integrating over height at each radius compare well with those derived from observations of DM Tau (Dutrey et al. 1997), although detailed comparison through radiative transfer (i.e. calculation of molecular line intensities from model disks) is still lacking. The freeze-out of molecules in the midplane explains the low average abundances of heavy-element-containing species relative to hydrogen, whereas the high abundance ratio of CN to HCN is caused by photodissociation in the surface layers.

To obtain this good agreement with observations, Aikawa & Herbst (1999a) were forced to use the simplifying assumption that the probability $`S`$ for sticking upon collision of a molecule with a grain is significantly smaller than unity. If the sticking probability is unity and if only thermal desorption is considered, the molecular column densities obtained in the Kyoto model are much smaller than observed, which suggests that either there is some efficient non-thermal desorption mechanism, or that the disk temperature is higher than assumed in the Kyoto model. Aikawa & Herbst (1999a) adopted an artificially low sticking probability $`S=0.03`$ in order to reproduce the observed CO spectra, without specifying the non-thermal desorption process or modifying the temperature distribution in the Kyoto model. Since adsorption is such a dominant process in the disk, the cause of the apparently low sticking probability should be considered more seriously.

In order to explain the strong mid-infrared emission from disks, several groups suggested the possibility of higher dust temperatures than assumed in the Kyoto model due to efficient reception of stellar radiation by flared disks (e.g., Kenyon & Hartmann 1987, Chiang & Goldreich 1997, D’Alessio et al. 1998, 1999). In the two-layer Chiang & Goldreich (1997) model (C-G model hereafter), the upper layer (the so called “super-heated” layer) is directly heated by stellar radiation from the central star to temperatures $`T>50`$ K at radii of $`100`$ AU. Also, recent observations of high-frequency lines (e.g. CO $`J=65`$) support the possibility that the disk temperature is higher than assumed in the Kyoto model (Thi et al. 2001, van Zadelhoff et al. 2001). Willacy & Langer (2000) investigated the molecular distributions in the C-G model to see if the super-heated layer can maintain enough gaseous organic molecules to account for the observations. It was found, however, that molecules in this layer are destroyed by the harsh ultraviolet radiation from the star, whereas they are frozen out onto the grains in the cold lower layer. These authors therefore had to adopt a very high photodesorption rate in the lower layer to keep the molecules off the grains.

In this paper we report an investigation of molecular distributions in another disk model with a vertical temperature gradient: the model of D’Alessio et al. (1998, 1999). These scientists obtained the temperature and density distribution in steady accretion disks around T Tauri stars by solving the equations for local 1-D energy transfer (including radiation, convection and turbulence) and hydrostatic equilibrium in the vertical direction. Whereas in the C-G model, the disk is divided into two discrete layers – super-heated and interior – the model of D’Alessio et al. gives continuous distributions of temperature and density. The differences in the temperature and density distributions between these two models have a significant effect on the gaseous molecular abundances in the disk. Since a vertical distribution of density in the super-heated layer is not given explicitly in the C-G model, Willacy & Langer (2000) assumed a Gaussian density distribution with a scale height determined by the mid-plane temperature. In the model of D’Alessio et al., the gas is more extended than in a Gaussian distribution owing to the high temperature in the surface region. The higher densities (and thus the higher column densities) at large $`Z`$ shield the lower layers from stellar ultraviolet radiation. In addition, the temperature variation in the vertical direction is more gradual in the model of D’Alessio et al. than the step function assumed in the C-G model. Therefore, compared with the C-G model, the model of D’Alessio et al. contains more gas in a warm and shielded layer, in which high abundances of gaseous molecules are expected.

The rest of the paper is organized as follows. In §2 we describe the adopted model for protoplanetary disks and the chemical reaction network. Numerical results on the distributions of molecular abundances and column densities are discussed in §3. In §4, molecular column densities and line intensities in the D’Alessio et al. model are compared with observations. Our conclusions and a discussion are given in §5.

## 2 Model

The disk model by D’Alessio et al. (1998, 1999) has been adopted in our work. In this model, the two-dimensional structure ($`R,Z`$) of the disk is obtained by solving for hydrostatic equilibrium and energy transport in the $`Z`$-direction. Various heating sources are considered, such as viscous dissipation of accretion energy, cosmic rays, and stellar radiation. Among them, stellar radiation is dominant at $`R>2`$ AU, assuming typical parameters of T Tauri stars: $`T_{}=4000`$ K, $`M_{}=0.5`$ M, and $`R_{}=2`$ R. We adopt a disk with an accretion rate $`\dot{M}=10^8`$ M yr-1 and viscosity parameter $`\alpha =0.01`$ as the fiducial, or standard, model, but also consider cases with the same $`\alpha `$ and differing accretion rates $`\dot{M}=10^7`$ M yr-1 and $`\dot{M}=10^9`$ M yr-1 (D’Alessio et al. 1999). The distributions of density and temperature in these three models are shown in Fig. 1. The surface density in the fiducial model is similar to that in the minimum-mass disk, which was adopted by Aikawa & Herbst (1999a), and is approximately an order of magnitude larger (smaller) in the model with the larger (smaller) $`\dot{M}`$. In the fiducial disk, the masses inside 100 AU and 373 AU are about 0.017 M and 0.06 M, respectively.

We have selected a few radial points from the model of D’Alessio et al. (viz. $`R=`$ 26, 49, 105, 198, 290, and 373 AU), divided the disk at each radius into several (30-40) layers, or slabs, depending on height $`Z`$, and calculated the molecular abundances in each layer as functions of time (Aikawa & Herbst 1999a). We have not included any hydrodynamic motions in the disk, such as accretion or turbulence. The main goal of this paper is to investigate the effect of a vertical temperature gradient on molecular abundances and line intensities. Although the model of D’Alessio et al. is an accretion disk model, the temperature distribution is determined by the irradiation from the central star and radiation transfer in the vertical direction, while the contribution of the accretion energy as a heat source is negligible in the region we are concerned with. In addition, we find that the chemical time scale is shorter than the accretion timescale, which is $`10^6`$ yr, in a large fraction of the molecular layers (section 3.2), so that molecular distributions obtained in this paper can be a reasonable approximation of reality.

The chemical model and chemical reaction equations adopted in this paper are almost the same as those described in Aikawa & Herbst (2001). We use the “new standard model” for the gas-phase chemistry (Terzieva & Herbst 1998, Osamura et al. 1999), extended to include deuterium chemistry (Aikawa & Herbst 1999b). The ionization rate by cosmic rays is assumed to be the “standard” value in molecular clouds, $`\zeta =\mathrm{1.3\hspace{0.33em}10}^{17}`$ s-1, because the attenuation length for cosmic-ray ionization is much larger than the column densities in the outer regions of the disks, with which we are concerned (Umebayashi & Nakano 1981). Photoprocesses induced by ultraviolet radiation from the interstellar radiation field and from the central star are included. The ultraviolet flux from the central star varies with time and object, and can reach a value $`10^4`$ times higher than the interstellar flux at $`R=100`$ AU (Herbig & Goodrich 1986, Imhoff & Appenzeller 1987, Montmerle et al. 1993). This maximum value is adopted in this paper, as in Aikawa & Herbst (1999a). We assume that the ultraviolet radiation from the central star is not energetic enough to dissociate CO and H2. Self- and mutual shielding of H2 and CO from interstellar UV is considered as in Aikawa & Herbst (1999a). Chemical processes induced by X-rays from the central star are not included in this paper.

Regarding gas-grain interactions, the surface formation of H2, the surface recombination of ions and electrons, and the accretion and thermal desorption of ice mantles are included, but all other grain-surface reactions are not included. The sticking probability for accretion is assumed to be 1.0, unless stated otherwise. We have not considered any modifications to the grain-surface rate equation for H2 formation (cf. Caselli et al. 1998), because, given our initial conditions, the rate of molecular hydrogen formation is not important to the model. The total number of species and reactions included in our network are 773 and 10446, respectively.

The adopted elemental abundances are the so-called “low-metal” values (e.g., Lee et al. 1998, Aikawa et al. 1999). The initial molecular abundances are obtained from a model of the precursor cloud with physical conditions $`n_\mathrm{H}=\mathrm{2\hspace{0.33em}10}^4`$ cm-3 and $`T=10`$ K at $`\mathrm{3\hspace{0.33em}10}^5`$ yr, at which time observed abundances in pre-stellar cores such as TMC-1 are reasonably well reproduced (Terzieva & Herbst 1998).

## 3 Results

### 3.1 Vertical Distribution

Figure 2 contains assorted vertical distributions at the outermost radius $`R=373`$ AU in our fiducial disk model. The solid lines in Fig. 2 ($`ac`$) show the physical parameters: temperature, density, and Av. The attenuation of interstellar radiation (A$`{}_{}{}^{\mathrm{IS}}{}_{\mathrm{v}}{}^{}`$) is obtained from the equation

$$A_\mathrm{v}=\frac{N_\mathrm{H}}{\mathrm{1.59\; 10}^{21}\mathrm{cm}^2}\mathrm{mag},$$
(1)
where $`N_\mathrm{H}`$ is the vertical column density of hydrogen nuclei from the disk surface to each point in the disk. The attenuation of stellar radiation (A$`{}_{}{}^{\mathrm{star}}{}_{\mathrm{v}}{}^{}`$) is obtained via the same equation, but with $`N_\mathrm{H}`$ replaced by the column density from the central star. Fig. 2 ($`de`$) shows assorted molecular abundances at a disk age of $`t=\mathrm{1.0\hspace{0.33em}10}^6`$ yr, which is typical for T Tauri stars. Significant amounts of molecules exist in the gas phase due to thermal desorption at $`Z>100`$ AU, while most molecules are adsorbed onto grains below this height, where $`T<20`$ K. As discussed by van Zadelhoff et al. (2001), this molecular layer covers the region in which the lines of the main isotopes of the observed species become optically thick and thus where most of the observed emission arises. The results for models with different $`\dot{M}`$ are similar except that the height of the molecular layer is shifted in accordance with the distribution of $`A_\mathrm{v}^{\mathrm{IS}}`$, which is the main determinant of the vertical temperature distribution.

The vertical distribution of the physical parameters (dashed lines) and molecular abundances in the Kyoto model are also shown in Fig. 2 ($`ac`$, and $`f`$) for comparison. The mass of the central star in the latter model is 0.5 M, as in the D’Alessio et al. model, so that the density distribution is modified from that adopted in Aikawa & Herbst (1999a). The sticking probability of neutral species onto grain surfaces is assumed to be 0.03, as in Aikawa & Herbst (1999a). It can be seen that the density distribution in the D’Alessio et al. model is more extended than the Gaussian profile assumed in the Kyoto model (as well as in Willacy & Langer 2000), causing more efficient ultraviolet shielding of the warm layers just below the surface. Although we have not calculated the molecular distributions in the C-G model, the width of its warm molecular layer at $`R=373`$ AU can be estimated. In the C-G model, the Gaussian density distribution is similar to that in the Kyoto model, but the disk surface with $`A_\mathrm{v}^{\mathrm{star}}<1`$ mag is “super-heated” by stellar radiation. At $`R=373`$ AU, the boundary between the super-heated layer and the interior region is located at $`Z280`$ AU, a value estimated from the $`Z`$A$`{}_{}{}^{\mathrm{star}}{}_{\mathrm{v}}{}^{}`$ relation in the Kyoto model (Fig. 2 $`c`$). Since the temperature in the interior region is lower than 20 K at this radius in the C-G model, warm gaseous molecules exist only at heights larger than 280 AU. The upper boundary of the molecular layer, which is determined by photoprocesses, is estimated to be $`Z350`$ AU, again from the Kyoto model (Fig. 2 $`f`$). Therefore in the C-G model, the warm molecular layer, if any, is much narrower than in the D’Alessio et al. model. This estimate is consistent with the conclusion of Willacy & Langer (2000) that they need non-thermal desorption in the cold, more shielded layer to account for the observed molecular abundances within the C-G model.

### 3.2 Radial Distribution of Column Densities

Column densities are obtained by integrating the molecular abundances in the vertical direction. The column densities of assorted species as functions of disk radius are shown in Fig. 3 for three accretion disk models at $`t=\mathrm{1.0\hspace{0.33em}10}^6`$ yr with different mass accretion rates. Stable neutrals such as H2CO and HCN show little dependence on radius, because these molecules are abundant only in regions with certain physical conditions and the mass contained in the layer with these physical conditions does not vary much with radius. For example, the HCN abundance is high ($`n`$(HCN)/$`n_\mathrm{H}10^{10}10^9`$) only when $`n_\mathrm{H}<\mathrm{3\hspace{0.33em}10}^7`$ cm-3, $`T>`$ 20 K, and $`A_\mathrm{v}^{\mathrm{IS}}>0.2`$ mag (see Fig. 2). The critical temperature of $`20`$ K is not the sublimation temperature of HCN, but that of CO, which is the dominant form of carbon in the gas phase. For HCN, attenuation of interstellar radiation is more important than that of stellar radiation because the interstellar radiation penetrates deeper into the disk due to the effect of geometry. Radicals such as CN and C2H increase in column density with radius because of the lower density and lower flux of the destructive stellar UV in the outer regions. Radical column densities are more sensitive than HCN to stellar UV, because their abundances peak at greater heights. Carbon monoxide is abundant ($`n`$(CO)/$`n_\mathrm{H}10^4`$) in regions with $`T>20`$ K and $`A_\mathrm{v}^{\mathrm{IS}}>0.1`$ mag. The column densities of CO and HCO+ abruptly change at $`R100`$ AU, inside of which the temperature in the midplane is higher than 20 K. These characteristics of the radial distribution are similar to those in the Kyoto model (Aikawa et al. 1996, Aikawa & Herbst 1999a).

The amount of gas existing under physical conditions conducive for large abundances of molecules does not vary significantly among the three disk models with different accretion rates (and thus with different disk mass), either. Therefore, most (but not all) molecular column densities vary only by a very small factor among the three disk models, even though the total (H2) column density varies by two orders of magnitude. In the model with a mass accretion rate of $`\mathrm{1.0\hspace{0.33em}10}^9`$ M yr-1, the region with density $`n_\mathrm{H}10^510^6`$ cm-3, at which CN is abundant, is more shielded from stellar UV (Fig. 1), and thus the CN column density is higher than in the other two models.

Molecular column densities in our fiducial disk model at an earlier time of $`t=\mathrm{1.0\hspace{0.33em}10}^5`$ yr are shown with thick dot-dashed lines in Fig. 3. The variation in column density for most molecules during $`10^510^6`$ yr is less than a factor of 2; two exceptions are CS and H2O. In regions with $`A_\mathrm{v}^{\mathrm{IS}}`$ smaller than a few mag, which covers a large fraction of the molecular layer, the chemical timescales are short ($`<10^5`$ yr) because of photoprocesses. In the more shielded portion of the molecular layer at smaller $`Z`$, S-bearing gaseous molecules decrease in abundance after $`10^5`$ yr, since most sulfur is adsorbed onto grains in the form of CS, SO and OCS. Similarly, H2O gas decreases at $`10^6`$ yr, because most oxygen which is not in CO is adsorbed as H2O ice. On the other hand, abundances of other C-bearing molecules reach pseudo-steady-state values in a relatively short timescale ($`<10^4`$ yr), and do not show significant time variation during $`10^510^6`$ yr in the more shielded region because CO gas, their chemical precursor, remains the dominant component of carbon for more than $`\mathrm{1\hspace{0.33em}10}^6`$ yr. The pseudo-steady-state gas-phase chemistry of non-volatile carbon-containing species is balanced for a considerable period by formation reactions starting from CO and depletion onto the dust particles.

Values of the sticking probability $`S`$ are estimated to lie in the range $`0.1<S<1.0`$ (Williams 1993, and references therein). In order to check the dependence of the molecular column densities on $`S`$, we have performed calculations with $`S=0.1`$ in addition to our fiducial value of unity. Column densities of radicals such as CN and C2H do not depend on $`S`$, because they are abundant in the surface layer, in which adsorption is not the dominant process. Among the more stable species, some show significant dependence on $`S`$; the column densities of CO2 and OCS are larger by an order of magnitude in the model with the lower sticking probability at $`R=373`$ AU and $`t=\mathrm{1\hspace{0.33em}10}^6`$ yr. But the effect of $`S`$ is smaller for many other species; at $`R=373`$ AU and $`t=\mathrm{1\hspace{0.33em}10}^6`$: for example, the column densities of CS, HCN, and H2CO are larger only by factors of 3.2, 1.7, and 1.3, respectively, in the case of lower $`S`$. There are two reasons for this small dependence on $`S`$. First, adsorption is not always the dominant process in the molecular layer, depending on species and height from the midplane. Second, the higher abundance of gaseous O2 in the case of lower $`S`$ reduces the abundance of the carbon atom, and thus reduces the formation rate of organic molecules, which counteracts the lower adsorption rate.

### 3.3 D/H Ratios in Molecules

In the disks around LkCa15 and TW Hya (Qi 2000, Thi et al. in preparation), the deuterated species DCN and DCO+ have been detected, and the ratios of deuterated to normal species (DCN/HCN and DCO+/HCO+) are estimated to be $`0.01`$. The deuterated species HDO has been detected in LkCa15, but there is no estimated ratio of HDO to normal water. Although the estimated ratios might include uncertainties as discussed in §4, they would not alter the important conclusion that the molecular D/H ratios are higher than the cosmic elemental abundance ratio (D/H$`\mathrm{1.5\hspace{0.33em}10}^5`$) by orders of magnitude. It is interesting to see if our model reproduces the high D/H ratios and if these ratios can be used as probes of the chemical processes and physical conditions in disks. Fig. 4 shows some column density ratios of deuterated species to normal species in our models as functions of radius. The ratios are similar to the observed values of 0.01 for both HCN and HCO+. The mechanism of deuterium enrichment in disks is similar to that in molecular clouds; due to energy differences between deuterated and normal species, molecules such as H$`{}_{}{}^{+}{}_{3}{}^{}`$ have a high D/H ratio, which propagates to other species through chemical reactions (Millar, Bennet & Herbst 1989, Aikawa & Herbst 1999b, Roberts & Millar 2000). The D/H ratios decrease as the temperature rises, because the energy differences are less significant at higher temperatures. Thus, the D/H ratio decreases inwards (Fig. 4), because the temperature in the gaseous molecular layer is higher at inner radii. Because the temperature in the molecular layer does not depend much on the disk mass (see Fig. 1), the differences between models with different $`\dot{M}`$ are small.

The ratio of HDO/H2O shows a more complicated behavior than described above. It has a local peak at $`R100`$ AU, and the peak is higher in the more massive disk model. At $`R100`$ AU, the midplane is at 14-16 K for the three disk models, at which temperature CO is almost completely frozen onto grains but O can marginally be kept in the gas phase to produce water vapor. The D/H ratio in the vapor is enhanced by the CO depletion in the region, since CO is one of the main destruction channels of H2D+ (Brown et al. 1989). The layer at this critical temperature (14 K$`<T<16`$ K) is thicker in the more massive disk. At smaller radii, the midplane temperature is higher, which lowers the D/H ratio. At larger radii for disks with $`\dot{M}=10^7`$ and $`10^8`$ M yr-1, the midplane temperature is so low that O cannot remain in the gas phase to produce water vapor. The abundance of HDO has a sharp peak in a thin layer of 14 K $`<T<16`$ K offset from the midplane (Fig. 2 $`d`$). In the outer disk of the $`\dot{M}=10^9`$ M yr-1 model, the temperature is slightly higher than 20 K even in the midplane (see Fig. 1), which causes a lower HDO/H2O ratio than for the other two models.

## 4 Comparison with observation

### 4.1 Column densities

In the subsequent paragraphs, we discuss comparisons of our calculated column densities for assorted species at a particular radius with both single dish and interferometric data. It must be recognized, however, that it is not really possible to derive reliable columm densities at a particular radius from unresolved single dish data, making much of the subsequent discussion less quantitative than would be desirable.

Table 1 compares calculated molecular column densities obtained with three different accretion rates at a time of $`t=\mathrm{1.0\hspace{0.33em}10}^6`$ yr and a radius of $`R=373`$ AU with those estimated from the observations of DM Tau (Dutrey et al. 1997) and LkCa15 (Qi 2000, Thi et al. in preparation). In general, it is difficult to estimate the total H2 column density in the disk, especially at the outer radius. Dust observations suffer from uncertainties in the dust opacity at millimeter wavelengths, which depends on the grain size distribution. Moreover, the dust continuum is very weak at the outer radii and difficult to detect. The H2 column density cannot be directly estimated from molecular lines, because the molecular abundances relative to hydrogen are not known. Therefore Dutrey et al. (1997) estimated the H2 column density and averaged molecular abundances for DM Tau in a different way, paying attention to the critical density for excitation of the molecular lines and their optical depth. For the DM Tau disk, interferometric observations of 12CO ($`J=21`$) show that the CO gas extends to $`800`$ AU from the central star. Combining the different constraints, they obtained an H2 density $`n(R,Z)=\mathrm{5\; 10}^5(R/500AU)^3`$ exp$`[(Z/H)^2]`$ cm-3 over this range, in which $`H`$ is the scale height of the disk, with $`H`$ 175 AU at $`R=500`$ AU. From the line intensities obtained in single dish telescopes, they subsequently derived molecular abundances with respect to hydrogen assuming that the abundances are constant over the entire disk. We obtained the molecular column densities of DM Tau in Table 1 by vertically integrating this disk model using the average abundances listed in Table 1 of Dutrey et al. (1997).

Although the DM Tau disk extends to about 800 AU from the central star, and the outer radius would have a larger contribution to the line intensities because of the larger surface area, we list the calculated column density at $`R=373`$ AU, since the D’Alessio et al. model extends only out to this radius (the H2 column density at $`R=800`$ AU is about 3 times less than that at $`R=373`$ AU in the DM Tau disk). The model with an accretion rate $`\dot{M}=10^9`$ M yr-1 reproduces reasonably well the vertical H2 column density of DM Tau, and agrees with other DM Tau observations to within a factor of 2, except for CO and C2H. The former is overestimated, while the latter is underestimated. It should be noted that we have assumed that the stellar UV is not hard enough to dissociate CO, so that our model may overestimate the CO column density. Inclusion of CO photodissociation by stellar UV might also increase C2H and CN abundances, since more carbon is released from CO in the photodissociation region, in which C2H and CN are mainly formed (van Zadelhoff et al. in preparation).

Qi (2000) performed interferometric observations on LkCa15 and estimated beam averaged molecular column densities based on the velocity-integrated intensity measured over a much smaller beam than used to determine the DM Tau abundances. The vertical column densities are lower than the listed values by a small factor, which is less than 2 if the inclination is $`<60`$ degrees. Since the beam size is about $`0.6^{\prime \prime }13^{\prime \prime }`$ depending on the frequency of the line, the estimated values do not necessarily refer to 373 AU. If the emission is resolved and if we assume that the molecular column densities do not vary much with radius (as indicated by our disk models for $`R>100`$ AU), we can take the values listed by Qi (2000) to apply to 373 AU. The typical beam size of Qi (2000) is $`3^{\prime \prime }4^{\prime \prime }`$ ($``$300 AU radius at the distance of LkCa15), so the majority of the results will be just resolved, making it reasonable if not perfect to compare the observations with our result at $`R=373`$ AU.

In addition to the interferometric work on LkCa15, Thi et al. (in preparation) derived molecular column densities for this disk from high-frequency single dish observations of various species. Their beam-averaged values obtained with their assumption that the disk has a radius of 100 AU are included in Table 1 and differ from the interferometric column densities by up to an order of magnitude. Both values are much higher than the column densities found for DM Tau, so the agreement between our model and observations is worse for LkCa15 than for DM Tau. For CO and HCO+ the difference is less than a factor of 2, but the column densities of other species are about 1–2 orders of magnitude higher than the model results for all three mass accretion rates. Methanol in LkCa15 appears to be a special case; its calculated column density is more than five orders of magnitude too low, presumably because it is produced by the hydrogenation of CO on grain surfaces, which is not included in our model, and then evaporated into the gas.

The use of 100 AU for LkCa15 as a reference point by Thi et al. is somewhat arbitrary. If the values obtained by Thi et al. are taken to refer to a 373 AU radius disk, their LkCa15 column densities need to be reduced by a factor (100/373)2, and become closer to those found for DM Tau. Thi et al. and Qi (2000) also present data for a few other disks (TW Hya, HD 163296, and MWC 480) and determine column densities and abundances in a consistent way. Indeed, the DM Tau column densities and abundances are generally lower than the values for other disks, whereas those for LkCa15 are among the highest. Thus, these two sources appear to bracket the range of observed values.

### 4.2 Line intensities and profiles

In addition to the comparison with estimated molecular column densities at a single radius, a more direct comparison via line intensities from the entire model disk has been made. Since the most complete single-dish and interferometer data set is available for LkCa15, we restrict our efforts to this source. In estimating the molecular column densities from the observations, Qi (2000) assumed local thermal equilibrium (LTE) with an excitation temperature of 40 K throughout the LkCa15 disk, thereby deriving a mean column over the beam. In this paper, however, we have shown that a temperature gradient is important for the characteristics of the gaseous molecular layer in the disk and that the abundances are strongly varying with $`R`$ and $`Z`$. Also, the excitation of the molecules and the line emission depend on the density structure; van Zadelhoff et al. (2001) have shown that the assumption of LTE is not always valid for high frequency lines of molecules with high critical densities. Hence, it is better to compare directly simulated line emission with observations.

We calculate here the excitation of the molecules using the 2-dimensional (2D) NLTE molecular line radiative transfer code of Hogerheijde & van der Tak (2000). From the resulting level populations, the line profiles can be computed, taking into account the inclination of the source. NLTE molecular line radiative transfer in more than one dimension has been used in star-formation research only recently (e.g., Park & Hong 1995, Juvela 1997). The need for a full treatment of the radiative transfer in disks follows from the non-locality of the problem. The level populations depend both on the local parameters (temperature, density and radiation field) and on the global radiation field, which in turn depends heavily on the optical depth of the medium. The main problem is the slow convergence of the level populations. The adopted code by Hogerheijde & van der Tak (2000) uses an accelerated Monte Carlo method to enhance convergence. A more elaborate discussion of the methods is given in van Zadelhoff et al. (2001).

The main input parameters to the radiative transfer calculations in addition to the molecular data (collision rates, Einstein A coefficients) are the 2D distributions of temperature, density and abundance of the molecule of interest at each position ($`R,Z`$). The latter is taken directly from the chemical models at a chosen age. Other important parameters are the turbulent width, taken to be 0.2 km s-1 (van Zadelhoff et al. 2001), the systematic Keplerian velocity field for an assumed stellar mass, the thermal line width, the size of the disk, and the inclination of the source. Qi (2000) estimates a disk inclination of $`58\pm 10`$ degrees, and an outer radius of the CO disk of 435 AU for LkCa15. In the following, we assume a disk inclination of 60 degrees and a size (outer radius) of 400 AU. The dependence of line emission and profiles to disk radius and inclination can be roughly estimated (Omodaka et al. 1992, Beckwith & Sargent 1993), e.g., for optically thick species the peak intensities of molecular lines are proportional to the square of the disk radius. Dependence of the integrated intensities on disk inclination is also discussed in van Zadelhoff et al. (2001), and is found to be less than a factor of two when the inclination is varied from 0 to 60 degree.

Fig. 5 compares simulated line profiles for CO ($`J=65`$, $`32`$, and $`21`$), HCO+ ($`43`$), HCN ($`43`$), and CN ($`32`$) from the fiducial D’Alessio et al. model (thick solid line) with lines from LkCa15 observed with single-dish telescopes by van Zadelhoff et al. (2001) (solid line). Line profiles based on the D’Alessio et al. model with other accretion rates – $`\dot{M}=10^9`$ (dot-dashed line) and $`\dot{M}=10^7`$ (dashed line) M yr-1 – are also shown. In spite of the fact that the column densities of CO and HCO+ in our model are slightly smaller than those estimated by Qi (2000), the calculated intensities of these species are higher than observed in LkCa15 by a factor of $`23`$, which is caused by the higher disk temperatures in our model. The vertical temperature gradient lowers the optical depth of the disk, which further enhances emission intensities (van Zadelhoff et al. 2001). As opposed to CO and HCO+, the model intensities of CN and HCN are much lower than those observed in LkCa15. Those lines are optically thin in our model, and about 10 times more HCN and at least 50 times more CN are needed to fit the observed profiles, which is consistent with the comparison of column densities. The dotted lines in the CN and HCN panels show model profiles in which the molecular abundances are artificially enhanced. We conclude that CN and HCN in LkCa15 are much more abundant than in our model.

### 4.3 Discussion

What are the causes of the disagreement between our model results and LkCa15, and the difference between DM Tau and LkCa15 ? Let us first discuss the uncertainties in the estimate of the molecular column densities and size of the emission region. Although there are some differences in observed intensities, the lines in LkCa15 are, in fact, not much stronger than those in DM Tau, at most a factor of 2–3 after scaling the data to the same beam. The main reason for the different column densities in DM Tau and LkCa15 then seems to be the adopted size of the emission region, because the column densities are inversely proportional to the square of the size (radius) of the emission region. For CO, where interferometric observations have been performed both for DM Tau and LkCa15, the DM Tau disk is found to be somewhat larger than the LkCa15 disk — DM Tau is about 800 AU in radius, while LkCa15 is 435 AU. But whether this larger emission region in DM Tau holds for other species is not yet known without some degree of ambiguity because the DM Tau data come mainly from single-dish observations, and the size of the region and average molecular abundances are estimated from a model with constant abundances throughout the disk. For LkCa15, in contrast, the sizes of the emission regions are directly derived from interferometric observations and are at most a few arcseconds. Assuming that the interferometer resolves the emission, these column densities should therefore not be affected by uncertainties in the size of the emission region. Although interferometry is more direct, and thus seems more reliable for obtaining the size of the emission region, it should be noted that the derived size depends on the signal to noise ratio and dynamic range achieved in the observations. In fact, the size of the emission region obtained via interferometry seems to contain uncertainties; the size of the CO ($`J=21`$) emission region around LkCa15 is estimated to be $`600`$ AU by Duvert et al. (2000), which is larger than the 435 AU obtained by Qi (2000).

Since our 2D modeling procedure convolves the calculated molecular emission with the actual beam of the observations, a larger disk size cannot explain the discrepancy with the LkCa15 interferometer observations, but it does affect the calculated single dish intensities. Such an increase would not be sufficient to explain the discrepancies for LkCa15, however. For example, if we assume an outer disk radius of 600 AU, the calculated HCN and CN line intensities from our model disk (Fig. 5) are increased by a factor of 2 at most, while the calculated CO and HCO+ lines become even stronger compared with the observed profiles. There are a few possible explanations for these discrepancies. First, our model might overestimate the CO column density and underestimate the radical column density because we do not consider dissociation of CO via stellar UV radiation, as mentioned above. Detailed consideration of stellar UV radiation, including the CO and H2 dissociation, will be reported in a forthcoming paper. Indeed, CN is enhanced by an order of magnitude depending on the treatment of the radiation field, although HCN is not much changed. The inclusion of X-rays might be another solution. X-rays cause ionization and dissociation, which enhance chemical activity, and hence increase the transformation of CO to other organic molecules such as CN and HCN (Aikawa & Herbst 1999a, 2001). Different X-ray fluxes might also account for the differing molecular column densities in DM Tau and LkCa15, if they are intrinsic. X-rays also cause non-thermal desorption, which might enhance the CN and HCN abundances in the gas phase (Najita, Bergin & Ullom 2001). Finally, of the disks around T Tauri stars surveyed so far, LkCa15 stands out as the disk with the strongest molecular lines and richest chemistry in the interferometer and single-dish data (see §4.1, Qi 2000, Thi et al. in preparation).

## 5 Conclusion and Discussion

We have investigated the molecular distributions in protoplanetary disks by combining the “new standard” chemical model with the physical model of D’Alessio et al., which has a temperature gradient in the vertical direction. The calculated molecular column densities are in reasonable agreement with those estimated from single-dish observations of the DM Tau disk without the assumptions in previous calculations of non-thermal desorption and/or an artificially low sticking probability. In the warmer intermediate layers of our disk models, there are large amounts of gaseous molecules owing to thermal desorption and efficient UV shielding caused by large gas densities at large heights from the midplane, a phenomenon known as flaring. Gaseous molecules are abundant only in regions with certain physical conditions. The volume of the layers with these conditions, and thus the column densities of gaseous molecules, are not proportional to the total (H2) column density. Column densities of abundant molecules such as CN, HCN and HCO+ do not vary by more than a factor of three during the period $`t10^510^6`$ yr. Sulfur-bearing molecules and H2O show larger temporal variations.

Comparison of our model results with those of Willacy & Langer (2000), who adopted the C-G disk model with a Gaussian density distribution, indicates that gaseous molecular abundances are sensitive to the vertical structure of the disk model; in their model, molecules in the super-heated upper layer are destroyed by the harsh ultraviolet radiation from the star, while sufficient UV shielding is available in the warm upper layers of the D’Alessio et al. model. Deuterated species are also included in our chemical model. The molecular D/H ratios we obtain are in reasonable agreement with those observed in protoplanetary disks.

Despite our agreement with observations of DM Tau, the molecular column densities obtained in our models are smaller than those observed around LkCa15, except for CO and HCO+. The estimated column densities of all observed molecules around LkCa15 are higher than those around DM Tau by roughly an order of magnitude. This difference in derived molecular column densities in the two objects seems to derive, at least partially, from different and/or uncertain sizes of the emission regions. Comparison with other sources shows that some of the difference is likely to be intrinsic, and other physical parameters or processes, such as X-ray ionization and dissociation, are needed to account for the high column densities in the LkCa15 disk.

In addition to the calculation of molecular abundances and column densities, we have solved the equation of radiation transfer to obtain line profiles from our model disks, which can be directly compared with observations. Such a comparison is a much more detailed test of theory than is a comparison of column densities, since line intensities depend not only on the molecular column densities but also on the density and temperature of the molecular layer and the variation of the abundance of the molecule with $`R`$ and $`Z`$. The line intensities of HCN and CN obtained from the theoretical models are lower than the observed intensities in LkCa15, as expected from the comparison of column densities.

There are still several uncertainties in the vertical structure of protoplanetary disks, which might affect our results. Firstly, Chiang & Goldreich (1997) argue that the gas temperature could be lower than the dust temperature in the upper layers; although an important heat source of the gas is collisions with super-heated dust particles, gas-dust collisions are not frequent enough to equilibrate the gas temperature and grain temperature because of the low density. With a lower gas temperature, the disk would be less flared than in the model of D’Alessio et al., which assumes equal temperatures of gas and dust. Glassgold & Najita (2001) have pointed out, however, that gases in the upper layer can be heated by X-rays, which were not considered by Chiang & Goldreich (1997). Since Glassgold & Najita (2001) have listed the temperature only in the inner radius ($`R=1`$ AU), we made a rough estimate of the disk surface temperature for the outer radii $`R100300`$ AU based on the work of Maloney et al. (1996), which suggests that the surface temperature of the X-ray irradiated disk is indeed higher than the midplane (interior) temperature of the C-G and Kyoto models. Moreover, the UV photons can heat the gas through the photoelectric effect on grains and PAHs, as in models of photon-dominated regions. Hence we can at least expect disks to fall off more slowly in density with increasing height than in a simple Gaussian distribution, although detailed studies on the heating and cooling balance between gas and dust are desirable in order to obtain an accurate vertical structure for the disk. Another uncertainty lies in the size and distribution of dust particles. D’Alessio et al. (2001) and Chiang et al. (2001) have noted that their original models are geometrically too thick compared with the observations of edge-on disks, which suggests dust sedimentation and/or growth in the disk. Because the molecular abundances in our model depend on the efficiency of UV shielding by “small” (i.e. interstellar) dust grains (Aikawa & Herbst 1999a), we might have overestimated these abundances. Although a more detailed approach with dust sedimentation is beyond the scope of this paper, we emphasize that molecular abundances can help to resolve uncertainties in dust evolution and disk structure.

###### Acknowledgements.

The authors are grateful to P. D’Alessio for providing numerical tables of her models, to G. Blake, C. Qi and W.F. Thi for results of their observations prior to publication, and to M. Hogerheijde and F. van der Tak for use of their 2D Monte Carlo code. Astrochemistry in Leiden is supported through a Spinoza grant from the Netherlands Organization for Scientific Research (NWO). Astrochemistry at Ohio State is supported through a grant from the National Science Foundation. Numerical calculations were carried out at the Astronomical Data Analysis Center of National Astronomical Observatory in Japan.


format: text

text:
# SUB-MILLIMETER TESTS OF THE GRAVITATIONAL INVERSE SQUARE LAW

## 1 What might be special about gravity at length scales below 1 mm?

Very little is known about gravity at length scales below a few mm . Recently theorists, using several different arguments, have suggested that the unexplored short-range regime of gravitation may hold profound surprises , i.e. that the gravitational interaction could display fundamentally new behavior in the mm regime.

Many of these arguments are based on the notion, inherent in string or M theory, of more than 3 spatial dimensions. To maintain consistency with a vast body of observations the extra dimensions must be “curled up” in very small regions, usually assumed to be comparable to $`R_\mathrm{P}=\sqrt{G\mathrm{}/c^3}=1.6\times 10^{33}`$ cm, or else hidden in some other way . It has recently been noted  that the enormous discrepancy between natural mass scales of the Standard Model of particle physics ($`M_{\mathrm{SM}}1`$ TeV) and of gravity (the Planck mass $`M_\mathrm{P}=\sqrt{\mathrm{}c/G}=1.2\times 10^{16}`$ TeV) could be eliminated if gravity propagates in all the space dimensions while the other fundamental interactions are constrained to the three familiar dimensions. This unification scenario requires that some of the extra dimensions have radii $`R^{}`$ that are large compared to $`R_\mathrm{P}`$ with

$$R^{}=\frac{\mathrm{}c}{M^{}c^2}\left(\frac{M_\mathrm{P}}{M^{}}\right)^{2/n},$$
(1)
where $`M^{}`$ is the unification scale (usually taken as $`M_{\mathrm{SM}}`$) and $`n`$ is the number of large extra dimensions. The scenario with $`n=1`$ is ruled out by astronomical data. If there are 2 large extra dimensions, $`R^{}`$ must be about 1 mm, and the gravitational inverse-square law (which follows from Gauss’s Law in 3 spatial dimensions) will turn into a $`1/r^4`$-law (Gauss’s Law in 5 dimensions) at distances much smaller than $`R^{}`$.

Completely independent theoretical considerations also suggest that new effects may appear at short distances; string theories predict scalar particles (dilatons and moduli) that generate Yukawa interactions which could be seen in tests of the $`1/r^2`$ law. If supersymmetry is broken at low energies these scalar particles would produce mm-scale effects . Finally, there may be some significance to the observation  that the gravitational cosmological constant, $`\mathrm{\Lambda }3`$ keV/cm3, deduced from distant Type 1A supernovae  corresponds to a length scale $`\sqrt[4]{\mathrm{}c/\mathrm{\Lambda }}0.1`$ mm. These, and other, considerations suggest that the Newtonian gravitational potential should be replaced by a more general expression

$$V(r)=G\frac{m_1m_2}{r}(1+\alpha e^{r/\lambda }).$$
(2)
The simplest scenario with 2 large extra dimensions predicts $`\lambda =R^{}`$ and $`\alpha =3`$ or $`\alpha =4`$ for compactification on an 2-sphere or 2-torus, respectively , while dilaton and moduli exchange could produce forces  with $`\alpha `$ as large as $`10^5`$ for Yukawa ranges $`\lambda 0.1`$ mm.

## 2 Experimental Results

In February 2001 we published results of an inverse-square law test  obtained with a novel torsion pendulum/rotating attractor instrument. The active component of the pendulum was an aluminum ring with 10 equally-spaced holes bored into it. The pendulum was suspended just above a disk-shaped copper attractor that had 10 similar holes bored into it. As the attractor rotated slowly and uniformly underneath the pendulum, it produced a torque on the pendulum that varied back and forth 10 times for every revolution of the attractor. The attractor actually consisted to two concentric disks each with 10 holes: a thinner upper disk and a thicker lower disk. The holes in the lower disk were rotated by 18 degrees with respect to those in the upper disk so that, if inverse-square-law were correct, the torque on the ring from the lower disk canceled the torque from the upper disk. However, the torque from a short-range interaction could not be canceled simply because the lower disk was too far away to produce a short-range torque on the pendulum. We greatly reduced any electrostatic torques on the pendulum by placing a stationary, tightly stretched 20 $`\mu `$m thick Be/Cu membrane between the pendulum and attractor. Our design had several attractive features:

1. the signal occured at a different frequency than the disturbance (the revolution of the attractor). In this case the signals were at $`10\omega `$, $`20\omega `$, and $`30\omega `$ where $`\omega `$ is the attractor rotation frequency.
2. our test bodies were the “missing masses” of the holes in cylindrical rings and disks. This gave us accurately positioned test bodies with planar geometry (optimum because it maximizes the mass that can be placed in close proximity) that could be characterized very precisely.
3. the lower attractor disk that essentially canceled the Newtonian torque greatly reduced our sensitivity to nonlinearities and scale-factor uncertainties in our instrument.

This experiment, which constituted the PhD thesis work of C.D. Hoyle , is described in Ref. . The constraint on short-range Yukawa interactions from Ref.  is shown in Fig. 1. We encountered a surprising problem in the course of this measurement; for a while looked as if we were observing a substantial violation of the inverse-square law. Despite much effort, we could not account for the apparent violation. So we constructed a second 10-hole torsion pendulum and attractor having holes with different diameters and thicknesses to check the original result and again saw an apparent violation of the $`1/r^2`$ law. Blayne Heckel finally identified the problem: the commercial computer-controlled micropositioning stage from which the torsion fiber was suspended had a scale factor error–it actually moved only $`98\%`$ as far as it indicated. So of course we did not find that $`\stackrel{}{}\stackrel{}{g}=0`$; we were using correct distances along $`\widehat{x}`$ and $`\widehat{y}`$ and in incorrect distance along $`\widehat{z}`$! Reference  was based on the results from the first (Mark II) instrument. Hoyle has recently reanalyzed the data from the Mark II instrument as well as that from the second (Mark III) instrument. Figure 1 shows the improved constraint from the new analysis of the combined data.

## 3 Our second-generation instrument

Since the publication of our original results , we have made several improvements to our instrument. These were motivated by the recognition that the torque from a very short-range Yukawa interaction with $`\lambda s`$ ($`s`$ is the closest attainable pendulum-to-attractor separation) scales as

$$T=\frac{\mathrm{\Delta }E}{\mathrm{\Delta }\theta }\rho _p\rho _aA\lambda ^4e^{(s/\lambda )},$$
(3)
where $`E`$ is the interaction energy between the pendulum and attractor, $`\theta `$ is the twist of the pendulum, $`\rho _p`$ and $`\rho _a`$ are the densities of the pendulum and attractor, and $`A`$ is the area of the holes.

1. We increased the torque from a given short-range Yukawa interaction and reduced the Newtonian torque appreciably. This was done by:

* a new pendulum/attractor design that has two rows of 22 holes (to increase $`A`$) with thinner (1 mm thick) pendulum ring and upper attractor disk. The relative sizes of the penduolum and attractor holes was “tweaked” to put most of the power of a Yukawa torque into the fundamental $`22\omega `$ signal.
    * the pendulum and attractor are both made from molybdenum; this increases the $`\rho _p\rho _a`$ product by a factor of 4.5.
2. The minimum attainable spacing $`s`$ is less by at least a factor of two:

* we reduced the thickness of the conducting membrane to $`10\mu `$m.
    * we installed a passive damper that reduced the mean amplitude of the pendulum bounce mode by a factor of 6.
3. We cancelled Newtonian gravity to a higher degree by having thinner active components with more and smaller holes (Newtonian gravity, being long range, tends to average over several holes).
4. We reduced our torque noise by a factor of 6 by improving the autocollimator performance and increasing the pendulum’s quality factor $`Q`$ to $`4000`$.

The pendulum/attractor of our current Mark V instrument is shown in Fig. 2.

## 4 Conclusions

We have tested the gravitational inverse-square law at length scales well below 1 mm. We find no evidence for violations. For the proposed scenario with 2 “large” extra dimensions, our negative result corresponds to a 95% confidence upper limit of 150 $`\mu `$m on the size of the 2 dimensions; this implies a unification mass $`M^{}`$ greater than 4.0 TeV. Alternatively, we can use our results to set an upper limit on the size of the largest extra dimension, regardless of how many there are (subject only to the assumption that one of the extra dimensions is appreciably larger than the rest). In this case $`\alpha =1`$ , so that our data imply with 95% confidence that $`\lambda <200`$ $`\mu `$m. This upper limit constrains models that try to explain phenomena such as neutrino oscillations in terms of large extra dimensions.

We expect that, in the next year or so, our torsion-balance scheme for testing the gravitational $`1/r^2`$ law will provide good results for length scales down to 50 $`\mu `$m. This is less than the diameter of a typical human hair! If we do find evidence for violation of the $`1/r^2`$ law, we would then develop an instrument to check if the $`1/r^2`$ violating interaction also violates the weak equivalence principle. This would distinguish between exotic space-time scenarios of extra-dimensions etc. that do not violate the equivalence principle, from exotic particle exchange scenarios that must violate the principle.

However, until evidence for new physics is found, it is clearly better to work on tests of the inverse-square law than on equivalence-principle tests: the $`1/r^2`$ tests are more general (probing all finite-range effects), and more sensitive (in particle-exchange scenarios the composition-dependence is expected to be a relatively small fractional effect). But testing the gravitational $`1/r^2`$ law for length scales less than 50 $`\mu `$m will probably require a somewhat different technology. In a planar geometry (optimum because one gets the maximum amount of mass in close proximity) the signal of a short-range Yukawa interaction drops as roughly the 4th power of the Yukawa range while extraneous disturbances stay roughly the same size. This will present an interesting challenge for future experimental work.

## Acknowledgments

Jens Gundlach, Blayne Heckel, CD Hoyle, Dan Kapner, Ulrich Schmidt and Erik Swanson all made substantial contributions to the work reported here. We are grateful to Profs. David Kaplan and Keith Dienes for illuminating conversations. This work was supported by the NSF (grant PHY 9970987), and by DOE funding of the Center for Experimental Nuclear Physics and Astrophysics at the University of Washington.

## References


text:
# References

Universality and Scale Invariance in Hourly Rainfall

Pankaj Jain, Suman Jain and Gauher Shaheen

Physics Department, I.I.T. Kanpur, India, 208016

Abstract: We show that the hourly rainfall rate distribution can be described by a simple power law to a good approximation. We show that the exponent of the distribution in tropics is universal and is equal to $`1.13\pm 0.11`$. At higher latitudes the exponent increases and is found to lie in the range 1.3-1.6.

It is well known that daily rainfall over a localized region displays an almost chaotic behavior . This implies that it is very difficult to make even short term prediction of daily rainfall. The distribution of daily rainfall amount over a localized region is also well known to have a very large tail which implies that there is very high probability to have very large rainfall. The distribution shows considerable deviation from the normal distribution.

In the present paper we analyze the distribution of hourly rainfall amount. We are particularly interested in determining whether we can identify some universal features in the rainfall distribution which are observed at all spatial locations. The nature of distribution can also reveal some underlying features of the dynamical system. For example, it has been argued that dissipative systems with large number of variable tend to reach a critical state under the influence of an external perturbation. In the critical state the system displays an almost chaotic behavior with fluctuations displaying scale invariant power law distribution. This has been demonstrated numerically by postulating simplified models of the behavior of sand piles, earthquakes etc . The mechanism by which the dynamical systems achieve the critical state under the influence of external perturbation is called the self organized criticality . In the present paper we shall determine if the rainfall rate also displays this behavior by studying its distribution. Some preliminary results of this study have been presented in Ref. .

The data is taken from the http://daac.gsfc.nasa.gov web site and consists of SSM/I satellite data. The data consists of rainfall rate in units of mm/hour with the minimum rate equal to 0.1 mm/hour. It is available over the entire globe on a grid size of half a degree both in latitude and longitude. We extracted the data files from this website for the year 1997 for three different time periods. At different times the data is available at different locations. We randomly selected a total of 97 different locations on the globe out of the extracted data. In order that a reasonable amount of data is available for a given region, we group the rainfall data in intervals of 5 degrees both in latitude and longitude. The typical distribution for several different regions are shown in Fig. 1 and 2. We represent the power law distribution as

$$f(x)=a/x^b$$
(1)
where $`x`$ is equal to the hourly rainfall amount in mm and $`a`$ and $`b`$ are parameters. If the underlying dynamical system is in a critical state then we expect that certain observables in the system will follow a power law distribution. The rainfall amount per hour represents an interesting physical observable which we study in this paper.

We point out that while testing for a power law behaviour we are not testing the hypothesis that the distribution is a power law over the entire interval of rainfall per hour. The power law is expected to be valid only in some intermediate range of rainfall. At low as well as large rainfall we expect that the distribution will be distorted in comparison to a pure power behavior. At very large rainfall, for example, we expect the distribution to decay very rapidly since rainfall amounts larger than a certain value are physically not possible. We, therefore, cut off the distribution at some large value as well as at a low value of the rainfall. At the low rainfall value we make the distribution constant below a certain value $`x_{\mathrm{min}}`$, which is treated as an additional parameter of the distribution. At the large rainfall end we put a cut on the data so as to exclude all data with $`x>x_{\mathrm{max}}=10`$ mm/hr. The number of data points in the excluded region are typically found to be very small. The distribution is set equal to zero beyond $`x=x_{\mathrm{max}}`$. The power distribution, therefore, contains two adjustable parameters $`b`$ and $`x_{\mathrm{min}}`$. The parameter $`a`$ is obtained by normalization of the distribution.

We compare the fit obtained by the power law distribution with an alternative exponential distribution

$$g(x)=\alpha \mathrm{exp}(\beta x)$$
(2)
In analogy with the power distribution, $`g(x)`$ is also set equal to zero for $`x>x_{\mathrm{max}}`$. Hence the normalization $`\alpha `$ for this distribution is equal to $`\beta /(1\beta x_{\mathrm{max}})`$. We point out that we are not interested in finding the best possible distribution that describes the data. We are only interesting in demonstrating that a power law provides a good description of the data over a wide range of rainfall amounts. Hence we have not made extensive comparisons of various distributions and compare the power law only to a simple exponential fit.

In Fig. 1 and 2 we show the fits in several representative cases. The figure shows the power as well as the exponential fit to the rainfall distribution. It is clear from these plots that a power distribution fits the data reasonably well. It can also be seen that the fit is much better at lower latitudes in comparison to the higher latitudes. The log likelihood difference, defined to be the difference of the log likelihood for the power fit and the exponential fit, is generally found to lie in the range 50 to 150 in tropics and between 10 to 30 at the higher latitudes. Hence the power distribution gives an overall much better description of the data in comparison to an exponential distribution. In Fig. 3 and 4 we show a scatter plot of the exponent of the power distribution as a function of the latitude and longitude respectively. The error in the exponent values ranges from 0.02 to 0.1 in most cases. Only rarely does the error exceeds 0.1. The precise value depends on the number of data points obtained in each region which typically range from 100 to 2000. We also point out that the uneven distribution of the data points as a function of latitude and longitude is caused by our selection of times for which the data is extracted. However this does not have any influence on our results or conclusions since the data is spread over the entire globe.

We clearly see that the values of the exponent $`b`$ are different at lower latitudes in comparison to higher latitudes. The exponent close to equator is found to be close to unity whereas at higher latitudes the exponent is much larger. The exponent does not show very significant dependence on the longitude. Only a very marginal relationship is found after eliminating one outlier with a very large value of the exponent. In order to quantify the dependence of the exponent on the latitude we evaluate the correlation between the absolute value of the latitude and the exponent. For the entire set of 97 data points the correlation $`\rho =0.54`$. If we delete two outliers the correlation goes up to $`\rho =0.608`$. The probability $`p`$ that we can get this correlation from a random sample is very small. We find that for the entire set $`p=10^5`$% and after eliminating the two outliers $`p=3\times 10^7`$%.

The fact that at lower latitude the exponent values are clustered around unity irrespective of longitude is a clear indication of a universal behaviour in the rainfall distributions in the tropics. By putting a cut on the latitude in order to select only the region that lies between 20 N and 20 S we find that the mean value of the exponent $`b`$, defined in Eq. 1, is given by,

$$b=1.13\pm 0.11$$
(3)
The total number of data points contained within the latitude 20 N and 20 S are equal to 61. Two of these points were found to give an anomalously large values of the exponent compared to the rest and were treated as outliers. The mean value given in Eq. 3 is obtained after eliminating these outliers. The median value of the exponent is equal to 1.12 and remains unchanged with or without the inclusion of the outliers.

As this work was near completion we became aware of Ref.  where the authors have analyzed high-resolution rainfall data at Baltic coast Zingst in order to investigate its distribution. The authors find that the distribution is well described by a power law with exponent 1.36. Our results are in agreement with their findings since at higher latitudes we also find similar exponents.

In conclusion, we summarize the main results of the paper. We have shown that the hourly rainfall distribution is well described by a scale invariant power law. This is particularly true in the tropical region where the exponent is found to be $`1.13\pm 0.11`$ independent of the longitude. As we go towards the higher latitudes the exponent increases and generally is found to lie in a range 1.3-1.6. The power dependence and its universal character in the tropics indicates that the underlying dynamical system may be best describable in terms of self organized criticality.

Acknowledgements: We thank Ralf Bennartz, Ashu Jain and Mahendra Verma for useful discussions. We also thank the Distributed Active Archive Center at the Goddard Space Flight Center, Greenbelt, MD, 20771, for producing and distributing the data which was used in this work. Partial funding for this work is provided by DST.

All correspondence should be directed to Pankaj Jain (pkjain@iitk.ac.in).


text:
# 1 Transverse mass spectra of identified hadrons for minimum bias gold-gold collisions as measured for √𝑠=130 AGeV at RHIC (preliminary data taken from []). The data is normalized to a power law fit. The proton and kaon spectra are rescaled by a factor of 1/2 and 2, respectively.

SCALING PROPERTIES OF THE TRANSVERSE MASS SPECTRA
J. SCHAFFNER-BIELICHa, D. KHARZEEVb, L. MCLERRANb and R. VENUGOPALANb,c

aDepartment of Physics, Columbia University, New York, NY 10027, USA
bDepartment of Physics, Brookhaven National Laboratory, Upton, NY 11973, USA
cRIKEN BNL Research Center, Brookhaven National Laboratory,
Upton, NY 11973, USA

## Abstract

Motivated from the formation of an initial state of gluon-saturated matter, we discuss scaling relations for the transverse mass spectra at BNL’s Relativistic Heavy-Ion Collider (RHIC). We show on linear plots, that the transverse mass spectra for various hadrons can be described by an universal function in $`m_t`$. The transverse mass spectra for different centralities can be rescaled into each other. Finally, we demonstrate that $`m_t`$-scaling is also present in proton-antiproton collider data and compare it to $`m_t`$-scaling at RHIC.

In this work, we report on applying scaling features of the initial state of the color glass condensate to data at RHIC, in particular to the transverse momentum distribution of hadrons (see ).

It has been suggested that in collisions of heavy-ions at ultrarelativistic bombarding energies gluons are so densely packed in phase-space that they saturate  and form a color glass condensate . The initial momentum distribution of gluons for a color-glass condensate can be described by the scaling relation :

$$\frac{1}{\sigma }\frac{dN_g}{d\eta d^2p_t}=\frac{1}{\alpha _s(Q_s^2)}f_g\left(\frac{p_t^2}{Q_s^2}\right).$$
(1)
Here, $`\sigma `$ is the transverse area, $`Q_s`$ the saturation scale and $`f_g`$ an universal dimensionless function for the produced gluons. The saturation scale $`Q_s`$ depends on the energy, the centrality and atomic number of the two colliding nuclei.

We have to make several assumptions about the underlying dynamical picture to compare the above ansatz to actual experimental data: the initial gluon distribution is characterized by saturated gluons, there is a free streaming evolution which produces additional partons and finally, there is a freeze-out to hadrons. At all stages, the scaling properties of the initial state are preserved. Hence, we replace the scaling function for gluons and the saturation momentum with those for hadrons, i.e. $`f_gf`$ and $`Q_sp_s`$, respectively. Note, that $`p_s`$ is not equal to the saturation momentum for gluons but should have similar dependencies on energy and centrality. This simple picture works to explain the centrality dependence and the pseudorapidity dependence of the charged multiplicity at RHIC . The question we are asking is, whether the above scaling relation eq. (1) can be utilized to describe the transverse momentum spectra at RHIC. We point out that our assumptions are completely contrary to the notion of radial flow or of thermalized and equilibrated matter.

Guided by the scaling relation for saturated gluons we propose that the distribution of produced hadrons can be cast in the following form:

$$\frac{1}{\sigma }\frac{dN_h}{dyd^2m_t}=\frac{1}{\alpha _s(p_s)}\kappa _hf\left(\frac{m_t}{p_s}\right).$$
(2)
The appropriate replacements, $`f_gf`$ and $`Q_sp_s`$ are made from eq. (1) to eq. (2). In addition, we assume that the transverse momentum has to be replaced by the transverse mass for identified hadrons, in absence of any other boost invariant quantity. The new factors $`\kappa _h`$ take into account effects of conserved quantum numbers, as strangeness and baryon number and are constants. The parameters $`\sigma `$ and $`p_s`$ have to be determined from the data.

We take the minimum bias data of the transverse mass spectra from the PHENIX collaboration as measured for gold-gold collisions at $`\sqrt{s}=130`$ GeV  and plot it as a function of the transverse mass. The curves for different hadron species are very close to each other without any adjustments in the absolute normalization. The spectra is closer to a power law than to an exponential. Shifting the protons up a factor $`\kappa _h=2`$ (there is definitely a finite baryon number at midrapidity) and the kaons down by a factor $`\kappa _h=1/2`$ (there seems to be a suppression of strange particles) puts all curves on top of each other (see  for the corresponding plots). That means that the slope of the curve is the same for a given $`m_t`$ irrespective of the hadron. The $`\kappa _h`$ factors are indeed constants, they do not depend on $`m_t`$. The scaling pattern is present even for quite large values of $`m_t=3`$–4 GeV.

To better see any deviations from $`m_t`$-scaling, we divide the data by a power law fit $`1/(p_s+m_t)^n`$ with $`p_s=2.71`$ GeV and $`n=16.3`$. The value for $`n`$ seems to be high. Nevertheless, we remind the reader that the two parameters of the power law fit are tightly constrained by the mean $`p_t=2p_s/(n3)`$, so that lower values of $`n`$ can be achieved by lowering $`p_s`$ correspondingly with basically the same fit quality. The linear plot of the transverse mass spectra relative to that power law fit is shown in Figure 1. One sees, that the deviations from scaling are only about 30% including systematic errors from the fit!

Next, we consider the centrality dependence of the particle spectra using the transverse momentum distributions of $`\pi ^{}`$ and antiprotons from PHENIX . If there is indeed an universal function, the observed $`m_t`$-scaling should also hold for different centrality bins. In addition the centrality dependence is controlled solely by the transverse area $`\sigma `$ and the saturation momentum $`p_s`$. To test the universality, the distribution for different centralities is rescaled by

$$\frac{1}{\sigma }\frac{dN_h}{dyd^2m_t}\frac{1}{\lambda }\frac{1}{\sigma }\frac{dN_h}{dyd^2m_t}\text{ and }m_t\frac{m_t}{\lambda ^{}}.$$
(3)
The data points are then on top of each for all centrality bins (see  for a plot). Again, we make a linear plot by dividing the rescaled data by a power law. We take as reference curve a fit to the most central bin with the parameters $`p_s=1.65`$ GeV and $`n=11.8`$. The data relative to this power law is plotted in Figure 2. As one sees, scaling works within 30% up to say $`m_t2`$ GeV for all centrality classes.

The extracted scaling parameters $`p_s`$ and $`\sigma `$ are now expected to scale like $`N_{\mathrm{part}}^{1/6}`$ and $`N_{\mathrm{part}}^{2/3}`$, respectively. The saturation scale $`p_s`$ turns out to change as a function of centrality like $`p_s^2/p_{s,c}^2=0.61+0.39(N_{\mathrm{part}}/347)^{1/3}`$. This implies that there is an additional constant scale involved which is already present in pp or p$`\overline{\mathrm{p}}`$ collisions and can be associated with the Hagedorn temperature or equivalently the constant and finite mean transverse momentum seen at the ISR of about 300 MeV. The additional vacuum scale is compatible with the mean $`p_t`$ measured by the UA1 collaboration for p$`\overline{\mathrm{p}}`$ collisions at $`\sqrt{s}=200`$ GeV  of $`p_t=392\pm 3`$ MeV and the one by the STAR collaboration of $`p_t=508\pm 12`$ MeV for central gold-gold collisions at $`\sqrt{s}=130`$ GeV . The other scaling parameter, $`\sigma /\alpha _s(p_s)`$, follows the expected pattern, i.e. $`\sigma /\alpha _s(p_s)N_{\mathrm{part}}^{2/3}\mathrm{ln}\left(p_s^2/\mathrm{\Lambda }_{\mathrm{QCD}}^2\right)`$ .

Finally, we discuss $`m_t`$-scaling for p$`\overline{\mathrm{p}}`$-collisions. Identified hadron spectra were measured for energies of $`\sqrt{s}=540`$ and $`\sqrt{s}=630`$ GeV by the UA1, UA2 and UA5 collaborations  over a large range of transverse momenta. The transverse mass spectra of various hadrons is depicted in Figure 3. Without any change in the normalization, all curves are close to each other and have similar shapes. The proton data (black squares) are on top of the charged pion data, i.e. the $`\kappa _h`$ is one for the case of p$`\overline{\mathrm{p}}`$ collisions as the net baryon number is zero. Kaons are suppressed compared to pions and protons, similarly as seen in the RHIC data. The strange $`\mathrm{\Lambda }+\overline{\mathrm{\Lambda }}`$ data points are on the pion and proton data, but note that there is a factor two difference compared to the proton data, so that $`\mathrm{\Lambda }`$’s are suppressed by 1/2 like the kaons. We expect then a similar pattern for the hyperon data in heavy-ion collisions in RHIC, if corrected for the baryon number. So, $`m_t`$-scaling seems to work also for p$`\overline{\mathrm{p}}`$ collisions up to several GeV in transverse mass. The $`\kappa _h`$ factors, i.e. corrections to scaling due to strangeness and baryon number, follow the expected patterns. One notes, however, in Figure 3 that the $`\mathrm{\Lambda }`$ data start to have a slightly different slope above $`m_t2.5`$ GeV than the pions, so that it approaches the kaon data. This does not seem to be the case for the other mesons, $`\eta `$’s and kaons, though.

In summary, we showed that the transverse mass spectra of hadrons at RHIC follows one universal function of $`m_t`$. Radial flow is not necessary to describe the data. The same universal behavior seems to present in proton-antiproton collider data. The spectra at different centralities can be rescaled into each other up to $`m_t2`$ GeV. The scaling parameters follow the pattern expected from the formation of a gluon-saturated state (color-glass condensate) if an additional vacuum scale is introduced.

We thank Barbara Jacak, Julia Velkovska, and Nu Xu for many helpful discussions. JSB thanks RIKEN BNL Research Center and the Nuclear Theory Group at BNL for their kind hospitality. This manuscript has been authorized with the U.S. Department of Energy under Contracts No. DE-AC02-98CH10886 and No. DE-FG-02-93ER-40764.