text:
# Parameter estimation for Galactic binaries by LISA

## 1 Introduction

The Laser Interferometer space Antenna (LISA), a joint project of NASA and ESA (European Space Agency) would establish gravitational wave astronomy at low frequency band $`(10^4\text{H}zf10^1\text{H}z)`$. It would bring us essentially new information of the Universe (Bender et al. 1998). For example gravitational waves from merging super massive black holes (SMBHs) would be detected with significant signal-to-noise ratio ($`SNR`$) $`10^3`$, though event rate of such merging is highly unknown (e.g. Haehnelt 1994, Vecchio 1997). Galactic binaries are promising sources of LISA (Mironowskii 1965, Evans, Iben & Smarr 1987, Hils, Bender & Webbink 1991, Webbink & Han 1998). Gravitational waves from some known binaries (e.g. X-ray binary 4U1820-30) would be detected with $`SNR>5`$ by one year integration (Bender et al. 1998). In addition more than thousands of close white dwarf binaries (CWDBs) are expected to exist in LISA band (for recent studies see Yungelson et al 2001, Nelemans et al. 2001, Napiwotzki et al. 2002). Our target in this article is these Galactic binaries. We examine how accurately information of binaries can be extracted from gravitational waves observed by LISA.

Cutler (1998) studied the estimation errors for binary parameters with special attention to angular variables such as direction and orientation of binaries (see also Peterseim et al. 1997; Vecchio & Cutler 1998; Cutler & Vecchio 1998; Hughes 2002; Moore & Hellings 2002). He used approximation that emitted gravitational waves would be monochromatic, namely neglected the effects of the chirp signal $`\dot{f}`$. But the wave frequency or the chirp signal are fundamental quantities for gravitational wave astronomy. From the measured chirp signal $`\dot{f}`$ we can obtain the so called chirp mass $`M_c=M_1^{3/5}M_2^{3/5}(M_1+M_2)^{1/5}`$ ($`M_1,M_2`$: masses of two stars) for a binary whose orbital evolution is determined by gravitational radiation reaction. Furthermore the distance to the binary could be estimated from the chirp signal $`\dot{f}`$ and the amplitude of the wave signal (Schutz 1986). The frequency $`f`$ itself also contains important information. One of the authors (Seto 2001) pointed out that signature of the periastron advance could be detected in gravitational waves from an eccentric binary by measuring its wave frequencies preciously. If this method works well, we can estimate the total mass $`M_{total}=M_1+M_2`$ of the binary beside the chirp mass, and consequently each mass of the binary is obtained separately.

Estimation errors for fitting parameters correlate complicatedly to each other and depend largely on observational situations. For example longer observational periods would improve not only signal to noise ratio but also resolution of the frequency space. Note that the latter is crucial for reducing Galactic binary confusion noise, as the number of resolved binaries increases with decrease of the frequency bin $`T_{obs}^1`$ (see Seto 2002 for details). At LISA age it would become an interesting observational challenge to optically identify the binaries whose gravitational waves are detected by LISA. We also discuss impacts of these observational efforts on estimation of binary parameters.

The magnitude of the chirp signal $`\dot{f}`$ has strong dependence on the wave frequency $`f`$ as $`\dot{f}M_c^{5/3}f^{11/3}`$. At lower frequencies the estimation error $`\mathrm{\Delta }\dot{f}`$ would be larger than the signal $`\dot{f}`$ itself. Then it would be better to remove $`\dot{f}`$ from fitting parameters and simply put $`\dot{f}=0`$ from the beginning. By evaluating the estimation error $`\mathrm{\Delta }\dot{f}`$ we can quantitatively discuss these prescriptions for signal analysis.

In this article we only study parameter estimation errors on the assumption that (i) signal has been detected with high SNR, and (ii) noises are Gaussian distributed. We evaluate the estimation errors using the Fisher information matrix for maximum likelihood method. For low SNR our simple analysis would not be valid, as the probability distribution function of the fitting parameters could become highly complicated (e.g. multimodal) (Balasubramanian & Dhurandhar 1998). Thus our result should be regarded as a lower bound of the estimation errors. Non-Gaussian nature of noises would farther increases errors. This is a very important problem, but detailed quantitative analysis would be a formidable task.

This article is organized as follows. In §2 we briefly discuss the gravitational waveforms of chirping binaries and the data stream obtained by LISA. Then we mention the parameter estimation based on the Matched filtering analysis. In §3 we numerically evaluate the parameter estimation errors, and discuss its dependence on observational period or wave frequency in detail. §4 is devoted to summary and discussions.

## 2 Gravitational Waveforms and Parameter Extraction

### 2.1 Gravitational Wave Measurement with LISA

LISA consists of three spacecrafts forming an equilateral triangle, and orbits around the Sun, trailing $`20^{}`$ behind the Earth. The sides of the triangle are $`L=5\times 10^6`$ km in length, and the plane of triangle is inclined at $`60^{}`$ with respect to the ecliptic. The triangle rotates annually. The gravitational wave signal is reconstructed from the three data streams that effectively correspond to three time-varying armlength data, ($`\delta L_1,\delta L_2,\delta L_3`$) for gravitational waves with $`\lambda L`$. We basically analyze two data streams given by $`s_I(t)=(\delta L_1(t)\delta L_2(t))/L`$ and $`s_{II}(t)=(\delta L_1(t)+\delta L_2(t)2\delta L_3(t))/\sqrt{3}L`$. These data can be regarded as the response of two $`90^{}`$-interferometers rotated by $`45^{}`$ to one another (Cutler 1998). The data $`s_{I,II}(t)`$ contain both gravitational waves signal $`h_{I,II}(t)`$ to be fitted by matched filtering and the noise $`n_{I,II}(t)`$. The latter is constituted by the detectors noise and the binary confusion noise. As in Cutler (1998) we assume that the noises are stationary, Gaussian and uncorrelated with each other.

The gravitational wave signals $`h_{I,II}(t)`$ from a binary are written as

$$h_{I,II}(t)=\frac{\sqrt{3}}{2}\left[F_{I,II}^+(t)h_+(t)+F_{I,II}^\times (t)h_\times (t)\right],$$
(1)
where $`F_{I,II}^{+,\times }(t)`$ are the pattern functions which depend on the source’s angular position of the binary ($`\overline{\theta }_S,\overline{\varphi }_S`$) , its orientation ($`\overline{\theta }_L,\overline{\varphi }_L`$) and detector’s configuration. The angular variables with bars are defined in a fixed barycenter frame of the solar system, and the direction and orientation of the binary are assumed to be constant during the observation in this frame. The quantities $`h_{+,\times }(t)`$ are the two polarization modes of gravitational radiation from the binary. We can estimate both ($`\overline{\theta }_S,\overline{\varphi }_S`$) and ($`\overline{\theta }_L,\overline{\varphi }_L`$) from the time profiles of the two signals due to LISA’s annual rotation and revolution. Further discussion and details about the pattern functions are seen in Cutler (1998). We basically use his formulation but newly include the effects of the chirp signal $`\dot{f}`$.

### 2.2 Gravitational Waveforms

We study short-period ($`10^4\mathrm{Hz}f10^2\mathrm{Hz}`$) binaries such as the close white dwarf binaries (CWDBs) or the neutron star binaries (NBs) in our Galaxy (Hils, Bender & Webbink 1991). We only discuss binaries with circular orbits. This is an excellent approximation for CWDBs, as their orbits are circularized by strong tidal interaction in their earlier evolutionary stages. At LISA band some NBs or neutron star-white dwarf binaries would have non-negligible eccentricities $`e0.1`$ induced at their supernova explosions (see e.g. Brown, Lee, Zwart, & Bethe 2001). But extension to these eccentric binaries is straightforward.

The chirping gravitational waveform is given by the quadrupole approximation Peters (1964) as

$`h_+(t)`$ $`=`$ $`A\mathrm{cos}\left[2\pi \left(f+{\displaystyle \frac{1}{2}}\dot{f}t\right)t+\varphi _D(t)+\varphi _0\right]\times \left[1+\left(\widehat{L}\widehat{n}\right)^2\right],`$
$`h_\times (t)`$ $`=`$ $`2A\mathrm{sin}\left[2\pi \left(f+{\displaystyle \frac{1}{2}}\dot{f}t\right)t+\varphi _D(t)+\varphi _0\right]\times \left(\widehat{L}\widehat{n}\right),`$ (2)
where $`\widehat{L}`$ (given by $`\overline{\theta }_L,\overline{\varphi }_L`$) is the unit vector in the direction of the binary’s orbital angular momentum, $`\widehat{n}`$ (given by $`\overline{\theta }_S,\overline{\varphi }_S`$) is the unit vector toward the binary and $`\varphi _0`$ is an integral constant. We regard the frequency $`f`$ and its time variation $`\dot{f}`$ as constants in the above equations. A purely monochromatic waveform has $`\dot{f}=0`$. This is the case studied by Cutler (1998). When gravitational radiation reaction dominates evolution of the binary as in the case of CWDBs or NBs, the chirp signal $`\dot{f}`$ is given as $`\dot{f}=(96\pi ^{8/3}/5)f^{11/3}M_c^{5/3}`$ with the chirp mass $`M_c`$. The perturbative expansion for the intrinsic evolution of the wave frequency in equation (2) are valid under the condition $`\dot{f}T_{obs}f`$, where $`T_{obs}`$ is the observational period. This condition is expressed as

$$f0.14\left(\frac{T_{obs}}{10\text{yr}}\right)^{3/8}\left(\frac{M_c}{1M_{}}\right)^{5/8}\text{Hz}.$$
(3)
The amplitude $`A`$ in equation (2) is given in terms of the wave frequency $`f`$, chirp signal $`\dot{f}`$ and the distance $`D`$ as

$$A=\frac{5}{96\pi ^2}\frac{\dot{f}}{f^3D}.$$
(4)
Thus we could determine the distance $`D`$, if we could measure three observables $`f`$, $`\dot{f}`$ and $`A`$ (Schutz 1986). This is an important aspect of gravitational wave astronomy.

The term $`\varphi _D(t)`$ in equation (2) is caused by revolution of LISA around the Sun and called the Doppler phase. Its explicit form is given by

$$\varphi _D(t)=2\pi fR\mathrm{sin}\overline{\theta }_S\mathrm{cos}\left[\overline{\varphi }(t)\overline{\varphi }_S\right],$$
(5)
where $`R=1`$ AU and $`\overline{\varphi }(t)=2\pi t/T(T=1\text{yr})`$ is the direction of LISA in the fixed barycenter frame.

### 2.3 Parameter Extraction

Let us briefly discuss the matched filtering analysis and the parameter estimation errors Finn (1992); Cutler & Flanagan (1994). We assume that the signal $`h_\alpha (t)`$ is characterized by some unknown parameters $`\gamma _i`$ (eight parameters in the present case: $`(A,f,\dot{f},\varphi _0,\overline{\theta }_S,\overline{\varphi }_S,\overline{\theta }_L,\overline{\varphi }_L)`$). In the matched filtering analysis the variance-covariance matrix of the parameter estimation error $`\mathrm{\Delta }\gamma _i`$ is given by inverse of the Fisher information matrix $`\mathrm{\Gamma }_{ij}`$ as $`\mathrm{\Delta }\gamma _i\mathrm{\Delta }\gamma _j=\left(\mathrm{\Gamma }^1\right)_{ij}`$. For a quasi-monochromatic binary ($`\dot{f}T_{obs}f`$) the noise spectrum $`S_n(f)`$ is nearly constant in the frequency region swept by the binary and the Fisher matrix simply becomes Cutler (1998)

$$\mathrm{\Gamma }_{ij}=\frac{2}{S_n(f)}\underset{\alpha =I,II}{}_0^{T_{obs}}𝑑t\frac{h_\alpha (t)}{\gamma _i}\frac{h_\alpha (t)}{\gamma _j}.$$
(6)
The error boxes for the angular parameters $`\widehat{L}`$ and $`\widehat{n}`$ become ellipses in the celestial sphere due to the correlation of two parameters $`\theta `$ and $`\varphi `$. In this article we represent the estimation errors for direction and orientation of binaries in the form defined in Cutler (1998) as follows

$$\mathrm{\Delta }\mathrm{\Omega }=2\pi \sqrt{\mathrm{\Delta }\mu ^2\mathrm{\Delta }\varphi ^2\mathrm{\Delta }\mu \mathrm{\Delta }\varphi ^2},$$
(7)
where we have defined $`\mu =\mathrm{cos}\theta `$. In the same manner the signal to noise ratio ($`SNR`$) is given by

$$(SNR)^2=\frac{2}{S_n(f)}\underset{\alpha =I,II}{}_0^{T_{obs}}𝑑th_\alpha (t)h_\alpha (t).$$
(8)
From equations (6) and (8) it is apparent that the expressions for the estimation errors $`\mathrm{\Delta }\gamma _i\mathrm{\Delta }\gamma _j`$ do not depend on the noise spectrum $`S_n(f)`$ when they are normalized by the signal-to-noise ratio (Cutler, 1998). In this article we extensively use this normalization method. For example the parameter estimation errors for a simple wave form $`h(t)=A\mathrm{sin}[2\pi (f+\dot{f}t/2)t+\varphi _0]`$ with four fitting parameters $`(A,f,\dot{f},\varphi _0)`$ is easily evaluated as in Seto (2002). The explicit forms for $`\mathrm{\Delta }A,\mathrm{\Delta }f`$ and $`\mathrm{\Delta }\dot{f}`$ are given by

$`{\displaystyle \frac{\mathrm{\Delta }A}{A}}`$ $`=`$ $`{\displaystyle \frac{1}{SNR}}=0.1\left({\displaystyle \frac{SNR}{10}}\right)^1,`$ (9)
$`\mathrm{\Delta }f`$ $`=`$ $`{\displaystyle \frac{4\sqrt{3}}{\pi }}{\displaystyle \frac{T_{obs}^1}{SNR}}=0.22\left({\displaystyle \frac{SNR}{10}}\right)^1T_{obs}^1,`$ (10)
$`\mathrm{\Delta }\dot{f}`$ $`=`$ $`{\displaystyle \frac{6\sqrt{5}}{\pi }}{\displaystyle \frac{T_{obs}^2}{SNR}}=0.43\left({\displaystyle \frac{SNR}{10}}\right)^1T_{obs}^2.`$ (11)
This simple analysis does not include information of the angular parameters, but would be helpful to understand more detailed numerical analysis in the following section.

## 3 Results

### 3.1 General Behavior

We have numerically evaluated the uncertainties of the estimated parameters for various quasi-monochromatic binaries. In this subsection we show results for a typical example with a fixed set of angular parameters at $`\mathrm{cos}\overline{\theta }_S=0.3,\overline{\varphi }_S=5.0,\mathrm{cos}\overline{\theta }_L=0.2`$ and $`\overline{\varphi }_L=4.0`$. In Table 1 we show LISA’s measurement accuracy for parameters $`(A,f,\dot{f},\mathrm{\Omega }_S,\mathrm{\Omega }_L)`$ at frequencies $`f=10^4,10^3`$ and $`10^2`$Hz. We present our results for two observational periods, $`T_{obs}=1`$ and 10yr. All results are normalized by $`SNR=10`$ after integration over each observational period $`T_{obs}`$. These results simply scale as $`(SNR/10)^1`$ for errors $`\mathrm{\Delta }A,\mathrm{\Delta }f`$ and $`\mathrm{\Delta }\dot{f}`$ and $`(SNR/10)^2`$ for error ellipses $`\mathrm{\Delta }\mathrm{\Omega }_{S,L}`$. The first row for each observational period $`T_{obs}`$ represents the case when all the eight parameters $`(A,f,\dot{f},\varphi _0,\overline{\theta }_S,\overline{\varphi }_S,\overline{\theta }_L,\overline{\varphi }_L)`$ are included in the matched filtering analysis. The second row corresponds to the case without the chirp signal $`\dot{f}`$. Results for $`T_{obs}=1`$yr are obtained under the same condition with Table 1 (case A) in Cutler (1998). Our numerical values completely coincide with his ones. The third row is the case when direction of the binary ($`\overline{\theta }_S,\overline{\varphi }_S`$) is given exactly by other method (e.g. optical identification of the binary) and removed from the fitting parameters.

Figs 1 and 2 show clearly dependence of the orbital parameters on wave frequency and observational period $`T_{obs}`$. In Fig.1 we show LISA’s measurement accuracy as a function of observational period $`T_{obs}`$ at given frequencies $`f=10^4,10^3`$ and $`10^2`$ Hz. All results are normalized by $`SNR=10`$ at integration period $`T_{obs}=1`$ yr. The solid lines are results for fitting all the eight parameters. The dotted lines represent the case when the angular position $`(\overline{\theta }_S,\overline{\varphi }_S)`$ are removed from the fitting parameters (see also Hughes 2002). The dashed lines are results with fitting only the angular position $`(\overline{\theta }_S,\overline{\varphi }_S)`$. For observational period $`T_{obs}2`$ yr, the difference between the solid and the dotted lines is very small especially for $`\mathrm{\Delta }A`$, $`\mathrm{\Delta }f`$ and $`\mathrm{\Delta }\dot{f}`$ irrespective of the frequency. Thus optical determination of the source direction $`\mathrm{\Omega }_S`$ would only slightly reduce the estimation errors of other parameters. Asymptotic behaviors of errors $`\mathrm{\Delta }A`$, $`\mathrm{\Delta }f`$ and $`\mathrm{\Delta }\dot{f}`$ are given by

$`{\displaystyle \frac{\mathrm{\Delta }A}{A}}`$ $`=`$ $`0.20\left({\displaystyle \frac{SNR}{10}}\right)^1,`$ (12)
$`\mathrm{\Delta }f`$ $`=`$ $`0.22\left({\displaystyle \frac{SNR}{10}}\right)^1T_{obs}^1,`$ (13)
$`\mathrm{\Delta }\dot{f}`$ $`=`$ $`0.43\left({\displaystyle \frac{SNR}{10}}\right)^1T_{obs}^2,`$ (14)
(see also Table.1), and the asymptotic time dependence is given as $`\mathrm{\Delta }AT_{obs}^{1/2},\mathrm{\Delta }fT_{obs}^{3/2},`$ and $`\mathrm{\Delta }\dot{f}T_{obs}^{5/2}`$, since we have $`SNRT_{obs}^{1/2}`$ from equation (8). We have also $`\mathrm{\Delta }\mathrm{\Omega }_{S,L}SNR^2T_{obs}^1`$ for directions and orientations of binaries. Numerical results for $`\mathrm{\Delta }f`$ and $`\mathrm{\Delta }\dot{f}`$ in equations (13) and (14) are almost identical to analytical ones in equations (10) and (11). Hence we can expect that the information of the angular parameters (the direction and orientation) does not affect the accuracy of the estimation for the frequency $`f`$ and the chirp signal $`\dot{f}`$ for observational period $`T_{obs}2`$ yr. But for the amplitude $`A`$, the above result (in ) is two times as large as equation (9). This is due to the fact that the amplitude $`A`$ is tied with the angular parameters (the inclination $`\widehat{L}\widehat{n}`$) in equation (2) and the estimation error $`\mathrm{\Delta }A`$ strongly depends on the error in $`\widehat{L}\widehat{n}`$ (see the detailed discussion in the next subsection).

Fig.2 is same as Fig.1 but given as a function of frequency $`f`$ with fixed observational period $`T_{obs}=1`$ and 10 yr. We find that for $`T_{obs}=10`$ yr the errors $`\mathrm{\Delta }A,\mathrm{\Delta }f`$ and $`\mathrm{\Delta }\dot{f}`$ do not depend on the frequency, as given in equations (12)-(14). But the angular resolution $`\mathrm{\Delta }\mathrm{\Omega }_S`$ depends on the frequency and is given by

$`\mathrm{\Delta }\mathrm{\Omega }_S`$ $`=`$ $`4.8\times 10^4\left({\displaystyle \frac{f}{10^2\text{Hz}}}\right)^2\left({\displaystyle \frac{SNR}{10}}\right)^2\text{sr},`$ (15)
for higher frequencies $`f10^3`$ Hz, and nearly constant for lower frequency $`f10^3`$ Hz.

Now we discuss correlation between Fisher matrix elements for $`T2`$yr. As expected from Figs 1 and 2, the source direction $`\mathrm{\Omega }_S`$ has almost no correlation with other parameters at higher frequencies $`f>10^3`$Hz 111In other words, Fisher information matrix $`\mathrm{\Gamma }_{ij}`$ in equation (6) becomes diagonalized to the two parts; $`\mathrm{\Gamma }\left(\begin{array}{cc}\mathrm{\Gamma }_A& 0\\ 0& \mathrm{\Gamma }_S\end{array}\right),`$ where $`\mathrm{\Gamma }_A`$ represents a $`6\times 6`$ matrix which corresponds to $`(A,f,\dot{f},\varphi _0,\overline{\theta }_L,\overline{\varphi }_L)`$ components and $`\mathrm{\Gamma }_S`$ represent a $`2\times 2`$ matrix which corresponds to $`(\overline{\theta }_S,\overline{\varphi }_S)`$ components. and weakly correlates with source orientation $`\mathrm{\Omega }_L`$ at lower frequencies $`f<10^3`$Hz. This dependence seems reasonable considering the information used to determine the source direction $`\mathrm{\Omega }_S`$, namely, the Doppler phase at $`f>10^3`$Hz and amplitude modulation at $`f<10^3`$Hz. The orientation $`\mathrm{\Omega }_L`$ is strongly related to the latter. As a result angular variables $`\mathrm{\Omega }_S`$ and $`\mathrm{\Omega }_L`$ correlate at lower frequencies. The detected wave form $`h_\alpha `$ is modulated by annual revolution and rotation of the detectors, and its phase has oscillating part due to the modulation that is more prominent at higher frequencies. Amplitudes of the derivatives $`h/\mathrm{\Omega }_s`$ also show this oscillation. In the Fisher matrix their cross terms with other derivatives are significantly canceled with long time integration. As a result the estimation errors for the angular direction $`\mathrm{\Omega }_s`$ show very weak correlation with other errors.

In the same manner we can also understand the asymptotic frequency dependence of the angular resolution $`\mathrm{\Omega }_S`$. In the derivatives $`h_\alpha /\gamma _i`$ for the Fisher matrix elements in equation (6) the frequency $`f`$ appears only through the Doppler phase $`\varphi _D(t)`$ in equation (5). At lower frequencies $`f<10^3`$Hz terms from the Doppler phase ($`f^1`$) is smaller than the terms from the amplitude modulation ($`f^0`$), and consequently the estimation error $`\mathrm{\Delta }\mathrm{\Omega }_S`$ does not depend on the frequency $`f`$. At higher frequencies $`f>10^3`$Hz the Doppler phase term becomes dominant and the Fisher matrix are diagonalized. Thus the error $`\mathrm{\Delta }\mathrm{\Omega }_S`$ depends on the frequency $`f`$ as $`\mathrm{\Delta }\mathrm{\Omega }_Sf^2`$ (Cutler & Vecchio 1998, Moore & Hellings 2002).

Finally we analytically investigate how the estimation errors improve with the observational period $`T_{obs}`$. We reanalyze the simple toy model $`h(t)=A\mathrm{sin}[2\pi (f+\dot{f}t/2)t+\varphi _0]`$ (see the sentences before equation (9)) adding the information of the source direction $`(\overline{\theta }_S,\overline{\varphi }_S)`$ by the Doppler phase $`\varphi _D(t)`$. This waveform is given as

$$h(t)=A\mathrm{sin}[2\pi (f+\dot{f}t/2)t+\varphi _D(t)+\varphi _0],$$
(16)
with six fitting parameters $`\gamma _i=(A,f,\dot{f},\varphi _0,\overline{\theta }_S,\overline{\varphi }_S)`$. We evaluate the magnitude of the variance-covariance matrix $`\mathrm{\Delta }\gamma _i\mathrm{\Delta }\gamma _j`$ by its determinant as $`det\mathrm{\Delta }\gamma _i\mathrm{\Delta }\gamma _j=(det\mathrm{\Gamma })^1`$. After some algebra the time dependence of $`det\mathrm{\Gamma }`$ is given analytically from equations (6) and (16) as,

$$det\mathrm{\Gamma }T_{obs}^{12}\left(1\frac{6}{\pi ^2\left(T_{obs}/1\text{yr}\right)^2}\right)\left(1\frac{90}{\pi ^4\left(T_{obs}/1\text{yr}\right)^4}\right),$$
(17)
where we assume that the observational period $`T_{obs}(f^1)`$ is a integer in units of year. In the above expression the quantity $`(det\mathrm{\Gamma })^1`$ formally diverges at $`T_{obs}=\sqrt[4]{90}/\pi =0.98`$ \[yr\] which is very close to 1, thus the estimation errors have large value at $`T_{obs}=1`$ yr. Fig.3 shows the inverse of the determinant of the Fisher matrix in the form $`T_{obs}^{12}(det\mathrm{\Gamma })^1`$ as a function of the observational period $`T_{obs}`$. From Fig.3, $`T_{obs}^{12}(det\mathrm{\Gamma })^1`$ rapidly converges at $`T_{obs}2`$ yr. In this case we also obtain the estimation errors as follows;

$`{\displaystyle \frac{\mathrm{\Delta }A}{A}}`$ $`=`$ $`{\displaystyle \frac{1}{SNR}},`$ (18)
$`\mathrm{\Delta }f`$ $`=`$ $`{\displaystyle \frac{4\sqrt{3}}{\pi }}{\displaystyle \frac{T_{obs}^1}{SNR}}\times \sqrt{{\displaystyle \frac{1(45/8)x^2(1+x^2)}{(16x^2)(190x^4)}}},`$ (19)
$`\mathrm{\Delta }\dot{f}`$ $`=`$ $`{\displaystyle \frac{6\sqrt{5}}{\pi }}{\displaystyle \frac{T_{obs}^2}{SNR}}\times {\displaystyle \frac{1}{\sqrt{190x^4}}},`$ (20)
$`\mathrm{\Delta }\mathrm{\Omega }_S`$ $`=`$ $`{\displaystyle \frac{1}{\pi f^2R^2\mathrm{cos}\overline{\theta }_SSNR^2}}\times {\displaystyle \frac{1}{\sqrt{(16x^2)(190x^4)}}}`$ (21)
$``$ $`4.3\times 10^4\left({\displaystyle \frac{\mathrm{cos}\overline{\theta }_S}{0.3}}\right)^1\left({\displaystyle \frac{f}{10^2\text{Hz}}}\right)^2\left({\displaystyle \frac{SNR}{10}}\right)^2\text{sr},`$ (22)
where $`x=\pi ^1(T_{obs}/1\text{yr})^1`$ and equation (22) is valid for $`T_{obs}2`$ yr and is similar to equation (15). The above results for $`\mathrm{\Delta }f`$, $`\mathrm{\Delta }\dot{f}`$ and $`\mathrm{\Delta }\mathrm{\Omega }_S`$ are about $`40\%`$ smaller than the the results in Fig.1 at $`T_{obs}=1`$ yr, because we do not include the information of the source orientation. $`\mathrm{\Delta }f`$ and $`\mathrm{\Delta }\dot{f}`$ are asymptotically same as equations (10) and (11) for $`T_{obs}2`$ yr.

### 3.2 Statistical Analysis

So far we have studied a specific set of the angular parameters $`\widehat{n}`$ and $`\widehat{L}`$. In this subsection we present statistical results for their various combinations at the asymptotic region $`T_{obs}2`$yr.

We have made 100 realizations of $`\widehat{n}`$ and $`\widehat{L}`$ that are distributed randomly on celestial spheres. Then we calculate the estimation errors $`\mathrm{\Delta }f`$ and $`\mathrm{\Delta }\dot{f}`$ for each binary normalized by $`SNR=10`$ with $`T_{obs}=10`$yr. We find that theses errors depend very weakly (less than 10% scatter) on the directions $`\widehat{n}`$ and $`\widehat{L}`$ in contrast to the results for $`T_{obs}=1`$yr (scattering typically factor $`2^{\pm 1}`$ around the average). We calculate their mean values and obtain results given in the following forms

$`\mathrm{\Delta }f`$ $`=`$ $`0.22\left({\displaystyle \frac{SNR}{10}}\right)^1T_{obs}^1,`$ (23)
$`\mathrm{\Delta }\dot{f}`$ $`=`$ $`0.43\left({\displaystyle \frac{SNR}{10}}\right)^1T_{obs}^2.`$ (24)
Note that these results do not depend on the frequency $`f`$ and would be useful for quantitative analysis of quasi-monochromatic binaries with $`T_{obs}2`$yr. As expected from the previous subsection the simple analytical estimations given in equations (10) and (11) are very close to equations (23) and (24).

We have also studied the estimation errors $`\mathrm{\Delta }\mathrm{\Omega }_S`$, $`\mathrm{\Delta }\mathrm{\Omega }_L`$ and $`\mathrm{\Delta }A/A`$. Distributions of the latter two have very large scatters. This is because (i) the Fisher matrix elements relating to the orbital orientation $`\widehat{L}`$ become singular at the highly symmetric face-on configuration $`\widehat{n}=\pm \widehat{L}`$ and (ii) the estimation of the amplitude $`A`$ is closely related to the inclination. For majority of realizations with $`|\widehat{n}\widehat{L}|0.8`$ we have typically $`\mathrm{\Delta }A/A0.2(SNR/10)^1`$. To determine the source direction $`\widehat{n}`$ we can use the information of the Doppler phase in addition to the annual amplitude modulation caused by LISA’s rotation. Thus the error $`\mathrm{\Delta }\mathrm{\Omega }_S`$ does not show such a bad behavior. For the above 100 realization at $`f=0.01`$Hz with $`SNR=10`$ and $`T_{obs}=10`$ yr the errors $`\mathrm{\Delta }\mathrm{\Omega }_S`$ are distributed as $`1.3\times 10^4\mathrm{sr}\mathrm{\Delta }\mathrm{\Omega }_S3.7\times 10^3`$sr with the mean value $`\mathrm{\Delta }\mathrm{\Omega }_S=7.1\times 10^4`$sr. Therefore the following relation roughly gives the estimation error for the source direction with $`T_{obs}2`$yr

$$\mathrm{\Delta }\mathrm{\Omega }_S7.1\times 10^4\left(\frac{SNR}{10}\right)^2\left(\frac{f}{10^2\mathrm{Hz}}\right)^2\mathrm{sr},$$
(25)
at $`f2\times 10^3`$Hz where the Doppler phase becomes more important than the annual amplitude due to the rotation of LISA.

Now let us discuss issues concerned with the chirp signal $`\dot{f}`$ based on equation (24). We assume that the chirp signal is dominated by gravitational radiation reaction and introduce a parameter $`R\mathrm{\Delta }\dot{f}/\dot{f}`$ that represents relative accuracy of the measured signal $`\dot{f}`$. For a given threshold $`R`$ we can solve the corresponding frequency $`f_R`$ as

$$f_R=9.2\times 10^4\left(\frac{M_c}{1M_{}}\right)^{5/11}\left(\frac{SNR}{10}\right)^{3/11}\left(\frac{T_{obs}}{10\mathrm{y}\mathrm{r}}\right)^{6/11}\left(\frac{R}{1.0}\right)^{3/11}\mathrm{Hz}$$
(26)
for observational period $`T_{obs}2`$yr. The frequency $`f_{R=1}`$ can be regarded as the critical frequency for treatment of the chirp signal $`\dot{f}`$. At higher frequencies $`ff_{R=1}`$ the parameter $`\dot{f}`$ should be included in the matched filtering, but the simple prescription $`\dot{f}=0`$ would be better at $`ff_{R=1}`$ since the expected signal would be completely buried in error. If the chirp signal is measured with accuracy $`R(1)`$, the chirp mass can be estimated with relative accuracy $`\mathrm{\Delta }M_c/M_c3R/5`$ from relation $`\dot{f}M_c^{5/3}`$. The chirp signal $`\dot{f}`$ is also essential to determine the distance $`D`$ to the binary. From the simple relation $`A=5\dot{f}/96\pi ^2f^3D`$, the estimation error for distance $`D`$ is roughly evaluated as

$`{\displaystyle \frac{\mathrm{\Delta }D}{D}}`$ $``$ $`{\displaystyle \frac{\mathrm{\Delta }A}{A}}+{\displaystyle \frac{\mathrm{\Delta }\dot{f}}{\dot{f}}},`$ (27)
$``$ $`\mathrm{Max}\{{\displaystyle \frac{\mathrm{\Delta }A}{A}},{\displaystyle \frac{\mathrm{\Delta }\dot{f}}{\dot{f}}}\},`$ (28)
$``$ $`0.2\left({\displaystyle \frac{SNR}{10}}\right)^1\mathrm{Max}\{1,\left({\displaystyle \frac{f}{1.4\times 10^3\mathrm{Hz}}}\right)^{11/3}\left({\displaystyle \frac{M_c}{1M_{}}}\right)^{5/3}\left({\displaystyle \frac{T_{obs}}{10\mathrm{y}\mathrm{r}}}\right)^2\},`$ (29)
where we have used the typical error for the amplitude estimation $`\mathrm{\Delta }A/A0.2(SNR/10)^1`$. For compact binaries with chirp mass $`M_c1M_{}`$ such as NBs or CWDBs the chirp signal $`\dot{f}`$ is the dominant source of the error at lower frequencies $`f10^3`$Hz and observational period $`T_{obs}10`$yr. The amplitude $`A`$ becomes dominant one at $`f10^3`$Hz.

Finally, we calculate the estimation error for the three dimensional position of the Galactic binaries. The signal to noise ratio is calculated from equation (8) as $`SNR380(f/5\times 10^3\text{Hz})^{2/3}`$ $`(D/10\text{kpc})^1(T_{obs}/10\text{yr})^{1/2}(M_c/1M_{})^{5/3}(\sqrt{S_n}/8\times 10^{21}\text{Hz}^{1/2})^{1/2}`$. The noise spectrum $`S_n(f)`$ (in units of $`\mathrm{Hz}^1`$) is nearly constant for $`3\times 10^3\text{Hz}f10^2`$ Hz (see Fig.5 in Cutler 1998). Hence one could determine both the angular position and the distance to the binary with the accuracy of

$`\mathrm{\Delta }\overline{\theta }_S`$ $``$ $`2\left({\displaystyle \frac{f}{5\times 10^3\text{Hz}}}\right)^{5/3}\left({\displaystyle \frac{D}{10\text{kpc}}}\right)\left({\displaystyle \frac{T_{obs}}{10\text{yr}}}\right)^{1/2}`$ (30)
$`\times \left({\displaystyle \frac{M_c}{1M_{}}}\right)^{5/3}\left({\displaystyle \frac{\sqrt{S_n}}{8\times 10^{21}\text{Hz}^{1/2}}}\right)^{1/2}\text{arcmin},`$
$`{\displaystyle \frac{\mathrm{\Delta }D}{D}}`$ $``$ $`5\times 10^3\left({\displaystyle \frac{f}{5\times 10^3\text{Hz}}}\right)^{2/3}\left({\displaystyle \frac{D}{10\text{kpc}}}\right)\left({\displaystyle \frac{T_{obs}}{10\text{yr}}}\right)^{1/2}`$ (31)
$`\times \left({\displaystyle \frac{M_c}{1M_{}}}\right)^{5/3}\left({\displaystyle \frac{\sqrt{S_n}}{8\times 10^{21}\text{Hz}^{1/2}}}\right)^{1/2}.`$
## 4 Conclusion

We have calculated LISA’s measurement accuracy for short-period binaries $`(10^4\text{H}z`$ $`f`$ $`10^1\text{H}z)`$ in our Galaxy, including the effects of chirp signal $`\dot{f}`$ and dependence on observational period $`T_{obs}`$. We find that the measurement accuracy rapidly improves for $`T_{obs}2`$yr comparing with $`T_{obs}1`$yr. This might be an important element for discussing operation period of LISA.

At observational period $`T_{obs}2`$ the errors for quantities $`f`$ and $`\dot{f}`$ is independent on the information of the positions and orientations of binaries in contrast to $`\mathrm{\Delta }A`$ and $`\mathrm{\Delta }\mathrm{\Omega }_{S,L}`$. It is also found that the estimation errors $`\mathrm{\Delta }A`$, $`\mathrm{\Delta }f`$ and $`\mathrm{\Delta }\dot{f}`$ are almost independent on the frequency $`f`$. The fitting formulae of the estimation errors for the source parameters (such as frequency $`f`$, chirp signal $`\dot{f}`$, amplitude $`A`$, angular position $`\mathrm{\Omega }_S`$ and distance $`D`$) are given as functions of $`f`$ and $`T_{obs}`$. We expect these would be powerful tools for quantitatively studying possibility gravitational wave astronomy.

The relative motion between the Sun and a binary could affect the frequency $`f`$ by the Doppler factor, $`(1+v_R)`$ where $`v_R`$ is the relative radial velocity. The tangential velocity $`v_T`$ or accelerating motion $`\dot{v}`$ would affect the time variation of the frequency, $`(\dot{f}/f)_{gal}v_T^2/cD`$ or $`\dot{v}/cv_{rot}^2/cR_g`$ respectively where $`v_{rot}`$ is the rotation velocity ($`200`$km/s) of the Galaxy, $`R_g`$ is the Galactic radius $`10`$kpc and $`D`$ is the distance to the binary (see e.g. Damour & Taylor 1991). For most Galactic binaries in LISA band this should be negligibly small compared with $`(\dot{f}/f)_{GW}`$ due to gravitational radiation reaction. This condition is expressed as $`f3\times 10^5(M_c/1M_{})^{5/8}(v_{rot}/200\text{km})^{3/4}(D/10\text{kpc})^{3/8}`$ Hz. For binaries very close to us $`D100`$ pc these effects would be important. But note that we might measure the proper motion of such binaries by optical observation.

We would like to thank an anonymous referee for helpful comments to improve the manuscript, and Takeshi Chiba, Hideyuki Tagoshi and Hirotaka Takahashi for useful comments and discussions. This work was supported in part by Grant-in-Aid of Scientific Research of the Ministry of Education, Culture, Sports, Science and Technology No. 0001416.


format: text

text:
# Red Giant Branch Stars as Probes of Stellar Populations I: 2MASS Calibration and Application to 2MASS GC01

## 1 Introduction

The red giant branch (RGB hereafter) stars are among the brightest red stars in stellar systems, older than a few Gyrs. These stars appear in almost all galaxies, including II Zw 40 (Östlin ostlin (2000)), considered until recently as the best candidate for a primeval galaxy. Therefore the red giants are a promising tool for probing the parameters of old populations and the history of star formation in any galaxy.

Galactic globular clusters, with their single age and metallicity, are the ideal sites for calibrating the RGB parameters. Since Da Costa & Armandroff (daC90 (1990)) provided the first reliable calibration of the position of the RGB versus metallicity, there has been a significant advancement, both because of the improvement of the astronomical instrumentation and the development of the corresponding theory.

The infrared waveband is particularly compelling for such studies, in comparison with the optical one, because of the potential to probe the stellar populations of systems with high foreground and/or intrinsic extinction. The relatively small size and field of view of the IR arrays have made it more difficult to carry out photometry of large areas, and to compile uniform samples, necessary to calibrate the RGB parameters reliably, in comparison with the optical region. The Two Micron All Sky Survey111See http://www.ipac.caltech.edu/2mass/releases/second/ doc/explsup.html (2MASS hereafter) offered for the first time such an opportunity.

Previous calibrations suffered a number of drawbacks. Kuchinski, Frogel, & Terndrup (kft95 (1995)) and Kuchinski & Frogel (kf95 (1995)) studied only metal rich globular clusters, with \[Fe/H\]$``$1. Later on, Ivanov et al. (2000a ) added to their sample three metal poor globular clusters. Ferraro et al. (ferr00 (2000)) used exceptional quality data but in a photometric system, based on unpublished standards by Ian Glass (South African Astronomical Observatory) with no available transformations to any of the other systems. The similarity of their final results to those of Ivanov et al., who used the CIT system, leads to a conclusion that the two photometric systems are not radically different.

To expand the basis of the RGB studies, we used high quality uniform photometry of Galactic globular clusters from the 2MASS point source catalog, and assembled a sample of RGBs of clusters with well known distances and reddening. We calibrated the behavior of the RGB slope, zero point, and tip with metallicity in the 2MASS photometric system. For the first time we offer such calibration in a well defined photometric system with all-sky coverage. We present a “fan” diagram, suitable for abundance distribution analysis. Our results complement the recent work of Cho & Lee (cho01 (2001)) and Grocholski & Sarajedini (gro01 (2001)). who explored the properties of the RGB bump.

## 2 2MASS Calibration of the RGB Parameters

### 2.1 Sample

Infrared photometry of about 80 of the 147 globular cluster, listed in Harris (harr96 (1996), revision June 22, 1999) is currently available from the First and the Second 2MASS incremental releases. However, many of the observed clusters are rendered unsuitable for this project because they suffer from one or more of the following drawbacks:

1. Severe foreground contamination. This refers mostly to globular clusters in the general direction of the Galactic center.
2. The photometry of some distant halo clusters or highly reddened clusters is not sufficiently deep.
3. Compact globular clusters have unresolved cores, limiting the number statistics, and leading to crowding problems.

After a visual inspection of the luminosity functions and color-magnitude diagrams (CMD hereafter) of all available globular clusters, we selected a subset of 27 objects (Table 1) for which reliable estimates of at least one RGB parameters could be obtained. To calculate the absolute magnitudes, and the intrinsic colors of stars, we adopted the reddening and distance estimates from the compilation of Harris (harr96 (1996)), and the reddening law of Rieke & Lebofsky (rl85 (1985); used throughout this paper).

A special care was taken to assure the uniformity of the metallicity data. Harris (harr96 (1996)) based his compilation on the system established by Zinn & Wsst (zw84 (1984)), Zinn (zinn85 (1985)), and Armandroff & Zinn (az88 (1988)). Carretta & Gratton (cg97 (1997)) suggested that this scale may overestimate the abundance of metal-rich globular clusters, and developed a new one. We list their measurements in Table 1 as well. Values in brackets are estimates based on the transformation they derive, to the Zinn & West scale (see Eq. 7 of Carretta & Gratton).

Ferraro et al. (ferr99 (1999)) argued that the total amount of heavy elements is a better metallicity measurement, because it accounts naturally for the opacity variations, which depend on the total metallicity, not just the iron abundance (Salaris, Chieffi, & Straniero sal93 (1993)). They developed a news scale, designed to measure the that. For the clusters, absent in their sample, we adopted their prescription for $`\alpha `$-element enhancement (Sect. 3.4 in Ferraro et al.; estimates bracketed in Table 1).

Since the metallicities originate from different sources, and the corresponding uncertainties are often not quoted, we adopted an uniform error of 0.20 dex for all measurements. This is a typical value of the accuracy of the abundance estimates.

The first step before estimating the RGB parameters was to eliminate statistically the foreground contamination. As many stars were removed randomly from the CMD of the globular clusters, as the number of stars, present in a nearby field with the same area as the globular cluster field. The statistical removal is not reliable in case of severely contaminated globular clusters, where the number of foreground stars was comparable to the number of cluster members. A subset of CMDs for four clusters is shown in Figure 1 and 2. In addition, we removed the stars within 15 arcsec from the clusters centers, to minimize the crowding effects.

### 2.2 RGB Slope and Zero Point

Da Costa & Armandroff (daC90 (1990)) prepared a set of standard RGBs on the V versus V-I CMD, for metallicity estimates of globular clusters in the optical. The technique has been expanded toward the near infrared wavebands, where the giant branch is linear, and therefore, more obvious and easier to define compared with the optical or optical-infrared CMDs. The slope is related to the effective temperature of the stars along the RGB, and $`\mathrm{T}_{\mathrm{eff}}`$ in turn depends on the opacity, and the heavy element abundance. The slope is a reddening and distance-free metallicity indicator. It is significantly less demanding in terms of observing time and telescope collecting area than the spectroscopic methods, but unlike them it can be applied reliably only to uniform groups of stars.

We followed the procedure to determine the RGB slope, described in Ivanov et al. (2000a ), fitting the RGB on the $`[K_s,(JK_s)]`$ CMD with a linear equation: $`JK_s=a\times K_s+b`$. Only the stars above the horizontal branch were included. We applied a least square method, taking into account the uncertainties along both axis. Two iterations were carried, and all the stars outside $`5\sigma `$ from the first fit were removed from the calculations.

For the purpose of fitting the RGB of an individual cluster, the zero point is equivalent to previous calibrations of the RGB colors at fixed $`K`$-band magnitude levels (e.g. Ivanov et al. 2000a ; Ferraro et al. ferr00 (2000)), although more straightforward. To derive the zero-point relation to the abundance, we used only the clusters which suffer minimal extinction $`(\mathrm{E}_{BV}`$ 0.3 mag).

Notably, the fit to the zero point shows a very small variation with metallicity: 0.11 mag for \[Fe/H\]H96 varying from -2 to -0.5. For comparison, $`(JK_S)_0`$ at $`K_S=5.5`$ mag varies by 0.32 mag for the same metallicity range (Ferraro et al. ferr00 (2000)). This result offers the possibility for a more reliable reddening estimates than in case of RGB color calibrations.

The behavior with metallicity of the RGB slope and zero point are demonstrated in Figure 3. Fitting coefficients for individual cluster RGBs are listed in Table 2.

### 2.3 RGB Tip

The tip of the RGB is a well known distance indicator. Mould & Kristian (mou86 (1986)), Lee (l93 (1993)) and Lee et al. lfm93 (1993) pioneered the method in the optical, measuring the corresponding jump in the luminosity function in external galaxies. Later it was applied to number of resolved stellar systems, among the most challenging of which were some Virgo cluster members (Harris et al. harr98 (1998)) and NGC 3379 (Sakai et al. sak97 (1997)).

This is a statistically demanding technique, requiring 50-100 stars per bin. The limited size of infrared arrays explains the difficulty to apply it in the near infrared. Although it does cover a large area, the 2MASS photometry can not alleviate the intrinsic problem of the small number of giants in globular clusters. Thus, we assume that the brightest cluster member represents the RGB tip.

In most of the cases the tip is obvious, but sometimes additional criteria had to be applied to determine the brightest stars. We took advantage of the linearity of the RGB in the near infrared, and excluded from the considerations bright stars that deviated from the color of the RGB at a given magnitude level, predicted by the RGB slope fit by more than 0.5 mag. We also excluded some extremely bright stars, with luminosity higher than the rest of the RGB by 2-3 mag. They were obvious foreground contamination. Finally, the red variables from Clement et al. (cle01 (2001)) were excluded .

The formal uncertainties of the stellar magnitudes given in 2MASS were discarded, since they do not represent well the uncertainty in the tip magnitude. Instead, we adopted the difference in the magnitudes of the two brightest stars, accounting for the possibility that the brightest star may be a non-member. This led to typical error values of 0.2-0.4 mag, much larger than the 2MASS errors.

The RGB tip magnitudes for 20 globular clusters are given in Table 3. The behavior of the RGB tip with metal abundance is shown in Figure 4.

### 2.4 Results

The behavior of the derived RGB parameters with metallicity was fitted with linear equations, taking into account the errors along both axes. The coefficients are given in Table 4. Figure 5 shows a “fan”-like grid of RGBs in a reddening and distance corrected CMD, for the three metallicity scales discussed in Section 2.1. It demonstrates that the RGB behavior can be reduced to a simple rotation around a nearly-fixed point. This comes to no surprise, since the zero point of the RGB fit is almost independent of the abundance (Section 2.2).

The grid allows to obtain the metallicities of individual sdtars in resolved systems in the infrared, and to obtain the metallicity distributions as it has been done before in the optical (e.g. Saviane et al. sav00 (2000)). The linear representation of the RGB in absolute $`K_S`$ magnitude versus the intrinsic $`(JK_S)_0`$ color is:

$$M_{K_S}=RGB_{ZP}+RGB_{Sl}\times (JK_S)_0$$
(1)
Both the slope $`RGB_{Sl}`$ and the zero point $`RGB_{ZP}`$ were calibrated as functions of the abundance:

$$RGB_{Sl}=a_0^{Sl}+a_1^{Sl}\times [\mathrm{Fe}/\mathrm{H}]$$
(2)
$$RGB_{ZP}=a_0^{ZP}+a_1^{ZP}\times [\mathrm{Fe}/\mathrm{H}]$$
(3)
After substituting Eq. 2 and 3 in 1, and solving for the stellar metallicity:

$$[\mathrm{Fe}/\mathrm{H}]=[M_{K_S}a_0^{Sl}\times (JK_S)_0a_0^{ZP}]/[a_1^{Sl}\times (JK_S)_0+a_1^{ZP}]$$
(4)
This equation should be used with caution. Small photometric errors are crucial for accurate \[Fe/H\] estimates. The upper part of the infrared CMD is better suited for determining individual metallicities because near the root of the RGB the abundance relation degenerates. The exact faint limit is determined by the quality of the photometry, and by the presence of the horizontal branch at $`M_{K_S}1.5`$ mag.

## 3 2MASS GC01

Hurt et al. (hur99 (1999), hur00 (2000)) reported the serendipitous discovery of 2MASS GC01 (hereafter GC01) in 2MASS data. It is a heavily obscured globular cluster, laying in the Milky Way disk, in the general direction of the Galactic center. Ivanov, Borissova, & Vanzi (2000b ) determined $`(\mathrm{m}\mathrm{M})_0=12.414.0,`$ and $`\mathrm{A}_V=20.918.8,`$ assuming $`[\mathrm{Fe}/\mathrm{H}]=0.5`$ and $`2.0,`$ respectively. The main source of uncertainty in these estimates was the unknown metallicity of the cluster, although the location of GC01 suggested that it might be a metal rich object.

The first step toward a metallicity estimate of GC01 was to remove the foreground stars contamination. Unlike the clusters we used to derive the RGB parameters calibrations, the contamination here is severe, reaching $`30`$% in the RGB region. We performed 2000 foreground substractions, and estimated the RGB slope for each realization separately. This method yields distributions of the RGB slope and tip, and the widths of these distributions measure the respective uncertainties.

To carry out this procedure we defined the CMD area, encompassing the RGB: $`8.0K_S14.0`$ mag, and $`2.5JK_S6.0`$ mag. Then we divided it into 0.2 mag square bins. Experiments with different bin sizes indicated that any value between 0.2 and 1.0 mag leads to the same conclusions.

Next, we constructed CMDs for the cluster field, and for a surrounding field with an equal area. To minimize the crowding effects we omitted the stars within 15 arcsec from the cluster center. The outer limit of the cluster field was constrained by the cluster diameter (3.3$`\pm `$0.2 arcmin, Hurt et al. hur99 (1999), hur00 (2000)). We carried out our calculations to two values of the outer radii: 1.0 and 1.5 arcmin. Smaller values limit the number statistics, and larger ones increase the fraction of the foreground contamination. The foreground field was defined as a circular annulus with $`10.4910.54`$ or $`10.4310.54`$ arcmin size, respectively for 1.0 and 1.5 armin cluster fields.

Finally we counted the stars in each bin, and subtracted randomly from the “cluster” bins as many stars as were present in the “field” bins. If the latter bin had more stars than the former one, we subtracted stars from the nearby bins, again in a random way.

We carried out a linear fit on the RGB stars in the foreground subtracted CMD in the same manner as for the calibration clusters (see Section 2.2). For stars without error measurements in 2MASS we adopted $`\sigma =0.20`$ mag. Using the linearity of the RGB, we imposed a faint limit of the stars, included in the fit, just above the horizontal branch level. The luminosity function of GC01 (Ivanov et al. 2000b , Figure 5) indicates that the horizontal branch is at $`K_S`$=13.0-13.2 mag. To minimize the uncertainties of the RGB slope, we also imposed color limits on the stars we used in the fit. The red one was set to $`JK_S`$=5.5 mag, and has no effect on the slopes because of the negligible number of stars to the red of the RGB. The results are somewhat more sensitive to the blue limit. We choose to impose $`JK_S`$=3.50 and 3.80 mag, because of the well-defined limit of the RGB at this colors (Figure 6).

A summary of the results for the RGB slope and zero point GC01 is presented in Figure 7. The RGB tip is omitted because it is sensitive only to the adopted cluster radius. Clearly, the effects from the assumed parameters are smaller or comparable with the uncertainties originating due to the foreground contamination. We calculated the error-weighted averages for the realizations with $`K_S^{lim}=12.513.0,`$ to avoid the possible influence of the cluster horizontal brunch. The determined RGB parameters for GC01 are:

1. $`<\mathrm{Slope}>=0.0965\pm 0.0171`$
2. $`<\mathrm{ZeroPoint}>=5.40\pm 0.18`$
3. $`<\mathrm{Tip}>=8.87\pm 0.13`$

Using the equations, given in Table 4, we obtain: \[Fe/H\]H96= $`1.19\pm `$0.38, \[Fe/H\]CG97= $`1.06\pm `$0.34, and \[M/H\]F99= $`0.88\pm `$0.32. The errors given here include both the uncertainties in the slope, and in the calibrations. GC01 appears to resemble closely NGC 104.

The estimated absolute magnitude for the RGB tip, for a cluster with such abundance is $`M_{K_S}=6.30\pm 0.35`$ mag, using our new calibration for the Zinn metallicity scale. This leads to a distance modulus of $`(mM)_0+A_{K_S}=15.17\pm 0.38`$ mag consistent with Ivanov et al. (2000b ). We refrain from further considerations based on the RGB tip because of the poor statistics at the brighter end of the RGB. Instead, we will adopt $`(mM)_0+A_{K_S}=15.0\pm 0.4`$ mag, a result of interpolation between the values for \[Fe/H\]=$`1.0`$ and $`2.0,`$ in Table 2 of Ivanov et al. (2000b ).

We can also verify if the cluster reddening is consistent with the previous estimates. First, we subtract the RGB slope equations, written for the GC01 in apparent, and in absolute magnitudes. Respectively:

$$(m_Jm_{K_S})=\mathrm{RGB}_{\mathrm{ZP}}^{\mathrm{app}}+\mathrm{RGB}_{\mathrm{Sl}}\times m_{K_S}$$
(5)
$$(M_JM_{K_S})_0=\mathrm{RGB}_{\mathrm{ZP}}^{\mathrm{abs}}+\mathrm{RGB}_{\mathrm{Sl}}\times M_{K_S}$$
(6)
Here the notation is the same as in Table 4, and $`\mathrm{RGB}_{\mathrm{ZP}}^{\mathrm{abs}}=0.36\pm 0.20`$ mag, for the \[Fe/H\]H96 given above. The result of the subtraction is:

$$E(JK_S)=\mathrm{RGB}_{\mathrm{ZP}}^{\mathrm{app}}\mathrm{RGB}_{\mathrm{ZP}}^{\mathrm{abs}}+\mathrm{RGB}_{\mathrm{Sl}}\times [(mM)_0+A_{K_S}]$$
(7)
Substituting, we obtain: $`E(JK_S)=3.58\pm 0.37`$ mag, close to the values of 3.49 mag, interpolated as above from Ivanov et al. (2000b ). Our new estimate corresponds to $`A_V=21.07\pm 2.20`$ mag, and $`A_{K_S}=2.36\pm 0.25`$ mag.

## 4 Summary

The behavior of the RGB in the infrared was quantified based on an unprecedented large sample of 2MASS photometry of Milky Way globular clusters. The RGBs were fitted by straight lines. We produced new calibrations of the RGB slope, tip, and - for the first time - zero point, as functions of abundance. The introduction of the zero point streamlines greatly the RGB analyses in comparison with the traditional approach where RGB colors at given levels were used. Notably, the zero point is fairly insensitive to the abundance, varying by only 0.11 mag over a range from \[Fe/H\]$`{}_{H96}{}^{}=2`$ to $`0.5.`$ We present a “fan”-like diagram, suitable for analyses of the metallicity spread in resolved stellar systems.

The derived calibrations were applied to estimate the metal abundance of the recently discovered globular cluster GC01. It is a particularly challenging object because of the severe foreground contamination. We removed it with a random procedure, and used the RGB slope of the remaining pure cluster population to derive \[Fe/H\]H96=$`1.19\pm `$0.38. The uncertainty is dominated by the foreground contamination, and albeit large, allows to exclude extremely high abundances, expected from the cluster location. GC01 is likely to be moderately metal poor. The RGB tip and zero point yield distance modulus and extinction, consistent with our previous estimates.

###### Acknowledgements.

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, funded by the National Aeronautics and Space Administration and the National Science Foundation. The authors thank Dr. Ivo Saviane for the useful discussions, and the referee Dr. M. G. Lee for the comments that halped to improve the paper.


text:
# Localization, Coulomb interactions and electrical heating in single-wall carbon nanotubes/polymer composites

## Abstract

Low field and high field transport properties of carbon nanotubes/polymer composites are investigated for different tube fractions. Above the percolation threshold $`f_c`$0.33%, transport is due to hopping of localized charge carriers with a localization length $`\xi `$10-30 nm. Coulomb interactions associated with a soft gap $`\mathrm{\Delta }_{CG}`$2.5 meV are present at low temperature close to $`f_c`$. We argue that it originates from the Coulomb charging energy effect which is partly screened by adjacent bundles. The high field conductivity is described within an electrical heating scheme. All the results suggest that using composites close to the percolation threshold may be a way to access intrinsic properties of the nanotubes by experiments at a macroscopic scale.

Because of their structure, single-wall carbon nanotubes (SWNT) are ideal mesoscopic 1D systems. As such, their electronic transport properties are controlled by quantum size and charging effects. Most of these properties expected for 1D mesoscopic conductors have been observed on individual SWNT or bundles: Coulomb blockade and level quantizationtans , Luttinger liquid characteristicsbockrath , ballistic transportfrank …. On the other hand, measurements on macroscopic SWNT matfisher ; gaal ; yosida  show classical transport properties. In most of the cases, the conductivity is hopping like at low temperature suggesting localized charge carriers without any contribution of Coulomb interactions. Indeed, the large number of tubes involved along the conduction path should be partly responsible of the disappearance of the mesoscopic and metallic characteristics. However the reasons for the loss of the Coulomb interaction effects are unclear.

In this paper, we propose to investigate these differences by studying the transport properties of SWNT/polymer composites thin films where we limit the number of tube-tube contacts by varying the SWNT content. In all samples, the conductivity is hopping like with a localization length which seems to be independent of the SWNT content. In contrast with other worksfisher ; gaal ; yosida  Coulomb interactions associated with a soft Coulomb gap opening play a dominant role in the low temperature transport. We suggest that it originates from the Coulomb charging energy involved in the tube-tube transfer. In addition, large non linearities of the conductivity are observed at high electric field that we discuss in terms of electrical heating of the conduction gas. All our results indicate that close to the percolation threshold, some intrinsic properties of the tubes can be reached at a macroscopic scale. Coulomb interactions are essential in these systems. It corroborates the use of the Luttinger Liquid model for individual nanotubes. Still the mesoscopic character is lost because of the large number of tubes involved.

The SWNT are produced by arc dischargejournet  with a typical purity of about 70-90 %. Their characteristic diameter is in the range 1.3-1.5 nm and most of them are embedded in bundles with a typical size of 7-12 nm as shown by transmission electron microscopy (TEM). SWNT/PMMA composites are formed by drop casting a sonicated solution of SWNT, polymethylmetacrylate (PMMA) and toluene. Thin films are obtained with a thickness $``$10 $`\mu `$m. The SWNT content varies between a volume fraction $`f`$=0.1% to $`f`$=8%. Electron microscopies show that the films do not present large heterogeneities at their surface or in the bulknote1 . Transport measurements are performed within a planar four probe configuration with gold evaporated contacts $``$ 20 $`\mu `$m apart.

As already discussed elsewherebenoit , the evolution of the room temperature conductivity $`\sigma `$ of the composite thin films with $`f`$ is well described within the standard percolation theory. The critical behavior $`\sigma (ff_c)^\beta `$ with a percolation threshold $`f_c`$=0.33$`\pm `$0.03% and a critical exponent $`\beta `$=2.1$`\pm `$0.1 is shown in Fig.1. It is obeyed over two orders of magnitude in reduced fraction at room temperature and at 100 K. Indeed, the low value of the threshold $`f_c`$ arises from the very high aspect ratio (length over radius ratio) of the tubes which is found to be close to some hundredsbenoit . The critical exponent $`\beta `$=2.1 is in good agreement with the conventional 1.94 exponent found for random 3D connectivity. This quasi-ideal percolation behavior confirms that the dispersion of the conducting fillers in the matrix is good. Indeed, conduction is due to metallic SWNTsigma . Since a percolation behavior is obeyed, we expect that close to the threshold, very few metallic tubes contribute to the conduction pathway.

The temperature dependence of the resistivity $`\rho `$ of our samples is plotted in the lower panel of Fig.2. The resistivity is not metallic like in the investigated temperature range. It varies with temperature $`T`$ according to

$`\rho =\rho _0\mathrm{exp}(T_\alpha /T)^\alpha `$ (1)
with $`\alpha `$=1/2 below 30-50 K. This is the signature of localized carriers hopping. The exponent becomes $`\alpha `$=1/4 at higher temperatures. The composite films can be compared to a pressed pellet i.e. an ideal $`f`$=100 % sample. As shown in the upper panel, the pellet resistivity perfectly obeys a $`\alpha `$=1/4 law even at low temperature, in agreement with published datagaal ; yosida . The different $`\alpha `$ exponents as well as the slopes evolution in Fig.2 mean that some additional mechanisms set in at low temperature in addition to the network connectivity effect. It is reflected by the shift of the conductivity from the critical dependence shown in Fig.1 for $`f`$ 1% at $`T`$=20 K. It contrasts with more conventional percolating networksputten .

Hopping transport is frequently observed in SWNT matgaal ; yosida . In Mott’s modelmott  ($`\alpha `$=1/4), $`T_{\alpha =1/4}`$ is related to the localization length $`\xi `$ and to the density of states $`n(E_F)`$ at the Fermi level. Assuming that the charge carriers originate from the metallic tubes leads to a localization length $`\xi `$8.2 nm for the pressed pellet. Indeed, the use of Mott’s model here is questionable. It requires a constant DOS at the Fermi level which is the case for metallic SWNT. It is also restricted to a single phonon process. This condition usually limits the validity of the model to low temperature. However, the Debye temperature of SWNT should be quite highdebye  and thus the single phonon process is expected to extend to rather high temperatures.

As already noticed, the composite films show a different temperature dependence of the resistivity, described by Eq.1 with $`\alpha `$=1/2. This behavior indicates the presence of Coulomb interactionsefros  which open a soft Coulomb gapgap  $`\mathrm{\Delta }_{CG}`$ at the Fermi level related to the slope of the resistivity curves $`k_BT_{1/2}=2.8e^2/\kappa \xi `$ where $`\kappa =\kappa _{host}+4\pi e^2n(E_F)\xi ^2`$ is the effective dielectric constant ($`\kappa _{host}`$3 here). The evolution of the experimental $`T_{1/2}`$ with $`f`$ is given in inset of Fig.3. Below $`f`$=1%, $`T_{1/2}`$1000 K is almost constant. When $`f`$ increases up to 8%, $`T_{1/2}`$ decreases down to 200 K as 1/$`f`$. This scaling arises from the dilution of the SWNT in the matrix where $`n(E_F)=n^0(E_F)\times f`$ has to be considered instead of the pure SWNT DOS $`n^0(E_F)`$. The scaling parameter gives a localization length $`\xi `$7.4 nm in excellent agreement with the previous result. It indicates that the localization mechanism is intrinsic to the SWNT system independently of the $`f`$ fraction.

To our knowledge, while $`\alpha `$=1/4 hopping is widely reported in the literature on SWNT mat, this is the first time that a $`\alpha `$=1/2 regime is observed. It suggests that Coulomb interactions are effective in dilute SWNT systems but not in non dilute mat. Below $`f`$=1% the constant $`T_{1/2}`$ means a constant (soft) Coulomb gap $`\mathrm{\Delta }_{CG}`$2.5 meV. It decreases down to 0.5 meV when $`f`$ increases up to 8%. We argue that this Coulomb gap is due to the Coulomb charging energy involved in the transport process from one bundle (or tube) to the neighboring one. The charging energy is responsible for the Coulomb blockade phenomenon reported on isolated SWNT or bundles. Its order of magnitudetans ; bockrath  is 1-3 meV, in agreement with our value for $`\mathrm{\Delta }_{CG}`$. When increasing $`f`$ in the composite, the bundles (tubes) come closer to each other and a mutual screening occurs which reduces the charging energy. A crude estimate can be given using basic electrostatic theory. We compare the charging energy $`E_{c1}`$ of an isolated cylindric wire (an isolated bundle) with the charging energy $`E_{c2}`$ of a coax (an isolated bundle surrounded by its neighbors): $`E_{c1}/E_{c2}=\mathrm{ln}(L/R)/\mathrm{ln}((s+R)/R)`$ where $`L`$ is the length of the bundle, $`R`$ its radius and $`s`$ the bundle separation. A charging energy of $`E_{c1}`$3 meV corresponds to an isolated bundle with $`R`$=5 nm and $`L`$=2 $`\mu `$m. For $`f`$=10%, the average separation between bundles is 18 nm assuming that they run parallel to each other. Within this simple model, a reduction of the charging energy by a factor 4 is thus expected in rough agreement with our results. Below $`f`$=1%, the bundles are sufficiently far from each other for such a screening to be inefficient. Conversely increasing $`f`$ further leads to fully screened interactions in the case of the pellet ($`f`$40% gives already a screening effect by a factor of 10) and thus restores the Mott’s model. These results suggest that even in macroscopic composite films, intrinsic (nanoscopic) charging effects are still detectable close to the threshold.

Still localization is observed. The fact that $`\xi `$7-8 nm is much larger than the tube diameter (but much lower than the mean free path) suggests that localization does not occur at the tube scale. It means also that inter-tube transfers are not limiting in the investigated temperature rangestahl . Possible origins of localization are. i- The charges are localized at the bundle boundaries (transverse localization). Then $`\xi `$ gives an average diameter of the bundles. Our result is of the correct order of magnitude. ii- Localization occurs along the bundles by bundle-bundle (or tube-tube) contacts (longitudinal localization). Then the localization length is expected to depend on the filling of the composites. This is not observed. iii- Localization is due to defects along the tubes or into the bundles (bad matching of chirality, structural defects…). We cannot disregard this possibility. However, Raman and TEM characterizations show that the SWNT are of high quality.

Intrinsic localization is also reflected in the high field transport. Figure 4 shows also that large non linearities are observed in the $`I`$($`V`$) characteristics of the films (shown here for $`f`$=0.9 %) for moderate electric field strength $`E`$.

The experiments were carried out using standard pulse current techniques. Non linearities are observed in all samplesf1 . The non linear effects can reach up to 6 orders of magnitude at low temperature for the smallest SWNT contents. Increasing the temperature or $`f`$ leads to smaller effects. Indeed, we have checked that it is not due to self-heating. In localized systems, such non linear effects can have different originsmott ; jonscher : injection of excess carriers, modification of the emission probability, modification of the conduction pathway… We will not discuss these models however they do not apply here mainly because our system consists of conducting samples with a small DOS and a large localization length. Alternatively we suggest that the non linearities are due to an energy redistribution of the carriers caused by the high electric field jonscher ; shklovskii . Such a redistribution requires that the density of states is low and that empty states exist close to the filled states which is indeed the case in metallic SWNT. It is equivalent to a thermal effect and it can be associated with an ”electric temperature”shklovskii  $`T_{elec}`$:

$`k_BT_{elec}=\gamma eEa`$ (2)
where $`e`$ is the electron charge, $`a`$ the localization length and $`\gamma `$ a numerical factor close to 1. When applying high electric field at finite temperature $`T`$, both thermal and electrical heating play a role. The question of the effective temperature $`T_{eff}`$ seen by the carriers is not easy to answercleve . The simplest assumption is that $`T_{eff}=T+T_{elec}`$. Within this model, the conduction mechanism remains the same at low or high field but the thermostat temperature $`T`$ has to be replaced by the effective one $`T_{eff}`$. The curves of Fig.4 have been fitted according to $`J/E=\sigma _0\mathrm{exp}\{(T_{1/2}/[T+\beta E])^{1/2}\}`$. We adjust $`\beta `$ at each temperature $`T`$. As can be seen, the agreement is quite good. Unexpectedly, we find that $`\beta `$0.035 K.cm/V is almost constant for all the curves. It suggests that electrical heating really occurs. Typically $`T_{elec}`$35 K when $`E`$=100 V/cm. Obviously at high fields, $`T_{elec}T`$ and $`J/E`$ does not depend any more on $`T`$. According to Eq.2, $`\beta `$ is related to a localization length. Using $`\gamma `$=1 gives a length $`a`$28 nm. This value is higher but still of same order of magnitude as $`\xi `$8 nm found from the temperature dependence of the ”zero field” conductivity. It supports the use of the electric temperature concept.

We have suggested previously that localization is intrinsic to the SWNT system and it does not depend on $`f`$. Following this argument, we definitively fix $`\beta `$=0.035 K.cm/V (i.e. the localization length) for all samples and we plot the high field conductivity $`J/E`$ as well as the ohmic (low field) one $`\sigma `$ versus $`T_{eff}^{}{}_{}{}^{1/2}`$ on the same figure. The results are shown in Fig.5 for $`f`$=0.4, 0.9 and 4% respectively.

For a given $`f`$ fraction, all the curves (isothermal $`J(E)`$ and $`\sigma (T)`$) merge into a master plot characteristic of the $`f`$ fraction. It definitively confirms that localization is intrinsic and does not depend on the way the SWNT are dispersed into the insulating matrix. It shows also that a unique transport mechanism is responsible of both the low field and the high field conductivity.

In summary, this work shows that in SWNT/PMMA composites: i- Transport is due to localized carriers originating from the metallic SWNT with a characteristic localization length of 8-28 nm due to bundle boundaries. ii- Coulomb interactions arising from the charging energy limit the transport at low temperature. These interactions are partly screened by adjacent bundles. For a sufficiently large dispersion of the SWNT, screening is ineffective and the charging energy regime found in isolated SWNT transport is recovered at low temperature. iii- The connectivity of the nanotubes network into the matrix dominates the transport properties as soon as the temperature is high enough to wash out the Coulomb gap. iv- Because of the specific DOS of metallic SWNT, high electric fields can redistribute the energy of the carriers which leads to large non linearities of the $`I(V)`$ characteristics. Finally this work suggests that using composites close to the percolation threshold may be a way to access intrinsic tube properties at a macroscopic level. This approach may be developed for other transport properties such as thermal conductivity.
We acknowledge P. Bernier for providing us with the SWNT. This work is partly supported by the EEC COMELCAN HPRN-CT-2000-00128 contract