text:
# WIYN Open Cluster Study. X. The K-Band Magnitude of the Red Clump as a Distance Indicator

## 1 Introduction

During the past few years, the helium burning red clump (RC) has gained considerable attention for its potential as a standard candle. The primary advantage of the RC is the ease with which it can be recognized in the color-magnitude diagram. However, there is currently a great deal of controversy in the literature regarding the appropriate treatment of possible metallicity and age effects on the I-band absolute magnitude of the RC ($`M_I(RC)`$). There are two schools of thought on this issue; the first assumes a constant value for $`M_I(RC)`$ which is then used to facilitate a single-step distance determination via knowledge of the apparent RC magnitude and the extinction (e.g. Paczyński & Stanek 1998; Stanek & Garnavich 1998). The second approach is founded on the claim that both age and metal abundance have a significant influence on the luminosity of RC stars (e.g. Cole 1998; Sarajedini 1999) and must be accounted for in determining $`M_I(RC)`$ and therefore the distance.

Both Paczyński & Stanek (1998) and Stanek & Garnavich (1998) use HIPPARCOS RC stars with parallax errors of less than 10% to calculate the I-band absolute magnitude of the solar neighborhood red clump. In their analysis, Paczyński & Stanek (1998) find that $`M_I(RC)`$ shows no variation with color over the range $`0.8<(VI)_0<1.4`$ and, from a Gaussian fit to the RC luminosity function, find $`M_I(RC)=0.28\pm 0.09`$. Following the same methodology and building upon the earlier work, Stanek & Garnavich (1998) find a similar result with $`M_I(RC)=0.23\pm 0.03`$. With this calibration, a single step calculation is then used to determine the distance to the Galactic center (Paczyński & Stanek 1998) and M31 (Stanek & Garnavich 1998). Both of these investigations found little or no variation in $`M_I`$ of the RC stars with color; this was taken to imply that $`M_I(RC)`$ does not vary significantly with metallicity.

In contrast, theoretical models from Girardi & Salaris (2001) and the earlier models of Seidel, Demarque, & Weinberg (1987; see also Cole 1998) show that $`M_I(RC)`$ is dependent on both age and metallicity, becoming fainter as both increase. These models are in good agreement with the observations presented by Sarajedini (1999, hereafter S99). Using published photometry for 8 open clusters, S99’s most important result is that while $`M_I(RC)`$ is less sensitive to metal abundance than $`M_V(RC)`$, both still retain a considerable dependence on the age and metallicity of the stellar population. As a result, the single-step method of applying the solar-neighborhood $`M_I(RC)`$ to populations with a different age-metallicity mix could be problematic.

Alves (2000) also uses the HIPPARCOS RC for his calibration; however, he relies upon the K-band luminosity ($`M_K`$) of the RC stars in the hope that, since the K-band is less sensitive to extinction (and possibly metallicity as well) than the I-band, it might make a better choice as a standard candle. Alves (2000) restricts his RC stars to those that have metallicities from high resolution spectroscopic data. For this group of 238 RC stars, he finds a peak value of $`M_K(RC)=1.61\pm 0.03`$ with no correlation between $`[Fe/H]`$ and $`M_K`$. However, he is not able to explore the effect of age on $`M_K(RC)`$ due to the lack of such information for the individual stars in his sample.

These previous works have prompted us to combine the approaches of S99 and Alves (2000) and to investigate the influence of age and metal abundance on $`M_K(RC)`$ for a number of open clusters with well-known distances and metallicities. In section 2 we discuss the observational data. Section 3 compares our data with the results of theoretical models and presents a discussion of the results; our conclusions are summarized in Section 4.

## 2 The Data

### 2.1 Open Clusters

In the present study, the most important criterion that the observational data must fulfill is that of internal consistency. For example, we must ensure that the distance moduli, reddenings, ages, and metallicities of all of the clusters in our sample have been determined using the same techniques. In addition, it is imperative that the IR photometry we rely upon be measured and calibrated in a consistent manner. For the former, we use the database of open cluster properties as measured by Twarog et al. (1997), supplemented by cluster ages from the WEBDA111http://obswww.unige.ch/webda/webda.html database, and, for the latter, we utilize the 2nd Incremental Data Release of the 2MASS Point Source Catalog222http://irsa.ipac.caltech.edu/. We now discuss each of these in more detail.

Twarog et al. (1997) have compiled a list of 76 open clusters for which they provide reddenings, distance moduli and metallicities. For the purposes of the present paper, we limit ourselves to distance moduli derived via the technique of main sequence fitting so as to remain independent of methods that rely on the luminosity of the RC. Their metallicities have all been measured on the same system and the reddenings have been determined using an internally consistent method. The vast majority of these values are consistent with those found in the literature, except for the reddening of NGC 6819 for which the Twarog et al. (1997) value is much higher than other published values. As a result, we have decided to adopt the S99 reddening for NGC 6819 instead of the apparently discrepant value tabulated by Twarog et al. (1997). Additionally, because the determination of the reddening and distance modulus are coupled, we also adopt the S99 distance modulus for NGC 6819.

The ages of the open clusters have been obtained from WEBDA, which is a compilation of open cluster data from various sources. To check the reliability of the WEBDA ages, we compare them with the cluster ages determined by S99 in Fig. 1. S99 presents isochrone-fitting ages for 8 open clusters that have been determined in the same manner using the Bertelli et al. (1994) theoretical isochrones. The left panel of Fig. 1 plots the ages from S99 versus the ages in WEBDA, where both axes are in log space and the dashed line represents a zero age difference between the systems. From this plot, it is evident that there is a systematic offset between the two systems with the WEBDA ages being younger than those of S99. The average difference, $`\mathrm{\Delta }`$$`log(Age)=0.191`$, is used to shift the ages given in WEBDA onto the S99 system. The right panel of Fig. 1 shows the ages from S99 plotted against the shifted WEBDA ages; it is clear from Fig. 1 that the shifted ages are in better agreement with those of S99; as a result, we will apply this shift to nine of the clusters in our study and use the S99 ages for the five clusters that are common to both studies.

For the open clusters in the Twarog et al. (1997) study that possess main sequence fitting distances, we extracted $`JHK`$ photometry from the 2nd Incremental Data Release of the 2MASS Point Source Catalog. As noted above, these data have been obtained using similar instruments and reduced with the same pipeline techniques. For each cluster, we have utilized the same criteria for the 2MASS data retrieval. The field size is originally set to 30 arcmin in radius and then reduced to fields as small as 5 arcmin in radius in an attempt to isolate the cluster stars. The sources are limited to a brightness of 6th magnitude or fainter due to saturation effects at the bright end (Carpenter 2000). Lastly, we have extracted only the highest quality photometry from the 2MASS catalog, which provides a read flag (rd\_flg) indicating how the photometry of each star was measured. We have chosen to exclude any source that has a read flag of zero in any band since this implies that the source was not detected in that band and the magnitude given is an upper limit. We note that the vast majority of the stellar magnitudes used in this study are based on point-spread-function fitting (i.e. rd\_flg = 2); however, in order to include the brighter red clumps of more nearby clusters, we have had to use the aperture photometry in a minority of cases.

The 2MASS program uses a K-short ($`K_S`$) filter for their observations. We have chosen to convert these magnitudes to the K-band adopting the Bessell & Brett (1988, hereafter BB) system, which is also used in the theoretical models of Girardi et al. (2000) and Girardi & Salaris (2001). The transformation equations are derived by Carpenter (2001) and are adopted as follows:

$$(JK_{BB})=[(JK_S)(0.011\pm 0.005)]/(0.972\pm 0.006)$$
(1)
and

$$K_{BB}=[K_S(0.044\pm 0.003)](0.000\pm 0.005)(JK_{BB}).$$
(2)
We note, however, that transformation to the Koornneef (1983) K-band (used by Alves 2000, Sec. 3.2) would have a negligible effect on our results. To correct for the interstellar reddening, we adopt the extinction law determined by Cardelli et al. (1989), which, using their value of $`R_V=3.1`$, gives $`A_K=0.11A_V`$ and $`A_J=0.28A_V`$. From this, it is a simple matter to calculate the absolute K-band magnitude and dereddened $`JK`$ color of the open cluster stars.

We have determined the RC luminosity for our clusters by taking the median value of $`M_K`$ for all stars within a standard sized box placed around the RC. We use the median value of $`M_K`$ along with a constant box size in an attempt to eliminate any selection effects that may occur in choosing the location of the RC and to limit the effect of outliers on $`M_K(RC)`$. Figure 2 shows the color-magnitude diagrams (CMDs) for all 14 clusters focused on the RC and main sequence turnoff. The box used to select the RC stars for each cluster is also shown. We note that, where available, we have used published optical CMDs for our clusters to help isolate the approximate RC location. The uncertainty in $`M_K(RC)`$ is calculated by combining the standard error about the mean K-magnitude for all stars inside the red clump boxes along with the errors in $`E(BV)`$ and $`(mM)_V`$, all added in quadrature. Except where otherwise noted, we adopt 20% of the value as the error in $`E(BV)`$ and 10% of the value as the error in $`(mM)_V`$.

### 2.2 Globular Clusters

It is difficult to ensure that the distances, ages, and metallicities of globular clusters with RCs are on the same system as those of the open clusters. Fortunately, the ages and metallicities of the two globulars in our sample - 47 Tuc and NGC 362 - are sufficiently different from the bulk of the open clusters that small systematic discrepancies in these quantities should not be a significant hindrance to the interpretation of the results. In any case, we have decided to adopt literature values for the basic cluster parameters and use the globular cluster RCs as a consistency check.

In the case of 47 Tuc, we adopt the metallicity quoted by Carretta & Gratton (1997) of $`[Fe/H]=0.70\pm 0.07`$, which happens to be very close to the Zinn & West (1984) value. For the distance modulus and reddening, we average the published values listed in Table 2 of Zoccali et al. (2001) to obtain $`(mM)_V=13.45\pm 0.21`$ and $`E(BV)=0.044\pm 0.008`$, where the errors represent half of the range of tabulated values. Lastly, for the age of 47 Tuc, we adopt the oldest age for which the models predict the presence of a red clump at its metallicity - 12 Gyr.

For NGC 362, we adopt a similar approach. The metal abundance of $`[Fe/H]=1.15\pm 0.06`$ is taken from Carretta & Gratton (1997), which is approximately 0.1 dex more metal-rich than the Zinn & West (1984) value. Our search of the literature has revealed distance moduli that range from $`(mM)_V`$ = 14.49 (Zinn 1985) to 14.95 (Burki & Meylan 1986; see also Bolte 1987) and reddenings in the range $`E(BV)`$ = 0.032 (Vandenberg 2000) to 0.08 (Alcaino 1976) leading to adopted values of $`14.70\pm 0.23`$ and $`0.048\pm 0.024`$ for the apparent distance modulus and reddening of NGC 362, respectively. Once again, we adopt an age of 12 Gyr.

The red clumps of these globulars have been isolated in the 2MASS point source catalog in the same way as for the open clusters. The ($`M_K,(JK)_0`$) CMDs for 47 Tuc and NGC 362 are shown in Fig. 3 along with the box used to define their RCs. All of the relevant observational parameters for the open and globular clusters are listed in Table 1.

## 3 Results and Discussion

### 3.1 Cluster Data

As described in Sec. 1, we are interested in exploring the dependence of $`M_K(RC)`$ on $`[Fe/H]`$ and age. Plotted in Fig. 4 are the $`M_K(RC)`$ values for the 14 open clusters (open circles) and the two globulars (filled circles) in our sample versus the logarithm of the age (top right panel) and the metallicity (top left panel). The red clumps start out very bright at young ages and decrease in brightness by almost one magnitude as the cluster ages approach $`10^9`$ years after which they brighten by $``$ 0.5 magnitude. Then, the red clumps become slightly fainter as the cluster ages increase up to $`10^{10}`$ years. The top two panels of Figure 4 also include $`M_I(RC)`$ (open squares) from S99. Keeping in mind that the numbers of clusters is small, over the age and metallicity range common to both studies, the K-band absolute magnitude of the RC exhibits less sensitivity to age and metal abundance than does $`M_I(RC)`$.

### 3.2 Field Star Data

Alves (2000) reports K-band absolute magnitudes of solar-neighborhood stars in the HIPPARCOS catalog along with parallaxes and proper motions. Following the analysis of Alves (2000), we have limited ourselves to stars with $`2.2<(VK)_0<2.5`$ and $`2.5<M_K<0.8`$ in Table 1 of Alves’ paper in order to isolate a sample of nearby red clump stars. We have plotted $`M_K`$ vs $`[Fe/H]`$ for these stars in the bottom panel of Figure 4 (open circles). For comparison, the open and globular cluster data from the present work (filled circles) are also shown. Keeping in mind that there are far fewer open clusters than field stars in Fig. 4, we find good consistency in the locations of the two samples. Alves (2000) finds $``$$`M_K(RC)`$$``$$`=1.61\pm 0.03`$, while the open clusters in our sample give $``$$`M_K(RC)`$$``$$`=1.62\pm 0.06`$. This is remarkable given the fact that the open cluster distances are based on the main sequence fitting results of Twarog et al. (1997) and the field stars are on the HIPPARCOS distance scale. Both show mean K-band magnitudes of $``$–1.6 and very little if any dependence on metal abundance over the same range.

### 3.3 Comparison With Theoretical Models

Leo Girardi has kindly provided us with theoretical models that represent the median magnitude of the red clump as a function of age and metal abundance (Girardi & Salaris 2001; Crowl et al. 2001). Figures 5 and 6 show these theoretical models in the K-band compared with our open and globular cluster data. The nine panels of Fig. 5 display the K-band luminosity of the red clump as a function of metallicity for a range of ages. The five panels in Fig. 6 illustrate the variation of the K-band luminosity with age for a range of metal abundances. In both of these figures, clusters that are similar in metallicity or age to the model plotted in each panel are represented by filled circles while the remaining clusters are denoted by open circles. Figures 5 and 6 suggest that at ages younger than $``$2 Gyr, the red clump luminosity is greatly dependent on the age of the cluster and shows little effect from the metallicity whereas clusters older than $``$2 Gyr show the exact opposite, having little age dependence while still showing the effects of metallicity.

To facilitate a more detailed comparison between the models and the observations, we utilize an interpolation routine based on low order polynomials to compare the theoretical $`M_K`$ values with the observed ones. As a test of the interpolation, we have applied it to the observational data alone comparing the interpolated $`M_K(RC)`$ values to the actual values at the age and abundance of each cluster. We find that the root-mean-square of the residuals is negligible in $`M_K(RC)`$ with no systematic trends as a function of age or abundance. We have also tested the interpolation routine on the theoretical models with similar encouraging results. The accuracy of the interpolation allows us to compare the $`M_K(RC)`$ values predicted by the models for a given age and metallicity to the observed $`M_K(RC)`$ for each cluster. From this comparison, we find that the root-mean-square deviation of the theoretical models from the observations is 0.16 mag, with no systematic variation in the residuals as a function of age or metallicity. This deviation is slightly larger than the mean error in the $`M_K(RC)`$ values of the 16 clusters, which we find to be 0.13 mag. Given that the mean deviation of the models from the observations is roughly consistent with the errors inherent in the latter, it is reasonable to conclude that the models are generally consistent with the observational data.

Figure 7 shows the Girardi models plotted along with the Alves (2000) field red clump star data. The models reinforce the conclusion drawn by Alves (2000) that $`M_K`$ is insensitive to metallicity for nearby stars in this abundance range. Furthermore, given that a typical $`M_K`$ error in Alves (2000) data is 0.11 mag, this figure suggests that the vertical spread in the $`M_K`$ values is mainly the result of age effects among the field stars. Both the $`10^{8.8}`$ and $`10^{9.29.4}`$ year isochrones agree with the majority of the data. However, given the expectation that stars in the Solar neighborhood are likely to be near Solar-age, most of the HIPPARCOS stars in Fig. 7 probably have Log ages between 9.2 and 9.6 (1.6 to 4.0 Gyr).

It is interesting to note that the ages of the solar neighborhood stars as predicted by the models show a lack of stars around $`10^9`$ years (Fig. 7). In contrast, using a model of the solar neighborhood RC that assumes a constant star formation rate, Girardi & Salaris (2001) expect an age distribution for the RC stars that peaks at 1 Gyr with approximately 60% of the stars having this age. The discrepancy in this result with the apparent ages of the HIPPARCOS RC stars likely indicates a non-constant star formation rate in the solar neighborhood; this is not surprising if the formation of stars is triggered by density waves traveling through the solar neighborhood, which is an intrinsically episodic process.

## 4 Application As a Distance Indicator

An important aspect of this study is the application of the K-band red clump absolute magnitude as a distance indicator. To optimize this application in the present work, we seek a range of age and abundance over which variations in $`M_K(RC)`$ are minimized. Inspecting Fig. 4, we see that if the age of the stellar population is in the range 2$`<`$Age$`<`$6 Gyr and the metal abundance is between –0.5$`<`$\[Fe/H\]$`<`$0.0, then the intrinsic variation in $`M_K(RC)`$ is minimized suggesting that uncertainties in our knowledge of these properties are inconsequential in the determination of the distance. Based on these considerations, we have selected the open cluster NGC 2158. This cluster possesses 2MASS photometry, and it is included in the study of Twarog et al. (1997), so we have a metallicity value ($`[Fe/H]=0.24\pm 0.06`$) which is on the same system as the other clusters in this study. The age shift described in section 2.1 is also applied to NGC 2158 giving us an age of $`1.6\pm 0.2`$ Gyr. We note in passing that NGC 2158 was not included as part of our $`M_K(RC)`$ calibration because the distance given in Twarog et al. (1997) was determined using the magnitude of the RC and not main sequence fitting.

For the reddening toward NGC 2158 we can utilize the data in Table 1 to parameterize the intrinsic color of the red clump ($`(JK)_0`$) in terms of the metal abundance and age. Figure 8 shows $`(JK)_0`$ versus $`[Fe/H]`$ (left panel) and age (right panel) for the clusters in our sample. Using the interpolation discussed in section 3.3, we can determine the intrinsic color of NGC 2158 given its metallicity and age, for which we find $`(JK)_0=0.597\pm 0.003`$. We calculate the error in $`(JK)_0`$ by determining the uncertainty resulting from $`\sigma _{age}`$ and $`\sigma _{[Fe/H]}`$ and adding these in quadrature. Comparing the implied intrinsic color of the red clump with the apparent color, $`(JK)=0.821\pm 0.006`$, we find $`E(JK)=0.224\pm 0.007`$. Converting this to a color excess in the optical regime, we find $`E(BV)=0.43\pm 0.013`$, which is in good agreement with published values (e.g. Christian et al. 1985; Twarog et al. 1997). The preceding method represents an internally consistent formalism which can be utilized to estimate the reddening of a cluster.

The interpolation on $`M_K(RC)`$ using only the open cluster data predicts $`M_K=1.64\pm 0.08`$ for NGC 2158. Along with $`E(BV)=0.43`$ and the apparent RC K-band magnitude, $`K(RC)=11.55\pm 0.02`$, we find $`(mM)_V=14.38\pm 0.09`$. Our distance modulus for NGC 2158 agrees within the errors with the main sequence fitting modulus of $`(mM)_V=14.4\pm 0.2`$ found by Christian et al. (1985), but is slightly lower than that determined by Twarog et al. (1997) of $`(mM)_V=14.5`$.

## 5 Conclusions

In this paper, we have sought to establish the K-band absolute magnitude of the helium burning red clump stars ($`M_K(RC)`$) as a distance indicator. To facilitate this, we have utilized infrared photometry from the 2MASS catalog along with distances, metallicities, and ages for 14 open clusters and 2 globular clusters. Our sample encompasses an age range from 0.63 Gyr to 12 Gyr and metallicities from –1.15 to 0.15 dex. Based on an analysis of these data, we draw the following conclusions.

1. There is a statistically significant range of $`M_K(RC)`$ values among the star clusters in our sample. In particular, for the 14 open clusters, we calculate $``$$`M_K(RC)`$$``$$`=1.62`$ with a standard deviation of 0.21 mag. In contrast, the mean error in these $`M_K(RC)`$ values is 0.13 mag.

2. Upon inspection of figures 5 and 6, we find that for clusters younger than $``$2 Gyr, $`M_K(RC)`$ is insensitive to metallicity but shows a dependence on age. In contrast, for clusters older than $``$2 Gyr, $`M_K(RC)`$ is influenced primarily by the metallicity of the population and shows little or no dependence on the age.

3. In general, $`M_K(RC)`$ is less sensitive to age and metallicity than $`M_I(RC)`$ over the parameter range common to both this paper and Sarajedini (1999) from which the $`M_I(RC)`$ values are taken.

4. Over comparable metallicity and age ranges, our average $`M_K(RC)`$ value of –1.62 mag is consistent with that of Alves (2000) which is based on solar-neighborhood red clump stars with HIPPARCOS parallaxes. We also suggest that the significant scatter in the Alves (2000) $`M_K`$ data is likely due to a range of ages between $``$1.6 and $``$4 Gyr among these stars.

5. The theoretical red clump models based on the formalism of Girardi et al. (2000) agree reasonably well with our observational data, indicating that age plays an important role in determining $`M_K(RC)`$ for younger populations while metallicity mainly affects older populations.

6. Using the K-band absolute magnitude of the red clump, we are able to compute the distance to the open cluster NGC 2158. Adopting an age of $`1.6\pm 0.2`$ Gyr and $`[Fe/H]=0.24\pm 0.06`$, our calibration yields a distance of $`(mM)_V=14.38\pm 0.09`$.

7. When determining distances for star clusters having $`0.5[Fe/H]0.0`$ and $`10^{9.2}age10^{9.9}`$, one can ignore the interpolation discussed in section 3.3 and simply use $`<M_K(RC)>`$ = $`1.61\pm 0.04`$.

We would like to thank Leo Girardi and David Alves for providing us with electronic copies of their data. We would also like to thank Glenn Tiede, Ted von Hippel, Barbara Anthony-Twarog, and Bruce Twarog for making useful comments that greatly improved the clarity of an early manuscript. The comments of an anonymous referee are also greatly appreciated. A.S. is grateful to the National Science Foundation for support through grant No. AST-9819768. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.


format: text

text:
# Spectral flow in superconducting point contacts

## I Introduction

Effect of impurities on transport in superconducting nanostructures is one of the key issues in the physics of mesoscopic superconductors, among which point contacts are the simplest devices of interest. One can distinguish contacts with only a few transverse eigenmodes (or conducting channels) relevant for a given temperature (quantum contacts) and those where the number of channels $`N_{\mathrm{ch}}`$ is relatively large (mesoscopic contacts). In quantum contacts, the effect of impurities on both static and transport properties can be satisfactorily described by introducing certain transmission and reflection probability $`D+R=1`$ for each conducting channel Beenakker/rev . Mesoscopic weak links with disorder are sometimes also treated within this scheme. To this end, a quantum eigenstate problem for a charge carrier is to be solved for the particular realization of disorder, each state being considered as an individual channel with certain transmission probability. The problem of a multichannel weak link is then mapped into that for a set of tunnel contacts having an effective distribution of transmission probabilities. This scheme has been shown recently Kopnin/pointcont  to be equivalent to the semi-classical theory of superconductivity for calculating some static properties of multichannel contacts.

The individual channels found by solving the quantum-mechanical problem are independent in the sense that transitions between them are absent under static conditions. One might now raise a question how useful is the model of independent channels for transport properties of multichannel contacts when time-dependent perturbations induce transitions between the individual channels. This is of a crucial importance for further theoretical studies of transport in superconducting mesoscopic multichannel devices in presence of disorder. We address this problem using a specific example of a ballistic point contact between two clean superconductors assuming that the thickness $`d`$ of the insulating layer and the size $`a`$ of the orifice in it are much larger than atomic dimensions but shorter than the coherence length $`\xi `$ and the mean free path $`\mathrm{}`$, i.e., $`ad\xi \mathrm{}`$. The contact thus has a very large number of conducting channels, $`N_{\mathrm{ch}}\left(p_Fa\right)^21`$. The spectrum of the midgap states localized within the coherence length $`\xi `$ near the point contact was first found in Kulik

$$ϵ_\pm =|\mathrm{\Delta }|\mathrm{cos}(\varphi /2)$$
(1)
for a phase difference $`\varphi `$ between the two superconductors. The spectrum is degenerate: all particles flying from one superconductor into another have the same energy independently of their momentum directions within the corresponding hemispheres. Another example is a contact having a tunnel barrier with a transmission probability $`D`$. The spectrum barrierspect

$$ϵ=\pm |\mathrm{\Delta }|\sqrt{1D\mathrm{sin}^2(\varphi /2)}$$
(2)
has a (mini)gap $`\sqrt{R}|\mathrm{\Delta }|`$ near the phase difference $`\varphi =\pi `$, which separates the two branches of the spectrum. The mid-gap states have a profound effect on dynamics of the contacts. When the applied voltage is larger than the inelastic relaxation rate, $`eV\tau _ϵ^1`$, a nonlinear regime sets in which originates from multiple Andreev reflections (MAR) Zaikin ; Cuevas ; AverinBardas2 ; AverinBardas1 ; AverinBardas3 ; Bratus  of localized quasiparticles. For a contact without both impurities and a barrier, excitations from below the gap edge with $`ϵ=|\mathrm{\Delta }|`$ are accelerated through the midgap states Eq. (1) and re-appear above the gap with $`ϵ=|\mathrm{\Delta }|`$ and vice versa, thus producing a highly nonequilibrium situation. On the contrary, if a contact has a barrier with a reflectivity $`R`$, the localized excitations cannot accelerate due to the minigap and retain their original states. It is only Zener tunnelling through the minigap that provides nonequilibrium current through the contact AverinBardas1 ; AverinBardas3 ; Bratus .

The model of independent channels was applied to kinetic effects in multichannel point contacts in Ref. AverinBardas3 . The key ingredients of that model were (i) the MAR processes in separate channels with (ii) Zener transitions through a minigap between two spectral branches in Eq. (2) within each individual channel, and (iii) a standard distribution Dorokhov  of reflection coefficients over the channels. Here we study the transport in mesoscopic multichannel contacts with impurities and show that the above model of individual channels consisting of processes (i)–(iii) is not suitable for dynamic phenomena when the number of channels is large, $`N_{\mathrm{ch}}1`$, because transitions between different channels are excited in time-dependent conditions. As a result, the system finds the optimal way to relax through channels with small reflection probabilities avoiding weak Zener tunneling. We further demonstrate that the concept of MAR is equivalent to the spectral flow appearing in a wide variety of quantum systems specflow ; KopVol/specflow ; Stone . Using the quasiclassical theory of nonstationary superconductivity we show that the impurity scattering provides an effective momentum exchange between particles moving from the opposite sides of the contact and impedes the MAR acceleration. For low inelastic relaxation rate, the I–V curve displays a non-trivial behavior: At low voltages, impurity relaxation slows the acceleration down and results in a new linear current-voltage dependence with a large conductance. This differs from a square-root dependence obtained within the independent-channel model (i)–(iii).

## II Kinetic equation

The density of states (DOS) for the midgap states in a classical point contact was calculated in Refs. Kulik ; Kopnin/pointcont . It is $`\nu (s)=\nu (0)e^{k|s|}`$ with $`k=2v_F^1\sqrt{|\mathrm{\Delta }|^2ϵ^2}`$ where $`s`$ is the distance along the particle trajectory; $`\nu `$ is normalized to the normal single-spin DOS at the Fermi level, $`N_F`$. We use the standard kinetic equation approach LO ; KL  where the distribution function is split into two components $`f^{\left(1\right)}`$ and $`f^{\left(2\right)}`$ that are respectively odd and even in the variables $`ϵ,𝐩`$. They are coupled by the Boltzmann equation

$$\nu v_F\frac{1}{2}\frac{\dot{\chi }}{s}\frac{f^{(1)}}{ϵ}+\nu \frac{f^{(1)}}{t}+v_F\frac{}{s}\left(\nu f^{(2)}\right)=J.$$
(3)
The second equation ensures that the distribution function $`f^{\left(1\right)}`$ is constant along the trajectory KL . In Eq. (3), the first term in the l.h.s. originates from the Doppler energy and $`\dot{\chi }`$ is the time derivative of the phase, and $`J`$ is the impurity collision integral. A detailed discussion of the kinetic equations is given in Refs. LO ; KL . The standard technique of averaging over impurities applies since the number of impurities within the volume of the orifice is large. Indeed, the mean free time is $`\tau ^1N_Fn_{\mathrm{imp}}|u|^2`$ where $`|u|Up_F^3`$ is the Fourier transform of the impurity potential $`U`$. The number of impurities $`n_{\mathrm{imp}}da^2(d/\mathrm{})(E_F/U)^2N_{\mathrm{ch}}`$ can be very large even for $`UE_F`$ because of a macroscopic number of quantum channels in the orifice $`N_{\mathrm{ch}}1`$. The large number of impurities is the essential parameter in our theory.

We integrate Eq. (3) along the trajectory across the contact over distances longer than $`\xi `$. The first term in the l.h.s. of Eq. (3) has a sharp maximum at $`sd`$ where $`\nu (s)`$ is nearly constant. The third term in the l.h.s. disappears due to the decay of $`\nu `$. The collision integral contains Kopnin/pointcont  a combination of angular averages $`\nu \left[\nu f^{(1)}\nu f^{(1)}\right]`$. We write $`\nu =\left[\nu _++\nu _{}\right]/2`$ where $`\nu _\pm =(2\pi )^1_{\pm p_x>0}\nu (𝐩,𝐫)𝑑\mathrm{\Omega }_𝐩`$. The upper (lower) sign is for right-moving (left-moving) particles with $`p_x>0`$ ($`p_x<0`$). We choose the direction of the $`x`$ axis from the region with the phase $`\chi _1`$ into that with $`\chi _2`$. Since $`f^{(1)}`$ depends on the sign of $`p_x`$ only, being independent of the direction of $`𝐩`$ within each hemisphere $`p_x>0`$ or $`p_x<0`$,

$$J_\pm =\tau ^1\nu _\pm \nu _{}\left(f_\pm ^{(1)}f_{}^{(1)}\right).$$
(4)
The collision integral Eq. (4) is only nonvanishing if $`\nu _+`$ and $`\nu _{}`$ are nonzero simultaneously. It is localized at distances $`sd`$: The averages $`\nu `$ decay at $`sa`$ away from the contact, being proportional to the solid angle at which the orifice is visible from the position point. Since $`\nu `$ is constant within the orifice, $`\mathrm{}^1𝑑s\nu _{}=\gamma \nu _{}(0)`$ where $`\gamma d/\mathrm{}`$ is a contact-dependent geometric factor; we assume a constant $`\gamma `$ for simplicity. Finally, putting $`\varphi =\chi _2\chi _1`$ we find

$$\nu _\pm (0)\left[\pm \frac{1}{2}\frac{\varphi }{t}\frac{f_\pm ^{(1)}}{ϵ}+\frac{1}{\sqrt{|\mathrm{\Delta }|^2ϵ^2}}\frac{f_\pm ^{(1)}}{t}\right]=\gamma \nu _\pm (0)\nu _{}(0)\left(f_\pm ^{(1)}f_{}^{(1)}\right).$$
(5)
## III Spectral flow in a ballistic contact

In the collisionless limit Kulik ,

$$\nu _\pm \left(0\right)=\pi \sqrt{|\mathrm{\Delta }|^2ϵ^2}\delta \left(ϵϵ_\pm \right).$$
(6)
Since $`ϵ_\pm /\varphi =\pm \frac{1}{2}|\mathrm{\Delta }|\mathrm{sin}(\varphi /2)`$ we can write Eq. (5) as $`\nu _\pm (0)df^{(1)}/dt=0`$ where the total time derivative of the distribution function is taken along the spectrum Eq. (1). The total time derivative vanishes: the distribution is constant in time and thus it is independent of energy. To find it we note that, as $`\varphi `$ varies from $`0`$ to $`2\pi `$, the distribution of particles with $`p_x>0`$ remains the same as it was at $`t=0,\varphi =0`$ for an energy $`ϵ=\left|\mathrm{\Delta }\right|`$, i.e., $`f^{(1)}=f_\mathrm{\Delta }`$ where $`f_\mathrm{\Delta }=\mathrm{tanh}\left(|\mathrm{\Delta }|/2T\right)`$. For particles with $`p_x<0`$, the distribution remains as it was at $`t=0,\varphi =0`$, and $`ϵ=\left|\mathrm{\Delta }\right|`$, i.e., $`f^{(1)}=f_\mathrm{\Delta }`$. If $`\varphi `$ decreases from $`2\pi `$ to $`0`$, the signs are opposite,

$$f_\pm ^{(1)}=f_\mathrm{\Delta }\mathrm{sign}(\dot{\varphi }).$$
(7)
Excitations are accelerated by the applied voltage undergoing multiple Andreev reflections Zaikin ; Cuevas ; AverinBardas1  before their energies grow by $`2|\mathrm{\Delta }|`$ to let them escape from the potential well.

The same process can also be viewed differently. Under a small voltage bias $`eV|\mathrm{\Delta }|`$, quasiparticles are accelerated along the spectrum Eq. (1) with $`\varphi =2eVt`$. As the phase varies from $`0`$ to $`2\pi `$, the quasiparticle energy changes from $`|\mathrm{\Delta }|`$ to $`\pm |\mathrm{\Delta }|`$. A quasiparticle that has been captured below the gap edge emerges above the gap after one period of Josephson oscillations (see Fig. 1). Another particle gets captured and then released during the next period, and so on AverinBardas2 . This is known as the spectral flow. It plays a crucial part in vortex dynamics KopVol/specflow ; Stone . In a context of weak links, it was first discussed in Ref. Volovik . In the absence of collisions, the distribution of excitations localized near the contact Eq. (7) is determined by distribution of excitations at the gap edges, which are delocalized and are thus nearly in equilibrium due to the relaxation in the bulk. However, emerging quasiparticles have a “wrong” distribution with respect to the “native” ones; they relax producing dissipation.

The current through the contact with the Sharvin conductance $`R_N^1=e^2N_Fv_FS/2`$ is

$$I=\frac{1}{eR_N}\frac{dϵ}{2}\left[\nu _+f_+^{(1)}\nu _{}f_{}^{(1)}\right].$$
(8)
where $`Sa^2`$ is the area of the contact. With Eq. (7) the current due to Andreev states is

$$|e|R_NI=\pi |\mathrm{\Delta }|\left|\mathrm{sin}(\varphi /2)\right|f_\mathrm{\Delta }\mathrm{sign}(V).$$
(9)
Excitations with energies $`|ϵ|>|\mathrm{\Delta }|`$ give a much smaller current, which can be neglected. Equation (9) agrees with the results of Refs. Zaikin ; Cuevas ; AverinBardas1 .

## IV  Spectral flow and impurities

A small amount of impurities in the area of the contact lifts the degeneracy of Eq. (1) and modifies the spectrum in a way that an impurity band appears near the point $`\varphi =\pi `$ where the two spectral branches cross Kopnin/pointcont . In the presence of impurities the DOS is $`\nu _{ϵ\pm }(0)=\left(E_{}/2\gamma \right)\sqrt{\left(4\gamma /E_+E_{}\right)1}`$ where $`E_\pm =ϵ/|\mathrm{\Delta }|\pm \mathrm{cos}(\varphi /2)`$. Both $`\nu _+`$ and $`\nu _{}`$ are nonzero near the lines $`E_\pm =0`$ within the impurity band $`0<E_+E_{}<4\gamma `$ shown in Fig. 1. To simplify the partial differential equation (5) we note that, in the region far from the crossing point, $`|ϵ||\mathrm{\Delta }|\sqrt{\gamma }`$, $`|\varphi \pi |\sqrt{\gamma }`$, the ratio $`\nu _{}/\nu _+=E_+/E_{}`$ is of the order of $`\gamma |\mathrm{\Delta }|^2/ϵ^21`$ near the line $`E_+=0`$, and it is $`\nu _{}/\nu _+ϵ^2/\gamma |\mathrm{\Delta }|^21`$ near $`E_{}=0`$. Therefore, far from the crossing point, $`\left|\varphi \pi \right|\beta `$ where $`\sqrt{\gamma }\beta \pi `$, the DOS has the form of Eq. (6), and the l.h.s. of Eq. (5) reduces to the full derivative $`df_\pm ^{(1)}/dt`$ along the respective lines $`E_\pm =0`$. Since the collision integral also vanishes for $`\left|\varphi \pi \right|\beta `$, the distribution is constant in time far from the crossing point, in the same way as it was without impurities.

Equation (5) can be solved by perturbations in the limit of small voltages. The zero-order term corresponds to vanishing of the l.h.s. and results in the condition $`f_+^{(1)}=f_{}^{(1)}`$ that holds through the crossing point for a given energy. We find that the distribution functions at the four lines in Fig. 1 entering the crossing region satisfy $`f_\pm ^{(1)}(\pi +\beta )=f_{}^{(1)}(\pi \beta )=\pm f_\mathrm{\Delta }`$. To the zero approximation, particles emerging in the continuum with an energy $`ϵ=\pm |\mathrm{\Delta }|`$ have the distribution $`\pm f_\mathrm{\Delta }`$ that corresponds to the local equilibrium at the gap edges

$$f_\pm ^{(\mathrm{eq})}=\mathrm{sign}[\mathrm{cos}(\varphi /2)]f_\mathrm{\Delta }.$$
(10)
In this approximation, no relaxation in the continuum occurs. The deviation $`f_{1\pm }`$ from the local equilibrium is proportional to the time derivative, i.e., to the voltage. We write $`f_\pm ^{(1)}=f_\pm ^{(\mathrm{eq})}+f_{1\pm }`$. The l.h.s. of Eq. (5) can be estimated as $`\nu eV\left(f_\mathrm{\Delta }/ϵ^{}\right)`$ where $`ϵ^{}|\mathrm{\Delta }|\sqrt{\gamma }`$ is the characteristic scale of energy near the crossing point. Using $`\nu 1/\sqrt{\gamma }`$ we find $`f_{1\pm }\left(eV/|\mathrm{\Delta }|\gamma \right)f_\mathrm{\Delta }`$. The deviation $`f_{1\pm }`$ is thus linear in $`V`$ for low voltages.

In general, $`f_\pm ^{(1)}`$ can be found by matching the four constants $`f_{}^{(1)}(\pi \pm \beta )`$ while integrating Eq. (5) across $`\varphi =\pi `$. We note that for a contact with a barrier of reflectivity $`R=1D`$, the excitation spectrum is given by Eq. (2) with a gap $`\sqrt{R}|\mathrm{\Delta }|`$ that separates states with positive energies from those with negative energies. As a result, the distribution for excitations with $`ϵ>0`$ is decoupled from that for $`ϵ<0`$ unless the Zener tunneling mixes them AverinBardas1 ; AverinBardas3 . On the contrary, in macroscopic contacts with impurities, the states near $`\varphi =\pi `$ form a continuum where the regions $`ϵ>0`$ and $`ϵ<0`$ are connected along the lines of a nonzero DOS passing through $`\varphi =\pi `$, see Fig. 1. According to Eq. (5), the distribution of excitations is thus continuous within the entire region of nonzero DOS. Therefore, the distributions on all four lines entering the crossing region in Fig. 1 are coupled to each other.

Since we are interested the asymptotic values of $`f_\pm ^{(1)}`$ for $`\varphi =\pi \beta `$ and $`\varphi =\pi +\beta `$, we neglect the fine structure of $`\nu `$ and approximate it by Eq. (6). The collision integral is then nonzero only at the crossing point of $`ϵ_+`$ and $`ϵ_{}`$ defined by Eq. (1). Equation (5) yields

$$\frac{df_\pm ^{(1)}}{dt}=\pi \gamma |\mathrm{\Delta }|^2\mathrm{\Gamma }\delta \left(ϵϵ_{}\right).$$
(11)
$`\mathrm{\Gamma }`$ is a combination $`f_+^{(1)}(ϵ,\varphi )f_{}^{(1)}(ϵ,\varphi )`$ taken in the vicinity of $`ϵ=0`$, $`\varphi =\pi `$, which couples the four values $`f_\pm ^{(1)}(\pi \pm \beta )`$. We shall see that $`\mathrm{\Gamma }`$ varies slowly near $`\varphi =\pi `$ thus we approximate

$$\mathrm{\Gamma }=\frac{1}{2}\left[f_+^{(1)}(\pi +\beta )+f_+^{(1)}(\pi \beta )f_{}^{(1)}(\pi +\beta )f_{}^{(1)}(\pi \beta )\right].$$

We integrate Eq. (11) over time along the dependence $`E_+=0`$ for particles with $`p_x>0`$ or along $`E_{}=0`$ for $`p_x<0`$ and obtain

$$f_\pm ^{(1)}(\pi +\beta )f_\pm ^{(1)}(\pi \beta )=\alpha \mathrm{\Gamma }$$
(12)
where $`\alpha =\pi |\mathrm{\Delta }|\gamma /2eV`$. Adding the two Eqs. (12) we find that the quantity

$$A(\beta )=f_+^{(1)}(\pi +\beta )f_{}^{(1)}(\pi \beta )$$
(13)
satisfies $`A(\beta )=A(\beta )A`$, i.e., it is a slowly varying function near $`\beta =0`$. This proves that $`\mathrm{\Gamma }=[A(\beta )+A(\beta )]/2`$ can indeed be treated as a constant, $`\mathrm{\Gamma }=A`$.

Consider first the case $`\dot{\varphi }=2eV>0`$. The initial condition at $`t=0`$, $`\varphi =0`$ is the equilibrium distribution at an energy $`ϵ=|\mathrm{\Delta }|`$, i.e., $`f_\pm ^{(1)}=f_\mathrm{\Delta }`$. Since the distribution is constant outside the crossing region, we have also $`f_\pm ^{(1)}(\pi \beta )=f_\mathrm{\Delta }`$. To find the functions $`f_\pm ^{(1)}(\pi +\beta )`$ we use Eqs. (12, 13) whence $`A=2f_\mathrm{\Delta }/\left(1+\alpha \right)`$ and

$$f_\pm ^{(1)}(\pi +\beta )=\left[\left(1\alpha \right)/\left(1+\alpha \right)\right]f_\mathrm{\Delta }.$$

For weak scattering $`\alpha 1`$, the distribution remains unchanged as compared to its initial value, $`f_\pm ^{(1)}(\pi +\beta )=f_\pm ^{(1)}(\pi \beta )=f_\mathrm{\Delta }`$. Quasiparticles relax towards the local equilibrium $`\pm f_\mathrm{\Delta }`$ as they emerge in the continuum from inside the gap at an energy $`\pm |\mathrm{\Delta }|`$. If scattering is strong $`\alpha 1`$, the distribution changes considerably at the crossing point $`\varphi =\pi `$ due to exchange between particles moving from the opposite sides of the contact. Emerging particles with an energy $`ϵ=\pm |\mathrm{\Delta }|`$ have a distribution close to the “local equilibrium” Eq. (10) such that no relaxation occurs. This is not the true equilibrium for $`|ϵ|<|\mathrm{\Delta }|`$: the latter can only be established when the spectral flow rate $`eV`$ is much smaller than the inelastic rate. The local equilibrium is achieved through a very effective branch mixing due to the momentum exchange on impurities when particles pass through the crossing region of the spectral branches. The deviation from the local equilibrium is zero as long as $`\varphi <\pi `$ while it appears for $`\varphi >\pi `$:

$$f_{1\pm }f_\pm ^{(1)}(\pi +\beta )\left[\pm f_\mathrm{\Delta }\right]=\left[2/\left(1+\alpha \right)\right]f_\mathrm{\Delta }.$$
(14)
This agrees with our earlier perturbative estimate for the low-voltage limit.

Calculating the current for $`\varphi <\pi `$ from Eq. (8) we find $`I=I_s`$ where the “supercurrent” is defined through the “equilibrium distribution” Eq. (10):

$$eR_NI_s=\pi |\mathrm{\Delta }|\mathrm{sin}(\varphi /2)f_\mathrm{\Delta }\mathrm{sign}\left[\mathrm{cos}(\varphi /2)\right].$$

This expression differs from the well-known result for a steady-state supercurrent Kulik . The reason is that Eq. (10) is not the true equilibrium, see the discussion above. For $`\varphi >\pi `$ we have $`I=I_s+I_n`$; the “normal” current contains the deviation from equilibrium Eq. (14). Considering also the case $`\dot{\varphi }=2eV<0`$ we find for the normal current

$$|e|R_NI_n=2\pi |\mathrm{\Delta }|\left(1+|\alpha |\right)^1\left|\mathrm{sin}(\varphi /2)\right|f_\mathrm{\Delta }\mathrm{sign}(V).$$

As in Ref. AverinBardas3 , it exists only when $`\pi <\varphi <2\pi `$ for $`\dot{\varphi }>0`$ or when $`0<\varphi <\pi `$ for $`\dot{\varphi }<0`$.

For a voltage-biased contact, the supercurrent vanishes after averaging over time:

$$|e|R_N\overline{I}=2|\mathrm{\Delta }|\left(1+|\alpha |\right)^1f_\mathrm{\Delta }\mathrm{sign}(V).$$
(15)
For high voltages, $`|\alpha ||\mathrm{\Delta }|\gamma /|eV|1`$, the MAR result Eq. (9) is reproduced. For low voltages $`|\mathrm{\Delta }|\gamma |eV|\tau _ϵ^1`$, the current is linear, $`\overline{I}=GV`$ where $`G=(4/\pi \gamma R_N)f_\mathrm{\Delta }`$. This conductance is much larger than the contribution $`R_N^1`$ from the states with $`|ϵ|>|\mathrm{\Delta }|`$. To our knowledge, this low-voltage region has not been studied experimentally though the contacts with suitable parameters are now under intensive investigation contacts .

## V Discussion

The linearity of the I-V curve for low voltages disagrees with a square-root behavior predicted in AverinBardas3  using the model based on Zener transitions within independent channels averaged over the universal distribution of reflection coefficients Dorokhov  valid for diffusive conductors. This distribution cannot be applied directly for nearly ballistic contacts. The proper distribution for contacts with a small disorder was found Kopnin/pointcont  to be $`P(R)=(2\pi \gamma )^1\sqrt{\gamma R}/\sqrt{R}`$. Its denominator does also contain $`\sqrt{R}`$ like the universal distribution, which, at the first glance, should have also resulted in a square-root current-voltage dependence. However, the kinetic-equation approach applicable for a multichannel contact leads to a different result. The reason is that, in a macroscopic contact, the quantum channels are very close to each other in energies with the level spacings inversely proportional to the number of channels $`N_{\mathrm{ch}}`$. Thus, in addition to the transitions within each channel (which are the only ones considered in the independent-channel model), time dependent processes excite even more effectively transitions between different channels. As a result, the distribution function is spread over the entire energy interval of nonzero DOS or over the range $`R\gamma `$ of reflection coefficients. The role of each particular channel within the range $`\delta R`$ should thus be weighted with $`\delta R/\gamma `$. If we now average the low-voltage Zener processes having the rate $`\mathrm{exp}(|\mathrm{\Delta }|R/eV)`$ with $`P(R)`$ we would obtain the essential reflection coefficients $`ReV/|\mathrm{\Delta }|`$. Thus the fraction of essential channels contributing to the current in multichannel contacts is decreased by a factor $`eV/|\mathrm{\Delta }|\gamma `$ as compared to the independent channel model. This yields a much smaller current compared to what is predicted by the kinetic-equation approach. Therefore, the Zener processes can be neglected in favor of the momentum relaxation on impurities. To summarize, the model based on Zener transitions within independent channels does not seem to be suitable for kinetic processes in multichannel contacts. We believe, however, that calculations using the Landauer approach with the full account of interchannel transitions followed by the appropriate averaging would be equivalent to the kinetic equation technique.

We thank D. Averin, M. Feigel’man, A. Larkin, and G. Volovik for valuable discussions. This work was supported by the U.S. Department of Energy Office of Science through contract No. W-31-109-ENG-38 and by the Russian Foundation for Basic Research.


text:
# A Non-parametric Analysis of the CMB Power Spectrum

## 1. Introduction

There has been growing evidence for the existence of peaks and valleys in the temperature power spectrum of the CMB. From a theoretical standpoint, such features are a direct result of the physics in the primordial photon-electron plasma, predicted by gravitational instability models of structure formation (Peebles & Yu 1970; Hu and Sugiyama 1995). These features are important for constraining the cosmology of our Universe. For instance, in many models, the ratio of the height of the first peak to the second peak is dependent on the spectral tilt, $`n_s`$ and the baryon fraction, $`\mathrm{\Omega }_{baryons}/\mathrm{\Omega }_{matter}`$. The ratio of the third peak to the second peak is dependent on $`\mathrm{\Omega }_{matter}h^2`$ and $`n_s`$ (see Hu et al. 2001 for further discussion).

Most often in the literature, the CMB power spectra are fit to a suite of cosmological models (Tegmark et al., 1999,2000; Jaffe et al. 2001). These physical models are well-motivated and sophisticated, but they contain many free parameters (*e.g.*, eleven in the work of Wang, Tegmark, and Zaldarriaga 2001– WTZ), some of which are unknown (ionization depth, contribution from gravity waves) or degenerate (e.g. see Efstathiou 2001). Typically, some sort of likelihood analysis is performed to determine which cosmological model best fits the data.

There is however, another approach: place constraints on the features of the power spectrum and use these features to determine the cosmological parameters. The assumptions here are that the peaks and valleys are best described by the broad range of cosmological models (as in Hu et al.) or by parabolas or some other chosen function (as in Knox & Page 2000 and de Bernardis et al. 2001). A potential problem in all of these approaches is that it is difficult to get valid statistical confidence intervals (see Section 2). There is also the concern that the fitted features may be artifacts from the multitude of assumptions.

In this paper, we take what may be considered a more conservative approach: we make no assumptions whatsoever about the true underlying function. Our new statistical technique is non-parametric and allows for valid confidence intervals to be measured for peak characteristics. One theme of our work is that confidence intervals for any quantity of interest can be extracted from a confidence set for the unknown spectrum. These techniques for fitting and inference are applicable to a wide variety of astrophysical data-analysis problems.

## 2. Overview of Nonparametric Analysis

In general, non-parametric statistical methods estimate functions without imposing a finite-dimensional parametric form. The resulting estimates are obtained by carefully smoothing the data to balance bias and variance. See Hastie and Tibshirani (1990) for details and examples.

The CMB data, after suitable preprocessing, take the form $`(X_1,Y_1),\mathrm{},(X_n,Y_n)`$ where $`X_i`$ is the multipole moment (usually denoted $`\mathrm{}`$) ordered according to increasing $`X`$, and $`Y_i`$ is the estimated power spectrum at $`X_i`$ (usually denoted $`C_{\mathrm{}}`$ with some constants). Let $`f(X_i)`$ be the true power spectrum at $`X_i`$. Then,

$$Y_i=f(X_i)+ϵ_i$$
(1)
where $`ϵ_i`$ is the error in estimating $`f(X_i)`$. We require that the $`ϵ_i`$ are uncorrelated, zero-mean Gaussians. Therefore, we take the uncorrelated statistical errors as given by the experiments. There are additional, correlated noise terms in all of the experiments due to calibration, beam width and pointing uncertainties, which were not used in our analysis. The magnitude of these correlated errors is typically 5% - 20%, although the Boomerang effective beam uncertainties can be higher for $`\mathrm{}>600`$. The method we use in this work can be generalized to include correlated errors (Genovese et al., in preparation). We assume that $`f^2`$ is integrable; otherwise, we make no further assumptions. Our non-parametric technique yields a vector $`\widehat{𝐟}`$ that estimates the vector $`𝐟=(f(X_1),\mathrm{},f(X_n))`$ of the true spectrum’s values at the $`X_i`$s (see Figure 1 top).

After we perform the fit, we need to quantify the uncertainty to make any inferences. We begin by constructing a set of vectors $`C_n`$ from the data that traps the true power spectrum, $`𝐟`$, with probability $`1\alpha `$, where $`1\alpha `$ is a pre-specified confidence level. With $`C_n`$ we can derive a confidence interval for any quantity of interest. Consider, for example, the number of peaks. For any vector $`𝐟C_n`$, define $`P(𝐟)`$ to be the number of peaks in $`𝐟`$. Then, because $`C_n`$ contains the true spectrum with confidence $`1\alpha `$, the range on the measured peaks is $`(\mathrm{min}_{𝐟C_n}P(𝐟),\mathrm{max}_{𝐟C_n}P(𝐟))`$. We can find confidence intervals for the locations and heights of peaks in a similar way. There are important advantages for these confidence measurements over “standard” $`\chi ^2`$ or even Bayesian techniques, which we discuss further in the next section.

## 3. Methodology

Refer to the model in equation (1). Let the functions $`\varphi _1,\varphi _2,\mathrm{}`$ be an orthonormal basis over the range of the observed $`X_i`$s. The choice of basis is somewhat important. For instance, if we were fitting to a galaxy spectrum, with highly peaked emission lines on top of a smooth, broad continuum, a wavelet basis would allow for the simultaneous fitting of broad and narrow features. On the other hand, a cosine basis would require large amplitude high frequency terms to match the emission lines. This might cause the continuum fit to be wiggly. In our work, we use a discrete cosine basis, $`\varphi _1(x)1`$, $`\varphi _2(x)=\sqrt{2}cos(\pi x)`$, $`\varphi _3(x)=\sqrt{2}cos(2\pi x),\mathrm{}`$, since it has well determined properties for confidence limits. For simplicity in the derivation, we assume that the variance (second moment about the mean) of each $`ϵ_i`$ is the same value $`\sigma ^2`$. This is not necessary in practice, nor was it assumed in the full analysis.

Any square integrable function $`f`$ can be expanded as $`f(x)=_{j=1}^{\mathrm{}}\beta _j\varphi _j(x)`$. Estimating $`f`$ amounts to estimating the $`\beta _j`$’s. If we have chosen a good basis for representing $`f`$, the higher-order terms in this series will tend to decay rapidly. Hence, we can approximate the infinite sum with finite sum $`f(x)_{j=1}^N\beta _j\varphi _j(x)`$. So, to estimate the underlying function, we need to find the $`\beta `$s.

Let $`Z_j=N^{1/2}_{i=1}^NY_i\varphi _j(X_i)`$, for $`j=1,\mathrm{},N`$. Based on the theory in Beran (2000), we take $`N=n`$. This choice ensures that the estimate of $`f`$ is optimal and that the resulting confidence intervals remain valid. It can be shown that each statistic $`Z_j`$ has approximately a Gaussian distribution with mean $`\theta _j=\sqrt{n}\beta _j`$ and variance $`\sigma ^2`$. This re-parameterization means that we need to estimate $`\theta =(\theta _1,\mathrm{},\theta _n)`$. Given an estimate $`\widehat{\theta }`$ of $`\theta `$, we can estimate $`f(X_i)`$ via $`\widehat{f}(X_i)=(1/\sqrt{n})_{j=1}^n\widehat{\theta }_j\varphi _j(X_i)`$.

We could use $`Z_j`$ as an estimate of $`\theta _j`$, but this yields a poor estimate of $`f`$ because it is too variable. A smoother estimate can be obtained by damping the higher frequency terms in the expansion. In statistics, this is called a “shrinkage estimator”. We consider shrinkage estimators of the form $`\widehat{\theta }=(\gamma _1Z_1,\gamma _2Z_2,\mathrm{},\gamma _nZ_n)`$, where $`1\gamma _1\gamma _2\mathrm{}\gamma _n0`$ are called the shrinkage coefficients. The smaller $`\gamma _j`$, the smaller the contribution of $`\varphi _j`$ in the estimating expansion for $`f`$. With the cosine basis, for example, such shrinkage estimators damp down the contribution from high-frequency terms.

Every choice of $`\gamma =(\gamma _1,\mathrm{},\gamma _n)`$ gives an estimate $`\widehat{\theta }^\gamma `$ which then yields an estimate $`\widehat{f}^\gamma `$ of the function $`f`$. We would like to choose $`\gamma `$ to minimize the mean squared error, $`\mathrm{MSE}(\gamma )=(f(x)\widehat{f}^\gamma (x))^2𝑑x`$. Unfortunately $`\mathrm{MSE}(\gamma )`$ is unknown, because it depends on the true $`f`$, but it can be estimated by Stein’s Unbiased Risk Estimator (Stein 1981):

$$\widehat{\mathrm{MSE}}(\gamma )=\underset{j}{}\left[\sigma ^2\gamma _j^2+\left(Z_j^2\sigma ^2\right)(1\gamma _j)^2\right].$$
(2)
We use the Pooled Adjacent Values (PAV) algorithm (Robertson, Wright, and Dykstra, 1988) to minimize $`\widehat{\mathrm{MSE}}(\gamma )`$ as a function of $`\gamma `$ while maintaining the ordering constraint on $`\gamma `$. The minimizer is denoted $`\widehat{\gamma }`$ and the final estimate is therefore $`\widehat{\theta }=(\widehat{\gamma }_1Z_1,\widehat{\gamma }_2Z_2,\mathrm{},\widehat{\gamma }_nZ_n)`$.

Next we construct a confidence set, $`C_n`$ for the vector of function values at the observed data, $`𝐟_n=(f(X_1),f(X_2),\mathrm{},f(X_n))`$. Throughout the paper we have said that the confidence sets are “valid.” Formally, what this means is that for any $`c>0`$,

$$\underset{n\mathrm{}}{lim}\underset{𝐟_nc}{sup}\left|\mathrm{Pr}(𝐟_nC_n)(1\alpha )\right|0$$
(3)
as $`n\mathrm{}`$, where the operation $`a`$ denotes $`\sqrt{n^1_ia_i^2}`$. This means that for large $`n`$, the confidence set traps the values of the true function with probability very close to $`1\alpha `$. The confidence set is an ellipse of the form

$$C_n=\{\theta :n^1\underset{j}{}(\theta _j\widehat{\theta }_j)^2\widehat{MSE}(\widehat{\gamma })+n^{1/2}\tau z_\alpha \}$$
(4)
where $`z_\alpha `$ is the number that has probability $`\alpha `$ to the right under a standard Gaussian. For instance, if we choose 95% confidence ($`\alpha =0.05`$), then $`z_\alpha =1.645`$, while for 67% confidence ($`\alpha =0.33`$), the confidence “radius” is smaller with $`z_\alpha =0.44`$. Here, $`\tau `$ is defined similar to Beran (2000):

$$\tau ^2=4\sigma ^2\underset{j}{}(Z_j^2\sigma ^2)(1\widehat{\gamma }_j)^2+2\sigma ^4\underset{j}{}(2\widehat{\gamma }_j1)^2.$$
(5)
We can in turn write the confidence set for $`𝐟`$ as

$$𝒟_n=\{𝐟:𝐟=n^{1/2}\underset{j=1}{\overset{n}{}}\theta _j\varphi _j\text{ and }\theta 𝒞_n\}.$$
(6)
Here, the true power spectrum $`𝐟=(f(X_1),\mathrm{},f(X_n))`$ is estimated at each of the original data points. Prior information about $`𝐟`$ of the form $`𝐟𝒫_n`$ allows us to replace $`𝒟_n`$ with $`𝒟_n𝒫_n`$ while maintaining a $`1\alpha `$ confidence level. In particular, we take $`𝒫_n`$ to be the set of vectors $`𝐟`$ corresponding to spectra with zero to three peaks for $`\mathrm{}1100`$. From the confidence set $`𝒟_n𝒫_n`$, we can derive confidence sets for any interesting feature of $`𝐟`$.

A key point is that the resulting intervals on any measured quantity are *simultaneously valid*, meaning that all of the intervals contain the corresponding true quantity with probability $`1\alpha `$. In contrast, deriving $`1\alpha `$ confidence intervals from a collection of individual chi-squares does not obtain $`1\alpha `$ simultaneous coverage, but often substantially lower coverage. For example, a common technique to determine the 95% confidence range for a specific parameter is to “marginalize” over the other parameters (see Tegmark et al. 1999,2000). Such a technique will provide full coverage for that one parameter. However, when parameter ranges are combined, the confidence is lower than 95%. Bayesian intervals derived from a posterior distribution suffer from a similar problem in that the long-run frequency that the interval contains the true quantity may be much less than $`1\alpha `$.

## 4. Results

Figure 1 shows the combined data from BOOMERANG, MAXIMA, and DASI (Halverson et al. 2001; Lee et al. 2001; and Netterfield et al. 2001). The bottom panel compares our non-parametric fit to the the fit of WTZ who use more experiments than the three examined here. Recall, our fit requires no assumptions about the data or the underlying cosmology. The WTZ fit, on the other hand, requires an 11 dimensional parameter space, with numerous prior assumptions placed on those parameters (for example, the Hubble constant is constrained to $`72\pm 8`$ km s-1Mpc-1). The agreement between the two fits is very good, considering the difference in methodology.

The power of this technique lies in the ability to make quantitative statements about the true function with some specified confidence. First, we checked the 95% confidence set and find that every function within this set has at least one peak. Specifically, we set $`\alpha =0.05`$ and determined the “radius” of this confidence ellipse (ie. the right side of the inequality in Eq 4). We then searched all possible functions with zero to three peaks over the specified range in $`\mathrm{}`$ to see if the condition in Eq 4 was met. Our definition of one peak requires the sorted data (according to increasing $`\mathrm{}`$) to have a section with increasing temperature followed by a section with decreasing temperature. For $`\alpha =0.05`$, no zero-peaked functions met this condition. Our set of zero-peaked functions includes those with constant temperature as well as those with either increasing temperature or decreasing temperature, but not both. We rule out functional forms with zero peaks at the level 95% confidence level. We also compared specifically against the best fit constant function for the power spectrum. We rule out a flat temperature spectrum (based on the weighted average of the data) at $`>8\sigma `$ confidence. We perform the same analysis at the 68% confidence level (e.g. $`\alpha =0.32`$ in Eq 4). The 68% confidence set rules out all single-peak functions. So at the one sigma level, we have found at least two peaks in the data. Only in the 40% ($`\alpha =0.6`$) confidence set can we rule out two-peak functions. From our confidence sets, the data supports two peaks out of three peaks in the CMB power spectrum (for $`0\mathrm{}1100`$).

We calculate ranges for the peak heights and locations in Table 1. Figure 2 shows the non-parametric fit with confidence intervals for peak heights and locations at 95% confidence. Finally, we computed confidence intervals for the ratios of successive peaks under a three-peak model. The 95% confidence interval for the ratio of the height of the first peak to the height of the second peak is $`(1.06,4.27)`$. The 95% confidence interval for the ratio of the height of the second peak to the height of the third peak is $`(0.41,2.5)`$. This rules out equal heights for the first two peaks at the 95% level. These results are consistent with Hu et al. (2001) who find much stronger constraints on the height-height ratios by fitting to cosmological models.

## 5. Discussion and Conclusions

We present an application of a new and powerful non-parametric technique to CMB temperature data. Past approaches were based on complicated cosmological models or parameterized forms. There is superb visual agreement between the non-parametric fit and the best fitting cosmological model. Quantitatively, we provide constraints on the peak locations, heights, and height ratios of the power spectrum. These constraints can be used to place corresponding limits on the cosmological parameters that they describe. For instance, Hu et al. (2001) derive relationships which could in principle, be used for this purpose.

At the $`2\sigma `$ confidence level, we find at least one peak in the current CMB power spectra data, while at the $`1\sigma `$ level, we find two or more peaks. Only for a very low confidence, 40%, can we rule out two peak functions. Therefore, the data do not yet show the three expected peaks for $`\mathrm{}1100`$ (in the three CMB datasets examined here). There are two explanations for this: the model is right, but there is insufficient precision in the current data, or the model is wrong. If in fact the errors on the current measurements are simply too large, then these standard errors would have to be reduced to one-third of their current values to rule out a two-peak spectrum at the 95% confidence level. This suggests a range for the maximum required errors for future CMB experiments (via MAP and Planck) to “discover” three peaks in the CMB spectrum.

We point out that the lack of assumptions used to arrive at our best fit is conservative. On the other hand, results from fitting assumed cosmological models are optimistic, since those models all have a multi-peaked spectra (e.g. Hu et al. 2001). While the physical underpinnings for cosmological models are well founded, the last 50 years (or even five) have seen radical changes in those models which best fit the data. Therefore, a method to describe the CMB that is “cosmology free” has scientific value. Finally, we note that the methods described here can be applied to the many astrophysical problems that are not well suited for standard parametric techniques.

This work was done in collaboration with the Pittsburgh Computational Astrostatistics Group (`www.picagroup.org`). The authors would like to thank the referee for suggestions which improved the readability and usefulness of this work. RCN, LW, and CG were partially supported by NSF KDI grant DMS-9873442.

During the refereeing process of our paper, two related papers came to our attention (Durrer, Novosyadlyj & Apunevych 2001; Douspis & Ferreira 2001). These papers perform a model–independent measurement of the CMB power spectrum but they are not non-parametric estimates of the CMB acoustic peaks, as discussed herein, since they use phenomenological models to describe the underlying power spectrum. It is interesting to note however, that all three analyses find low statistical significances for the detection of the second and third peaks. We await higher precision measurements of the CMB power spectrum to secure the detection, location, shape and amplitude of these peaks.


text:
# Dusty Acoustic Turbulence in the Nuclear Disks of two LINER Galaxies NGC 4450 and NGC 4736

## 1 Introduction

Dust spirals in the inner kpc regions of galaxies reveal a source of compression that can affect the angular momentum distribution of the gas and possibly drive accretion to an AGN. Barred galaxies with an inner Lindblad resonance tend to have two long and symmetric dust spirals near the resonance that are a continuation of the leading-edge dust lanes in the bar (Athanassoula 1992). Many barred galaxies have ILR rings too (Buta & Crocker 1993; Pérez-Ramírez et al. 2000; Knapen et al. 2000). Non-barred galaxies, galaxies without an ILR (e.g., late Hubble types), and regions of barred galaxies inside their ILRs can have more irregular dust spirals.

Whether the nuclear spirals are regular or irregular, their presence in optical images suggests a density variation that is at least a factor of 2 and therefore likely to involve shocks. These shocks are oblique for azimuthal flows, so the gas will experience a torque when it enters a spiral and an opposite torque when it leaves. In cases where the gas moves faster than the spirals, the net torque is negative and the gas loses angular momentum and energy at the shock, causing it to spiral inward. If the spirals move faster than the gas, as might be the case outside a fast nuclear bar, then the gas gains angular momentum and may gradually move out.

We are interested in the source of compression for the irregular nuclear dust spirals that appear in non-barred galaxies or inside the ILRs of some barred galaxies. We have proposed that some of these are caused by random sonic noise that amplifies weakly as it propagates toward the center (Elmegreen et al. 1998; Montenegro, Yuan, & Elmegreen 1999; Englmaier & Shlosman 2000). A signature of this process is an irregularity of structure, a wide range of pitch angles reflecting different times of origin and different wave propagation directions, a tendency to fill the volume available with a spiral separation comparable to the epicyclic radius for motions at the sound speed, and a general trend of increasing density with decreasing radius for each spiral. Because the velocities inside these spirals cannot yet be measured, there is no way to be sure they are driving inflow. Nevertheless, if any of the spirals come close to the nucleus, then such accretion would seem to be likely.

Nuclear dust spirals and clouds have been studied previously using Hubble Space Telescope (HST) data. Van Dokkum & Franx (1995) found them on WFPC V-band images of early type galaxies. Malkan, Gorjian, & Tam (1998) did a large snapshot survey of active galaxies and classified the nuclear morphology including the dust. Elmegreen et al. (1998) noted dust spirals in WFPC2 images of the interacting galaxy NGC 2207 and proposed they might drive accretion to the nucleus. Regan & Mulchaey (1999) found nuclear dust spirals in 6 AGN galaxies using WFPC2 and NICMOS images, and Martini & Pogge (1999) found them with HST data in 20 Seyferts; both studies concluded that the spirals could drive accretion. Ferruit, Wilson & Mulchaey (2000) studied 12 early type Seyfert galaxies in more detail and noted the dust spirals too. Tomita et al. (2000) used HST archival images to study dust features in E and S0 galaxies, but did not comment on spirals specifically. Tran et al. (2001) studied nuclear dust in elliptical galaxies.

Here we measure the extinction and structural properties of irregular dust spirals in the central kpc of two LINER galaxies, NGC 4450 and NGC 4736. Their Hubble types are about the same, SA(s)ab and RSA(r)ab (de Vaucouleurs et al. 1991), and their distances are taken to be 16.8 and 4.3 Mpc (Tully 1988).

## 2 Observations

Hubble Space Telescope preview images of the galaxies NGC 4450 and NGC 4736 are shown in Figure 1. These images were obtained from the HST archive web site and are degraded in resolution by a factor of 2 from the full HST data; they are from proposals 5375 by V. Rubin and 5741 by J. Westphal, respectively. For our analysis of dust structure, we used the original Planetary Camera (PC) images with $`0.0455^{\prime \prime }`$ per pixel, corresponding to 3.7 pc and 0.95 pc per pixel, respectively. The PC fields of view cover 3000 pc and 760 pc. Multiple images in the V passband (F555W filter) were combined to remove cosmic ray tracks and bad pixels. Figure 1 shows the full WFPC2 image to illustrate the outer disk morphologies, their position angles, and the orientation of the HST images relative to north.

### 2.1 Dust Absorption

Intensity decrements for dust cloud absorption were measured for many features using radial and azimuthal profiles on the PC images. Figure 2 shows a sample profile from the V-band image of NGC 4736 taken along the strip shown in Figure 3 (which uses the G4-G12 enhanced image from Fig. 6 below). Most of the intensity dips deeper than a few hundredths of a magnitude are dust features, as can be determined by comparing Figures 2 and 3. The thinner features may be verified as dust also by comparing Figure 2 with the G1-G3 and G2-G6 images in Figure 6.

Intensity decrements do not measure the cloud extinction because there is likely to be a smooth intensity contribution from foreground starlight. If dust with opacity $`\tau `$ in the midplane is next to an intercloud medium with opacity $`\tau _0`$, then the cloud to intercloud intensity ratio is $`\left(1+e^\tau \right)/\left(1+e^{\tau _0}\right)`$. For $`\tau >\tau _0`$, this ratio is between 0.5 and 1, making the maximum cloud contrast 0.75 mag. Thus the magnitude difference is not a good measure of intrinsic extinction unless $`\tau `$ is small and $`\tau _0`$ is either independently known or close to zero. An intensity ratio lower than 0.5 corresponds to an opaque dust cloud on the near side of the galaxy.

Radiative transfer models of the intensity distributions in several passbands are necessary to determine the intrinsic opacities of these clouds, but the present HST data were not acquired for this purpose and the wavelength coverage and integration times are not adequate for a proper analysis. Dust spirals in another, similar, nuclear region were analyzed by radiative transfer previously using HST data in V and I passbands (Elmegreen et al. 1998). We found that typical extinctions are several magnitudes inside the darkest clouds, and that the average extinction in the inner disk is so low that the gas is not likely to be self-gravitating. We expect the same situation here because of a lack of star formation associated with the dust. Radiative transfer models of dust features for main galaxy disks are more common than for nuclear disks (e.g., Howk & Savage 1997; Kuchinski et al. 1998; Regan 2000). Galaxy extinctions have also been measured using background sources (e.g., Pizagno & Rix 1998; White, Keel, & Conselice 2000; Keel & White 2001; Elmegreen et al. 2001)

The typical absorption trough for a dust feature in Figure 2 is several hundredths of a magnitude in V band. The deep feature indicated is part of a long dust spiral that covers a factor of two in radius and then branches into two thinner spirals in the inner part (see Fig. 3). Figure 4 shows the magnitude difference for this spiral and several others in the two galaxies as a function of radius. These differences were measured by comparing the intensities of the dust lanes to the adjacent interdust regions at several radii. Estimated measuring uncertainties are 0.002 mag. Each line in Figure 4 is a different spiral; the top dashed line corresponds to the dust feature shown in Figure 2. The darkness of the dust generally increases inward, suggesting that either the density of the spirals or their thickness relative to the stellar disk increases inward.

### 2.2 Image Enhancement, Fractal Structure, and Dust Opacity

The spiral dust structures are very faint. They require enhancement with unsharp masks or other filtering techniques. Figures 5 and 6 show composites of four filtered images for each galaxy. Each filter isolates structures within a range of scales spanning a factor of 3. The top left panels contain unsharp-masked images made by subtracting a smoothed version from the original image, using a smoothing function that is a Gaussian with a dispersion of $`\sigma =3`$ pixels. This image shows the smallest structures, from 1 to 3 pixels. In the top right, a smoothed image with a $`\sigma =6`$ px Gaussian was subtracted from a smoothed image with a $`\sigma =2`$ px Gaussian (this will be called a G2-G6 image). This shows features with scales between 2 px and 6 px. The bottom left and right images subtract $`\sigma =12`$ px from $`\sigma =4`$ px convolutions (G4-G12) and $`\sigma =24`$ px from $`\sigma =8`$ px convolutions (G8-G24), respectively. Each image has considerable structure, indicating that clouds span a wide range of scales.

The fractal dimension of multiscale structure can be determined from box counting or cloud counting techniques, plotting the number of structures as a function of scale in log-log coordinates (Mandelbrot 1982; Westpfahl et al. 1999). Because we have already isolated certain scales with each filter, the cloud counting results on the subtracted images will give the same slope as the derivative with respect to length of the box counting results on the original images. Recall that box counting includes all structures larger than the scale of the box. The fractal dimension $`D`$ is the slope of the function $`N(>S)`$ versus $`S`$ on log-log coordinates for number of structures larger than $`S`$ determined by the box counting technique. If we write $`N(>S)=_S^{\mathrm{}}n(S)𝑑S`$ for number of structures $`n(S)`$ between $`S`$ and $`S+dS`$, then $`n(S)dSS^{D1}dS`$ and $`n(S)d\mathrm{log}SS^Dd\mathrm{log}S`$ for linear and log intervals of size, respectively. Thus the slope on a log-log plot of the number of clouds in each interval of filter size, for logarithmically spaced intervals like we have in Figures 5 and 6, gives the fractal dimension of the cloud structures. We choose the filter technique rather than the box technique so we can see the nature of the structure on each scale, whether it is filamentary or shell-like, for example. Figures 5 and 6 indicate that it is filamentary on all the resolved scales.

Figure 7 shows the number of spirals intercepted by deprojected azimuthal scans on the G2-G6 image plotted as a function of radius. The radii used for the measurements are 22, 45, 67, 90, and 112 pixels. A deprojected azimuthal scan of this image showing average counts for the radial interval from 66 to 68 pixels is shown in Figure 8. The number of clouds plotted in Figure 7 at this radius, 22 clouds, is shown again in Figure 8, with each cloud indicated by an arrow. The full circle of the azimuth in Figure 8 corresponds to 421 pixels. Generally the cloud count can be done easily by eye on a computer screen down to the brightness limit of the image. The number of dark features in an azimuthal scan is always about one-third of the radius, which makes the azimuthal spacing between them $`6\pi `$ pixels. Because the average cloud size on this image is $`4`$ px, the clouds fill $`4/\left(6\pi \right)=21`$% of the circumference.

This spiral filling factor makes sense given the extinction contrast. If the spirals formed by the compression of a fraction $`f`$ of the gas and the spiral spacing is $`C5`$ times the spiral width ($`C`$ is the inverse of the filling factor given above), then the column density of each spiral is $`(C1)f+1`$ times the ambient average. This column density has a contribution from uncompressed material that was at the spiral position already, which is the fraction $`1f`$ of the ambient column density, plus the compressed material from the whole volume out to the next spiral, which is the amount $`fC`$ times the ambient column density. The sum is $`1f+fC=(C1)f+1`$ times the ambient column density, as above. The assumption here is that not all of the interstellar medium will be compressed by the turbulence. A warm or hot component, with a thermal speed larger than the mean turbulent speed, should flow around the moving clouds and not compress easily. Thus $`f`$ may be considered the mass fraction of the gas that is compressible by the prevailing turbulence.

In this situation, suppose the spiral opacity is $`\tau `$, the interspiral opacity is $`\tau _0`$, and the ambient average before the compression was $`\tau _a`$, all from gas concentrated in the midplane. Then the compression gives $`\tau =\left[\left(C1\right)f+1\right]\tau _a`$, the conservation of mass gives $`\tau +\left(C1\right)\tau _0=C\tau _a`$, and the density contrast in magnitudes, $`\mathrm{\Delta }m`$, gives

$$\frac{1+e^\tau }{1+e^{\tau _0}}=10^{0.4\mathrm{\Delta }m}$$
(1)
These equations can be solved for the three opacities given $`C`$, $`f`$, and $`\mathrm{\Delta }m`$. This solution has a simple expression in the limit of small $`\mathrm{\Delta }m`$. In that case, there are two solutions for $`\tau `$, $`\tau _0`$ and $`\tau _a`$, one where the opacities are very large so hardly any background light gets through, and another where the opacities are very small. The relatively smooth continuation of the outer disk starlight into the inner region implies the latter. Then in the limit of small $`\mathrm{\Delta }m`$,

$$\tau _a0.8\mathrm{ln}10\mathrm{\Delta }m/\left(fC\right).$$
(2)
For $`C5`$ and $`f1/2`$, this gives $`\tau _00.7\mathrm{\Delta }m`$ on average. Considering $`\mathrm{\Delta }m0.04`$ from Figure 4, the ambient opacity of $`\tau _a0.03`$ corresponds to a mass column density of $`0.5`$ M pc-2. This result is $`10`$ times smaller than the average extinction in the center of NGC 2207, which was about 0.6 mag (Elmegreen et al. 1998), but the dust features here are ten times fainter than in NGC 2207, where a prominent dust lane used for the opacity measurement had a V-band intensity contrast of 0.3 mag. Evidently, the column densities in NGC 4450 and NGC 4736 are too small for the inner gas disks to be significantly self-gravitating; for reasonable epicyclic frequency and velocity dispersion, the Toomre stability parameter $`Q`$ exceeds 100.

We can do this same exercise for spirals of different sizes. The top part of Figure 9 shows the number of dust spirals counted in different azimuthal scans versus the size of the Gaussian filter, plotting the G1-G3 filter combination as a point with an average Gaussian radius of 2, for example. The radii used for the azimuthal scans are 45, 67, 90 and 112 pixels. The number of dust spirals decreases approximately inversely with Gaussian filter radius for all galactic radii, suggesting again that the spirals are packed in tightly, even on a variety of scales. If the spacing between spirals in the azimuthal direction is $`\lambda `$, and this spacing is the factor $`C`$ times the spiral width, $`W`$, then the number $`N`$ of spirals at radius $`R`$ is $`2\pi R/\lambda `$ and the inverse of the filling factor is $`C=2\pi R/\left(NW\right)`$. This inverse filling factor, $`C`$, is shown in the bottom of Figure 9. For all but the smallest features, which are probably undercounted, $`C35`$ for all sizes and radii. This is the same value of $`C`$ that was used in the previous discussion of $`\tau _a`$. Thus the spirals are likely to be made on a wide variety of scales by compressions of a factor of $`C35`$ from an ambient medium with an opacity slightly less than 0.1 mag. We suggest below that the mean gas density could be $`1.5`$ cm-3.

The top part of Figure 9 may also be used to determine a fractal dimension, which equals the slope according to the discussion at the beginning of this subsection. The number of clouds decreases a little less rapidly than the inverse of the smoothing scale $`\sigma `$, giving an apparent fractal dimension slightly less than 1. If we write the dependence of the number of counted clouds on the cloud size as a power law, $`\sigma ^D`$, for logarithmic intervals of $`\sigma `$, then $`D`$ is the fractal dimension. From Figure 9, $`D`$ equals $`0.79\pm 0.16`$ for NGC 4450 and $`0.51\pm 0.04`$ for NGC 4736, considering Gaussian radii from 4 to 32.

### 2.3 Fourier Transform Power Spectra

Fourier transform power spectra of the azimuthal intensity profiles at ten radii in each galaxy are shown in Figure 10. The radii $`R`$ were chosen so that $`2\pi R`$ equals an allowable number of pixels in a Fast Fourier Transform done with the SRCFT subroutine in the IBM ESSL software package. These numbers give radii $`R=`$ 33, 57, 81, 107, 134, 153, 183, 204, 229, and 244 pixels.

The power spectra are the sums of the squared real and imaginary parts of the FFTs. These are one-dimensional FFTs, taken from the azimuthal intensity traces along the deprojected images. Deprojection was made with a position angle of $`7^{}`$ and an inclination of $`42.2^{}`$ for NGC 4450 on the HST image, and with a position angle and inclination of $`2^{}`$ and $`35.6^{}`$ for NGC 4736. The inclinations were taken from the Third Reference Catalogue of Bright Galaxies (de Vaucouleurs et al., 1991). The position angles were measured from the HST images of the nuclear regions and differ from the position angles of the outer disks. In terms of the number of degrees counter-clockwise from North, the position angles of the nuclear-disk major axes are $`1^{}`$ and $`21^{}`$ for NGC 4450 and NGC 4736, respectively.

The projected ellipses used for the intensity traces are shown on the G4-G12 images in Figure 11. Each intensity trace is an average over 9 deprojected azimuthal scans spaced by one pixel and centered on one of the radii $`R`$ given above. Each pixel used for an azimuthal average is blackened in Figure 11.

The power spectra are illustrated in Figure 10 by stacking the different radii on top of each other with the smallest radius at the bottom. The vertical scale corresponds to a factor of 10 in the power spectrum for each tic mark; the absolute scale is arbitrary. The solid, dashed, and dot-dashed lines have slopes of $`5/3`$, $`1`$, and $`8/3`$. Solid lines with a slope of $`5/3`$ are drawn through the low frequency portions of each power spectrum to illustrate the general trend.

The power spectra all have slopes of about $`5/3`$ at low frequency. The slopes flatten at high frequency to a slope of about $`1`$ or flatter. The low frequency power spectrum could be steeper than $`5/3`$ at the smallest radius, possibly as steep as $`8/3`$, as shown by the dot-dashed line.

Turbulence gives a one-dimensional power spectrum with a slope of $`5/3`$ for structures that are thinner on the line-of-sight than the transverse wavelength, and a slope of $`8/3`$ for structures that are thicker on the line-of-sight than the wavelength (Lazarian & Pogosyan 2000). Two-dimensional power spectra have slopes for these two cases that are steeper by one, i.e., $`8/3`$ and $`11/3`$. The azimuthal scan results here, giving slopes of $`5/3`$ at small spatial frequencies, are consistent with turbulent structures that are wider than the galaxy disk scale height.

The same slope of $`5/3`$ was found for power spectra of azimuthal scans of the HI intensity of the Large Magellanic Cloud (Elmegreen, Kim, & Staveley-Smith 2001). Thus the origin of the structure is probably the same. Power law power spectra have also been found in other interstellar emission line surveys, using HI (Crovisier & Dickey 1983; Green 1993) and CO (Stutzki et al. 1998) in the Milky Way and HI in the Small Magellanic Cloud (Stanimirovic et al. 1999). The slopes are always consistent with turbulence as the origin of the structure.

The LMC showed a turnover in the power spectrum at high spatial frequency in both the azimuthal and linear scans, going from slopes of $`5/3`$ to $`8/3`$ at high spatial frequency, and it showed a similar turnover in the two-dimensional power spectrum of the whole disk, going from $`8/3`$ to $`11/3`$ at high frequency. These turnovers were proposed to correspond to the transition from structures that are larger in the transverse direction than the galaxy thickness to structures that are smaller than the galaxy thickness (Elmegreen, Kim, & Staveley-Smith 2001), as predicted in the general case by Lazarian & Pogosyan (2000). Elmegreen et al. also found that LMC thickness increased with radius. This result was confirmed by Padoan et al. (2001), who found a similar thickness-dependent feature in the spectral correlation function of HI emission from the LMC, and used this feature to map the inferred disk thickness around the plane. Analogous turnovers in the power spectrum are not seen for the nuclear dust spirals studied here, presumably because the noise level is too high. This problem is suggested by the flattening of the power spectra at intermediate to high frequencies – before the turnover can occur. Longer integration times on these sources could possibly get the noise down to a level where the power spectrum turns over. Then the line-of-sight disk thicknesses could be measured. All we can say here is that the thickness is less than $`50`$ pc in NGC 4450 and less than $`10`$ pc in NGC 4736. These limits are not unreasonable because they can still give the inner gas disks about the same aspect ratios as a typical main galaxy gas disk. When combined with the average extinction estimate given in the previous subsection, $`A_V0.03`$ mag, a disk thickness of $`10`$ pc corresponds to an average gas density of $`1.5`$ cm-3, which is about the same as that in the Solar neighborhood.

The morphologies of the dust structures in the nuclear regions of NGC 4450 and NGC 4736 differ from the morphologies of the gas emission structures in the Milky Way and Magellanic Clouds surveys. In the nuclear regions, the dust is mostly in the form of spirals, which show a power-law power spectrum in the azimuthal direction and a gradient of intensity in the radial direction (cf. Fig. 4). For the HI and CO emission line surveys, the structures are more blob-like or shell-like, with essentially no spirals, even in the whole LMC and SMC. Evidently, the origin of these gaseous structures is the same, i.e., random and persistent compressions from turbulence, but the nuclear regions of NGC 4450 and NGC 4736 have a lot of shear also, and this shear distorts the compressed regions into spirals. The LMC and SMC have relatively little shear, as do local interstellar clouds.

Numerical simulations of this shearing effect for gravitating media were in Toomre & Kaljnas (1991), Huber & Pfenniger (1999), Wada & Norman (1999), and Semelin & Combes (2000). Pressures from star formation were also considered by Wada, Spaans, & Kim (2000) and Wada & Norman (2001). The present case differs from these others slightly because the nuclear dust structures here are mostly without star formation and they are probably also without significant self-gravity. The origin of the turbulence is unknown, although a mild, pressure-driven instability was suggested by Montenegro et al. (1999).

## 3 Discussion

The nuclear spirals in NGC 4450 and NGC 4736 vaguely resemble the outer spirals, which are somewhat flocculent in each case, but the nuclear spirals do not continue smoothly from the outer spirals and there are important structural differences. The B-band image of NGC 4450 in the Carnegie Atlas of Galaxies (Sandage & Bedke 1994) shows two long dust spirals in the main disk, along with some flocculent structure; the stellar spiral arms are smooth. In contrast, the nuclear region has no stellar arms and at least 7 prominent dust arms, some with pitch angles as high as $`45^{}`$ and some crossing each other. The eastern side of the nuclear region shows more dust than the western side because of the galaxy’s inclination. Some small dust feathers extend nearly radially from the center toward the south, reminiscent of jets. HST spectral observations by Ho et al. (2000) reveal double-peaked line profiles with high velocity wings, characteristic of accretion disk activity observed in other LINERs.

NGC 4736 is an early-type galaxy with an outer ring and a circumnuclear starburst ring. Its main disk structure is flocculent and defined primarily by the dust. Sandage & Bedke’s B-band print shows the inner disk structure as composed of many tightly-wrapped arms, but the central region is saturated in the reproduction. Waller et al. (2001) present UIT UV and ground-based R-band images of the central regions, including an unsharp-masked image showing the complicated flocculent structure of the main disk. They also show an HST FOC image of the main nuclear dust arms. Ground-based NIR observations by Mollenhoff, Matthias, & Gerhard (1995) suggested a weak stellar bar with a length of $`20^{\prime \prime }`$, which was also noted by Maoz et al. (1995) from an HST FOC image. The bar was observed in CO by Sakamoto et al. (1999) and Wong & Blitz (2000). In the HST image, the region corresponding to the bar shows up as an elongated disk with a position angle nearly perpendicular to the major axis of the galaxy. The structure inside the circumnuclear ring, which is really two tightly wrapped arms, consists of a dozen dust arms within a radius of 50 pc from the center, branching to dozens more dust spirals out to 200 pc. The nuclear dust spirals are not attached to the main inner disk dust spirals.

The nuclear dust in NGC 4450 and NGC 4736 has several characteristics that differ from spiral arms and dust clouds in main galaxy disks. These are:

* The nuclear dust spirals shown here have no associated star formation. Other nuclear spirals in different galaxies have star formation (e.g., Coma D15 in Caldwell, Rose, & Dendy 1999), so the gaseous nature here is not universal. The lack of star formation suggests that the inner gas disks in NGC 4450 and NGC 4736 are not strongly self-gravitating. The same was true for NGC 2207 (Elmegreen et al. 1998) and for several other inner disks in the study by Martini & Pogge (1999). Our opacity estimate in Section 2.2 is also consistent with this.
* The nuclear dust is in the form of spiral arms of various pitch angles, widths, and lengths. Some of the arms are trailing, a few are leading, and many cross each other. This pattern is generally more irregular than main disk flocculent arms (see atlas in Elmegreen 1981). Main disk flocculent arms are rarely leading. They generally do not cross each other; if they branch into spurs, then this branching is toward larger radii (Elmegreen 1980). They also have star formation that gives them a thicker, more patchy quality, rather than a filamentary quality.
* The nuclear dust spirals in NGC 4450 and NGC 4736 have decreasing contrast with increasing radius. Ambient dust extinctions generally decrease with galactocentric distance because of the exponential distribution of gas column density. The flocculent galaxy NGC 5055 has such a decrease, for example, as measured by the extinctions of OB associations (Acarreta et al. 1996). Nuclear dust spirals are not just ambient extinctions, however. They are morphologically more similar to main disk spiral arms than diffuse cloud extinctions because they are organized and most likely formed by compressive processes in the presence of shear. From this point of view, nuclear dust spirals should be compared to main disk spirals, and then the radial decrease in nuclear spiral amplitude is unusual. Density wave spirals in non-barred galaxies tend to get stronger with increasing radius, out to at least the corotation zone (Elmegreen & Elmegreen 1984; Elmegreen et al. 1996). The unusual result that nuclear dust spirals get weaker with radius is presumably the result of crowding near the center for waves that move inward, as predicted for solutions to Bessel’s wave equation (Elmegreen et al. 1998; Montenegro et al. 1999). This is a different dynamical situation than for main disk spiral arms, for which the curvature terms ($`1/kr`$ for wavenumber $`k`$) in the wave equations can usually be ignored (e.g., Bertin et al. 1989).
* The number density of dust features is about constant with radius, indicating that the inner disk is completely filled with structure. This is unlike the situation for main galaxy disks which often have a small number of arms (e.g., 2-5) that get further apart with radius. The nuclear spirals also have some indication of a hierarchical or fractal structure because of a non-integer slope of the size distribution, examined with unsharp masks. This size distribution is approximately a power law with a slope in the range from 0.5 to 0.8. This power law is reminiscent of other properties of interstellar clouds formed by turbulence, but complicated in this case by the effects of shear, which make spirals rather than clumps, and by the Coriolis force, which resists turbulent motions on large scales. The dust features are also very weak, and the smaller clouds, as well as those further from the center, are difficult to see above the pixel noise. Fractal structure in the dust of another galaxy was also found by Keel & White (2001) using a background illumination source. Fractal structure in galactic clouds is well known (e.g., Falgarone, Phillips, & Walker 1991).
* Fourier transform power spectra in the azimuthal direction show the characteristic signature of turbulence compression, which is a power law slope of $`5/3`$ for one-dimensional structures that are larger than the line-of-sight thickness of the galaxy disk. This makes the nuclear dust features studied here resemble the HI clouds in the LMC, with the important difference that the nuclear clouds are spiral filaments, presumably affected by shear, and the LMC clouds are globular and shell-like in a low shear environment.

There is no direct evidence in our observations for accretion driven by the dust spirals. Radial velocities will have to be measured to determine this. However, the increase in dust opacity with decreasing radius for the main spirals is consistent with the amplification that is expected for inward motions. In that case, the spirals could drive nuclear accretion.

D.M.E. and B.G.E. gratefully acknowledge HST archival grant HST-AR-09197. We also acknowledge useful comments by the referee.