text:
# Atomic lines in infrared spectra for ultracool dwarfs

## 1 Introduction

Ultracool dwarfs extend from the cool M dwarfs into the brown dwarf regime and include the spectral types L and T. It is likely that most of the ultracool dwarfs are not massive enough to undergo nuclear fusion in their cores to burn hydrogen and are thus brown dwarfs rather than stars. Here we are concerned in general terms with the spectroscopy of cool objects and use the broad temperature classification “ultracool dwarf” in preference to more specific terms with a physical meaning such as M dwarf and brown dwarf (Jones & Steele 2001).

The astrophysical importance of ultracool dwarf stars derives from the large-scale cosmological significance of their mass density, to their extremely slow chemical evolution and to their formation process and their distinction from planets. Ultracool dwarfs span the mass range from the coolest stars through two orders of magnitude in mass to giant planets. By determining their fundamental parameters, such as effective temperature, metallicity and surface gravity, their position in the Hertzsprung-Russell diagram can be determined. So far these parameters, in particular surface gravity and metallicity, are poorly determined. Yet to derive a reliable mass function across the brown dwarf regime and to discern the abundance patterns of low-metallicity brown dwarfs, it is essential to design much finer tools of analysis. Changes in effective temperatures, surface gravity and metallicity can lead to similar spectral changes. To disentangle these effects requires observations of several spectral diagnostics that respond differently to changes in these parameters (e.g., Reid & Hawley 2000).

Colour information has proved crucial in the identification of most ultracool dwarfs. However, their complex energy distributions mean that colours do not change monotonically with changing luminosity and temperature. This makes differentiation of second order effects such as metallicity and gravity very difficult to discern. Spectroscopic observations are crucial and low-mass objects have been successfully classified both in the optical (Kirkpatrick et al. 2000; Basri et al. 2000) and infrared (McLean et al. 2001; Geballe et al. 2001; Leggett et al. 2001) regimes across the M and L spectral classes. However by the T spectral class there is insufficient flux at short wavelengths to use optical wavelengths. Furthermore even by the early L spectral class the energy distribution at optical wavelengths is almost completely dominated by broad alkali resonance lines and dust opacities (Pavlenko 1998; Pavlenko et al. 2000; Burrows et al. 2000; Tsuji 2002). The infrared region is also difficult to model with strong molecular opacities of water, iron hydride and at cooler temperatures methane as well. However, there are also a number of atomic lines observable in the infrared spectra of ultracool dwarfs (e.g. McLean et al. 2003). These infrared atomic lines lie across the peak in energy distribution for these stars and are thus comparatively easier to observe than lines of the same atoms in the visible.

In recent years there has been a focus on improvements to the atomic data in the optical region, e.g., NIST (Kelleher et al. 1999), Opacity Project (Seaton et al. 1992). However, currently the accuracy and completeness of the atomic data, in particular oscillator strengths, in the IR region are far poorer. In many cases either the oscillator strengths are highly inaccurate or there are simply no measurements available. This becomes evident when synthetic spectra are generated to match infrared spectra, e.g., Jones et al. 1996, Jones & Viti 2000. Here we identify suitable sensitive infrared atomic lines that may one day be used for reliable temperature, gravity and abundance analyses for a wide range of ultracool dwarfs, from M dwarf stars through brown dwarfs as well as extra-solar giant planets. Our paper is the first step in the development of a database of infrared atomic lines with reliable oscillator strengths suitable for the study of low temperature astrophysical objects. Our specific focus is the interpretation of ultracool dwarfs covering the temperature space from the M dwarfs through the brown dwarfs to extra-solar giant planets.

Our focus is the prioritisation of atomic lines whose oscillator strengths we are either already measuring or planning to measure. Section 2 describes the observational material, atomic line lists and theoretical spectra necessary for this study. Our procedures for identifying important infrared atomic transitions are discussed in Section 3. In Section 4 we discuss our results and conclusions. In Section 5 we outline our plans for experimental measurements of the atomic data.

## 2 Input datasets

### 2.1 Observed spectra of Arcuturus and Sunspot umbra

For comparisons and the robust identification of atomic lines we used both electronic and hard copies of the “Infrared Atlas of the Arcturus spectrum, 0.9 – 5.3 $`\mu `$m ” (hereafter “Arcturus atlas”) (Hinkle et al. 1995) and “An Atlas of a Dark Sunspot Umbral spectrum from 1970 to 8640 cm-1 (1.16 to 5.1 $`\mu `$m )” (hereafter “Sunspot umbra atlas”) (Wallace & Livingston 1992). We used the electronic versions of both atlases ( ftp://ftp.noao.edu/catalogs/arcturusatlas/ir/ and ftp://ftp.noao.edu/fts/spot1atl/ ) to compare with our synthetic spectra.

Arcturus ($`\alpha `$ Boo) is a K1 III star of effective temperature Teff =4300$`\pm `$30 K and log g=1.5$`\pm `$0.15 (Peterson et al. 1993). These parameters are far from adequate to describe ultracool dwarf atmospheres ( Teff $``$2600 K, log g$``$4.5 – 5.0 ) though the Arcturus atlas provides careful identification of atomic lines and gives an opportunity to assess the behaviour of atomic lines when decreasing effective temperature from 4000 K to 2000 K. The atlas contains winter and summer observations. Since only the telluric spectrum is affected by seasons, in this paper we only use the telluric free summer data. The resolution of this atlas is of the order of R$``$100000. The presence of a magnetic field in Arcturus is indirectly confirmed by strong emission lines in the ultraviolet spectral region formed in chromosphere of Arcturus (Ayres et al. 1986). However the average value of the magnetic field is, probably, similar to the solar one, since atomic line splitting is not seen in the observed sprectrum.

The Sunspot umbra atlas has a spectral type of M2 – M5 V, which corresponds to an effective temperature range of around Teff =3170 – 3520 K (Wallace & Livingston 1992; Viti et al. 1998). The magnetic field strength is estimated as $``$ 3360 Gauss. Due to the position of the umbra near the centre of the solar disk we see only $`\sigma `$-components of Zeeman split lines of large Lande factor $`g`$ in the Sunspot umbra atlas. Due to the strong magnetic field in the spot most of the atomic lines in the Sunspot umbral spectrum have doublet-like profiles. Molecular lines, with few exceptions, are magnetically insensitive and show single profiles. The Sunspot umbra atlas provides a spectrum corrected for atmospheric absorption. Its resolution is R$``$200000 between 10000 and 25000 Å.

### 2.2 Synthetic spectra

Computations of synthetic spectra were carried out using the WITA6 program (Pavlenko 2000) assuming Local Thermodynamic Equilibrium.

For the computation of synthetic specta we used the NextGen model atmospheres structures (Hauschildt, private communication, Hauschildt et al. 1999) for cool spectra, and model atmosphere structures of Kurucz (Kurucz 1994) for the higher temperature synthetic spectra of Arcturus.

The atomic linelist used for our spectral modeling and line identification was taken from the VALD database (Kupka et al. 1999). VALD is a compilation of several different lists of atomic line data which were obtained from experimental measurements and theoretical calculations by various authors.

Chemical equilibrium were computed for the mix of around 100 molecular species. Input molecular data and continuum opacity sources are described elsewhere (Pavlenko et al. 1995; Pavlenko 1997; Pavlenko 2000). The profiles of absorption lines are described by the Voigt function $`H(a,v)`$. The formulae of Unsold (1955) were used for calculation of damping constants. For these model atmospheres we adopt the conventional value of microturbulent velocity $`v_t`$ = 2 km/s. Theoretical spectra were computed with a wavelength step of 0.05 Å and then convolved with gaussians to match the instrumental broadening.

We should note, that we didn’t consider dusty effects in our synthetic spectra calculations, since the main purpose of the present work was not to model the real spectra of ultracool objects, but to reveal the prominent atomic features to be observed in such objects. The effects of dust will modify atomic line strengths in ultracool dwarf spectra. But dust will not affect the quantum- mechanical parameters which we intend to measure at the present time.

## 3 Procedure

Our procedure of line identification and selection is based on the comparison of synthetic spectra with the observed spectra of Arcturus and Sunspot umbra from 10000 to 25000 Å.

For ultracool dwarfs we used only synthetic spectra. The choice of model atmosphere parameters was aimed at obtaining a reliable atomic line list for ultracool dwarfs in the IR region. A model atmosphere of Teff /log g/\[M/H\] = 2000/5.0/0.0 was selected (hereafter “ultracool dwarf model”). This corresponds to an ultracool dwarf spectral class of around L0 (Kirkpatrick et al. 1999, Martin et al. 1999). The half width of the gaussian we used to account for instrumental broadening in the “ultracool dwarf” synthetic spectrum is 0.4Å which corresponds to resolution R$``$20000 – 40000, approximately the resolution available to modern infrared echelle spectrometers operating between 10000 – 25000Å. Another model atmosphere 4000/4.5/0.0 (hereafter “Arcturus-like model”) was chosen to select some observable spectral atomic features formed in hotter atmospheres. The half width of the gaussian for this synthetic spectrum is 0.13Å to match the line broadening by macroturbulent velocity $`v_{\mathrm{macro}}`$=3.5 km/s in the Arcturus spectrum (Peterson et al. 1993). The Teff = 4000 K is slightly lower than Arcturus Teff =4300 K (Peterson et al. 1993), and higher than the temperature of the Sunspot umbra Teff $``$3170 – 3520 K (Wallace & Livingston 1992).

In the electronic versions of the Arcturus and Sunspot umbra atlases the data are given on a wavenumber scale. We transform the wavenumbers to wavelengths in air using the formula from Allen’s Astrophysical Quantaties (2000).

It is well known that the difficulties of line identification in IR spectra of ultracool dwarfs have a number of causes. There are many strong molecular lines of H2O (Jones et al. 2002), FeH (Cushing et al. 2003), CO (Pavlenko & Jones 2003), CH4 (Noll et al. 2000) native to ultracool dwarfs as well as telluric molecular lines, particulary strong between the J, H and K bands. Nevertheless, we were able to identify many absorption features in the atmospheres of Arcturus and Sunspot umbra spectra likely to be present in ultracool dwarfs. To do this, we developed a procedure based on the inter-comparison of computed and observed spectra designed to minimise possible errors in line identification. The line identification files were made using
a) “ultracool dwarf” synthetic spectrum (2000/5.0/0.0),
b) “Arcturus-like” synthetic spectrum (4000/4.5/0.0).
In addition we used line identification files from hard copy Arcturus and Sunspot umbra atlases.

In our identification procedure we assume that lines absorb only at a given wavelength. This assumption provides an upper limit of absorption. In the real spectra the residual intensity may be lower due to contribution of other lines. For every line from the VALD linelist in the range of 10000 – 25000Å the residual flux $`r_\nu =\frac{(\mathrm{Flux}\mathrm{in}\mathrm{line})_\nu }{(\mathrm{Flux}\mathrm{in}\mathrm{continuum})_\nu }`$ was computed. Atomic lines from VALD with central intensities deeper than 0.8 of the residual flux in the spectrum of 2000/5.0/0.0 atmosphere were selected to restrict the number of lines to those likely to be observable in cool spectra. These lines are given in the Table 1. The set of the selected lines is rather large and thus we assign different levels of priority to the lines. The main criteria for selecting lines of the highest priority is that the central intensity of the identified line ought to be greater than 0.6 of the residual flux (column 4 in the Table 1) in the “ultracool dwarf” synthetic spectrum.

Plots for all spectral regions containing atomic lines of interest are shown in Figure 1 ( they are also available at http:// www.astro.livjm.ac.uk/ hraj/ spectralatlas/ index.html ). The identifications of atomic features are labelled in the spectra. Arrows and labels at the top of the plots show identification in the “ultracool dwarf” synthetic spectrum, the same labels in the bottom of the plots are for the ”Arcturus-like” spectrum. The Sunspot umbral spectrum is shown by bold line, the Arcturus spectrum by a dotted line, the “ultracool dwarf” synthetic spectrum as a solid line and “Arcturus-like” synthetic spectrum as a dashed line. In order to simplify our plots molecular line identifications are not included. Usually molecular lines have narrow deep profiles allowing them to be distinguished from atomic lines. Their identification in the Sunspot umbra and Arcturus can be found in the hard copy atlases.

## 4 Results

Our results are shown in two tables. Table 1 contains all identified lines with central intensities deeper than 0.8 of residual flux using “ultracool dwarf model” identifications. Table 2 investigates the sensitivity of priority 1 lines identified in Table 1 to tempereture, gravity and metallicity.

In the first column of Table 1 we estimate the priority of lines for measurement of atomic data.
a) Priority 1 indicates that atomic line is deeper than 0.6 in the 2000/5.0/0.0 synthetic spectrum. Also we include as first priority three lines of Ti (22211.229, 22232.838, 22274.012Å) which are shallower than 0.6. They are located far from strong atomic and molecular lines and thus are good for identification purposes in the spectra of ultracool dwarfs.
b) Priority 2 indicates that the residual flux of an atomic line is between 0.8 – 0.6.
c) Priority 3 indicates that the atomic line is located in the wing of a nearby stronger line or in a wavelength region with strong telluric absorption.

Atomic lines without an assigned priority might be relevant in the analysis of ultracool dwarf spectra, but they are predicted to be weak in the synthetic spectrum.

The second column shows the wavelength of the line in Å in air . These values were taken from VALD. We check that wavelengths from Arcturus and Sunspot umbra atlases identification files are in a good agreement with VALD. The wavelengths of a few lines differ from VALD values by less than 0.1Å in the region 10000 – 20000 Å and by less than 0.3Å beyond 20000 Å.

The third column shows the atomic identifications of the line.

Columns 4 – 7 show the linedepths of identified lines. Column 4 – the residual flux in the line obtained from the “ultracool dwarf” synthetic spectrum; column 5 – from the ”Arcturus-like” synthetic spectrum; column 6 – measured from the hard copy of Arcturus atlas depth of identified atomic lines; column 7 – the same as column 6 for the hard copy of the Sunspot umbra atlas. As was mentioned above we used the hard copy versions of the Arcturus and Sunspot umbra atlases to measure linedepths. The accuracy of our line depth measurements from the hardcopy atlases is $``$0.01 – 0.02 for the Arcturus atlas (column 6) and $``$0.02 – 0.03 for the Sunspot umbra atlas (column 7). We note that linedepth information given in columns 6 and 7 is only used to check the presence of a line in the observed spectra. Thus we do not need high accuracy linedepth measurements for the observable spectra of the Sunspot umbra and Arcturus.

We note that the data in the Sunspot umbra atlas provide a spectrum from 11600Å, rather than 10000Å where our figures and tables start from. Some lines are labelled (see column 7) by the symbol ”$``$” which means that part of the line in the hard copy of the atlas is obscured by a gap in the data and the value given in the table is for the visible part of the line. In columns 6 and 7 there are also some lines with ”bl” suffices. This means that these lines are blended with stronger nearby lines.

There are two regions with strong telluric absorption ( marked in columns 6 and 7 by $``$): 13600–14770Å and 18260–19450Å , regions between the J and H, and H and K photometric bands. Lines identified from VALD in these regions are difficult to observe in ground based observations but can be used in probable future space observations of ultracool dwarfs.

Finally we note, that in Table 1 we italicise atomic lines of Rb, Y, Ba and Lu, identified only using theoretical computations. Although these lines are not seen in the observations of Arcturus and Sunspot umbra, they are strong in synthetic spectrum computed for the 2000/5.0/0.0 model atmosphere, and therefore also likely to be strong below 2000 K.

In Table 2 we show only lines which were marked in Table 1 as first priority lines. The first column are wavelengths in Å; the second column are names of elements. In the third column we give the central line intensity for the “ultracool dwarf” (2000/5.0/0.0) spectrum. We compute some synthetic spectra for ultracool objects varying parameters of Teff , log g and metallicity to show the sensitivity of the central intensities of identified lines to these parameters. For our computations we chose the minimal step in model atmospheres grid: $`\mathrm{\Delta }`$Teff = 100 K, $`\mathrm{\Delta }`$log g = 0.5, $`\mathrm{\Delta }`$\[M/H\]=-0.5. In columns 4 – 6 we show the difference of central intensity ( in % ) computations with “new” models from central intensity computations for our “ultracool dwarf model”. Negative values in columns 4 – 6 indicate that a line is stronger in the new model relative to the reference model 2000/5.0/0.0.

Column 4 shows the dependence of residual fluxes on effective temperature. For comparison a 2100/5.0/0.0 synthetic spectrum was used. One can see that only a few lines show strong temperature dependence for $`\mathrm{\Delta }`$Teff =100 K. The NextGen 2000/4.5/0.0 model was used with 2000/5.0/0.0 to test gravity sensitivity. The dependence for $`\mathrm{\Delta }`$(log g)=0.5 is also relatively weak for most lines and is shown in column 5. $`\mathrm{\Delta }`$\[M/H\]=0.5 from computations using NextGen model atmospheres 2000/5.0/-0.5 and 2000/5.0/0.0 is shown in column 6. The metallicity dependence is relatively high for all lines. Thus Table 2 indicates that relatively small variations of model atmosphere parameters will not seriously affect our line selection and prioritisation.

## 5 Experimental work

The identification and prioritisation of atomic lines presented in this paper is only the first step in our project and is primarily based on a theoretical treatment. Our program of measurements of atomic data for the priority lines given in this paper is underway. We are currently measuring oscillator strengths by measurements of sets of relative line intensities using the high resolution Fourier transform spectrometers at Imperial College (Pickering 2002) and NIST (National Institute of Standards and Technology, US). We have already recorded spectra of Ti, Mn and Na, and measurements of Mg, K, Ca and Fe are in progress. Measurements of other species listed in Table 1 are underway in the near future. These line intensities will be used to obtain branching ratios which are then combined with level lifetimes to obtain f-values. We expect to achieve uncertainties in f-values of around 10-15 %.

## 6 Acknowledgements

We would like to thank the authors of the Arcturus and Sunspot umbra atlases for making their data available through an ftp site and Prof. Hauschildt and the PHOENIX team for NextGen model atmospheres. NSO/Kitt Peak FTS data used here were produced by NSF/NOAO. SV thanks the Italian Space Agency (ISA) for financial support. RBW thanks Gillian Nave at NIST for the assistance and facilities provided for the IR experimental measurements. All the authors thank the Royal Society and PPARC of the UK for support of travel and experimental work.


format: text

text:
# The study of 𝛿 Sct stars in the transition era from ground–based to space photometry

## 1. Introduction

The different research teams involved in the study of $`\delta `$ Sct stars have slightly changed their strategy in the past years. The observational effort to secure worldwide coverage of case studies is continued, but the specifications become more severe, especially about target characterization and frequency resolution.

After the successful launch of the Canadian satellite MOST, which will be the pioneer of asteroseismology from space, the future missions are programmed to properly take into account the need of an adequate frequency resolution: COROT will spent 30 and 150 days (additional and core programs, respectively) on the target, while EDDINGTON will arrive up to 1 year. Such a requirement is a direct consequence of the observational results on $`\delta `$ Sct, $`\gamma `$ Dor, SPB stars, … obtained from ground. It should be noted that without these results (see Poretti 2000 for a review about $`\delta `$ Sct stars) the scientific background of the space missions would be much less defined and the risks of incomplete results (owing to inaccurate selection of targets, insufficient resolution, underestimate of the influence of the rotation) much higher.

## 2. The increasing number of frequencies detected in FG Vir

López de Coca et al. (1984) easily detected the main period of FG Vir on the basis of few nights. It can be considered the discovery step. After that, Mantegazza, Poretti & Bossi (1994) detected seven frequencies on the basis of a 14–nights run at the European Southern Observatory. It can be considered the second step, i.e., the limit for a single–site campaign. Owing to its equatorial position, FG Vir was a suitable candidate for multisite campaign, being possible to observe it from both hemispheres. Indeed, Breger et al. (1999) confirmed the seven frequencies and increased the total number of detected frequencies to 24, thanks to a campaign involving 6 observatories during 40 days. Again, this can be considered the limit for a multisite campaign which can be organized with a moderate effort. The next step was to set up a campaign not only involving more than 6 sites, bust also spanning a longer time baseline. Moreover, the recent results obtained by Mantegazza & Poretti (2002) strongly indicated how spectroscopy can supply useful hints about the inclination of the rotational axis and to propose some mode identifications. Therefore, the 2002–03 campaign on FG Vir has been performed by combining spectroscopy and photometry over a 2–years time baseline. The aim is to detect close doublets and investigate in a more powerful way the eventual amplitude variations. Such a long–term campaign constitutes the new fronteer of the ground–based observations.

Figure 1 shows how the number of detected frequencies increased by refining the observational approach. ¿From Fig. 1 it can be extrapolated that the lowering of the threshold amplitude by a factor 10-1 (which should be easily attainable from space) will allow us to detect hundreds of excited modes.

## 3. The preparation of the COROT space mission

The COROT satellite will monitor targets which must be selected in two fields centered at $`\alpha =6^h50^m`$, $`\delta =0`$ (Anticenter direction), and $`\alpha =18^h50^m`$, $`\delta =0`$ (Center direction). The prototype of the class, $`\delta `$ Sct itself, is included in the Center direction field, but it is considered too evolved to match the goals of the mission, as asteroseismic inferences for these stars are greatly complicated by the richness of excited modes in a narrow frequency range (275 possible modes in a 4 cd-1 interval; Templeton et al. 1997). Therefore, we searched for new $`\delta `$ Sct stars located in the COROT field–of–view which are unevolved stars (high–priority for COROT) or little evolved ones (low–priority stars). Efforts have been concentrated especially in the galactic Center direction, as there are no many solar–like candidates in such a direction and target allocation is easier.

As a result, we actually scanned the solar neighbourhood; $`\delta `$ Sct stars are commonly found and the incidence of variability is around 23% (13 variables out of 57 candidates) for stars located in the lower part of the instability strip, i.e., the region between the ZAMS and the blue hooks of the evolutionary tracks (Poretti et al. 2003). With such a percentage, the possibility to find $`\delta `$ Sct stars in a limited sky region is reasonable and they can provide a suitable basis of asteroseismic targets for any space mission. Figure 2 shows the specific result obtained for COROT. The fine tuning of the candidates is now possible and the ground–based preparatory work can greatly take advantage from the expertise piled–up in the past decades on the case studies. It is suitable that the $`\delta `$ Sct stars which will be selected as COROT targets will became case studies for spectroscopy and/or preliminary photometry.

## 4. Colour information from space missions

When performing photometry in different colours it is possible to calculate solutions for each passband. The phase differences and the amplitude ratios will be then used to discriminate the different spherical harmonics of a non-radially pulsating star, therefore helping to identify the oscillation mode (Garrido 2000).

The explanetary camera of the COROT space mission consists on a CCD with a bi-prism in front of it giving coloured spots on the detector (see Vauclair in these proceedings for details). The three colours basically supplied by the camera are sufficient to distinguish the different spherical degrees $`\mathrm{}`$ of the corresponding radial or non-radial excited periods, making therefore possible the mode identification of several types of pulsating stars, including $`\delta `$ Scuti stars. Physical information derived from these coloured data also includes limb-darkening coefficients, non-adiabatic phase lags (see Garrido 2003 for details) and, for some specific cases, clues concerning stellar rotation (see Daszynska-Daszkiewicz et al. 2002 for details).

EDDINGTON is an ESA space mission devoted to perform asteroseismology and planet detection by the method of transits. The present status of the mission consists of three telescopes: one will receive white light, while the other are equipped with two coloured filters blocking only a small (5 to 10 %) fraction of the light from the red (Ultrawide Blue Filter) or from the blue (Ultrawide Red Filter). Differences between the white light and these two coloured filters give narrow enough photometric bands to be useful for mode identification, as shown in detail in Deeg & Garrido (2003). In this way no many photons are lost, as required by the original mission specifications. Obviously the system is valid only for stars showing relatively high S/N photometric amplitudes (see Garrido 2003 for details). $`\delta `$ Scuti will have sufficient S/N ratios as to easily derive phases and amplitudes for the different colours.

## References

Breger, M., Handler, G., Garrido, R., et al. 1999, A&A, 349, 225

López de Coca, P., Garrido, R., Costa, V., & Rolland, A. 1984, IBVS 2465

Daszynska-Daszkiewicz, J., Dziembowski, W.A. 2002 & Goupil, M.J. A&A, 392, 151.

Deeg, H., & Garrido, R., 2003, Eddington Internal Report.

Garrido, R., 2000, in “Delta Scuti and Related Stars”, Eds. M. Breger and M.H. Montgomery, ASP Conf. Ser., vol. 210, p. 67

Garrido, R., 2003, 2nd Eddington workshop, Palermo.

Mantegazza, L., & Poretti, E., 2002, A&A, 396, 911

Mantegazza, L., Poretti, E., & Bossi, M. 1994, A&A, 287, 95

Poretti, E., 2000, in “Delta Scuti and Related Stars”, Eds. M. Breger and M.H. Montgomery, ASP Conf. Ser., vol. 210, p. 45

Poretti, E., Garrido, R., Amado, P.J., et al., 2003, A&A, 406, 203

Templeton, M., McNamara, B.J., & Guzik, J.A., 2003, AJ, 114, 1592

Zima W., Heiter U., Cottrell P.L., Lehmann H., Mathias P., Poretti E., & Breger M., 2003, in “Asteroseismology Across the HR Diagram”, Eds. M.J. Thompson, M.S. Cunha, M. Monteiro, Kluwer Academic Publ., P489


text:
# Why Massive Black Holes are Small in Disk Galaxies ?

## 1 Introduction

Supermassive black holes (SMBHs) are now thought to reside in virtually all galactic bulges. The recent compilation of the kinematical data of galactic centers in both active and inactive galaxies has shown that the mass of SMBH does correlate with the mass of a galactic bulge; the mass ratio of BH to bulge is $`0.001`$ as a median value (Kormendy & Richstone 1995; Richstone et al. 1998; Magorrian et al. 1998; Gebhardt et al. 2000; Ferrarese & Merritt 2000; Merritt & Ferrarese 2001; McLure & Dunlop 2001, 2002; Wandel 2002). On the other hand, it appears that the SMBH does not correlate with a disk component (Kormendy & Gebhardt 2001). In disk galaxies, the mass ratio is significantly smaller than 0.001 if disk stars are included (Salucci et al.2000; Sarzi et al. 2001, Ferrarese 2002, Baes et al. 2003). These observations imply that the formation of the SMBH is related to the bulge components of galaxies.

Recently, as a possible mechanism for the formation of SMBHs, the radiation-hydrodynamic model is proposed, where the radiation drag (Poynting-Robertson effect) due to stellar radiation in a galactic bulge drives the mass accretion to grow a massive BH (Umemura 2001). This is a parallel model of the formation of massive BHs by the Compton drag in an early universe (Umemura, Loeb, & Turner 1997). In an optically thick regime, the efficiency of radiation drag is saturated due to the conservation of the photon number (Tsuribe & Umemura 1997). The timescale of radiation drag-induced mass accretion is $`t_{\mathrm{drag}}=8.6\times 10^7R_{\mathrm{kpc}}^2(L_{}/10^{12}L_{})^1(Z/Z_{})^1\mathrm{yr}`$, where $`R_{\mathrm{kpc}}`$ is the galaxy size in units of kpc, and $`Z`$ is the metallicity of gas. The mass accretion rate is given by $`\dot{M}=M_\mathrm{g}d\mathrm{ln}J/dtL_{}/c^2`$, where $`J`$, $`L_{}`$, and $`M_\mathrm{g}`$ are the total angular momentum of gaseous component, the total luminosity of the system, and the total mass of gas, respectively. Thus, the final BH mass is estimated by $`M_{\mathrm{BH}}L_{}(t)/c^2𝑑t.`$ Hence, the final BH mass is proportional to the total radiation energy from stars, and the theoretical upper limit of BH-to-bulge mass ratio is determined by the energy conversion efficiency of nuclear fusion from hydrogen to helium, i.e., 0.007 (Umemura 2001). In practice, the inhomogeneity of interstellar medium helps the radiation drag to sustain high efficiency (Kawakatu & Umemura 2002). By incorporating the realistic chemical evolution, the radiation drag model predicts $`M_{\mathrm{BH}}/M_{\mathrm{bulge}}0.001`$ (Kawakatu, Umemura, & Mori 2003). However, the radiation drag efficiency could be strongly subject to the effect of geometry (Umemura et al. 1997,1998; Ohsuga et al. 1999a). If the system is spherical, the emitted photons are effectively consumed within the system, whereas a large fraction of photons can escape from disk-like systems, and thus the drag efficiency is likely to be reduced considerably. It is suggested that such a geometrical effect may be the reason why the BH is smaller in disk galaxies (Umemura 2001; Umemura 2003), although the details have not been clear quantitatively.

In this Letter, we scrutinize the efficiency of radiation drag in a system composed of bulge and disk. For the purpose, we accurately solve the 3D radiative transfer in an inhomogeneous interstellar medium. The main goal is to elucidate the relationship between the final BH mass formed by the radiation drag and the bulge-to-disk ratio of host galaxy.

## 2 Model

We consider a system composed of bulge and disk, which consist of dark matter, stars, and dusty interstellar matter (ISM). The total baryonic mass of the system, $`M_\mathrm{b}`$, is set to be $`10^{11}M_{}`$. But, it is noted that the results are scaled with $`M_\mathrm{b}`$. The dark matter is distributed following the profile by Navarro, Frenk, $`\&`$ White (1997), $`\rho (r)(r/r_\mathrm{s})^1(1+r/r_\mathrm{s})^2`$, where $`r_\mathrm{s}`$ is the characteristic radius, which is $`r_\mathrm{s}=50\mathrm{k}\mathrm{p}\mathrm{c}`$ here. The mass of the dark matter component, $`M_{\mathrm{DM}}`$, is equal to $`M_\mathrm{b}`$ within the galaxy size of $`r_\mathrm{g}=`$10kpc.

We incorporate the evolution of galaxy from the formation stage, where a galactic system starts with a protogalactic cloud purely composed of barynic gas. By invoking the star formation in Schmidt’s law, stellar components are generated in this system. To include the recycling of ISM, we employ an evolutionary spectral synthesis code, PEGASE (Fioc & Rocca-Volmerange 1997), with assuming a Salpeter-type initial mass function (IMF) as $`\varphi =dn/d\mathrm{log}m_{}=\varphi _0(m_{}/M_{})^{1.35}`$ for the mass range of $`[0.1M_{},60M_{}]`$, where $`m_{}`$ is the stellar mass. The star formation rate (SFR) per unit mass at time $`t`$, $`C(t)`$, is assumed to be proportional to the gas mass fraction, $`C(t)=kf_\mathrm{g}`$ with $`f_\mathrm{g}M_\mathrm{g}(t)/M_\mathrm{b}`$, where $`k`$ is a rate coefficient which is $`8.6\mathrm{Gyr}^1`$ for the bulge component and $`0.86\mathrm{Gyr}^1`$ for the disk component, respectively.

The formed stars in the bulge are distributed following Hernquist’s profile (1990), which mimics de Vaucouleur’s profile: $`\rho (r)(r/r_{\mathrm{core}})^1(1+r/r_{\mathrm{core}})^3`$ , where $`r_{\mathrm{core}}`$ is the core radius. The bulge radius, $`r_{\mathrm{bulge}}`$, is given by a mass-to-size relation as $`r_{\mathrm{bulge}}/r_\mathrm{g}=(M_{\mathrm{bulge}}/M_{\mathrm{galaxy}})^{1/2}`$ (Kormendy 1985; Chiosi & Carraro 2002), where $`M_{\mathrm{bulge}}`$ is the bulge mass and $`M_{\mathrm{galaxy}}`$ is the final stellar mass of the galaxy, and $`r_{\mathrm{core}}`$ is set to $`r_{\mathrm{core}}=0.1r_{\mathrm{bulge}}`$. In the present analysis, the bulge mass fraction ($`M_{\mathrm{bulge}}/M_{\mathrm{galaxy}}`$) is altered as a key parameter. In addition, angular momentum is smeared, assuming rigid rotation, according to the spin parameter of $`\lambda =0.05`$ (Barnes & Efstathiou 1987; Heavens & Peacock 1988).

As for the disk component, formed stars are distributed exponentially in the radial direction, with assuming the vertical scale-height which is changed from $`0.01r_\mathrm{g}`$ to $`0.1r_\mathrm{g}`$. The disk radius, $`r_{\mathrm{disk}}`$, is given by a mass-to-size relation of $`r_{\mathrm{disk}}/r_\mathrm{g}=(M_{\mathrm{disk}}/M_{\mathrm{galaxy}})^{1/2}`$, where $`M_{\mathrm{disk}}`$ is the disk mass. The disk component is assumed to be in rotation balance.

As for ISM, we consider clumpy matter, since the ISM is observed to be highly inhomogeneous in active star-forming galaxies (Sanders et al. 1988; Gordon, Calzetti, & Witt 1997). Here, $`N_\mathrm{c}(=10^4)`$ identical clouds are distributed randomly in a system of the bulge and disk. The number of clouds is minimal enough to resolve the disk scal-height of $`0.01r_\mathrm{g}`$. (It is noted that simulations with a three times larger number of clouds did not lead to any fundamental difference.) The distribution and motion of the ISM is the same as those of stellar components. We assume that the cloud covering factor is unity according to the analysis by Kawakatu & Umemura (2002). The internal density in a cloud is assumed to be uniform, and the opacity is determined by a dust-to-gas ratio which is calculated by PEGASE. The optical depth of a cloud and therefore the overall optical depth of galaxy, $`\tau _\mathrm{T}`$, depend on the cloud size $`r_\mathrm{c}`$, where $`r_\mathrm{c}=100`$pc is assumed as a fiducial case. But we have confirmed that no essential difference in the final BH mass is found by changing $`r_\mathrm{c}`$ so as to enhance $`\tau _\mathrm{T}`$ by an order.

The radiation drag, which drives the mass accretion, originates in the relativistic effect in absorption and reemission of radiation. This effect is involved in relativistic radiation hydrodynamics (RHD) (Umemura, Fukue, & Mineshige 1997, 1998; Fukue, Umemura, & Mineshige 1997). The angular momentum transfer in RHD is given by the azimuthal equation of motion in cylindrical coordinates,

$$\frac{1}{r}\frac{d(rv_\varphi )}{dt}=\frac{\chi }{c}[F^\varphi (E+P^{\varphi \varphi })v_\varphi ],$$
(1)
where $`\chi `$ is the mass extinction coefficient, $`E`$ is the radiation energy density, $`F^\varphi `$ is the radiation flux, and $`P^{\varphi \varphi }`$ is the radiation stress tensor. All radiative quantities are determined by radiation from stars diluted by dusty ISM. Here, we solve radiative transfer including the dust extinction by ISM and obtain these radiative quantities. Then, we evaluate the total angular momentum loss rate (see Kawakatu & Umemura 2002 for the details of method). If an optically thick cloud is irradiated by the radiation, an optically-thin surface layer is stripped by the radiation drag (Tsuribe & Umemura 1997), and loses angular momentum. Since the bulge luminosity is lower than the Eddington luminosity for dusty gas, $`L_{Edd}=4\pi cGM_{\mathrm{dyn}}/\chi `$ with the dynamical mass $`M_{\mathrm{dyn}}`$ within $`r_\mathrm{g}`$, during the whole history of galaxy, the stripped gas is not blown away by the radiation pressure, but confined in the galaxy. Hence, the gas moves in a gravitational potential weakened by radiation pressure and loses angular momentum by radiation drag. The stripped gas losing angular momentum falls onto the galactic center to form a massive object, which is likely to undergo the inside-out viscous collapse (Mineshige, Tsuribe, & Umemura 1998; Tsuribe 1999) and eventually evolve into a SMBH by the general relativistic instability (Saijo et al. 2002; Shibata & Shapiro 2002). Using equation (1), we can calculate the mass accretion rate onto the central object by $`\dot{M}_\mathrm{g}/M_\mathrm{g}=\dot{J}/J`$, where $`J`$ is the total angular momentum of ISM. Then, the BH mass ($`M_{\mathrm{BH}}`$) is assessed as

$$M_{\mathrm{BH}}=_0^{t_0}\dot{M}_\mathrm{g}𝑑t=_0^{t_0}M_\mathrm{g}\frac{\dot{J}}{J}𝑑t,$$
(2)
where $`t_0`$ is the Hubble time.

## 3 Results

By changing the bulge mass fraction, $`f_{\mathrm{bulge}}M_{\mathrm{bulge}}/M_{\mathrm{galaxy}}`$, we investigate the relation between the final BH mass and the bulge-to-disk ratio of galaxy. The growth of BH depends upon the optical depth of galaxy, $`\tau _\mathrm{T}`$, but no essential difference in the final BH mass is found by enhancing $`\tau _\mathrm{T}`$ by an order, because the mass accretion rate is saturated when the system becomes optically-thick (Umemura 2001; Kawakatu & Umemura 2002). In Figure 1, the resultant BH-to-galaxy mass ratio, $`M_{\mathrm{BH}}/M_{\mathrm{galaxy}}`$, is shown against the bulge mass fraction ($`f_{\mathrm{bulge}}`$) by thick lines for three values of disk scale-height. We can see a clear trend that $`M_{\mathrm{BH}}/M_{\mathrm{galaxy}}`$ decreases with decreasing $`f_{\mathrm{bulge}}`$. This tendency comes from three effects: First, a larger number of photons escape from the disk surface, since the surface-to-volume ratio is larger for the disk component compared to the bulge. This effect is called the geometrical dilution here. Secondly, the radiation from disk stars is heavily diminished across the disk, because the edge-on optical depth becomes larger quickly. Thirdly, the velocity difference between a star and an absorbing cloud becomes closer to zero in an optically-thick disk, so that the radiation drag cannot work effectively. For instance, if $`f_{\mathrm{bulge}}=0.03`$ which corresponds to Sd galaxies, $`M_{\mathrm{BH}}/M_{\mathrm{galaxy}}`$ is reduced by a factor of $`20`$ to $`100`$, compared to $`M_{\mathrm{BH}}/M_{\mathrm{galaxy}}0.001`$ for elliptical galaxies ($`f_{\mathrm{bulge}}=1`$). Furthermore, it is predicted that a bulgeless (pure disk) system can also possess a massive black hole with the mass ratio of $`M_{\mathrm{BH}}/M_{\mathrm{galaxy}}0.001\times \left(h/2r_{\mathrm{disk}}\right)`$, where $`h`$ is the scale height of disk. This ranges from $`6\times 10^6`$ to $`5\times 10^5`$, corresponding to $`h=0.01r_{\mathrm{disk}}`$ and $`h=0.1r_{\mathrm{disk}}`$, respectively.

On the other hand, the bulge component itself is likely to be impervious to the above three geometrical effects. In Figure 1, by focusing on only the bulge component, the BH-to-bulge mass ratio, $`M_{\mathrm{BH}}/M_{\mathrm{bulge}}`$, is also shown. As seen in Figure 1, $`M_{\mathrm{BH}}/M_{\mathrm{bulge}}`$ is almost constant as expected. Hence, it is concluded that the BH-to-galaxy mass ratio in disk galaxies is substantially determined by bulge components, not by disk components: $`M_{\mathrm{BH}}/M_{\mathrm{galaxy}}0.001`$, where the constant value of 0.001 is basically determined by the energy conversion efficiency of nuclear fusion from hydrogen to helium, i.e., 0.007 (Umemura 2001). The BH accretion can trigger an active galactic nucleus (AGN) as well, and then the further mass accretion is induced via radiation drag by the AGN. The effect can enhance the final BH mass maximally by a factor of 1.7 (Umemura 2001). Also, a different level of ISM clumpiness can reduce the final BH mass by a factor of 2 (Kawakatu & Umemura 2002). As for the effect of IMF, if the slope and the mass range of IMF are changed to satisfy the spectrophotometric properties, then the final BH mass is altered by a factor of $`\pm 40\%`$ . Similar uncertainties attach also to the final BH mass in disk galaxies.

## 4 Comparison with Observations

The present results are compared with observational data for SMBHs in disk galaxies in Figure 2. The data we adopted are summarized in Table 1. In the upper panel of Figure 2, $`M_{\mathrm{BH}}/M_{\mathrm{galaxy}}`$ is shown against $`f_{\mathrm{bulge}}`$. The hatched area is the prediction shown in Figure 1. Two dotted lines are respectively the upper and lower bound of the prediction by including other theoretical uncertainties mentioned above. The observational data are categorized into four types: disk galaxies which do not possess AGNs, Seyfert 1 galaxies (Sy1s), narrow line Seyfert 1 galaxies (NLSy1s), and Seyfert 2 galaxies (Sy2s). Also, the objects accompanied with starburst events are specified. The observational data in Figure 2 shows the decrease of $`M_{\mathrm{BH}}/M_{\mathrm{galaxy}}`$ by more than an order of magnitude from $`f_{\mathrm{bulge}}=1`$ to $`f_{\mathrm{bulge}}=0.03`$. This trend is broadly consistent with the theoretical prediction. But, it is worth noting that several objects possess relatively small BHs compared to the prediction, and they are mostly Seyfert galaxies with starburst events (NGC 7469, NGC 1068, Circinus) or NLSy1s (Mrk 590, NGC4051). In recent years, an evolutionary model for AGNs beyond the AGN unified theory has been considered from both observational and theoretical points of view (Maiolino et al. 1995; Radovich, Rafanelli, & Barbon 1998; Ohsuga & Umemura 1999, 2001). In this model, AGNs with starbursts are in earlier evolutionary phase, and therefore the BH mass is likely to be smaller, compared to a post-starburst phase. Furthermore, NLSy1s are thought to be in a rapidly growing phase of an immature BH (Mathur et al. 2001). NGC1068 could be also in such a phase (Kawaguchi 2003). This point may be intriguing to analyze further with larger sample data in the future.

In the lower panel of Figure 2, $`M_{\mathrm{BH}}/M_{\mathrm{bulge}}`$ is shown against $`f_{\mathrm{bulge}}`$. When only the bulge component is focused, the observational data for $`M_{\mathrm{BH}}/M_{\mathrm{bulge}}`$ lies at a level of $`0.001`$, and they roughly agree with the prediction. But, again the ratios for Seyferts with starbursts and NLSy1s fall appreciably below $`0.001`$.

The present model also predicts a SMBH even in a pure disk system without bulge, and the mass of BH depends upon the disk scale-height. NGC4395 is an archetype of bulgeless galaxy in which a massive BH resides (Filipenko & Ho 2003). The BH fraction in this galaxy is estimated to be $`M_{\mathrm{BH}}/M_{\mathrm{galaxy}}1.5\times 10^5`$. In the present model, this level can be realized by a disk scale height of $`h0.04r_{\mathrm{disk}}`$.

We thank T.Nakamoto, T.Saitou and A.Yonehara for helpful discussions. We also thank the anonymous referee for precious comments. Numerical simulations were performed with facilities at the Center of Computational Physics, University of Tsukuba. This work was supported in part by the Grant-in-Aid of the JSPS, 11640225.


text:
# Diffusion on Complex Networks : A way to probe their large scale topological structures

## 1 Introduction

Complex networks are natural structures for representing many real-world systems. Their flexible and adaptive features make them ideal for interrelating various types of information and to update them as time goes by. Recently, there has been tremendous interest, from various fields of science, in the study of the statistical properties of such networks . By now, many interesting features of complex networks have been established. So far, the main attention in the study of complex networks has been on the properties of individual nodes and how they connect (link) to their nearest neighbors. However, exploring the local properties of the network outside the nearest neighbor level, has not been well studied. A noteworthy exception is the recent study by Girvan and Newman . Such studies will naturally have to address the cluster (modular) structure of the network. To know if a network is modular or not, is important if one tries to assess its robustness and stability, since only a few critical links are responsible for the inter-modular communication in clustered networks. Hence to make the network more robust, such critical links have to be identified and strengthen.

In this work we will address the large scale topological structures of networks. This problem will be approached by considering an auxiliary diffusion process on the underlying complex network. As it will turn out, it is the slowest decaying diffusive eigenmodes that will be of most interest to us, since such modes will contain, as we will see, information about the weakly interacting modules of the network. Hence, by studying how an ensemble of random walkers slowly reaches a state of equilibrium, one can learn something about the topology of the network.

To study spectral properties of networks is not new; Its variants has previously been applied to random graphs , to social networks (the correspondence analysis) , random and small-world networks (the Laplace equation analysis) , artificial scale-free networks , and community structures . A diffusion approach has also made it to practical applications. For instance, the analysis of a diffusion process lies at the heart of the popular search engine Google .

## 2 A Markovian diffusion process on complex networks

The original physical motivation behind the method to be outlined below, was that the relaxation of some arbitrary initial state (of walkers on the network) toward the steady state distribution, was expected to be fast in regions that are highly connected, while slow in regions that had low connectivity. By definition, a module is highly connected internally, but has only a small number of links to nodes outside the module. Hence, it was reasoned, that by identifying the slowly decaying eigenmodes of the diffusive network process, one should be able to obtain information about the large scale topological structures, like modules, of the underlying complex network. We will now formalize this idea and outline the method in some detail.

Let us start by assuming that we are dealing with a (fully connected) complex network consisting of $`N`$ nodes and characterized by some degree distribution $`p(k)`$ where $`k_i`$ denotes the connectivity of node $`i`$. Imagine now placing a large number ($`N`$) of random walkers onto this network. The fraction of walkers on node $`i`$ (relative the total number) at time $`t`$ we will denote by $`𝝆_i(t)`$. In each time step, the walkers are allowed to move randomly between nodes that are directly linked to each other. Since there is no way that a walker can vanish from the network, the total number of walkers must be conserved globally, i.e. one must have $`_i𝝆_i(t)=1`$ at all times. Furthermore, locally a continuity equation must be satisfied; If we consider, say node $`i`$, that directly links to other nodes $`j`$, one must require that (the balance equation):

$`𝝆_i(t+1)𝝆_i(t)`$ $`=`$ $`{\displaystyle \underset{j}{}}A_{ij}{\displaystyle \frac{𝝆_j(t)}{k_j}}{\displaystyle \underset{j}{}}A_{ji}{\displaystyle \frac{𝝆_i(t)}{k_i}}.`$ (1)
In writing this equation, we have introduced the so-called adjacency matrix $`A_{ij}=A_{ji}`$ defined to be $`1`$ if node $`i`$ and $`j`$ are directly linked to each other, and $`0`$ otherwise , and $`k_i`$, we recall, is the degree of node $`i`$, i.e. its number of nearest neighbors $`k_i=_jA_{ij}`$. The first term on the right hand side of Eq. (1) describes the flow of walkers into node $`i`$, while the last term is associated with the out-flow of walkers from the same node. As will be useful later, Eq. (1) can be casted into the following equivalent matrix form:

$`_t𝝆(t)`$ $`=`$ $`𝐃𝝆(t),`$ (2)
where $`_t𝝆(t)=𝝆(t+1)𝝆(t)`$, and $`𝐃`$ is the diffusion matrix to be defined below. Eq. (2) should be compared to the continuous diffusion equation for the particle density: $`_t𝝆(𝐫,t)=D\mathbf{}^2𝝆(𝐫,t)`$, where $`D`$ is the diffusion constant. A complex network obviously has an inherent discreetness, and hence no continuous limit can be taken. However, if the continuous diffusion equation is understood in its discreet form, one is lead to regard Eq. (2) as the master equation for the random-walk process taking place on the underlying network. By comparing Eqs. (1) and (2), as well as taking advantage of $`_jA_{ji}𝝆_i(t)/k_i=𝝆_i(t)`$ (because $`_iA_{ij}=k_j`$), one is lead to $`D_{ij}=A_{ij}/k_j𝜹_{ij}`$. Hence, an equivalent formulation of Eq. (2), is

$`𝝆(t+1)`$ $`=`$ $`𝐓𝝆(t),`$ (3)
where $`𝐓`$ is the transfer matrix, related to the diffusion matrix by $`𝐓=𝐃+\mathrm{𝟏}`$. In component form one thus has $`T_{ij}=A_{ij}/k_j`$. Since the transfer matrix, “transfers” the walker distribution one time step ahead, it can therefore be thought of as a time-propagator for the process. Formally the time development, from some arbitrarily chosen initial state $`𝝆(0)`$, can be obtained by iteration on Eq. (3) with the result $`𝝆(t)=𝐓^t𝝆(0)`$. Here $`𝐓^t`$ means the transition matrix to the $`t`$’th power.

In general, the transfer matrix (or equivalently, the diffusion matrix) will not be symmetric. However, $`𝐓`$ can be related to a symmetric matrix $`𝐒`$ by the following similarity transformation $`𝐒=\mathrm{𝐊𝐓𝐊}^1`$ with $`K_{ij}=𝜹_{ij}/\sqrt{k_i}`$. Thus, $`𝐓`$ and $`𝐒`$ will have the same eigenvalue spectrum, and all eigenvalues $`𝝀^{(𝜶)}`$ (of $`𝐓`$) will be real. It is this eigenvalue spectrum that will control the time-development of the diffusive process through $`𝚲^t`$ with $`\mathrm{\Lambda }_{ij}=𝜹_{ij}𝝀^{(𝜶)}`$. It should be noted that since the total number of walkers is conserved at all times, one must have $`1<𝝀^{(𝜶)}1`$, and that at least one eigenvalue should be one.111For the diffusion matrix, this means that its spectrum, $`\{𝝀_D^{(𝜶)}\}`$, must satisfy $`2<𝝀_D^{(𝜶)}0`$ since $`𝝀_D^{(𝜶)}=𝝀^{(𝜶)}1`$. We will here adopt the convention and sort the eigenvalues so that $`𝜶=1`$ corresponds to the largest eigenvalue $`𝝀^{(1)}=1`$, $`𝝀^{(2)}`$ to the next to largest one, and so on. Physically the principal eigenvalue $`𝝀^{(1)}=1`$, corresponds to a stationary state where $`𝝆(\mathrm{})𝝆^{(1)}`$, the diffusive current flowing from node $`i`$ to node $`j`$ is exactly balanced by that flowing from $`j`$ to $`i`$. The stationary state is unique for single component networks, and is fully determined by the connectivities of the network according to $`𝝆_i(\mathrm{})=k_i/(_ik_i)k_i`$. This can be easily checked by substituting this relation into Eq. (3). All modes corresponding to eigenvalues $`|𝝀^{(𝜶)}|1`$, are decaying modes since the time dependence of $`𝝆(t)`$ enters through $`𝚲^t`$ and one recalls that $`|𝝀^{(𝜶)}|1`$. Notice that $`𝝀^{(𝜶)}>0`$ represent non-oscillatory modes, while $`𝝀^{(𝜶)}<0`$ correspond to states where oscillations will take place with time, but this latter possibility will not be considered here.

The large scale topology of a given complex network reflects itself in the statistical properties of its diffusion eigenvectors $`𝝆^{(𝜶)}`$. One such property is the Participation Ratio (PR), that will be defined below. It quantifies the effective number of nodes participating in a given eigenvector with a significant weight. Since the stationary state of a network depends on the connectivities of its nodes, $`𝝆_i(\mathrm{})k_i`$, it is convenient to introduce a normalized eigenvector

$`c_i^{(𝜶)}(t)`$ $`=`$ $`{\displaystyle \frac{𝝆_i^{(𝜶)}(t)}{k_i}}.`$ (4)
Observe that $`c_i^{(𝜶)}`$ is nothing but the walker density per link of node $`i`$. Hence, in effect, $`c_i^{(𝜶)}`$ represents the outgoing currents flowing from node $`i`$, along each of its links, toward its neighbors. In the steady state, where $`𝝆_i(\mathrm{})k_i`$, it follows that these currents are all the same for any link in the network. Hence, highly connected nodes are not treated differently from less connected nodes. More formally, $`𝐜^{(𝜶)}`$ are the eigenvectors of the transposed transfer matrix $`𝐓^{}`$ corresponding to the same eigenvalue $`𝝀^{(𝜶)}`$ (of $`𝐓`$). Here one alternatively could have defined $`𝐓^{}`$ (instead of $`𝐓`$) as the transfer matrix from the very beginning. Doing so, would have resulted in a master equation of the form (3), but for the currents $`𝐜(t)`$ with corresponding eigenmodes $`𝐜^{(𝜶)}`$. The physical interpretation of such an alternative equation is as follows: Instead of walkers, think of a signal propagating on the underlying network. The signal at node $`i`$ at time $`t+1`$ is then the average of the signal at the nearest neighbors at time $`t`$.

If the link currents are normalized to unity, i.e. if $`𝐜^{(𝜶)}_2=1`$ in the $`L_2`$-norm, then the participation ratio (PR) is defined as

$`𝝌_𝜶`$ $`=`$ $`\left[{\displaystyle \underset{i=1}{\overset{N}{}}}\left(c_i^{(𝜶)}\right)^4\right]^1.`$ (5)
For the stationary state, where all the currents are equal to $`c_i^{(1)}=1/\sqrt{N}`$, one has $`𝝌_1=N`$ where we recall that $`N`$ is the total number of nodes in the network. Hence, when $`𝜶1`$, the ratio $`𝝌_𝜶`$ can be regarded as an effective number of nodes participating in the eigenvector $`𝐜^{(𝜶)}`$ with a significant weight . It should be noted that the participation ratio is a simple and crude measure of size. Strictly speaking, for this ratio to be able to say something with confidence about the size of a given module, the main contribution to $`𝝌_𝜶`$ should come from nodes within that module. For instance, if $`𝝌_𝜶`$ gets major contributions from current elements $`c_i^{(𝜶)}`$ of different signs, i.e. from different topological structures, it is not trivial to relate $`𝝌_𝜶`$ to the size of a single module.

For all types of complex networks where a modular structure is of interest, it is important to try to quantify the number of modules. Such information is of interest since this number may say something about the organizing principles being active in the generating process of the network. However, the measurement of this number is often hampered by statistical uncertainties due to temporary (non-robust) topological structures being counted as modules. Hence the challenge is to obtain the significant (or robust) number of modules that are due to the “rule of creation” and not just happened to be there by chance. For this purpose, it is useful to introduce a null model — a randomized version of the network at hand — but so that the degree distribution of the original network is not changes . A rough estimate of the number of different modules contained in a network could be given by the number of slowly decaying non-oscillatory modes that have a participation ratio, $`𝝌_𝜶`$, significantly exceeding the (ensemble) averaged participation ration of the corresponding randomized network, $`𝝌_𝜶^{\text{rd}}`$. Hence, it is suggested that a module is significant if $`𝝌_𝜶𝝌_𝜶^{\text{rd}}`$.

We will close this section by noting that the (outgoing) currents, $`c_i^𝜶`$, can also be used for vitalization purposes. The basis of this approach is the observation that the outgoing currents on links within a module are almost constant, and the size of this constant depends on to which extent the module participates in the given eigenmode. The constants thus varies from module to module and from eigenmode to eigenmode. Hence, by sorting $`c_i^{(𝜶)}`$ by size, nodes belonging to the same module will be located close to each.

## 3 The data sets

During a two year period in the mid 1970’s, W. Zachary studies the relations among $`34`$ members of a university karate club during a period of trouble . A serious controversy existed between the trainer of the club (node $`1`$ in Fig. 1a) and its administrator (node $`34`$). Ultimately this conflict resulted in the breakup of the club into two new clubs of roughly the same size. In his original study , W. Zachary mapped out the strength of friendship between the various members. In our study, however, we will be considering the unweighted version of his network  (see Fig. 1a). Recently, this same network was used by Girvan and Newman  in their study of community structures.

The second network that will be considered is a much larger network taken from the organization of the Internet. The Internet consists of a large number of individual computers that are identified by their so-called IP-address, and they are usually grouped together in local area networks (LAN). Such computer networks are connected to one another via routers — a complex network linking device. In addition to being able to direct pieces of information to its intended destination, a router also has the ability to determine the best path to a given destination (routing). This is done by keeping available an updated routing table telling the router how to reach certain destinations specified by the network administrator. In much the same way as computers are being organized into networks, the same is done for routers. If the router network becomes large and coherent enough, it may make out what is called an Autonomous System (AS). An AS is a connected segment of a network that consists of a collection of subnetworks (with hosts attached) interconnected by a set of routes. Usually it is required that the subnetworks and the routers should be controlled by a coherent organization — like, say, a university, or a medium-to-big enterprise. Importantly for the efficiency of the Internet, each Autonomous System, identified by a unique AS number, is expected to present to the other systems a consistent list of destinations reachable through that AS, as well as a coherent interior routing plan.

The particular network we will be considering is an Autonomous System network collected by the Oregon Views project on January 3, 2000 . In this network the AS will act as nodes, while their routing plans, i.e. with which other AS a given one directly shares information, will correspond to the links of the network.

## 4 Results and discussion

### 4.1 The Zachary Friendship Network

In Sec. 2, it was argued that from the currents corresponding to the slowly decaying eigenmodes, one should be able map out the large scale structure of the underlying network. Our hypothesis will now be put to the test, and we start with the small friendship network due to Zachary . This network is depicted in Fig. 1a, and it is small enough to enable us to see how it all comes about. By visual inspection of this figure, it is apparent that there are, at least, two large scale clusters — one corresponding to the trainer of the karate club (node $`1`$) and his supporters, and one to the administrator (node $`34`$). In Fig. 1a, the karate club members (according to Zachary) supporting the trainer in the ongoing conflict are marked with squares ($`16`$ nodes), while the followers of the administrator are marked with open circles ($`18`$ nodes). Notice that this subdivision was done by Zachary , and no sophisticated clustering techniques were applied to obtain these results.

We will now see if the same clustering structure can be obtained by the method outlined in Sec. 2 of this paper. Let us start by considering the two most slowly decaying modes of the network, namely $`𝜶=2`$ and $`3`$ (with our ordering), corresponding to the eigenvalues $`𝝀^{(2)}=0.87`$ and $`𝝀^{(3)}=0.71`$, respectively. If one sorts the elements of the current vector $`𝐜^{(2)}`$ by size, and group them according to their signs, one recovers two clusters (see the abscissa of Fig. 1b). This grouping, or topological structure, fits nicely with the original assignments made by Zachary (Fig. 1a) and indicated by the open squares and circles in Fig. 1b. Hence it is the trainer-administrator separation of the karate club members that is mapped out by the $`𝐜^{(2)}`$-currents. Observe that node $`3`$, which has an equal number of links to supporters of the trainer and administrator, has a current that is almost zero. Our identification of this node is the only one that differs from that made originally by Zachary. Interestingly, our classification, including node $`3`$, fits that previously obtained by Girvan and Newman  in their hierarchical tree clustering approach.

If the same procedure was repeated, but for the currents $`𝐜^{(3)}`$, one would map out other topological features of the network. In Fig. 1b, a $`c^{(2)}c^{(3)}`$-plot is presented. Such a plot groups the nodes along two axis; the $`c^{(2)}`$-axis — corresponding to the trainer-administrator axis — and the $`c^{(3)}`$-axis. The most striking feature of Fig. 1b is the group of nodes corresponding to $`c^{(2)}<0`$ and $`c^{(3)}>0`$, i.e. to the following group of nodes $`\{5,6,7,11,17\}`$. In Fig. 1a, these nodes are located in the upper left corner of the graph. They are in addition to being connected among themselves, only connected to the rest of the network via the trainer. Thus, they represent a sub-cluster, and this is indeed what is apparent from Fig. 1b.

The participation ratios for the two slowest decaying diffusive modes were found to be $`𝝌_2=14.0`$ and $`𝝌_3=13.1`$, respectively. As can be seen from Fig. 1b, these participation ratios receive substantial contributions from both positive and negative current elements, i.e. from more then one topological structure. Hence, we suspect that these numbers being close to the size of the administrator and trainer “clan” (of size $`18`$ and $`16`$ nodes, respectively) is somewhat accidental.

### 4.2 The Autonomous System Network

We will now focus our attention on the Autonomous System (AS) network introduced in the previous Section. This network, consisting of $`N=6,474`$ nodes and $`12,572`$ undirected links, is so big that plotting it for the purpose of study its modular structure is not a practical option. Motivated by what we found for the small Zachary friendship network, we will now go ahead and use somewhat similar techniques for its study.

In Fig. 2 the participation ratio, $`𝝌_𝜶`$, of eigenvectors $`c_i^{(𝜶)}`$ (top) and the eigenvalue density (bottom) are plotted as functions of the corresponding eigenvalues $`1<𝝀_i^{(𝜶)}<1`$. The data for the Internet, that is an example of a Scale-Free Network , are displayed together with the data for its randomized counterpart (the null model). From Fig. 2 it is apparent that while the density of states is rather similar for these two networks, the participation ratios of the slowly decaying modes, especially for $`𝝀^{(𝜶)}`$ close to $`1`$, are markedly higher in the Internet network than in the accompanying randomized network. These differences signal large scale topological structures  that are real and not accidental. For the most slowly decaying diffusive eigenmode, the participation ratio is $`𝝌_2100`$ (cf. Fig. 2a). From Fig. 3, that depicts the currents of the two slowest decaying diffusive eigenmodes, one observes that the main contribution to $`𝝌_2`$ comes from current elements $`c_i^{(2)}`$ of the same sign. Thus, a module has been detected, and $`𝝌_2`$ roughly measures its size. From, Fig. 3 one also makes the interesting observation that the main contributing nodes to $`𝝌_2`$ are Autonomous Systems located in Russia (denoted by open squares in Fig. 3). Hence, the diffusive mode $`𝜶=2`$ maps out Russia, as a Russian module ! In total there are $`174`$ Russian nodes in our data set. Moreover, from Fig. 3, modules corresponding to countries like the USA, France and Korea are also easily identified. The edges of the AS Internet network, defined as the nodes corresponding to the most distant values of the current elements $`c_i^{(2)}`$, are in this case represented by a Russian and a US military site located in the South Pacific. These nodes can thus be said to represent the extreme edges of the Internet. One has also studies the number of significant modules (where $`𝝌_𝜶𝝌_𝜶^{\text{rd}}`$), and found that this number for the AS network is roughly $`M=100`$ . The number of different nodes participating in one or more of these $`M`$ modules was found to be about $`N_M1800`$. Hence the modularity of the network is (at least) $`N_M/N30\%`$.

### 4.3 Generic features of the current-current plots

One common feature of the current-current plots, Figs. 1b and 3, is their line or star-like structures. Such structures are in fact rather generic, and we will here try to explain why. We have argued earlier that current elements $`c_i^{(𝜶)}`$ should almost be the same within a module. Hence, for two nodes $`i`$ and $`j`$ belonging to one and the same module, the fraction $`c_i^{(𝜶)}/c_j^{(𝜶)}`$ is expected to be a constant unique for that module, but independent of the diffusive mode $`𝜶`$. For two different (significant) diffusive modes, say $`𝜶`$ and $`𝜷`$, the fraction $`𝜸=c_i^{(𝜶)}/c_i^{(𝜷)}`$ (for node $`i`$) will therefore as a consequence also be constant. Thus, under the assumptions made above, one predicts a straight line in a current-current plot (that pases through the origin)  : $`c_i^{(𝜶)}=𝜸c_i^{(𝜷)}`$, for nodes $`i`$ belonging to the same module. From Fig. 3, this simple argument appears to be valid. That Fig. 1b seems to not pass through the origin, but still being straight lines, we believe is due to the diffusive modes being excited over a non-trivial background in this case.

### 4.4 Numerical implementation

Before presenting the conclusions of this paper, we will add a few closing remarks on the numerical implementation of the method. For this work to have any relevance for large real-world networks, a fast, memory saving, and optimized algorithm for the calculation of the largest eigenvalues and the corresponding eigenvectors of a sparse matrix is required. Fortunately, such an algorithm has already been implemented and made available e.g. through the TRLan software package . This software is optimized for handling large problems, and it can run on large scale parallel supercomputers.

## 5 Conclusions

We have generalized the normal diffusion process to diffusion on (discreet) complex networks. By considering such a process, it has been demonstrated that topological properties — like modular structures and edges — of the underlying network can be probed. This is achieved by focusing on the slowest decaying eigenmodes of the network. The use of the procedure was exemplified by considering a small friendship network, with known modular structure, as well as a routing network of the Internet, where the structure was not so well known. For the friendship network the known structure was well reproduced, and the Internet was indeed found to be modular. The detected modules of the Internet were consistent with the geographical location of the nodes, and the individual modules corresponded roughly to the national structure. Interestingly, it was observed that a political subdivision of the Internet was also one of the predictions of the algorithm presented in this paper; The two most poorly connected nodes of the Internet (extreme edges), where found to be represented by a Russian and a US military site located in the South Pacific.

## Acknowledgment

Work at Brookhaven National Laboratory was carried out under Contract No. DE-AC02-98CH10886, Division of Material Science, U.S. Department of Energy.


text:
# Geometrical and physical properties of maths-type 𝑞-deformed coherent states for 01

## Acknowledgments

C.Q. would like to thank J.R. Klauder for an interesting discussion. She is a Research Director of the National Fund for Scientific Research (FNRS), Belgium.

## Figure captions

Fig. 1. The Mandel parameter $`Q_q(t)`$ as a function of $`t=|z|^2`$ for $`q=0.98`$ (solid line), $`q=0.96`$ (dashed line) and $`q=0.94`$ (dot-dashed line).

Fig. 3. The variance ratio $`R_q(t)`$ as a function of $`t=|z|^2`$ for imaginary $`z`$ and $`q=0.98`$ (solid line), $`q=0.96`$ (dashed line) and $`q=0.94`$ (dot-dashed line).


text:
# The Effect of Metallicity on the Cepheid Period-Luminosity Relation

## 1. Introduction

Ever since the work of Edwin Hubble, the Cepheid Period-Luminosity ($`PL`$) relation is a fundamental tool in determining Galactic and extragalactic distances. In spite of its paramount importance, to this day we still lack firm theoretical and empirical assessment on whether or not chemical composition has any significant influence on the pulsational properties of Cepheids.

The debate on metallicity effects on the Period-Luminosity relation is an old one (see, for example, the reviews by Madore & Freedman 1991 and, more recently, Feast 2003), but the matter is far from being settled. As an illustration, different prescriptions on how to account for the metallicity effect lead to differences as large as 10% in the the distance to the Large Magellanic Cloud (e.g. Groenewegen & Oudmaijer 2000), which directly translates into a 10% uncertainty on $`H_0`$. To make things worse, this effect is *systematic* and only a careful calibration of the dependence of the $`PL`$ relation on metallicity can provide firmer ground on which to build the astronomical distance ladder. Regrettably, the current understanding of the subject is poor.

Theoretical pulsational models by different groups lead to markedly different results. On the one side computations based on *linear models* (e.g. Chiosi et al 1992, Sandage et al 1999, Alibert et al 1999, Baraffe & Alibert 2001) suggest a mild dependence of the $`PL`$ relation on chemical composition (less than 0.1 mag at $`\mathrm{log}(P)=1`$ for a change in metallicity between $`Z=0.004`$ and $`Z=0.02`$). This result is challenged by the outcome of the *non-linear convective models* (e.g. Bono et al 1999, Caputo et al 2000) which find that both the slope and the zeropoint of the Period-Luminosity relation depend significantly on the adopted chemical composition. Again for $`\mathrm{log}(P)=1`$ and the same variation in metallicity as above, they predict a change as large as 0.4 magnitudes in $`V`$, 0.3 magnitudes in $`I`$ and $`0.2`$ magnitudes in $`K`$. Moreover, the change is such that metal-rich Cepheids are fainter than metal-poor ones, again at variance with the results of Alibert et al (1999) and Baraffe & Alibert (2001). Finally, recent calculations by Fiorentino et al (2002) based on non-linear models indicate that the $`PL`$ relation also depends on the helium content of the stars.

From an observational standpoint, the subject was approached in essentially two ways: a direct measurement of the iron content in nearby, *i.e.* bright, Cepheids, or by measuring a secondary metallicity indicator in external galaxies known to contain Cepheids, under the assumption that they would have the same chemical composition. Regrettably, *direct metallicity determinations* for individual Cepheids exist only for a handful of stars. The most accurate studies are by Fry & Carney (1997, \[Fe/H\] and \[$`\alpha `$/H\] for 23 Galactic Cepheids) and Andrievsky et al (2002, complete abundance analysis of 77 Galactic Cepheids; see also Luck et al 2003). The former find a spread in \[Fe/H\] of 0.4 dex, which they argue is real, the latter focus on investigating Galactic abundance gradients ($`0.05`$ dex/kpc). Luck et al (1998) have measured the iron abundance of 16 Cepheids in the Large and Small Magellanic Cloud (LMC and SMC, respectively) with the echelle spectrograph at the CTIO 4 meter telescope. The relatively small size of the telescope forced them to choose only stars at the bright end of the Period-Luminosity relation, *i.e.* long periods, whereas the entire range needs to be sample in order to reach meaningful conclusions on the $`PL`$ relation (e.g. Kennicut et al 1998, Caputo et al 2000). As for the *indirect measurements* in external galaxies (e.g. Kochanek et al 1997, Kennicutt et al 1998), which probe a broader range in metallicity than the one available in the Galaxy, they find that metal-rich Cepheids are brighter in $`V`$ and $`I`$ than metal-poor ones by $`0.24`$ to $`0.40`$ magnitudes per dex in metallicity. These results have been widely used in the literature, including in the final paper of the hst Key Project on $`H_0`$ (Freedman et al 2001).

## 2. The Data Sample

With the scenario described above in mind, we have collected high quality spectra of an unprecedented dataset of 40 Galactic, 22 LMC and 14 SMC Cepheids with the aim of measuring metallicity directly on a statistically significant sample of stars covering about one order of magnitude in chemical composition. The Galactic stars were observed with FEROS at the ESO 1.5m telescope on La Silla at a spectral resolution of 48,000 with a signal-to-noise ratio of about 70 to 150. As for the spectra of the Magellanic Cloud objects, we used the UVES spectrograph on the Kueyen telescope on Cerro Paranal, yielding a resolution of 40,000 and signal-to-noise ratio of about 50 to 70.

Deriving metallicities for Cepheids is by no means a trivial task, but we have now carefully analyzed and derived \[Fe/H\] for an initial subsample of 14 Cepheids in the LMC, 12 in the SMC and 6 in the Galaxy. The data reduction procedures and analysis are described in detail in Mottini et al (2003).

## 3. Results: dependence of the PL relation on metallicity

The main result from the subset analyzed so far is summarized in Figure 1. There, we plot the V-band residuals from the standard Madore & Freedman (1991) $`PL`$ relation as a function of the iron abundance derived from FEROS and UVES spectra. A positive magnitude difference here means that the star is fainter than the mean $`PL`$ relation. For consistency with the adopted $`PL`$ relation, the LMC is assumed to have a distance modulus of 18.50. The SMC is considered 0.44 mag more distant (e.g. Cioni et al 2000) and the distances to Galactic stars were taken from Fouqué et al (2003), Barnes et al (2003) and Tammann et al (2003). Photometry and reddening are taken from Laney & Stobie (1994).

The data are binned in metallicity to reflect the typical observational error on our determination ($`\pm 0.1`$ dex). The Magellanic Cloud Cepheids populate the bins at lower metallicity, the higher \[Fe/H\] one containing only Galactic stars. The median magnitude residual in each bin is plotted as a square and the rms and error on the median are shown as boxes and errorbars, respectively. The magnitude difference plotted on the ordinates can be directly compared to the *“c”* quantity computed by Fiorentino et al (2002). This is the correction to be applied to the distance modulus derived using a universal, metallicity independent $`PL`$ relation to take into account the influence of metallicity as derived from their non-linear pulsation models.

The agreement between our data and the predictions by Fiorentino et al (2002), indicating a strong dependence of the $`PL`$ relation on metallicity, is tantalizing, even though the errorbars are still somewhat large. As it can be seen, the error on the inferred distance can be as large as 10% if the effect of metallicity is neglected.

## References

Alibert, Y. et al 1999, A&A344, 551

Andrievsky, S.M. et al 2002, A&A, 381, 32

Baraffe, I. & Alibert, Y. 2001, A&A, 371, 592

Barnes, T.G. et al 2003, ApJ, 592, 539

Bono, G. et al 1999, ApJ, 512, 711

Caputo, F. et al 2000, A&A, 359, 1059

Chiosi, C. et al 1992, ApJ, 387, 320

Cioni, M.-R. et al 2000, A&A, 359, 614

Feast, M. 2003, in Stellar Candles for the Extragalactic Distance Scale, Lect.Notes Phys., 635, 45

Fiorentino, G. et al 2002, ApJ, 576, 402

Fouqué, P. et al 2003 in Stellar Candles for the Extragalactic Distance Scale, Lect.Notes Phys., 635, 21

Fry, A.M. & Carney, B.W. 1997, AJ, 113, 1973

Groenewegen, M. & Oudmaijer, R.D. 2000, A&A, 356, 849

Kennicutt, R.C. et al 1998, ApJ, 498, 181

Kochanek, C.S. 1997, ApJ, 491, 13

Laney, C.D. & Stobie, R.S. 1994, MNRAS, 266, 441

Luck, R.E. et al 1998, AJ, 115, 605

Luck, R.E. et al 2003, A&A, 401, 939

Madore, B. & Freedman, W. 1991, PASP, 103, 933

Mottini, M. et al 2003, this conference

Sandage, A. et al 1999, ApJ, 522, 250

Tammann, G. et al 2003, A&A, 404, 423