text:
# Chandra X-ray Observations of NGC 4258: Iron Absorption Lines from the Nucleus

## 1 Introduction

The Seyfert 1.9 galaxy NGC 4258 (M 106) is a prime candidate for the study of low luminosity active galactic nuclei. The mass of the central black hole in NGC 4258 is well determined from position and velocity measurements of the H2O maser emission which show a thin, slightly warped, accurately-Keplerian disk of radius 0.12 pc to 0.25 pc in orbit around a $`(3.9\pm 0.1)\times 10^7M_{}`$ black hole (Nakai et al. 1993; Greenhill et al. 1995; Miyoshi et al. 1995; Herrnstein et al. 1999). Wilkes et al. (1995) have found strongly polarized, broad optical emission lines indicative of an obscured active nucleus. Early X-ray observations with Einstein (Cecil et al. 1992; Fabbiano et al. 1992) and ROSAT (Pietsch et al. 1994; Vogler & Pietsch 1999; Cecil et al. 1995) detected soft, extended X-ray emission, but could not penetrate the gas obscuring the nucleus. ASCA observations were the first to detect the active nucleus, which exhibited a power law spectrum of photon index $`\mathrm{\Gamma }1.8`$ and an absorption corrected 2 – 10 keV luminosity of $`L_\mathrm{x}=4\times 10^{40}`$ erg s-1 obscured by a column density of $`N_\mathrm{H}1.5\times 10^{23}`$ cm-2 (Makishima et al. 1994). The X-ray luminosity is a small fraction of the Eddington luminosity, $`L_\mathrm{X}/L_{\mathrm{Edd}}<3\times 10^5`$. Subsequent ASCA observations have shown the nuclear luminosity to be variable on year-long timescales with $`L_\mathrm{x}`$ in the range $`(0.41)\times 10^{41}`$ erg s-1 (Reynolds et al. 2000; Terashima et al. 2002; Ptak et al. 1999) while the photon index and column density have remained approximately constant. BeppoSAX observations revealed 100% continuum variability on a timescale of a $`\mathrm{few}\times 10^5`$ s, and 10 – 20% variability on timescales as short as one hour (Fiore et al. 2001). Interesting variability has also been seen by ASCA in the 6.4 – 6.9 keV iron lines. In 1993, a narrow Fe line was seen at $`6.66_{0.07}^{+0.20}`$ keV with an equivalent width (EW) of $`180_{64}^{+73}`$ eV (Terashima et al. 2002), while in 1996, the Fe line energy had changed to $`6.31_{0.10}^{+0.09}`$ keV with an EW of $`54_{27}^{+25}`$ eV (Terashima et al. 2002). This change in line energy, by 0.35 keV, is significantly larger than the systematic uncertainty in the ASCA calibration (Iwasawa 2003). In a 1999 ASCA observation, a narrow iron line was observed at $`6.45_{0.07}^{+0.10}`$ keV with an EW of $`107_{37}^{+42}`$ eV (Reynolds et al. 2000). An XMM-Newton observation in 2000 constrained the EW of a 6.45 keV iron line to be $`<45`$ eV (Pietsch & Read 2002), indicating that the iron line had weakened significantly since the 1999 ASCA observation (a $`100`$ eV EW narrow iron line would have easily been detected in the XMM-Newton observation).

We have obtained a series of Chandra observations to study various aspects of NGC 4258. At 7.2 Mpc (Herrnstein et al. 1999) $`1\mathrm{}=35`$ pc, and we take the Galactic hydrogen column towards NGC 4258 to be $`N_\mathrm{H}=1.19\times 10^{20}`$ cm-2 (Murphy et al. 1996).

## 2 Observations and Reduction

NGC 4258 was observed by Chandra on four separate occasions with ACIS-S (Table 1). The first observation on 2000-03-08 had a duration of 2.9 ksec, with a 128 row sub-array and alternating 0.4 sec and 0.1 sec frame-times, in order to obtain an approximate count rate and assess the effects of pile-up. The second observation on 2000-04-17 had a duration of 13 ksec with a full frame and 3.2 sec frame-time, primarily to study the “anomalous arms” of NGC 4258 (Wilson et al. 2001, hereafter WYC); the nucleus suffers from the effects of pile-up. The third observation on 2001-05-28 had a duration of 20 ksec with a full frame and 3.2 sec frame-time, again to study the anomalous arms, and the nucleus again suffers from the effects of pile-up. The fourth observation, taken immediately after the third, on 2001-05-29 had a duration of 7 ksec with a 128 row sub-array and a 0.4 sec frame-time to mitigate the effects of pile-up and to study the nucleus of NGC 4258.

The latest gain files were applied to the data and background flares were filtered out by excluding data from those times during which the background (i.e. source-free regions on the S3 chip) count rate was $`>3\sigma `$ above or below the mean background count rate. For the 2.9 ksec duration observation, we used a dead-time correction (DTCOR) of 0.552 to calculate the exposure time of the 0.4 sec frames. The resulting exposure times are listed in Table 1. We shall refer to the 4 observations as the 1.6 ksec, 13 ksec, 20 ksec and 7 ksec exposures respectively. Note that because the 0.4 sec frame-time 1.6 ksec exposure was taken in an alternating mode, the 1.6 ksec of exposure time is spread over the 2.9 ksec duration of the observation. CIAO 2.2.1 with CALDB 2.17 were used, and the spectra were analyzed using XSPEC 11.2 (Arnaud 1996). To obtain spectra of the nucleus, counts were extracted from circular apertures of radius $`1\stackrel{}{\mathrm{.}}5`$ centered on the brightest pixel, and background counts were extracted from an annulus of inner radius $`4\mathrm{}`$ and outer radius $`10\mathrm{}`$ also centered on the nucleus. The spectra were grouped to $`15`$ counts per bin and the $`\chi ^2`$ statistic was used to compare models to the data. The confidence region for a model parameter was determined by varying that parameter until the $`\chi ^2`$ value, minimized by allowing all of the remaining free parameters to vary, was equal to the best fit value plus $`\mathrm{\Delta }\chi ^2`$ where $`\mathrm{\Delta }\chi ^2=2.706`$, 6.635 and 7.879 for the 90% 99% and 99.5% confidence regions for a single interesting parameter, respectively. Unless stated otherwise 90% confidence regions have been used.

When analyzing spectra of the nucleus taken with a 3.2 sec frame-time it is important to correct for the effects of pile-up (the 0.4 sec frame-time observations do not require pile-up correction), and we used the pile-up correction routine of Davis (2001) that is incorporated into the XSPEC software package. We find the grade morphing parameter $`\alpha `$ (the fraction of “good grades” given by $`p`$ piled-up photons is $`\alpha ^{(p1)}`$) to be $`\alpha =0.41_{0.07}^{+0.06}`$ for the 20 ksec data set and $`\alpha =0.43_{0.03}^{+0.03}`$ for the 13 ksec data set.

## 3 X-ray Morphology

The 3.2 s frame-time observations of 2000-04-17 and 2001-05-28 have been combined to produce the grey scale image shown in Fig. 1. Here we concentrate on the nucleus and its immediate environment, and do not discuss the much larger scale “anomalous arms” (see WYC). As noted by WYC and by Pietsch & Read (2002), there is a weaker source $`2\stackrel{}{\mathrm{.}}5`$ southwest of the nucleus. Soft, diffuse emission is seen to extend to the northwest and southeast of the nucleus, and this emission is the beginning of the “anomalous arms”. Comparing the spatial distribution of the 3 – 8 keV counts from the nucleus in the 2001-05-29 observation, which does not suffer from pile-up, with a theoretical point spread function at 6.4 keV, we find that the hard X-ray nucleus is unresolved.

## 4 X-ray Light Curves

X-ray light curves of the nucleus of NGC 4258 in 450 s bins are shown in Fig. 2. The background count rate is negligible ($`<\mathrm{few}\times 10^4`$ cts s-1) and has not been subtracted. The 3.2 sec frame-time data sets have been corrected for the effects of pile-up, using the continuum models given in Section 5.1. The light curves show 10 – 14% rms variability on short timescales and $`>`$20% variability on year-long timescales between the first and last observations (panels a and d, neither of which is affected by pile-up). One of the effects of moderate to severe pile-up is to reduce the sensitivity of the observed count rate to changes in the actual count rate (i.e., as the actual count rate increases the observed count rate increases at a much slower rate). This effect adds a systematic error to our determinations of the un-piled mean count rates (panels b and c) as well as reducing the observed amplitude of short timescale variability. This systematic error is difficult to estimate, but may be as large as a $`\mathrm{few}\times 10`$%.

## 5 X-ray Spectra

### 5.1 Nucleus

#### 5.1.1 Continuum Spectrum

Initially we studied the 2001-05-29 observation (which has an exposure time of 7 ksec and a frame-time of 0.4 sec) because it is the highest S/N observation not affected by pile-up. Two components are required to describe the spectrum of the nucleus: i) a low energy component below $`2`$ keV, such as a 1 keV bremsstrahlung absorbed by the Galactic column, and ii) a heavily absorbed ($`N_\mathrm{H}7\times 10^{22}`$ cm-2) high energy component above $`2`$ keV, such as a hard power law of photon index $`\mathrm{\Gamma }1.4`$ or a very hot bremsstrahlung with a temperature $`>20`$ keV. We then modeled the other data sets in a similar manner, correcting for pile-up as necessary. Table 2 summarizes our models of the continuous spectra of the nucleus for the different observations, and Fig. 3 shows a model fit to the 2001-05-28 data set. For each of these observations, the soft X-ray component is the same while the hard X-ray component declined from 2000 to 2001. During this variation, the intrinsic absorbing column density and photon index of the hard component remained constant (Table 2). The models presented here assume absorption by a single column of neutral gas. If the absorber is ionized then a higher column density may be permitted (see discussion in Sections 6.1.1 and 6.1.2). If a partial covering model is used the best fit has a covering fraction of $`>`$ 99% and parameters very similar to our single absorber fits; the uncovered fraction, $`<`$ 1% of the nucleus, accounts for the observed soft emission. It is also possible for an additional absorption component to be added to our models; $``$ 40% of the nucleus may be absorbed by $`N_\mathrm{H}7\times 10^{22}`$ cm-2 while $``$ 60% of the nucleus may be absorbed by an additional $`N_\mathrm{H}8\times 10^{23}`$ cm-2, and this model requires a softer power law, with $`\mathrm{\Gamma }1.9`$. In addition, an arbitrary amount of complete obscuration may be present that we cannot constrain. Our decision to use a single column density is the most conservative approach and this column density is a lower limit to the true column density.

#### 5.1.2 Iron Emission Lines

We now focus our attention on the spectrum in the region of the Fe K$`\alpha `$ lines. The ratios of the data to the best fitting power law models above 5 keV are shown in Fig. 4, in which no strong iron K$`\alpha `$ emission lines are apparent. If a narrow, neutral iron K$`\alpha `$ emission line at a rest frame energy of 6.4 keV is added to the model, the upper limit to the EW of such a line is ranges between 94 and 887 eV (90% confidence) for the different observations (see Table 3). If the data are grouped to have $`10`$ or $`3`$ counts per bin, there is tentative evidence of iron K$`\alpha `$ emission lines in the first observation and the last two observations at the energies expected for neutral and H-like iron, although these features have low statistical significance.

#### 5.1.3 Iron Absorption Lines

During the 2000-04-17 observation two strong absorption features are seen at NGC 4258 rest frame energies of $`6.41_{0.06}^{+0.05}`$ keV and $`6.87_{0.06}^{+0.09}`$ keV, with equivalent widths of $`166`$ eV and $`215`$ eV, respectively (see Fig. 4 and Table 3; the heliocentric redshift of NGC 4258 is $`z=0.00155`$ (Cecil et al. 1992)). These absorption lines at 6.41 keV and 6.87 keV are both statistically significant at $`>`$ 99.5% confidence (Table 3). An F-test may be used to quantify the improvement to the simple power law model when the absorption features are added. The power law model has $`\chi ^2`$ / degrees of freedom (dof) of 115.8 / 101. Adding a single absorption line at 6.41 keV (the energies of all absorption features discussed here have been left as free parameters) to the power law model gives an improvement of $`\mathrm{\Delta }\chi ^2`$ / $`\mathrm{\Delta }`$ dof = 12.9 / 2, which is significant with 99.5% confidence. Adding a single absorption line at 6.87 keV to the power law model gives an improvement of $`\mathrm{\Delta }\chi ^2`$ / $`\mathrm{\Delta }`$ dof = 7.7 / 2, which is significant with 96.2% confidence, although the absence of the 6.41 keV line is so significant that it affects the continuum fit in the iron line band. If we add the 6.87 keV line to the power law plus 6.41 keV line model we find an improvement with the addition of the 6.9 keV line of $`\mathrm{\Delta }\chi ^2`$ / $`\mathrm{\Delta }`$ dof = 11.6 / 2, which is significant with 99.5% confidence. If we compare the power law model to the model with two absorption lines we find $`\mathrm{\Delta }\chi ^2`$ / $`\mathrm{\Delta }`$ dof = 23.5 / 4, significant at $`>`$ 99.9%. Only three of our observations have sufficient signal-to-noise ratio to detect absorption lines of this strength, and if we assume they are either “cold” lines at $`6.4`$ keV or “hot” lines at $`6.9`$ keV, this give us only six “chances” to see a strong iron absorption feature, so the presence of absorption lines at these energies in this data set is significant with at least 99.5% confidence. Continuing to assume the absorption lines are subject to only the redshift of NGC 4258, the 6.41 keV line is consistent with either of or a blend of the $`n=12`$ resonance lines of Fe XVIII and Fe XIX at 6.43 and 6.46 keV, respectively (the energy expected for Fe XX is 6.50 keV and inconsistent with the data, while ionization states below Fe XVIII do not have an L-shell vacancy). The feature at 6.87 keV is probably the $`n=12`$ resonance line of Fe XXVI at 6.97 keV (the energy expected of Fe XXV is 6.70 keV which is inconsistent with the measured value).

There is no evidence of a strong iron K absorption edge at 7.8 keV (from Fe XVIII; we use the edge energies given by Verner & Yakovlev 1995), 7.9 keV (from Fe XIX), with a limit of $`\tau _{\mathrm{edge}}<0.5`$. We are insensitive to an absorption edge at 9.3 keV (from Fe XXVI).

If the 2000-04-17 observation is divided into two consecutive segments, each 6 ksec in duration, we can attempt to investigate any short timescale variability of the absorption features. The spectra of the iron line region for each of these segments is shown in Fig. 5. A strong absorption line is seen at $`6.9`$ keV in the first segment, but has disappeared in the second segment 6000 sec later. In the first segment, the addition of a narrow absorption line at $`6.95_{0.12}^{+0.22}`$ keV to the continuum model is statistically significant with $`>`$ 99.5% confidence (see Table 3). In the second segment, the addition of an absorption feature at 6.9 keV to the continuum model does not improve the fit at all. If the model of the second segment includes an absorption line at 6.95 keV, identical to that observed in the first segment, we find that removing this line improves the quality of the fit with $`>`$ 95% confidence. These considerations indicate that the narrow absorption feature at 6.9 keV in the first segment is real and that it is significantly weaker in the second segment. The absorption feature at 6.4 keV is seen during both segments; in the first, it is significant at only the $`>`$ 90% confidence level, but is marginally stronger in the second segment where it is statistically significant at $`>`$ 99.5% confidence.

### 5.2 Off-nuclear Source

The spectrum of the off-nuclear source taken from the 20 ksec observation on 2001-05-28 is well described by an absorbed power law. The intrinsic column density is $`N_\mathrm{H}=4.3_{1.6}^{+2.8}\times 10^{21}`$ cm-2, the power law photon index is $`\mathrm{\Gamma }=1.8_{0.4}^{+0.5}`$, and the unabsorbed 2 – 10 keV luminosity is $`L_\mathrm{x}=6.1\times 10^{38}`$ erg s-1 (assuming the source is in NGC 4258). These parameters are consistent with those found by WYC from the observation on 2000-04-17.

## 6 The Fe Absorption Lines

### 6.1 Curve of Growth Analysis

In order to see a strong absorption line, the solid angle of the resonantly scattering material as seen by the central X-ray source must be small, otherwise the absorption line would be “filled in” by emission lines. We thus imagine the absorbing material as a “blob” of gas, and use the strength of an absorption line to estimate the column density of the iron ion responsible for the absorption. Following the method described by Kotani et al. (2000) (also see Spitzer 1978), the equivalent width (EW), $`W_\nu `$, of a resonance line at a frequency $`\nu `$ is given by

$$W_\nu =_0^{\mathrm{}}\left\{1e^{N_{\mathrm{Fe}}\sigma (\nu )}\right\}𝑑\nu [\mathrm{Hz}]$$
(1)
where $`N_{\mathrm{Fe}}`$ is the column density of the iron ion in question and $`\sigma (\nu )`$ is the cross section of that ion at a frequency $`\nu `$. If the ions obey a Maxwell-Boltzmann distribution, then $`\sigma (\nu )`$ has a Voigt profile that depends on the temperature of the blob, being broader at a higher temperature. If the energy distribution of the ions is dominated by turbulence or bulk motion, then the appropriate turbulent or bulk velocity dispersion is used in place of the temperature-dependent velocity dispersion. The oscillator strengths and Einstein coefficients required to determine the cross sections were taken from Verner et al. (1996) for Fe XXVI K$`\alpha `$ and Behar & Netzer (2002) for Fe XVIII and Fe XIX K$`\alpha `$. Equation (1) was integrated for different values of $`N_{\mathrm{Fe}}`$ giving the “curve of growth” ($`W_\nu `$ vs. $`N_{\mathrm{Fe}}`$), shown in Fig. 6 for Fe XXVI (cf. Kotani et al. 2000, their Fig. 3), Fe XVIII and Fe XIX K$`\alpha `$ (the absorption line at $`6.4`$ keV may be a blend of Fe XVIII and Fe XIX, and here we assume that the iron ions are equally divided between Fe XVIII and Fe XIX). From Fig. 6, we note that the large EW of the absorption lines we see ($`166`$$`215`$ eV; Table 3) can correspond to a broad range of iron column densities depending on the kinematic temperature of the blob. For example, the column density of Fe XXVI required to produce a 215 eV EW absorption line decreases from $`N_{\mathrm{Fe}\mathrm{XXVI}}=6\times 10^{21}`$ cm-2 at a temperature of $`kT100`$ keV (equivalent to a $`1\sigma `$ line of sight velocity dispersion of the iron ions of $`\mathrm{\Delta }v4.1\times 10^7`$ cm s-1) to $`N_{\mathrm{Fe}\mathrm{XXVI}}=5\times 10^{18}`$ cm-2 at $`kT=10,000`$ keV (equivalent to a $`1\sigma `$ velocity dispersion of the iron ions of $`\mathrm{\Delta }v4.1\times 10^8`$ cm s-1).

A complication is that, while we see strong absorption lines from a blend of Fe XVIII – Fe XIX and from Fe XXVI, we do not see a strong absorption line from Fe XXV. For a single temperature, collisionally ionized or photoionized plasma, such is not possible (see Arnaud & Raymond 1992; Krolik & Kallman 1984, respectively), so there must be two distinct blobs with different ionization states. Since the absorption line from Fe XVIII – Fe XIX did not vary during the 13 ksec exposure while the Fe XXVI absorption line did vary, it is likely that the $`6.4`$ keV absorption line originates in cooler gas at larger radii. We shall, therefore, assume that the $`6.4`$ keV and $`6.9`$ keV absorption lines originate in spatially distinct blobs. We also assumed that the ion velocities within each absorbing blob are well described by a Maxwell-Boltzmann distribution since the CCD-resolution spectra provide little information about the line profile. This need not be the case, and the absorption lines may be made up of a number of individual components some of which may be saturated while others are unsaturated (equivalently the blob may be made up of sub-blobs of different column densities and velocity dispersions). It is difficult to quantify how these differences would affect our analysis, but the velocity dispersion that we derive can be thought of as a column density weighted mean velocity dispersion of (and within) the sub-blobs.

In the following, we derive some physical parameters of the “blobs”. In particular, we are interested in their opacity in the 2 – 5 keV region, as this is the range of the spectrum from which the equivalent hydrogen absorbing column of $`N_\mathrm{H}=7\times 10^{22}\mathrm{cm}^2`$ was derived, assuming the gas is neutral and of solar abundance (Section 5.1.1).

#### 6.1.1 Photoionized Plasma

We first investigate the blob responsible for the Fe XXVI K$`\alpha `$ absorption line. For $`W_\nu =215`$ eV (Table 3), our curve of growth analysis indicates that $`N_{\mathrm{Fe}\mathrm{XXVI}}>6\times 10^{18}`$ cm-2 (see Fig. 6) and, if the iron abundance is equal to the solar value of $`3.5\times 10^5`$ Fe ions per hydrogen nucleus (Anders & Grevesse 1989), this column corresponds to $`N_\mathrm{H}(\mathrm{ionized})>2\times 10^{23}`$ cm-2. If the gas in this blob is photoionized so that the ionization state of iron is dominated by Fe XXVI, an ionization parameter $`\xi =L_\mathrm{x}/(n_eR^2)3000`$ erg cm s-1 is required (Kallman & McCray 1982), where $`n_e`$ is the electron density and $`R`$ is the distance of the blob from the ionizing source. The plasma temperature corresponding to $`\xi =3000`$ erg cm s-1 is $`T3\times 10^6`$ K, and $``$ 60% of the iron ions are in the form of Fe XXVI, with the remainder being evenly split between Fe XXV and Fe XXVII (Kallman & McCray 1982). The X-ray opacity of such a gas is reduced by a large factor below the Fe K edge (Krolik & Kallman 1984) and hence a column density of $`N_\mathrm{H}(\mathrm{ionized})>2\times 10^{23}`$ cm-2 is consistent with $`N_\mathrm{H}(\mathrm{obs})=7\times 10^{22}`$ cm-2. If we consider instead the stronger 6.9 keV line seen during the first segment of the 2000-04-17 observation the column density requirements are increased by a factor of a few, and are still consistent with $`N_\mathrm{H}(\mathrm{obs})`$.

We now turn our attention to the blob producing the Fe K$`\alpha `$ 6.4 keV absorption feature from a blend of Fe XVIII–Fe XIX K$`\alpha `$. For $`W_\nu =166`$ eV (Table 3) our curve of growth analysis indicates that $`N_{\mathrm{Fe}}>10^{19}`$ cm-2 (see Fig. 6) and, if the iron abundance is equal to the solar value, this corresponds to $`N_\mathrm{H}(\mathrm{ionized})>3\times 10^{23}`$ cm-2. A photoionized plasma with an ionization parameter of $`\xi 100`$ erg cm s-1 has almost all of its iron ions evenly distributed between Fe XVIII and Fe XIX (Kallman & McCray 1982, model 1). The plasma temperature corresponding to $`\xi =100`$ erg cm s-1 is $`T10^5`$ K. From 2 to 5 keV the opacity is $`3`$ times smaller than that of a neutral gas (Krolik & Kallman 1984). Given the approximations inherent to our analysis, the required column density of $`N_\mathrm{H}(\mathrm{ionized})>3\times 10^{23}`$ cm-2 is consistent with the equivalent hydrogen column density for neutral gas ($`N_\mathrm{H}=7\times 10^{22}`$ cm-2 – see Table 2) needed to account for the observed *photoelectric* absorption. This consistency, however, requires a very high kinematic temperature of $`kT>1000`$ keV, which presumably represents bulk motion with a $`1\sigma `$ line of sight velocity dispersion of the iron ions of $`\mathrm{\Delta }v>1.3\times 10^8`$ cm s-1. We conclude that the gas responsible for the observed photoelectric absorption may be the same as that responsible for the 6.4 keV Fe absorption line.

#### 6.1.2 Collisionally ionized Plasma

We first consider the blob responsible for the Fe XXVI K$`\alpha `$ absorption line. If the absorbing gas is collisionally ionized, the ionic fraction of Fe XXVI is maximal at a temperature $`T10^8\mathrm{K}`$ (Arnaud & Raymond 1992), and such a plasma will have negligible opacity in the 2 – 5 keV band (Krolik & Kallman 1984).

Considering next the blob responsible for the $`6.4`$ keV absorption line, the ionic species Fe XVIII – Fe XIX have maximal fractions at a temperature of $`T10^7`$ K (Arnaud & Raymond 1992) and such a plasma will have a negligible opacity at energies $`<2`$ keV, and roughly an order of magnitude lower opacity than neutral gas in the 2 – 5 keV band (Krolik & Kallman 1984). If such a blob is responsible for *both* the observed photoelectric absorption and the $`166`$ eV EW absorption line at 6.4 keV, a kinematic temperature of $`kT1000\mathrm{keV}`$ (equivalent to a $`1\sigma `$ line of sight velocity dispersion of the iron ions of $`\mathrm{\Delta }v1.3\times 10^8\mathrm{cm}\mathrm{s}^1`$) is needed. In this case, the absorbing gas could be in collisional equilibrium at $`10^7\mathrm{K}`$, with the above value of $`\mathrm{\Delta }v`$ representing the turbulent or bulk outflow velocity. Thus, again, it is possible that the gas responsible for the observed photoelectric absorption is the same as that which gives the Fe K$`\alpha `$ absorption line at 6.4 keV.

We emphasize that all other observations of NGC 4258 made by us (Table 3) and others (Section 1) have found either no Fe lines or Fe lines in emission. Thus the spectrum is clearly variable (cf. Terashima et al. 2002). If we adopt the conclusion that the gas responsible for the Fe K$`\alpha `$ 6.4 keV absorption seen on 2000-04-17 is indeed the same as that which causes the photoelectric absorption, then the variability at 6.4 keV would result from variability in the emission line, since the photoelectric absorbing column density appears to be constant (Table 2 and papers cited in Section 1).

#### 6.1.3 Nature of Absorbing Plasma

Since our line of sight to the nucleus grazes the surface of the thin, masing disk of NGC 4258 (inclination $`i=82^{}\pm 1^{}`$; Herrnstein et al. 1999), the absorbing blobs may be part of a disk-wind boundary layer or an outflow. The absence of strong iron K$`\alpha `$ emission is consistent with a hot, low radiative efficiency, advection dominated accretion flow (ADAF; e.g., Narayan & Yi 1994; Abramowicz et al. 1995) with a lower limit to the radius of an outer thin disk of approximately $`100R_G`$ (Reynolds et al. 2000). Such accretion flows may show strong outflows (Blandford & Begelman 1999) or convection (Quataert & Gruzinov 2000; Narayan et al. 2000) and the ionized absorbing “blobs” we have found in NGC 4258 may be inhomogeneities in such a flow with high turbulent velocities (to account for the large equivalent widths of the absorption lines – see Fig. 6). The short-timescale variability of the 6.9 keV absorption line (Section 5.1.3) might result from absorption by a blob at $`10R_G`$ which corresponds to a dynamical timescale of $`6000`$ sec.

The X-ray emission from the nucleus of NGC 4258 may originate from the base of the jet. As gas transitions from the ADAF to the jet it is shock heated, and the X-ray continuum emission may be synchrotron self-Compton emission from the accelerated electrons (Yuan et al. 2002). The mass ejection process may be unstable and the absorption lines we see could be from blobs of cooling thermal gas entrained in the jet outflow.

It is also very interesting to note that similar X-ray absorption features have been seen from the Galactic black hole or neutron star binaries GRO J1655-40 (Ueda et al. 1998; Yamaoka et al. 2001), Circinus X–1 (Brandt & Schulz 2000; Schulz & Brandt 2002), GRS 1915+105 (Kotani et al. 2000; Lee et al. 2002), GX 13+1 (Ueda et al. 2001), MXB 1659-298 (Sidoli et al. 2001), X 1624-490 (Parmar et al. 2002) and X 1254-690 (Boirin & Parmar 2003), all of which have high inclination accretion disks ($`i60^{}`$$`80^{}`$, with the exception of GX 13+1 for which the inclination is not known). If a high inclination is required to produce a strong absorption line, as is suggested by these observations, then NGC 4258 is an important galaxy to study the connection between accretion disks in AGN and Galactic systems.

## 7 Conclusions

Our four Chandra observations of the low luminosity active galactic nucleus of NGC 4258 have shown:

i) The continuum spectrum contains a component below 2 keV, such as 1 keV thermal emission absorbed by the Galactic column, plus a component above 2 keV, which is well described by a hard power law ($`\mathrm{\Gamma }=1.4`$) absorbed by an equivalent hydrogen column of $`N_\mathrm{H}=7\times 10^{22}`$ cm-2.

ii) The soft component of the continuum spectrum does not vary, while the flux of the hard component declined from March/April 2000 to May 2001. During this decline, the intrinsic absorbing column density and photon index of the hard component remained constant (Table 2).

iii) We do not detect an iron K$`\alpha `$ emission line, with the upper limit to the EW of a narrow, neutral iron line ranging between 94 and 887 eV (90% confidence) for the different observations (Table 3). In the second observation (on 2000-04-17) absorption lines are seen at 6.4 and 6.9 keV that are statistically significant with $`>`$ 99.5% confidence. There is evidence that the absorption line at 6.9 keV varied on a timescale of 6000 secs.

iv) We argue that the absorption lines are Fe K$`\alpha `$ resonant absorption, with that at 6.4 keV by Fe XVIII – Fe XIX, and that at 6.9 keV by Fe XXVI ions. A curve of growth analysis has been performed for each line in order to obtain the relationship between the kinematic temperature and the ionic column for the observed equivalent width.

v) Given that the maser disk is viewed at an inclination $`i=82^{}`$, we suggest that the 6.9 keV absorbing gas may be in an inner disk, a disk-wind boundary layer or be thermal gas entrained at the base of the jet.

vi) The gas responsible for the 6.4 keV absorption line may be photoionized or collisionally ionized. In each case, either the ionization temperature is lower than the kinematic temperature or (and more likely) the absorbing matter has a high bulk velocity dispersion ($`1\sigma `$ line of sight velocity dispersion $`\mathrm{\Delta }v1.3\times 10^8`$ cm s-1). We have shown that this gas may be the same as that which is responsible for the observed photoelectric absorption.

Further X-ray observations of the nucleus of NGC 4258 are required to monitor the variability of the iron absorption and emission lines. Such observations should have higher sensitivity and better spectral resolution than those presented here, with the goal of clean separation of absorption and emission lines.

We thank Yuichi Terashima and Chris Reynolds for discussions. This research was supported by NASA through grants NAG 8-1027 and NAG 8-1755.


format: text

text:
# Power Spectrum Normalisation and the Non-Gaussian Halo Model

## 1 Introduction

In the inflation-driven $`\mathrm{\Lambda }`$CDM paradigm, the initial fluctuations are assumed to obey Gaussian statistics. Current observations of the cosmic microwave background (CMB) of galaxy clustering and of galaxy clusters are consistent with this hypothesis (eg. Verde et al., 2001; Gaztañaga and Wagg, 2003). However, models with non-Gaussian initial conditions have been proposed and are not yet ruled out by observations (see Avila-Reese et al., 2003, and reference therein). These models include ones based on inflation theories with non-linearities of a single scalar field (Martin et al., 2000; Gangui et al., 2002) or with multi-fields (Peebles, 1999; Antoniadis et al., 1997; Linde and Mukhanov, 1997; Bernardeau and Uzan, 2002), or on cosmological defects (Turok and Spergel, 1990; Avelino et al., 1998).

In this paper, we study how non-Gaussian initial conditions can affect the determination of $`\sigma _8`$, the normalisation of the matter power spectrum on 8 $`h^1`$ Mpc scales. This normalisation is an essential ingredient to predict the phonemenology of the low redshift universe and has been measured using several methods. The recent CMB measurements with WMAP (Spergel et al., 2003) yield a constraint on the normalisation of $`\sigma _8=0.9\pm 0.1`$ (68%CL). The abundance of galaxy clusters depends strongly on $`\sigma _8`$ and yields values in the range of $`\sigma _8=0.7`$-0.8 (eg. Pierpaoli et al., 2003) for $`\mathrm{\Omega }_m0.3`$. Weak lensing by large-scale structure, or ’cosmic shear’, provides a direct measurement of mass fluctuations and has recently been detected and measured (see Refregier, 2003; Van Waerbeke and Mellier, 2003, for recent reviews). The different cosmic shear surveys yield values of $`\sigma _8`$ in the range of $`0.70.9`$ for the same value of $`\mathrm{\Omega }_m`$. The marginal discrepancies within and between these different techniques may be due to residual systematics, such as uncertainties in the mass-temperature of X-ray clusters, or the calibration of the shear in cosmic shear surveys. They could also be due to physics beyond the standard $`\mathrm{\Lambda }`$CDM model, which has different effects on the various methods and surveys. In particular, non-Gaussianity appears as a prime candidate to explain discrepancies between the different methods. Cluster abundance and cosmic shear indeed both probe non-linear structures at low redshifts, but are sensitive to matter fluctuations on different scales. Non-Gaussianity may thus affect these two methods differently, while leaving the CMB power spectrum essentially unaffected. Reciprocally, it is interesting to establish whether the comparison between the different methods with present and future surveys may constrain non-Gaussianity.

To study the impact of non-Gaussianity on measurements of $`\sigma _8`$, we present a generalisation of the halo model (Ma and Fry, 2000; Seljak, 2000). Building on the work of Robinson and Baker (2000) for the mass function, Ma and Fry (2000) and Koyama et al. (1999) for the halo bias, and Navarro et al. (1996) and Avila-Reese et al. (2003) for the halo profiles, we compute the non-linear power spectrum for arbitrary non-Gaussian models. This allows us to compute the cluster temperature function and the cosmic shear statistics for these models. We then study whether a discrepancy in the determination of $`\sigma _8`$ from the different observational methods can be explained by non-Gaussianity. We also discuss how future surveys can be compared to set constraints on non-Gaussianity. The paper is organised as follows. In §2, we present the different non-Gaussian models and discuss how they are constrainted with current observations. In §3, we describe the non-Gaussian halo model. In §4, we study the effect of non-Gaussianity on the determination of $`\sigma _8`$ with the different methods. Our conclusions are summarised in §5.

## 2 Non-Gaussian Models

The statistics of the post-recombination linear density field $`\rho `$ are usually described by first forming the density contrast $`\delta =(\rho \overline{\rho })/\overline{\rho }`$ with respect to the average matter density $`\overline{\rho }=\rho `$ in the universe. We consider the probability distribution function (PDF) $`p(\delta ,M)`$ of the density contrast $`\delta `$ smoothed with a spherical top-hat of radius $`R`$ corresponding to an enclosed mass $`M=4\pi \overline{\rho }R^3/3`$. It is defined so that

$$𝑑\delta p(\delta ,M)=1,𝑑\delta p(\delta ,M)\delta =\delta =0.$$
(1)
The variance of the distribution is given by

$$𝑑\delta p(\delta ,M)\delta ^2=\delta ^2=\sigma ^2(M)$$
(2)
and can be computed from the linear matter power spectrum $`P_{\mathrm{lin}}(k)`$ using

$$\sigma ^2(M)=\frac{d^3k}{(2\pi )^3}P_{\mathrm{lin}}(k)W^2(k,M),$$
(3)
where $`W(k,M)=3(\mathrm{sin}kR/kR\mathrm{cos}kR)/(kR)^2`$ is the top-hat window function. The linear power spectrum is calculated using the transfer function from Bond and Efstathiou (1984).

To quantify the level of non-Gaussianity of the distribution, Robinson et al. (2000) have considered the excess probability $`\alpha `$ in the $`3\sigma `$ tail over that of a Gaussian PDF, i.e.

$$\alpha (M)=\frac{_{3\sigma (M)}^{\mathrm{}}𝑑\delta p(\delta ,M)}{_{3\sigma (M)}^{\mathrm{}}𝑑\delta p_\mathrm{G}(\delta ,M)}.$$
(4)
Another way of quantifying non-Gaussianity is the skewness

$$\mu _3(M)=\delta ^3=𝑑\delta p(\delta ,M)\delta ^3$$
(5)
To quantify non-Gaussianity from CMB measurements, several authors (eg. Komatsu et al., 2003; Verde et al., 2001) have considered the non-linear coupling $`f_{NL}`$ of the primordial (pre-recombination) potential field $`\mathrm{\Phi }`$ to the Gaussian field $`\mathrm{\Phi }_G`$, defined as

$$\mathrm{\Phi }=\mathrm{\Phi }_G+f_{NL}[\mathrm{\Phi }_G^2\mathrm{\Phi }^2].$$
(6)
This parametrisation is motivated by a class of inflation models. A conversion between $`f_{NL}`$ and the skewness $`\mu _3`$ of the (post-combination) density field can be derived using the results of Matarrese et al. (2000) and Verde et al. (2001). Current constraints from galaxy clusters (Avila-Reese et al., 2003; Robinson et al., 2000) yield $`\alpha 4`$ at 95%CL. As we will see in §4.1, the recent WMAP CMB measurement (Komatsu et al., 2003; Robinson et al., 2000) translates into $`\alpha 1.6`$ at 95%CL for $`M10^{14}M_{}`$.

As a special case, we consider scale free models, where the shape of the distribution is independent of the scale R. In this case, the PDF can be written as

$$p(\delta ,M)=\sigma (M)^1g\left(\sigma (M)^1\delta \right),$$
(7)
where $`g`$ is a function which obeys

$$_0^{\mathrm{}}𝑑xg(x)=1,_0^{\mathrm{}}𝑑xxg(x)=0,_0^{\mathrm{}}𝑑xx^2g(x)=1,$$
(8)
so as to verify equations (1-2). For the Gaussian case,

$$g_\mathrm{G}(x)=(2\pi )^{\frac{1}{2}}e^{x^2/2}.$$
(9)
As an example of a scale invariant distribution, we consider the log-normal distribution which provides an adequate description of a range of non-Gaussian models based on non-linear or multi-field inflation and cosmic strings (Robinson and Baker, 2000). For this distribution,

$$g(x)=\frac{1}{\sqrt{2\pi }(x+x_0)\beta }\mathrm{exp}\left[\frac{\mathrm{ln}^2\left(\frac{x+x_0}{\gamma }\right)}{2\beta ^2}\right],$$
(10)
where

$$x_0=\left(e^{\beta ^2}1\right)^{\frac{1}{2}},\gamma =\left(e^{2\beta ^2}e^{\beta ^2}\right)^{\frac{1}{2}},$$
(11)
so as to enforce equation (8). The parameter $`\beta `$ controls the non-Gaussianity of the PDF and tends towards zero as the distribution becomes Gaussian. Figure 1 shows this distribution for several values of $`\beta `$ and corresponding values of $`\alpha `$.

## 3 Generalised Halo Model

A number of tools exist for modelling the growth of structure since inflation. On the largest scales, fluctuations in density are sufficiently small that their growth can be accurately modelled using linear theory. On small scales, however, the complex nonlinear interactions of gravitational in fall can, at present, only truly be modelled by using N-body simulation. These numerical simulations are however CPU intensive, making the exploration of parameter space prohibitively slow. To circumvent this, semi-analytic tools have been developed to provide fast and accurate fits to numerical simulations.

One such tool is the halo model (Ma and Fry, 2000; Seljak, 2000), which can be used to calculate the full non-linear power spectrum for a given cosmological model and has been shown to be in good agreement with high resolution numerical simulations. The halo model breaks the power spectrum into two parts: one that dominates on large scales and is determined by the clustering of halos of different masses, and a second term that is dominant on small scales and is determined by the density profile of the halos. The halo model thus requires three pieces of information: the mass-function of halos, the bias of halos compared to the background density field and the halo density profile. The halo model as described by Ma and Fry (2000) and Seljak (2000) calculates these three ingredients by assuming that the density perturbations are Gaussian. We show how this can be extended to include models that have non-Gaussian fluctuations.

### 3.1 Mass Function

Robinson and Baker (2000) have proposed a generalisation of the Press and Schechter (1974) formalism which provides a good fit to N-body simulations for a class of non-Gaussian models. In their formalism, the fraction of objects of a given mass $`M`$ or larger that collapse is calculated by integrating the PDF of the density fluctuations as

$$F(>M)=A_{\delta _c}^{\mathrm{}}P(\delta ,M)𝑑\delta ,$$
(12)
where $`A`$ is a normalisation factor and $`\delta _c`$ is a critical over density which can be computed for arbitrary cosmologies using the results of Lacey and Cole (1993) and Eke et al. (1996). The mass function, defined as the number density of halos per unit mass, is then given by

$$\frac{dn}{dM}=\frac{\overline{\rho }}{M}\frac{d\nu }{dM}f(\nu ),$$
(13)
where

$$f(\nu )=\left|\frac{dF}{d\nu }\right|$$
(14)
and

$$\nu =\frac{\delta _c}{\sigma (M)}.$$
(15)
The normalisation $`A`$ is chosen to ensure that all the matter is accounted for, so that

$$\frac{1}{\rho }_0^{\mathrm{}}𝑑MM\frac{dn}{dM}=_0^{\mathrm{}}𝑑\nu f(\nu )=1.$$
(16)
For a scale invariant PDF (Eq. ), it is easy to show that equations (12) and (14) reduce to

$$f(\nu )=Ag(\nu ),$$
(17)
with $`A=\left[_0^{\mathrm{}}𝑑xg(x)\right]^1`$. This simple form reveals that the mass function at mass M provides a measure of the normalised PDF $`g(x)`$ at $`x=\nu =\delta _c/\sigma (M)`$. This assumes of course that the generalised Press-Schechter formalism is valid.

For the Gaussian case (Eq. ), we obtain the usual result

$$f_\mathrm{G}(\nu )=\sqrt{\frac{2}{\pi }}e^{\nu ^2/2},$$
(18)
with the usual normalisation factor $`A=2`$ prescribed by Press and Schechter (1974).

For the log-normal distribution (Eq. ), we obtain

$$f(\nu )=\frac{A}{\sqrt{2\pi }(\nu +x_0)\beta }\mathrm{exp}\left[\frac{\mathrm{ln}^2(\frac{\nu +x_0}{\gamma })}{2\beta ^2}\right],$$
(19)
with

$$A=\left[\mathrm{erfc}\left(\frac{\mathrm{ln}(x_0/\gamma )}{\sqrt{2}\beta }\right)\right]^1,$$
(20)
where $`\mathrm{erfc}`$ is the complementary error function. Figure 2 shows the resulting differential mass function $`dn/dM`$ for this model at $`z=0`$, for several values of the non-Gaussian parameter $`\alpha `$. Non Gaussianity with $`\alpha >1`$ tends to increase the number of halos with $`M10^{14}M_{}`$ and, to a smaller degree, to deplete halos with masses smaller than this limit. The high sensitivity of high mass halos to non-Gaussianity is expected since these halos probe the positive tail of the distribution. Throughout this paper we assume a standard $`\mathrm{\Lambda }`$CDM model, with $`\mathrm{\Omega }_m=0.3,\mathrm{\Omega }_\mathrm{\Lambda }=0.7,h=0.7,n=1\mathrm{and}\mathrm{\Gamma }=0.21`$. For this figure we adopted $`\sigma _8=0.77`$, but we will choose other values for this parameter depending on the data set considered.

### 3.2 Halo Bias

The bias parameter in the halo model describes how halos cluster compared to the underlying linear density field. Generalising the formalism of Mo and White (1996), Koyama et al. (1999) showed that the bias parameter for a general scale invariant PDF is given by

$$b(\nu )=1\frac{1}{\delta _c}\left[1+\frac{d\mathrm{ln}f(\nu )}{d\mathrm{ln}\nu }\right].$$
(21)
Note that this expression satisfies the condition $`𝑑\nu f(\nu )b(\nu )=1`$, which is required to ensure that the linear power spectrum is recovered on large scales (see discussion in Seljak, 2000). For the scale invariant case, $`f(\nu )`$ is proportional to the PDF (see §3.1). The halo bias can therefore be thought of as a measure of the derivative of the PDF in this case. The measurement of halo bias, using the halo correlation function for instance, can thus be used to probe the PDF independently of the mass function.

For the Gaussian case (Eq. ), this reduces to the usual Mo & White result

$$b_\mathrm{G}(\nu )=1+\frac{\nu ^21}{\delta _c}.$$
(22)
More accurate bias functions (eg. Jing, 1998) have been devised to provide better fits to N-body simulations, but they cannot be readily generalised to non-Gaussian models. Since we are only interested in the effect of non-Gaussianity relative to the Gaussian case, the above expression suffices for our purposes.

For the log-normal distribution, the bias function is

$$b(\nu )=1+\frac{1}{\delta _c}\frac{\nu }{\nu +x_0}\left[\frac{1}{\beta ^2}\mathrm{ln}\left(\frac{\nu +x_0}{\gamma }\right)\frac{x_0}{\nu }\right],$$
(23)
which does indeed go to Eq.  as $`\beta `$ tends towards zero. Figure 3 shows the bias as a function of mass for several values of $`\alpha `$. For skew-positive PDFs ($`\alpha >1`$), objects with low (high) masses are more (less) clustered compared to the Gaussian case.

### 3.3 Halo profile

We model the density profile of the halos as

$$\rho (r)=\rho _su(r/r_s),u(x)=x^\mu (1+x)^{\mu 3}$$
(24)
where $`r_s`$ is a characteristic radius. The inner slope $`\mu `$ must be less than 3 for convergence, and is equal to 1 in the case of a Navarro, Frenk & White (1996, NFW) profile and 1.4 in the case of a Moore et al. (1998) profile. We fix the normalisation $`\rho _s`$ so that the density contrast within the viral radius $`r_v`$ is $`\delta _v=200`$, i.e. so that $`M=4\pi \overline{\rho }\delta _vr_v^3/3=4\pi _0^{r_v}𝑑rr^2\rho (r)`$. This yields $`\rho _s=\overline{\rho }\delta _vc^3U^1(c)/3`$ where we have defined the compactness parameter $`c=r_v/r_s`$, and

$$U(c)=_0^c𝑑xx^2u(x)=\frac{c^{3\mu }}{3\mu }{}_{2}{}^{}F_{1}^{}(3\mu ,3\mu ;4\mu ;c)$$
(25)
where $`{}_{2}{}^{}F_{1}^{}`$ is the hypergeometric function which can be evaluate numerically (Press et al., 1996). In the case of the NFW profile ($`\mu =1`$), this function reduces to $`U(c)=\mathrm{ln}(1+c)c/(1+c)`$.

For the $`\mathrm{\Lambda }`$CDM Gaussian model, a good fit to the non-linear power spectrum is given for an NFW profile ($`\mu =1`$) with the concentration parameter is chosen as (Cooray, 2000; Cooray et al., 2000)

$$c(M,z)=c_{}(z)\left[\frac{M}{M_{}(z)}\right]^{h(z)}$$
(26)
where $`c_{}(M,z)10.3(1+z)^{0.3}`$ and $`h(z)0.24(1+z)^{0.3}`$. Here $`M_{}(z)`$ is the non-linear mass scale defined by $`\nu (M_{},z)=1`$.

Avila-Reese et al. (2003) have used numerical simulations to measure the halo profiles in non-Gaussian models. For a model with $`\alpha 3`$, they find an increase in the compactness parameter $`c`$ of about 20% compared to the Gaussian case for halos with masses in the range $`10^{11}\mathrm{to}10^{15}\mathrm{h}^1\mathrm{M}_{}`$. They also find the inner slope to be steeper ($`\mu 1.6`$) compared to the Gaussian NFW case ($`\mu =1`$).

More simulations would be needed to derive new scaling relations for halos in general non-Gaussian models. For our purposes, we use the use the NFW profile as a reference, but we also study the effect of changing the compactness and inner slope to those found by Avila-Reese et al. (2003).

### 3.4 Non-linear Power Spectrum

Ma and Fry (2000) and Seljak (2000) showed that the halo model provides an accurate prescription to compute the non-linear power spectrum for Gaussian models. In this formalism, the non-linear power spectrum $`P(k)`$ at a given redshift is written as

$$P(k)=P_1(k)+P_2(k),$$
(27)
where the 1-halo term is

$$P_1(k)=_0^{\mathrm{}}𝑑M\frac{dn}{dM}\left[\frac{\stackrel{~}{\rho }(k,M)}{\overline{\rho }}\right]^2$$
(28)
and the 2-halo term is

$$P_2(k)=\left[_0^{\mathrm{}}𝑑M\frac{dn}{dM}b(M)\frac{\stackrel{~}{\rho }(k,M)}{\overline{\rho }}\right]^2P_{\mathrm{lin}}(k).$$
(29)
In these expressions, $`\stackrel{~}{\rho }(k,M)`$ is the radial Fourier transform of the density profile (Eq. ) of a halo of mass $`M`$. The details of our implementation of the halo model can be found in Refregier and Teyssier (2002).

Figure 4 shows the effect of non-Gaussian initial conditions on the non-linear power spectrum. Non-Gaussianity is seen to affect transition scales $`0.1k3h`$ Mpc-1 between the linear and the non-linear regime. The linear ($`k0.1h`$ Mpc-1) and strongly non-linear ($`k3h`$ Mpc-1) are not significantly affected. Figure 5 shows the effect of varying the inner slope $`\mu `$ and the concentration parameter $`c`$ as indicated by the non-Gaussian simulations of Avila-Reese et al. (2003). Their suggested alteration of the concentration parameter by 20% (which we applied to all masses) tends to only slightly increase the power spectrum on small scales ($`k3h`$ Mpc-1). On the other hand, the alteration of the inner profile $`\mu `$ by 60% increases the power spectrum significantly on these scales. Note however, that the mass dependence of $`c(M)`$ which we used as a reference (Eq. ) was chosen to reproduce the non-linear power spectrum in numerical simulations for an NFW profile ($`\mu =1`$). A selfconsistent fit to the non-linear power spectrum from non-Gaussian simulations would be needed to establish whether this effect is indeed so pronounced. In the following, we keep $`c`$ and $`\mu `$ unchanged from the Gaussian NFW case, but remain cautious in the interpretation of our model for the power spectrum on small scales.

## 4 Power Spectrum Normalisation

### 4.1 Cosmic Microwave Background

Using the recent measurements from the WMAP mission, Spergel et al. (2003) found $`\sigma _8=0.9\pm 0.1`$ (68%CL). When combined with other CMB missions, galaxy surveys and Lyman alpha surveys, the authors find $`\sigma _8=0.84\pm 0.04`$ (68%CL). The CMB power spectrum provides a 2-point statistics of the matter fluctuations in the linear regime and is thus not affected by non-Gaussianity111Strictly speaking, non-Gaussianity may slightly bias the algorithms used to measure the CMB power spectrum, which often assume exact Gaussian statistics. In this paper, we ignore this small effect. This is also true when it is combined with measurement of the galaxy power spectrum on linear scale, while assuming that the galaxy bias is scale independent. As a result, the determination of $`\sigma _8`$ based on fits to the CMB power spectrum (with and without galaxy surveys) is insensitive to non-Gaussianity.

The measurement of higher order statistics of CMB maps, however, does currently present tight constraints on the level of primordial non-Gaussianity. Using the WMAP measurements, Komatsu et al. (2003) were able to place limits on the non-Gaussian coupling defined in equation (6) of $`58<f_{NL}<134`$ (95%CL). According to Matarrese et al. (2000, figure 2), the upper limit corresponds to an upper limit on the skewness of the post-recombination density field (Eq. ) on scale of $`M=10^{14}M_{}`$ of $`\mu _38.0\times 10^2`$, for a $`\mathrm{\Lambda }`$CDM model with $`\sigma _8=1`$. For the log-normal distribution (Eq. ), this corresponds to $`\alpha 1.6`$ at the 95% confidence level. Note that this conversion assumes that the scale dependence of the non-Gaussianity is that of models given by equation (6).

### 4.2 X-ray Cluster Counts

Counting X-ray clusters is a direct probe of the mass function and provides a sensitive measure of $`\sigma _8`$. The main difficulty is to estimate the masses of X-ray selected clusters. The most reliable method to do so consists of measuring the X-ray temperature of the cluster and to use a mass-temperature ($`M`$-$`T`$) scaling relation. Using an M-T relation derived form hydrodynamical simulations, several groups have derived a value of $`\sigma _8`$ around 0.9 from measurements of the X-ray temperature function (see Pierpaoli et al., 2001, and reference therein). More recently, the use of an observational $`M`$-$`T`$ relation led to a revision of $`\sigma _8`$ in the range $`0.7`$-$`0.8`$ (see eg. Pierpaoli et al., 2003; Seljak, 2001). The uncertainty in the $`M`$-$`T`$ relation is, at present, the largest source of uncertainty in cosmological parameter estimates using this method.

In order to study how the determination of $`\sigma _8`$ using this method is affected by non-Gaussianity, we consider the compilation of cluster temperature measurements presented by Pierpaoli et al. (2001). The catalogue consists of temperature measurements for 38 clusters taken using ROSAT and ASCA. The cluster temperatures range from $`kT=3.5`$ to $`12`$ keV, with those with temperatures greater than $`7.5`$ keV forming a volume limited sample. The resulting cumulative X-ray temperature function $`n(>kT)`$ from this compilation is shown in figure 6. To compute the cluster masses, Pierpaoli et al. (2001) have used the $`M`$-$`T`$ relation,

$$\frac{M(T,z)}{10^{15}h^1M_{}}=\left(\frac{T}{T_{}}\right)^{\frac{3}{2}}\left(\mathrm{\Delta }_cE^2\right)^{\frac{1}{2}}\left(1\frac{2\mathrm{\Omega }_\mathrm{\Lambda }}{\mathrm{\Delta }_c}\right)^{\frac{3}{2}},$$
(30)
where $`T`$ is the temperature in keV, and $`E^2=\mathrm{\Omega }_m(1+z)^3+\mathrm{\Omega }_\mathrm{\Lambda }+\mathrm{\Omega }_k(1+z)^2`$. The parameter $`\mathrm{\Delta }_c`$ is the mean over density inside a virial radius and can be computed using the fitting functions given by Pierpaoli et al. (2001). In a later work, Pierpaoli et al. (2003) fit the observed temperature function using this relation and find $`\sigma _8=0.77_{0.04}^{+0.05}`$ (68%CL).

Figure 6 also shows the predicted temperature function for different values of the non-Gaussian parameter $`\alpha `$. For the masses corresponding to $`kT1`$ keV (corresponding to $`M3\times 10^{13}h^1M_{}`$), a skew-positive PDF ($`\alpha >1`$) enhances the cluster counts (see Figure 2), the effect being more pronounced for larger temperatures. As a result, the determination of $`\sigma _8`$ from the X-ray cluster counts will be affected by non-Gaussianity. Figure 7 shows the value of $`\sigma _8`$ derived from a fit to the temperature function as a function of $`\alpha `$. Following Seljak (2001), we considered a fit to $`n(>kT)`$ at a single temperature $`kT`$, which we choose to be 4, 7, and 10 keV, in turn. For $`kT7`$ keV and for the Gaussian case ($`\alpha =1`$), we find $`\sigma _80.77`$ in agreement with Pierpaoli et al. (2003) and Seljak (2001). The value of $`\sigma _8`$ required to fit the temperature function decreases as $`\alpha `$ increases. The effect is more pronounced for larger cluster temperatures, as expected.

### 4.3 Cosmic Shear

Weak lensing by large scale structure, or cosmic shear, provides a direct measure of the mass distribution in the universe (see Bartelmann and Schneider, 2001; Refregier, 2003; Van Waerbeke and Mellier, 2003, for recent reviews). This effect is based on the shear that intervening large-scale structures induce on the images of background galaxies. Several groups have now performed accurate measurements of cosmic shear (Jarvis et al., 2003; Hamana et al., 2002; Brown et al., 2003; Hoekstra et al., 2002; Refregier et al., 2002; Van Waerbeke et al., 2002; Bacon et al., 2002), which have been used to set constraints on cosmological parameters.

The 2-point cosmic shear statistics can be described by the weak lensing power spectrum (see Bacon et al. (2000) for conventions)

$$C_{\mathrm{}}=\frac{9}{16}\left(\frac{H_0}{c}\right)^4\mathrm{\Omega }_m^2𝑑\chi \left[\frac{h(\chi )}{ar(\chi )}\right]^2P(\frac{\mathrm{}}{r},\chi ),$$
(31)
where $`\mathrm{}`$ is the multipole moment, $`\chi `$ is the comoving radius, $`r`$ is the comoving angular-diameter distance, $`a`$ is the expansion parameter, and $`H_0`$ and $`\mathrm{\Omega }_m`$ are the present value of the Hubble constant and matter density parameter, respectively. It can be computed for non-Gaussian models by using the halo-model approximation for the non-linear power spectrum $`P(k,\chi )`$ at a given comoving radius $`\chi `$, given in equation (27). The radial weight function is given by

$$h(\chi )=2_\chi ^{\mathrm{}}𝑑\chi ^{}n(\chi ^{})\frac{r(\chi )r(\chi ^{}\chi )}{r(\chi ^{})},$$
(32)
where $`n(\chi )`$ is the red shift distribution of galaxies and is normalised as $`𝑑\chi n(\chi )=1`$.

Figure 8 shows the weak lensing power spectrum for several values of $`\alpha `$. The galaxy redshift distributions was taken to be $`n(z)z^2\mathrm{exp}((z/z_0)^{1.5})`$ with $`z_0=z_m/1.41`$ and the median redshift was chosen to $`z_m=1`$ as suitable for several cosmic shear surveys. The cosmological parameters are those for our fiducial model given in §3.1 with $`\sigma _8=1.0`$. The linear power spectrum is also plotted for comparison and reveals that non-linear corrections are very important on scales $`\mathrm{}10^3`$. Non-Gaussianity has the most pronounced effect on intermediate scales $`10^2\mathrm{}10^4`$, which is consistent with its effect on the 3D power spectrum (see figure 4). For $`\alpha >1`$, non-Gaussianity tends to increase the lensing power spectrum on these scales.

In practice, it is often more convenient to consider the shear variance $`\sigma _\gamma ^2=|\gamma |^2`$ in randomly place circular cells of radius $`\theta `$. This is related to the shear power spectrum by

$$\sigma _\gamma ^2=\frac{1}{2\pi }_0^{\mathrm{}}𝑑\mathrm{}C_{\mathrm{}}|W_{\mathrm{}}|^2$$
(33)
where $`W_{\mathrm{}}=2J_1(\mathrm{}\theta )/(\mathrm{}\theta )`$ is the Fourier transform of the cell aperture. The shear variance of the different non-Gaussian models is shown in figure 9. Also shown are recent measurements of the shear variance from different surveys with median galaxy redshifts in the range $`0.8z_m1`$ (Brown et al., 2003; Refregier et al., 2002; Van Waerbeke et al., 2002; Bacon et al., 2002). The measurements are in rough agreement with each other and with the Gaussian $`\mathrm{\Lambda }`$CDM case. Non-Gaussian models with $`\alpha >1`$ produce larger shear variances on scales $`\theta 10^{}`$, i.e. on scales to which most surveys are sensitive. This is not surprising since this limit corresponds to the transition between the linear ($`\theta 10^{}`$) and non-linear ($`\theta 10^{}`$) regimes for the shear statistics.

We therefore expect non-Gaussianity to affect the determination of $`\sigma _8`$ from current cosmic shear surveys. Figure 10 shows the value of $`\sigma _8`$ derived from the shear variance measured at different scales and with $`z_m=1`$, as a function of the non-Gaussian parameter $`\alpha `$. We see that if the universe is non-Gaussian with, say, $`\alpha 3`$, the cosmic shear surveys (which assumed $`\alpha `$=1) would have overestimated $`\sigma _8`$ by about 5%. This shift is approximately unchanged if the galaxy redshift is set to $`z_m`$=0.6, as appropriate for shallower cosmic shear surveys (Jarvis et al., 2003; Hoekstra et al., 2002), which tend to find lower values of $`\sigma _8`$ than deeper ones (Brown et al. 2003; Refregier et al. 2002; Van Waerbeke et al. 2002; Bacon et al. 2002; see however Hamana et al. 2002). The discrepancy between shallower and some of the deeper surveys can therefore not be explained by non-Gaussianity.

As discussed in §3.3, the effect of non-Gaussianity depends on our assumption regarding the halo profile. We find that the bias in $`\sigma _8`$ is twice as pronounced if we use the modified profile suggested by Avila-Reese et al. (2003) instead of the NFW profile. Further high-resolution numerical simulations would be needed to resolve this theoretical uncertainty.

As expected, the effect of the non-Gaussian parameter, $`\alpha `$, depends on the scale being studied. On large angular scales, $`\theta >10^{}`$, $`\alpha `$ has a very weak effect on the shear variance, which is as expected since these scales are still in the linear regime. The intermediate scales, $`1^{}\theta 10^{}`$, are the most affected by our non-Gaussian PDF, as can also be seen in figure 10. This is also consistent with the behaviour of the weak lensing power spectrum (figure 8) since these intermediate scales correspond to $`10^2\mathrm{}10^3`$. On scales smaller than $`\theta =1^{}`$, the variance from models with different $`\alpha `$ values no longer diverge from each other. This makes measurements on small angular scales less sensitive to variations in the PDF. On even smaller scales, corresponding to $`\mathrm{}10^4`$, the shift in the power spectrum changes sign, suggesting that on these scales the spread in shear variance of our models will decrease. Since lensing surveys do not probe these angular scales, we do not investigate this further.

### 4.4 Combined Constraints

As discussed in the previous sections, the determination of $`\sigma _8`$ from the CMB with and without galaxy surveys is insensitive to non-Gaussianity. On the other hand, cluster counts and cosmic shear are affected by non-Gaussianity, as can be seen in Figures 7 and 10. For $`\alpha 5`$, we find that the error in $`\sigma _8`$ made by ignoring non-Gaussianity can be approximated by

$$\mathrm{\Delta }\sigma _8=\sigma _8^G\sigma _8s(\alpha 1),$$
(34)
where $`\sigma _8^G`$ is the value derived by (wrongly) assuming a Gaussian model (i.e. corresponding to $`\alpha =1`$ in these figures). The slope $`s`$ and intercept values $`\sigma _8^G`$ are listed in table 1 for cluster counts and cosmic shear at different temperatures and angular scales. As is apparent from the figures and in the table, cluster counts are more sensitive ($`s`$ larger by a factor of about 2) than cosmic shear to non-Gaussianity. However, as noted above, the use of the Avila-Reese et al. (2003) halo profile instead of the NFW profile would tend to increase the cosmic shear slope $`s`$ by about a factor of 2, which would then nearly match that for cluster abundance. While this needs to be established by further numerical simulations, we will use the NFW slopes $`s`$ for cosmic shear as a reference but remember that they are lower limits.

If we assume that the PDF shape is scale invariant over all scales, we have seen (see §4.1) that the recent WMAP skewness measurement imposes a constraint of $`\alpha 1.3`$ (68%CL). This limit corresponds to errors in $`\sigma _8`$ of $`\mathrm{\Delta }\sigma _80.02`$ and 0.01, for cluster counts and cosmic shear respectively. These are significantly smaller to both the current uncertainties ($`\mathrm{\Delta }\sigma _8`$0.05 and 0.1, respectively) and the discrepancy between the different methods and measurements ($`\mathrm{\Delta }\sigma _80.3`$). Thus, with this assumption, non-Gaussianity can not explain the internal and external discrepancies between these methods.

Note however that the above conclusion rests on an extrapolation from scales probed by WMAP ($`\mathrm{}1000`$ corresponding to $`25h^1`$ Mpc comoving at decoupling) down to scales relevant for clusters and cosmic shear, 0.1 to 5.0 comoving $`h^1`$Mpc. If we drop the assumption that the PDF is invariant over this large range of scales, the WMAP measurement no longer constrains the PDF relevant for these two latter techniques. In this case, we can attempt to reconcile the cluster and cosmic shear normalisations using non-Gaussianity. We find that this would require $`\beta 0.4`$. This level of negatively skewed PDF does not have a corresponding $`\alpha `$ parameter since the distribution does not have a 3$`\sigma `$ tail in this case. Such a negative skewness would be inconsistent with the observed shape of the temperature function (see figure 6) and of the shear variance (see figure 9).

Finally, we look at the level of accuracy future surveys need to reach in order to improve upon the current constraints on non-Gaussianity. Assuming scale invariance, we have shown that the current $`1\sigma `$ constraints set by WMAP correspond to uncertainties of $`\mathrm{\Delta }\sigma _80.02`$ for cluster counts and $`\mathrm{\Delta }\sigma _80.01`$ for cosmic shear. At present, cluster count measurements have $`1\sigma `$ uncertainties of $`\mathrm{\Delta }\sigma _80.05`$ which are dominated by an uncertainty in the normalisation $`T_{}`$ of the $`M`$-$`T`$ relation of about 10%. This means that if the error on $`\sigma _8`$ from clusters is reduced by a factor of about 2 then they will be comparable to that arising from our lack of knowledge of the statistics of the density field. This will be achieved if the uncertainty in $`T_{}`$ is reduced to about 4% in future surveys (see eg. Pierpaoli et al., 2003). Similarly, current lensing surveys of about 10 deg2 have $`1\sigma `$ errors of $`\mathrm{\Delta }\sigma _80.06`$ (eg. Van Waerbeke et al., 2002) and scale as $`\sqrt{A}`$ where $`A`$ is the survey area. Lensing surveys will thus have to increase in sensitivity by a factor of about 6, and thus have areas of the order of 400 deg2, before their measurement of $`\sigma _8`$ becomes sensitive to non-Gaussianity. Note, however, that it may be possible to probe primordial non-Gaussianity by studying higher order statistics of cosmic shear.

## 5 Conclusions

We have investigated the impact of non-Gaussian conditions on the determination of $`\sigma _8`$ via various observational methods. For this purpose, we have generalised the halo model to compute cluster statistics and the non-linear matter power spectrum. In passing, we noted that, for scale free non-Gaussian models, the halo mass function and bias provide a measure of the linear density PDF and of its derivative, respectively. The measurement of both of these statistics, therefore, provides an independent measure of the PDF and a test of its scale invariance.

The determination of $`\sigma _8`$ from the CMB power spectrum and from galaxy surveys in the linear regime is not sensitive to non-Gaussianity. On the other hand, the cluster X-ray temperature function tends to be enhanced at large temperatures ($`T3`$ keV) for skew-positive ($`\alpha >1`$) non-Gaussianity. If the primordial fluctuations are non-Gaussian and this is ignored, the determination of $`\sigma _8`$ from current X-ray cluster surveys would be overestimated by $`\mathrm{\Delta }\sigma _80.04(\alpha 1)`$. This depends weakly on the typical temperature of the clusters. Cosmic shear statistics from current surveys probe the non-linear part of the matter power spectrum and are thus also sensitive to non-Gaussianity. Using the generalised halo model, we find that the matter power spectrum is enhanced by a factor of 2 on scales $`0.1k3h\mathrm{Mpc}^1`$ for skew-positive non-Gaussianity of $`\alpha =5`$. As a result, the cosmic shear power spectrum is also enhanced by a factor of 2 on scales $`10^2\mathrm{}10^4`$. The error in $`\sigma _8`$ from current cosmic shear surveys if non-Gaussianity is ignored is $`\mathrm{\Delta }\sigma _80.02(\alpha 1)`$. This behaviour depends weakly on the angular scale $`\theta `$ and on the galaxy redshift $`z_m`$. We also find that the nonlinear power spectrum, and therefore cosmic shear surveys may be sensitive to the halo profile. This, in turn, may be sensitive to the primordial density field and needs to be studied in more detail via high-resolution N-body simulations (see Avila-Reese et al., 2003).

The strongest upper limits on primordial non-Gaussianity are provided by the recent WMAP measurements and correspond to $`\alpha 10.6`$ (95% CL). Within this limit and assuming scale invariance, non-Gaussianity cannot explain the differences between the $`\sigma _8`$ values derived from different cosmic shear surveys and X-ray cluster catalogs. Dropping the assumption of scale invariance from the scales probed by WMAP to those probed by cluster counts and shear, non-Gaussianity would tend to decrease the derived values of $`\sigma _8`$ for both cosmic shear and cluster counts and thus would need to be very pronounced to explain the discrepancy. Specifically, a highly negatively skewed PDF with $`\beta 0.4`$ (which does not have a corresponding $`\alpha `$) would be needed to resolve the discrepancy, but would not be compatible with the observed shape of the cluster temperature function and of the cosmic shear 2-point function. A comparison between future X-ray surveys with a two fold improvement in cluster mass calibration and future cosmic shear surveys with 400 deg2 will be required to constrain non-Gaussianity on small scales with a precision matching that of the current CMB constraints on larger scales. Our results suggest that the discrepancies are either due to systematics in one or several of the methods or to non-standard physics other than non-Gaussianity.

## Acknowledgements

The authors are grateful to Jerry Ostriker, Tarun Saini, Jochen Weller, Sarah Bridle and Romain Teyssier for useful discussions. AA was supported by a PPARC student grant and AR was supported in Cambridge by an Advanced PPARC Fellowship and a Wolfson College Fellowship.


text:
# Large scale simulations of lattice gauge theory: key NSF project of China and recent developments

## 1 INTRODUCTION

Lattice gauge theory (LGT) has been accepted as the most powerful tool for investigating non-perturbative strong interactions. Chinese physicists at Institute of High Energy Physics (IHEP), Institute of Theoretical Physics (ITP) and Peking Univ. (PKU) in Beijing, Nankai Univ. in Tianjin, Chengdu Univ. of Science and Technology and Sichuan Univ. in Chengdu, Zhejiang Univ. (ZJU) in Hangzhou, and Zhongshan Univ. (ZSU) in Guangzhou have been involved in the study of LGT since early 80’s. However, earlier investigations were mainly analytical ones (for a review, see Refs.), due to limited computational facilities.

Thanks to the rapid development of high performance supercomputers in China in late 90’s, success of the Symanzik improvement program, and increasing support from NSF (National Science Foundation), more and more Chinese lattice physicists are doing numerical simulations.

## 2 COMPUTING FACILITIES

There have been many commercial supercomputers installed at research institutions and universities in China. The Beowulf clusters, however, are becoming more and more popular.

The ZSU group built a cluster of 10 PC type computers. The major difference in our computers from one likely to be found in a home or business is that each is equipped with two CPUs. This allows us to roughly double our processing power without the extra overhead cost for extra cases, power supplies, network cards, etc. Specifically, we have installed two Pentium 3 processors in each motherboard.

The Legend Group built a cluster, named DeepComp 1800/LSSC-II, which is installed at Academy of Mathematics and System Science, Beijing, and also open to the Chinese lattice community. It consists of Myrinet and 512CPU (Xeon-2GHz). Its $`Rmax`$ and $`Rpeak`$ are 1.046 Tflops and 2.048 Tflops respectively. It was the number 43 of top 500 supercomputers and the 3rd fastest in Asia in 2002.

## 3 NATIONAL SCIENCE FOUNDATION

The Chinese NSF Committee was established in 1986, and it is an organization directly affiliated to the State Council for the management of the National Natural Science Fund. As shown in Tab. 1, these years have seen increasing funding. Each General Project is funded for 100 $``$ 200K yuan/3 years (1 Yuan=1/8USD); Each Project for Distinguished Young Scientists is funded for 0.8M yuan/4 years; Each Key Project is funded for 1 $``$ 2M yuan/4 years.

Chinese lattice physicists have got a lot of support from the NSF, including the fund for Distinguished Young Scientists and the Key project.

## 4 PLAN

In the recently funded Key NSF Project, we have planned to do large scale simulations of lattice QCD, using the parallel supercomputing facilities in China. We plan to develop new numerical methods and study the following hot topics: new hadrons such as glueballs and hybrid mesons, scattering of hadrons, topology of QCD vacuum, transition from the quark confinement phase to quark-gluon plasma phase, quantum instantons and quantum chaos. These investigations will provide useful information for the frontier experiments in particle physics.

## 5 RECENT PROGRESS

### 5.1 Improved actions and Hamiltonians

The idea of Symanzik improvement is to reduce the lattice spacing errors by adding new terms to the Wilson actions. The pursuit of the Symanzik program, combined with tadpole improvement, has recently led to significant progress in LGT, opening the possibility of approaching continuum physics on coarse and small lattices. Chinese lattice physicists also contributed actively to the improvement program: the development of the next to nearest neighbor improved quark action and improved Hamiltonians for the quarks and gluons.

### 5.2 Glueballs and hybrid mesons

Lattice computation of the glueball masses has a long history. Early investigations used the unimproved Wilson gluon action and Kogut-Susskind gluon Hamiltonian. Recent calculations used the tadpole improved action on anisotropic lattices. In Ref., the IHEP group performed a detailed study of the continuum behavior of the glueball $`J=0`$, $`2`$, $`4`$ operators and computed the corresponding spectrum. In Ref., the ZSU group performed the calculation of the $`0^{++}`$, $`2^{++}`$ and $`1^+`$ glueball masses, at $`\beta `$ = 2.2, 2.4, 2.6, 3.0, and 3.2, using the improved gluon action on the anisotropic lattices. In comparison with Ref., we have data at weaker coupling and larger lattices, so that the extrapolation to the continuum limit might be more reliable. Our results for the glueball masses are $`M_G(0^{++})=1737(\pm 52)(\pm 107)`$, $`M_G(2^{++})=2409(\pm 35)(\pm 103)`$ and $`M_G(1^+)=2946(\pm 63)(\pm 142)`$.

There have been several quenched lattice calculations of hybrid meson masses using respectively (1) the Wilson quark action, (2) the SW improved quark action, (3) the NRQCD action, (4) the LBO action, and (5) the improved KS quark action. In those simulations, configurations were generated using the Wilson gluon action.

In Ref., the ZSU group employed both improved gluon and quark actions on the anisotropic lattice, which should have smaller systematic errors, and should be more efficient in reducing the lattice spacing and finite volume effects. We obtained $`2013\pm 26\pm 71`$ MeV for the mass of the $`1^+`$ hybrid meson in the light quark sector, and $`4369\pm 37\pm 99`$MeV in the charm quark sector; the mass splitting between the $`1^+`$ hybrid meson in the charm quark sector and the spin averaged S-wave charmonium was estimated to be $`1302\pm 37\pm 99`$ MeV.

### 5.3 QCD at finite chemical potential

The Hamiltonian formulation of lattice QCD doesn’t encounter the complex action problem. In Refs., the ZSU group developed a Hamiltonian approach at finite chemical potential $`\mu `$ and at zero temperature obtained reasonable results in the strong coupling regime. A first order chiral phase transition was observed at some $`\mu _C`$. Quantum chaos and quantum instantons were also investigated; These phenomena might be relevant for the rich phase structure.

### 5.4 Hadron scattering

The lattice collaboration in Beijing studied scattering of hadrons using tadpole improved clover action on coarse anisotropic lattices. The scattering lengths of $`\pi \pi `$ in the $`I`$=2 channel and $`KN`$ in the $`I`$=1 channel were calculated within quenched approximation. Comparison with the new experiment and Chiral Perturbation Theory was made and good agreements were found.

### 5.5 Chiral properties of overlap quarks

The ZJU group studied the chiral symmetry relation and scaling of the overlap fermions. Divergence of quenched chiral logarithm was found in very light quark regime.

### 5.6 Algorithms

In Ref., the ZSU group developed the molecular dynamics equations for doing HMC simulations with dynamical clover quarks and proposed an even-odd preconditioning method for the quark determinant. In Ref., the algorithm was implemented.

The PKU group developed an algorithm for Ginsparg-Wilson quarks. The ZSU group extended an algorithm for computing the thermodynamical quantities of LGT with Kogut-Susskind quarks, including the chiral limit, to QCD.

The ZSU group developed an algorithm, named MC Hamiltonian to study the excited states of LGT and tested in the scalar models.

Some works mentioned above were done with the Laval Univ. and Kentucky Univ. groups.


text:
# The visual orbits of the spectroscopic binaries HD 6118 and HD 27483 from the Palomar Testbed Interferometer

## 1 Introduction

Following Michelson’s suggestion and his success in measuring diameters of Jupiter satellites with an optical interferometer (Michelson, 1891), Schwarzschild managed to resolve a number of double stars with his own grating interferometer (Schwarzschild, 1896). Not long afterward, Anderson determined the first visual orbit of Capella (Anderson, 1920), later improved by Merrill (1922) with the same instrument — the rotating interferometer attached to the 100-inch telescope at the Mount Wilson Observatory and originally used by Michelson. In the same experiment, Merrill also attempted to resolve a number of spectroscopic binaries but without success.

The modern generation of interferometers developed in the 80’s and 90’s finally allows to resolve and subsequently determine visual orbits of spectroscopic binaries in a fairly routine manner. Up to date about twenty spectroscopic binaries have had their orbits determined with optical interferometers (for a review see Quirrenbach, 2001). The Palomar Testbed Interferometer (PTI) itself has been successfully used to determine a number of visual orbits of double-lined spectroscopic binaries including RS Cvn (Koresko et al., 1998), $`\iota `$ Peg (Boden et al., 1999), 64 Psc (Boden et al., 1999), 12 Boo (Boden, Creech-Eakman, & Queloz, 2000), BY Dra (Boden & Lane, 2001) and Gliese 793.1 (Torres et al., 2002). The other instruments that have been successfully used in binaries studies include the Mark-III (Shao et al., 1988) and the Navy Prototype Interferometer (Armstrong et al., 1998).

The capability to determine visual orbits of spectroscopic binaries is an instrumental achievement very important for stellar astronomy. It allows to derive all geometric orbital elements, precise magnitude differences between the components and their angular diameters. This information combined with spectroscopic and photometric data yields colors, luminosities, masses and distance to the system. For binaries with resolved components it also provides their radii and temperatures. It is noteworthy that in some cases the derived component masses are already at the $`12`$ percent precision level sufficient to constraint modern main and post-sequence stellar models (see e.g. Torres et al., 2002). Significant progress in this area is expected when the new generation of optical interferometers offering baselines longer than a hundred of meters, e.g. CHARA (McAlister et al., 2000), and/or large apertures e.g. the Keck Interferometer (Colavita & Wizinowich, 2000) or the Very Large Telescope Interferometer (Glindemann et al., 2000) become fully operational. Still, observing spectroscopic binaries with currently available optical interferometers is a sensible goal as the potential of these instruments and the reservoir of accessible targets has not been depleted. Interferometric observations of currently accessible spectroscopic binaries combined with other observables deliver a complete set of physical parameters for the components in a much shorter period of time than astrometric observations of visual binaries (which have orbital periods measured in years) and in contrast to eclipsing binaries can be successfully carried out for systems with arbitrary orbital inclinations.

In Spring 2001 we have undertaken at the PTI an observing program to determine the orbits of resolvable spectroscopic binaries from the Batten catalog (Batten, Fletcher, & MacCarthy, 1989, 1997). Here we report our first results — the visual orbits of two double-lined binaries, HD 6118 one of the most eccentric spectroscopic binaries and HD 27483 a spectroscopic binary in the Hyades open cluster, which were observed with the PTI in K band (2.2 $`\mu `$m) band in 2001 and 2002.

## 2 Interferometric observable

### 2.1 Theory

From the theoretical standpoint optical interferometry differs little if at all from its radio counterpart. The fundamental observable from a two aperture interferometer is the fringe visibility (also called coherence function or just visibility) represented by a (usually complex) number $`\mathrm{\Gamma }=V\mathrm{e}^{i\varphi }`$ where $`V`$ is its amplitude and $`\varphi `$ phase. The relation between visibility and the source intensity distribution $`I(\mathrm{\Delta }\alpha ,\mathrm{\Delta }\delta )`$ is given (under source incoherence and small-field approximations) by the van Cittert-Zernike Theorem (Born & Wolf, 1999; Thompson, Moran, & Swenson, 2001)

$$\mathrm{\Gamma }(u,v)=I(\mathrm{\Delta }\alpha ,\mathrm{\Delta }\delta )\mathrm{e}^{2\pi \mathrm{i}(u\mathrm{\Delta }\alpha +v\mathrm{\Delta }\delta )}\mathrm{\Delta }\alpha \mathrm{\Delta }\delta ,$$
(1)
where $`(\mathrm{\Delta }\alpha ,\mathrm{\Delta }\delta )`$ are the coordinates in the tangent plane of the sky and $`(u,v)=𝐁_{}/\lambda `$ are the coordinates of the projected baseline vector ($`𝐁`$ represents the spatial separation of the apertures) in the same plane and expressed in units of the observing wavelength $`\lambda `$.

The amplitude of the visibility has dimensions of power and corresponds to the amount of power that the interferometer measures. Real instruments typically produce the normalized (by the total power received from the source), dimensionless visibility amplitude, $`0V1`$. For our further considerations two model visibility amplitudes are particularly useful: that of a uniform disk (i.e. a single star), and that of a binary system comprised of two model disks. It easy to show that these are given by respectively (see Boden, 2000)

$$V_{disk}^2=\left(\frac{2J_1(\pi \theta B_{}\lambda )}{\pi \theta B_{}\lambda }\right)^2$$
(2)
where $`B_{}=𝐁_{}`$ is the length of the projected baseline vector and

$$V_{binary}^2=\frac{V_1^2+r^2V_2^2+2rV_1V_2\mathrm{cos}(2\pi 𝐁_{}\mathrm{\Delta }𝐬/\lambda )}{(1+r)^2}$$
(3)
where $`V_1,V_2`$ are the visibilities of individual stars given by (2), $`r`$ is the brightness ratio at the observing wavelength $`\lambda `$ ($`r=P_2/P_1`$ where $`P_1,P_2`$ are the total powers of the binary components at the given wavelength) and $`\mathrm{\Delta }𝐬=(\mathrm{\Delta }\alpha ,\mathrm{\Delta }\delta )`$ is the separation vector between the primary and the secondary in the plane tangent to the sky.

### 2.2 Practice

The Palomar Testbed Interferometer (PTI) developed by the Jet Propulsion Laboratory, California Institute of Technology, for NASA to test interferometric techniques, is a long-baseline infrared interferometer located at Palomar Observatory, California. It consists of three fixed 40 cm apertures that combined pairwise offer baselines between 86 and 110 m. It operates in two wavelength ranges K ($`2.2\mu `$m) and H ($`1.6\mu `$m). Full details of the instrument’s architecture and operation can be found in Colavita et al. (1999). The data reduction and calibration procedures for the instrument are described in Boden, Colavita, van Belle, & Shao (1998); Colavita (1999); Lane, Colavita, Boden, & Lawson (2000).

Optical interferometers typically produce only square of visibility amplitude, $`V^2`$. The reason for this is that the phase of the measured visibility is usually corrupted by the atmosphere. In addition, we are able to construct unbiased $`V^2`$ estimators (Colavita, 1999; Mozurkewich et al., 1991). In the case of PTI data, the square of visibility amplitude is estimated as the average of the visibilities from all the channels of a given band (K, 2.0-2.4 $`\mu `$m; H, 1.4-1.8 $`\mu `$m) and corresponds to the SNR-weighted mean wavelength, $`\lambda `$. The remaining step that finally results in a reliable $`V^2`$ for an astronomical source is the calibration process.

Ideally, an interferometer’s response to a point source would be $`V^2=1`$. However since generally speaking a real instrument is not perfect and does not operate in a perfect environment, the response to a point source, called system visibility, is $`V_{sys}^2<1`$ and this value changes temporally (varying seeing index being one of the contributors) and spatially (e.g. with apparent source position). Thus the calibration process involves interlacing target observations with calibrator observations and using $`V_{sys}^2`$ to calibrate the target observation in a simple but effective way (Boden, Colavita, van Belle, & Shao, 1998; Mozurkewich et al., 1991)

$$V_{calibrated}^2=V_{measured}^2/V_{sys}^2$$
(4)
Choosing an appropriate set of calibrators is an important part of the interferometric observation. The choice of calibrators is dictated by a few, mostly common sense requirements: (1) a calibrator should have similar brightness to the target (although it should not be too faint if the target happens to be faint for a given interferometer) (2) it should be located in the sky not far from the target (3) finally, an ideal calibrator should resemble a point source as closely as possible. The last requirement is dictated by the estimation of the system visibility through

$$V_{sys}^2=V_{calmeasured}^2/V_{calexpected}^2$$
(5)
where $`V_{calmeasured}^2`$ is a measured visibility and $`V_{calexpected}^2`$ is an expected (model) visibility of a calibrator. From the above equation, it is clear that a calibrator should have as simple model visibility as possible to avoid introducing any model dependent systematic errors in the estimated system visibility. It can be realized quite well using single stars. Their model visibility can be described using equation (2) where the only unknown is the diameter of the star. With few exceptions, the diameters of stars are not know from direct observations. Hence to minimize the errors in the estimated system visibility due to errors in the (estimated) diameters of the calibrators (note that $`\mathrm{d}V_{sys}/\mathrm{d}\theta 0`$ as $`\theta 0`$), it is best to use calibrators which are as unresolved as possible (or alternatively use stars which have diameters very well determined). For the purpose of the calibration process, one can determine angular diameters of the calibrators by fitting the black body model to their archival photometric data. Assuming that the calibrators are essentially unresolved, such a procedure (i.e. errors in the estimated angular diameters) will not corrupt the calibration process.

It is sufficient to have two calibrators. However since more than 1/2 of all stars are in binary or multiple systems, one will sooner or later find out that one of the calibrators is a previously unknown binary. Hence a set of three calibrators is a practical and safe choice. From experience with PTI it follows that it is best to observe targets and calibrators in interleaved, short period ($`10`$ minutes) intervals to properly address any sky position-dependent and temporal variations in the system visibility. For the purpose of the orbit determination of a spectroscopic binary, it is beneficial to observe a target over 2-3 hour time span and gather $`V^2`$ measurements corresponding to different projections of the separation vector $`\mathrm{\Delta }𝐬`$ onto the baseline vector $`𝐁_{}`$ (conf. equation (3)). Also a final set of $`V^2`$ measurements should correspond to a fairly complete orbital phase coverage. An incomplete coverage may cause significant problems in assessing if the best-fit solution is indeed the correct orbital solution (see next section).

In the end, the observer is served by the reduction pipeline with the following six numbers $`\{t_i,V_i^2,\sigma _i,\lambda _i,u_i,v_i\}`$ where $`t_i`$ is the time of observation, $`V_i^2`$ is the calibrated visibility amplitude squared, $`\sigma _i`$ its error, $`\lambda _i`$ is the mean wavelength for the observation, and $`(u_i,v_i)`$ are the components of the projected baseline vector. While a (very small) error in $`t_i`$ is completely unimportant for further considerations, $`\{\lambda _i,u_i,v_i\}`$ are known with finite precision and a proper error analysis should take their uncertainties into account. Also to perform the calibrations, the diameters of the calibration stars are necessary. Since these are known with finite precision, their uncertainties should also be taken into account.

## 3 Data modeling

The interferometric observations of spectroscopic binaries share the description of the binary motion with regular astrometry of visual binary stars. The motion of the primary with respect to the secondary is thus described by the following equations (see e.g. Kovalevsky, 1995; Kamp, 1967)

$$\mathrm{\Delta }𝐬(t)=(\mathrm{\Delta }\alpha (t),\mathrm{\Delta }\delta (t))=\kappa a\left(𝐏(\mathrm{cos}E(t)e)+𝐐\sqrt{1e^2}\mathrm{sin}E(t)\right),$$
(6)
$$𝐏=𝐥\mathrm{cos}\omega +𝐦\mathrm{sin}\omega ,𝐐=𝐥\mathrm{sin}\omega +𝐦\mathrm{cos}\omega ,$$

$$𝐥=(\mathrm{cos}\mathrm{\Omega },\mathrm{sin}\mathrm{\Omega }),𝐦=(\mathrm{cos}i\mathrm{sin}\mathrm{\Omega },\mathrm{cos}i\mathrm{cos}\mathrm{\Omega })$$

where $`E=E(t)`$ is the eccentric anomaly given by the Kepler equation $`Ee\mathrm{sin}E=2\pi (tT_\mathrm{p})/P`$, $`P`$ is the orbital period, $`a,e,i,\omega ,\mathrm{\Omega },T_p`$ are the standard Keplerian elements—the semi-major axis, the eccentricity, the inclination, the longitude of pericenter, the longitude of ascending node and the time of periastron, and $`\kappa `$ is the parallax. The above equations can be used together with the equation (3) to model the observed visibilities. The typical approach is the least-squares fit with respect to the following set of parameters $`\{r,\widehat{a},P,e,i,\omega ,\mathrm{\Omega },T_p\}`$ where $`\widehat{a}=\kappa a`$ is the apparent semi-major axis. Obtaining the set of parameters that truly correspond to the global minimum of $`\chi ^2`$ is more tedious than in the case of radial velocity or astrometric data. From the relative orbit alone there are a few possible solutions for the pair of angles $`\omega ,\mathrm{\Omega }`$. Namely $`\omega ,\mathrm{\Omega }`$ as well as $`\omega \pm \pi ,\mathrm{\Omega }\pm \pi `$. Also, since the visibility itself depends on the scalar product $`𝐁_{}\mathrm{\Delta }𝐬`$, it is invariant under the rotations of the separation vector $`\mathrm{\Delta }𝐬`$ by $`\pi `$. Effectively, $`\chi ^2`$ has many local well defined minima and in order to find the global minimum one needs to carry out the least-squares fit on a grid of initial values of the model parameters that covers the entire parameter space.

Spectroscopic binaries that are observed with optical interferometers have usually well determined orbits from radial velocity (RV) measurements. It substantially reduces the computational effort to obtain the best-fit parameters as the spectroscopic orbit already supplies most of them (e.g. it enables to resolve ambiguity in $`\omega `$ but not in $`\mathrm{\Omega }`$). In fact it is useful to combine the fit to radial velocities and visibilities into one procedure by defining $`\chi ^2`$ as

$$\chi ^2=\underset{i=1}{\overset{N}{}}(V_i^2\widehat{V}_i^2)^2/\sigma _{V_i^2}^2+\underset{i=1}{\overset{N}{}}(RV_i\widehat{RV}_i)^2/\sigma _{RV_i}^2$$
(7)
where $`\sigma _{V^2},\sigma _{RV}`$ are the measurement errors and $`\widehat{V}^2,\widehat{RV}`$ denote model observables, and performing the fit over the entire set of parameters including those intrinsic to the radial velocity model. To perform the fits, we employ the Levenberg-Marquard algorithm (see e.g. Press, Teukolsky, Vetterling, & Flannery, 1992). The formal errors of the fitted parameters can easily be computed within the least-squares fit formalism. However one cannot forget about the additional systematic errors that come from several sources (1) the uncertainty in the projected baseline $`(u,v)`$ (2) the uncertainty in the mean wavelength $`\lambda `$ (3) the uncertainties in the diameters, $`\theta _j`$, of the calibrators (which affect the modeled visibility through equations (4)) and the components of the binary (which affect the modeled visibility through equations (3)). For the purpose of this work we adopt that following conservative, based on our experience with the PTI, estimates of these uncertainties (1) 0.01 percent in $`(u,v)`$, (2) 1 percent in $`\lambda `$ and (3) 10 percent in the calibrator and binary components diameters.

## 4 HD 6118

HD 6118 ($`\sigma `$ Psc, HR 291, HIP 4889) is a double-lined spectroscopic binary with an orbital period of 81.13 days, a spectral type of both components of B9.5V, $`V=5.47`$ mag, $`(BV)=0.05`$ mag and $`(UB)=0.18`$ mag (Cowley, Cowley, Jaschek, & Jaschek, 1969). Its binary nature was determined by Campbell (1918). The spectroscopic orbit was derived by Belserene (1947) based on 28 radial velocity measurements obtained in the years 1944-45 with the Lick Observatory Mills 3-prism spectrograph. The average error of these velocities is 5 km/s and it is the only set of RVs published for this star. HD 6118 has an orbital eccentricity of 0.89. The distance to the star is $`127\pm 11`$ pc from the Hipparcos parallax measurement of $`7.87\pm 0.68`$ mas (Perryman et al., 1997).

We observed HD 6118 in 2002 and collected 69 $`V^2`$ measurements in K band (see Table 4). Th best-fit solution to archival radial velocity (RV) and our $`V^2`$ measurements is shown in Table 2. It is characterized by an rms of $`4.8`$ km/s and $`3.6`$ km/s in the RVs of the primary and secondary and $`0.066`$ in $`V^2`$ (see Fig. 1-2). The reduced $`\chi ^2`$ of the combined RV/astrometric solution is 1.06. The derived physical parameters of the primary and secondary are shown in Table 3. The masses of the components are reasonably well determined as $`2.65\pm 0.27M_{}`$ and $`2.36\pm 0.24M_{}`$. The accuracy of the parallax, $`8.86\pm 0.07`$ mas, is ten times better than from Hipparcos but consistent with Hipparcos determination at $`1.5`$ of its formal error. The masses of the components accurate at the $`10`$ percent level do not allow us to perform any challenging tests of stellar evolution models. Nevertheless, they can be used to estimate the age of HD 6118. To this end, we used the theoretical isochrones from Bertelli et al. (1994) and the K-band photometric data for HD 6118 from the 2MASS catalog (Cutri et al., 2003). Even though we do not have any abundance information for the star, the error in the age is dominated by the errors in the component masses. From Figure 4 it follows that the age is 160-200 Myr (depending on the assumed abundance) with a rather large error of about 100-150 Myr. Finally, as can be seen in Figure 3, the apparent separation between the components is always too large to allow for an eclipse (the estimated angular diameters of the components are $`0.15`$ mas while the smallest orbital separation is about 0.8 mas).

## 5 HD 27483

HD 27483 (HR 1358, HIP 20284) is a double-lined spectroscopic binary in the Hyades open cluster It has an orbital period of 3.06 days, a spectral type for both components of F6V, $`V=6.16`$ mag, $`(BV)=0.46`$ mag and $`(UB)=0.02`$ mag (Johnson, MacArthur, & Mitchell, 1968). The radial velocity orbit was determined by Northcott and Wright (1952) and improved by Mayor and Mazeh (1987). The distance to the system is $`45.9\pm 1.8`$ pc from the Hipparcos parallax measurement of $`21.8\pm 0.85`$ mas.

We observed HD 27483 in 2001 and 2002 and collected 81 $`V^2`$ measurements in K band (see Table 7). Unfortunately, HD 27483 is only partially resolved with PTI because its semi-major axis is only 1.2 mas. In such a case there is a strong correlation between the apparent semi-major axis, $`\widehat{a}`$, and the brightness ratio, $`r`$. In effect, it is not possible to determine both parameters reliably (see also Koresko et al., 1998). Therefore, in order to derive the inclination of the system we proceeded as follows. One can assume that the brightness ratio is close to 1, since the mass ratio of both stars is very close to 1. Also, one can use the parallax measurement from Hipparcos, $`\kappa `$ , to additionally constrain the least-squares fit (through $`\widehat{a}=a\kappa `$). With such a setup ($`r=1`$, $`\kappa =21.8`$ mas), we fitted for the remaining orbital elements. The best-fit solution (see Table 5) to the archival radial velocity (RV) from Mayor & Mazeh (1987) and our $`V^2`$ measurements is characterized by the rms of $`1.1`$ km/s and $`1.8`$ km/s in the RVs of the primary and secondary and $`0.045`$ in $`V^2`$ (see Fig. 5-6). The reduced $`\chi ^2`$ of the combined RV/astrometric solution is 1.3. The derived physical parameters of the primary and secondary are shown in Table 5. Their errors include systematic contribution from the assumed brightness ratio (10 percent in $`r`$) and the parallax ($`0.85`$ in $`21.8`$ mas). The masses of the components are $`1.38\pm 0.13M_{}`$ and $`1.39\pm 0.13M_{}`$. Unfortunately, we cannot independently determine the distance to Hyades using the current set of $`V^2`$ measurements for HD 27483. The star is a perfect candidate for future measurements with CHARA (that will offer a much longer baseline) since for a 3-day period binary it is easy to acquire a complete orbital phase coverage in a relatively short period of time. However, our mass determination of the components and the Hipparcos distance (corresponding to $`M_V=3.6`$ mag for HD 27483 A and B if $`r=1`$) place the stars in the mass-luminosity diagram very close to the prediction based on current theoretical isochrones for the Hyades cluster (see Fig. 4 from Lebreton, Fernandes, & Lejeune, 2001).

## 6 Summary

We have resolved two double-lined spectroscopic binaries, HD 6118 and HD 27483, with the Palomar Testbed Interferometer. The data collected with PTI in 2002-2003 in the K band allow us to determine astrometric orbits and when combined with the radial velocity measurements also to derive all physical parameters of the systems. The masses of the components of HD 6118 are $`2.65\pm 0.27M_{}`$ and $`2.36\pm 0.24M_{}`$ and the distance to the system is $`112.9\pm 0.9`$ pc. Using the theoretical isochrones from Bertelli et al. (1994), we determine the age of HD 6118 AB as approximately 160-200 Myr. HD 27483 is a double-lined spectroscopic binary in the Hyades open cluster. The masses of its components are $`1.38\pm 0.13M_{}`$ and $`1.39\pm 0.13M_{}`$. Unfortunately, the system is only partly resolved with PTI and hence we are unable to reliably determine its apparent semi-major axis and thus the distance. However, our measurement of the component masses and the distance to the star from the Hipparcos catalog place the binary in the mass-luminosity diagram very close to the theoretical prediction by the current models for the Hyades cluster (Lebreton, Fernandes, & Lejeune, 2001).

The data presented in this paper were obtained at the Palomar Observatory using the Palomar Testbed Interferometer, which is supported by NASA contracts to the Jet Propulsion Laboratory. Science operations with PTI are possible through the efforts of the PTI Collaboration and excellent observational work of Kevin Rykoski. This research has made use of software produced at the Michelson Science Center, California Institute of Technology and the SIMBAD database, operated at CDS, Strasbourg, France. M.K. gratefully acknowledges the support of NASA through the Michelson fellowship program.