text:
# FIRST YEAR WILKINSON MICROWAVE ANISOTROPY PROBE (WMAP) OBSERVATIONS: TESTS OF GAUSSIANITY

## 1 INTRODUCTION

The Gaussianity of the primordial fluctuations is a key assumption of modern cosmology, motivated by simple models of inflation. Statistical properties of the primordial fluctuations are closely related to those of the cosmic microwave background (CMB) radiation anisotropy; thus, a measurement of non-Gaussianity of the CMB is a direct test of the inflation paradigm. If CMB anisotropy is Gaussian, then the angular power spectrum fully specifies the statistical properties. Recently, Acquaviva et al. (2002) and Maldacena (2002) have calculated second-order perturbations during inflation to show that simple models based upon a slowly-rolling scalar field cannot generate detectable non-Gaussianity. Their conclusions are consistent with previous work (Salopek & Bond, 1990, 1991; Falk et al., 1993; Gangui et al., 1994). Inflation models that have significant non-Gaussianity may have some complexity such as non-Gaussian isocurvature fluctuations (Linde & Mukhanov, 1997; Peebles, 1997; Bucher & Zhu, 1997), a scalar-field potential with features (Kofman et al., 1991; Wang & Kamionkowski, 2000), or “curvatons” (Lyth & Wands, 2002; Lyth et al., 2002). Detection or nondetection of non-Gaussianity thus sheds light on physics of the early universe.

Many authors have tested Gaussianity of CMB anisotropy on large angular scales ($`7^{}`$) (Kogut et al., 1996; Heavens, 1998; Schmalzing & Gorski, 1998; Ferreira et al., 1998; Pando et al., 1998; Bromley & Tegmark, 1999; Banday et al., 2000; Contaldi et al., 2000; Mukherjee et al., 2000; Magueijo, 2000; Novikov et al., 2000; Sandvik & Magueijo, 2001; Barreiro et al., 2000; Phillips & Kogut, 2001; Komatsu et al., 2002; Komatsu, 2001; Kunz et al., 2001; Aghanim et al., 2001; Cayón et al., 2002), on intermediate scales ($`1^{}`$) (Park et al., 2001; Shandarin et al., 2002), and on small scales ($`10^{}`$) (Wu et al., 2001; Santos et al., 2002; Polenta et al., 2002). So far there is no evidence for significant cosmological non-Gaussianity.

Most of the previous work only tested the consistency between the CMB data and simulated Gaussian realizations without having physically motivated non-Gaussian models. They did not, therefore, consider quantitative constraints on the amplitude of possible non-Gaussian signals allowed by the data. On the other hand, Komatsu et al. (2002), Santos et al. (2002), and Cayón et al. (2002) derived constraints on a parameter characterizing the amplitude of primordial non-Gaussianity inspired by inflation models. The former and the latter approaches are conceptually different; the former does not address how Gaussian the CMB data are or the physical implication of the results. In this paper, we adopt the latter approach, and constrain the amplitude of primordial non-Gaussianity in the WMAP 1-year sky maps.

Some previous work all had roughly similar sensitivity to non-Gaussian CMB anisotropy at different angular scales, because the number of independent pixels in the maps are similar, i.e., $`40006000`$ for COBE (Bennett et al., 1996), QMASK (Xu et al., 2002), and MAXIMA (Hanany et al., 2000) sky maps. Polenta et al. (2002) used about $`4\times 10^4`$ pixels from the BOOMERanG map (de Bernardis et al., 2000), but found no evidence for non-Gaussianity. The WMAP provides about $`2.4\times 10^6`$ pixels (outside the Kp0 cut) uncontaminated by the Galactic emission (Bennett et al., 2003a), achieving more than one order of magnitude improvement in sensitivity to non-Gaussian CMB anisotropy.

This paper is organized as follows. In § 2, we describe our methods for measuring the primordial non-Gaussianity using the cubic (bispectrum) statistics and the Minkowski functionals, and present the results of the measurements of the WMAP 1-year sky maps. Implications of the results for inflation models and the high-redshift cluster abundance are then presented. In § 3, we apply the bispectrum to individual frequency bands to estimate the point-source contribution to the angular power spectrum. The results from the WMAP data are then presented, and also comparison among different methods. In § 4, we present summary of our results. In Appendix A, we test our cubic statistics for the primordial non-Gaussianity using non-Gaussian CMB sky maps directly simulated from primordial fluctuations. In Appendix B, we test our cubic statistic for the point sources using simulated point-source maps. In Appendix C, we calculate the CMB angular bispectrum generated from features in a scalar-field potential.

## 2 LIMITS ON PRIMORDIAL NON-GAUSSIANITY

### 2.1 The WMAP 1-year Sky Maps

We use a noise-weighted sum of the Q1, Q2, V1, V2, W1, W2, W3, and W4 maps. The maps are created in the HEALPix format with $`nside=512`$ (Górski et al., 1998), having the total number of pixels of $`12\times nside^2=3,145,728`$. We do not smooth the maps to a common resolution before forming the co-added sum. This preserves the independence of noise in neighboring pixels, at the cost of complicating the effective window function for the sky signal. We assess the results by comparing the WMAP data to Gaussian simulations processed in identical fashion. Each CMB realization draws a sample from the $`\mathrm{\Lambda }`$CDM cosmology with the power-law primordial power spectrum fit to the WMAP data (Hinshaw et al., 2003b; Spergel et al., 2003). The cosmological parameters are in Table 1 of Spergel et al. (2003) (we use the best-fit “WMAP data only” parameters). We copy the CMB realization and smooth each copy with the WMAP beam window functions of the Q1, Q2, V1, V2, W1, W2, W3, and W4 (Page et al., 2003a). We then add independent noise realizations to each simulated map, and co-add weighted by $`N_{\mathrm{obs}}/\sigma _0^2`$, where the effective number of observations $`N_{\mathrm{obs}}`$ varies across the sky. The values of the noise variance per $`N_{\mathrm{obs}}`$, $`\sigma _0^2`$, are tabulated in Table 1 of (Bennett et al., 2003b).

We use the conservative $`Kp0`$ mask to cut the Galactic plane and known point sources, as described in Bennett et al. (2003a), retaining 76.8% of the sky (2,414,705 pixels) for the subsequent analysis. In total 700 sources are masked on the 85% of the sky outside the Galactic plane in all bands; thus, the number density of masked sources is 65.5 $`\mathrm{sr}^1`$. The Galactic emission outside the mask has visible effects on the angular power spectrum (Hinshaw et al., 2003b). Since the Galactic emission is highly non-Gaussian, we need to reduce its contribution to our estimators of primordial non-Gaussianity. Without foreground correction, both the bispectrum and the Minkowski functionals find strong non-Gaussian signals. We thus use the forground template correction described in section 6 of Bennett et al. (2003c) to reduce foreground emission to negligible levels in Q, V, and W bands. The method is termed as an ”alternative fitting method”, which uses only the Q, V, and W band data. The dust component is separately fitted to each band without assuming spectrum of the dust emission (3 parameters). We assume that the free-free emission has $`\nu ^{2.15}`$ spectrum, and the synchrotron has $`\nu ^{2.7}`$ spectrum. The amplitude of each component in Q band is then fitted across three bands (2 parameters).

### 2.2 Methodology

#### 2.2.1 Model for Primordial Non-Gaussianity

We measure the amplitude of non-Gaussianity in primordial fluctuations parametrized by a non-linear coupling parameter, $`f_{\mathrm{NL}}`$ (Komatsu & Spergel, 2001). This parameter determines the amplitude of a quadratic term added to the Bardeen curvature perturbations $`\mathrm{\Phi }`$ ($`\mathrm{\Phi }_\mathrm{H}`$ in Bardeen (1980)), as

$$\mathrm{\Phi }(𝐱)=\mathrm{\Phi }_\mathrm{L}(𝐱)+f_{\mathrm{NL}}\left[\mathrm{\Phi }_\mathrm{L}^2(𝐱)\mathrm{\Phi }_\mathrm{L}^2(𝐱)\right],$$
(1)
where $`\mathrm{\Phi }_\mathrm{L}`$ are Gaussian linear perturbations with zero mean. Although the form in equation (1) is inspired by simple inflation models, the exact predictions from those inflation models are irrelevant to our analysis here because the predicted amplitude of $`f_{\mathrm{NL}}`$ is much smaller than our sensitivity; however, this parameterization is useful to find quantitative constraints on the amount of non-Gaussianity allowed by the CMB data. Equation (1) is general in that $`f_{\mathrm{NL}}`$ parameterizes the leading-order non-linear corrections to $`\mathrm{\Phi }`$. We discuss the possible scale-dependence in Appendix C.

Angular bispectrum analyses found $`\left|f_{\mathrm{NL}}\right|<1500`$ (68%) from the COBE DMR 53$`+`$90 GHz coadded map (Komatsu et al., 2002) and $`\left|f_{\mathrm{NL}}\right|<950`$ (68%) from the MAXIMA sky map (Santos et al., 2002). The skewness measured from the DMR map smoothed with filters, called the Spherical Mexican Hat wavelets, found $`\left|f_{\mathrm{NL}}\right|<1100`$ (68%) (Cayón et al., 2002), although they neglected the integrated Sachs–Wolfe effect in the analysis, and therefore underestimated the cosmic variance of $`f_{\mathrm{NL}}`$. BOOMERanG did not measure $`f_{\mathrm{NL}}`$ in their analysis of non-Gaussianity (Polenta et al., 2002). The r.m.s. amplitude of $`\mathrm{\Phi }`$ is given by $`\mathrm{\Phi }^2^{1/2}\mathrm{\Phi }_\mathrm{L}^2^{1/2}\left(1+f_{\mathrm{NL}}^2\mathrm{\Phi }_\mathrm{L}^2\right)`$. Since $`\mathrm{\Phi }^2^{1/2}`$ measured on the COBE scales through the Sachs–Wolfe effect is $`\mathrm{\Phi }^2^{1/2}=3\mathrm{\Delta }T^2^{1/2}/T3.3\times 10^5`$ (Bennett et al., 1996), one obtains $`f_{\mathrm{NL}}^2\mathrm{\Phi }_\mathrm{L}^2<2.5\times 10^3`$ from the COBE 68% constraints; thus, we already know that the contribution from the non-linear term to the r.m.s. amplitude is smaller than 0.25%, and that to the power spectrum is smaller than 0.5%. This amplitude is comparable to limits on systematic errors of the WMAP power spectrum (Hinshaw et al., 2003a), and needs to be constrained better in order to verify the analysis of the power spectrum.

#### 2.2.2 Method 1: The Angular Bispectrum

Our first method for measuring $`f_{\mathrm{NL}}`$ is a “cubic statistic” which combines nearly optimally all configurations of the angular bispectrum of the primordial non-Gaussianity (Komatsu et al., 2003). The bispectrum measures phase correlations of field fluctuations. We compute the spherical harmonic coefficients $`a_{lm}`$ of temperature fluctuations from

$$a_{lm}=d^2\widehat{𝐧}M(\widehat{𝐧})\frac{\mathrm{\Delta }T(\widehat{𝐧})}{T_0}Y_{lm}^{}(\widehat{𝐧}),$$
(2)
where $`M(\widehat{𝐧})`$ is a pixel-weighting function. Here, $`M(\widehat{𝐧})`$ is the $`Kp0`$ sky cut where $`M(\widehat{𝐧})`$ takes 0 in the cut region and 1 otherwise. We filter the measured $`a_{lm}`$ in $`l`$-space and transform it back to compute two new maps, $`A(r,\widehat{𝐧})`$ and $`B(r,\widehat{𝐧})`$, given by

$`A(r,\widehat{𝐧})`$ $``$ $`{\displaystyle \underset{l=2}{\overset{l_{\mathrm{max}}}{}}}{\displaystyle \underset{m=l}{\overset{l}{}}}{\displaystyle \frac{\alpha _l(r)b_l}{\stackrel{~}{C}_l}}a_{lm}Y_{lm}(\widehat{𝐧}),`$ (3)
$`B(r,\widehat{𝐧})`$ $``$ $`{\displaystyle \underset{l=2}{\overset{l_{\mathrm{max}}}{}}}{\displaystyle \underset{m=l}{\overset{l}{}}}{\displaystyle \frac{\beta _l(r)b_l}{\stackrel{~}{C}_l}}a_{lm}Y_{lm}(\widehat{𝐧}).`$ (4)
Here $`\stackrel{~}{C}_lC_lb_l^2+N`$, where $`C_l`$ is the CMB anisotropy, $`N`$ the noise bias, and $`b_l`$ the beam window function describing the combined smoothing effects of the beam (Page et al., 2003a) and the finite pixel size. The functions $`\alpha _l(r)`$ and $`\beta _l(r)`$ are defined by

$`\alpha _l(r)`$ $``$ $`{\displaystyle \frac{2}{\pi }}{\displaystyle k^2𝑑kg_{\mathrm{T}l}(k)j_l(kr)},`$ (5)
$`\beta _l(r)`$ $``$ $`{\displaystyle \frac{2}{\pi }}{\displaystyle k^2𝑑kP(k)g_{\mathrm{T}l}(k)j_l(kr)},`$ (6)
where $`r`$ is the comoving distance. These two functions constitute the primordial angular bispectrum and correspond to $`\alpha _l(r)=f_{\mathrm{NL}}^1b_l^{\mathrm{NL}}(r)`$ and $`\beta _l(r)=b_l^\mathrm{L}(r)`$ in the notation of Komatsu & Spergel (2001). We compute the radiation transfer function $`g_{\mathrm{T}l}(k)`$ with a code based upon CMBFAST (Seljak & Zaldarriaga, 1996) for the best-fit cosmological model of the WMAP 1-year data (Spergel et al., 2003). We also use the best-fit primordial power spectrum of $`\mathrm{\Phi }`$, $`P(k)`$. We then compute the cubic statistic for the primordial non-Gaussianity, $`S_{\mathrm{prim}}`$, by integrating the two filtered maps over $`r`$ as (Komatsu et al., 2003)

$$𝒮_{\mathrm{prim}}=m_3^14\pi r^2𝑑r\frac{d^2\widehat{𝐧}}{4\pi }A(r,\widehat{𝐧})B^2(r,\widehat{𝐧}),$$
(7)
where the angular average is done on the full sky regardless of sky cut, and $`m_3=(4\pi )^1d^2\widehat{𝐧}M^3(\widehat{𝐧})`$ is the third-order moment of the pixel-weighting function. When the weight is only from a sky cut, as is the case here, we have $`m_3=f_{\mathrm{sky}}`$, i.e., $`m_3`$ is the fraction of the sky covered by observations (Komatsu et al., 2002). Komatsu et al. (2003) show that $`B`$ is a Wiener-filtered map of the underlying primordial fluctuations, $`\mathrm{\Phi }`$. The other map $`A`$ combines the bispectrum configurations that are sensitive to non-linearity of the form in equation (1). Thus, $`𝒮_{\mathrm{prim}}`$ is optimized for measuring the skewness of $`\mathrm{\Phi }`$ and picking out the quadratic term in equation (1).

Finally, the non-linear coupling parameter $`f_{\mathrm{NL}}`$ is given by

$$f_{\mathrm{NL}}\left[\underset{l_1l_2l_3}{\overset{l_{\mathrm{max}}}{}}\frac{(_{l_1l_2l_3}^{\mathrm{prim}})^2}{𝒞_{l_1}𝒞_{l_2}𝒞_{l_3}}\right]^1𝒮_{\mathrm{prim}},$$
(8)
where $`_{l_1l_2l_3}^{\mathrm{prim}}`$ is the primordial bispectrum (Komatsu & Spergel, 2001) multiplied by $`b_{l_1}b_{l_2}b_{l_3}`$ and computed for $`f_{\mathrm{NL}}=1`$ and the best-fit cosmological model. This equation is used to measure $`f_{\mathrm{NL}}`$ as a function of the maximum multipole $`l_{\mathrm{max}}`$. The statistic $`𝒮_{\mathrm{prim}}`$ takes only $`N^{3/2}`$ operations to compute without loss of sensitivity whereas the full bispectrum analysis takes $`N^{5/2}`$ operations. It takes about 4 minutes on 16 processors of an SGI Origin 300 to compute $`f_{\mathrm{NL}}`$ from a sky map at the highest resolution level, $`nside=512`$. We measure $`f_{\mathrm{NL}}`$ as a function of $`l_{\mathrm{max}}`$. Since there is little CMB signal compared with instrumental noise at $`l>512`$, we shall use $`l_{\mathrm{max}}=512`$ at most; thus, $`nside=256`$ is sufficient, speeding up evaluations of $`f_{\mathrm{NL}}`$ by a factor of 8 as the computational time scales as $`(nside)^3`$. The computation takes only 30 seconds at $`nside=256`$. Note that since we are eventually fitting for two parameters, $`f_{\mathrm{NL}}`$ and $`b_{\mathrm{src}}`$ (see sec. 3), we include covariance between these two parameters in the analysis. The covariance is, however, small (see Figure 8 in Appendix A).

While we use uniform weighting for $`M(\widehat{𝐧})`$, we could instead weight by the inverse noise variance per pixel, $`M(\widehat{𝐧})=N^1(\widehat{𝐧})`$; however, this weighting scheme is sub-optimal at low $`l`$ where the CMB anisotropy dominates over noise so that the uniform weighing is more appropriate. For measuring $`b_{\mathrm{src}}`$, on the other hand, we shall use a slightly modified version of the $`N^1`$ weighting, as $`b_{\mathrm{src}}`$ comes mainly from small angular scales where instrumental noise dominates.

#### 2.2.3 Method 2: The Minkowski Functionals

Topology offers another test for non-Gaussian features in the maps, measuring morphological structures of fluctuation fields. The Minkowski functionals (Minkowski, 1903; Gott et al., 1990; Schmalzing & Gorski, 1998) describe the properties of regions spatially bounded by a set of contours. The contours may be specified in terms of fixed temperature thresholds, $`\nu =\mathrm{\Delta }T/\sigma `$, where $`\sigma `$ is the standard deviation of the map, or in terms of the area. Parameterization of contours by threshold is computationally simpler, while parameterization by area reduces correlations between the Minkowski functionals (Shandarin et al., 2002). We use a joint analysis of the three Minkowski functionals (area $`A(\nu )`$, contour length $`C(\nu )`$, and genus $`G(\nu )`$) explicitly including their covariance; consequently, we work in the simpler threshold parameterization.

The Minkowski functionals are additive for disjoint regions on the sky and are invariant under coordinate transformation and rotation. We approximate each Minkowski functional using the set of equal-area pixels hotter or colder than a set of fixed temperature thresholds. The fractional area

$$A(\nu )=\frac{1}{A}\underset{i}{}a_i=\frac{N_\nu }{N_{\mathrm{cut}}}$$
(9)
is thus the number of enclosed pixels, $`N_\nu `$, divided by the total number of pixels on the cut sky, $`N_{\mathrm{cut}}`$. Here $`a_i`$ is the area of an individual spot, and $`A`$ is the total area of the pixels outside the cut. The contour length

$$C(\nu )=\frac{1}{4A}\underset{i}{}P_i$$
(10)
is the total perimeter length of the enclosed regions $`P_i`$, while the genus

$$G(\nu )=\frac{1}{2\pi A}(N_{\mathrm{hot}}N_{\mathrm{cold}})$$
(11)
is the number of hot spots, $`N_{\mathrm{hot}}`$, minus the number of cold spots, $`N_{\mathrm{cold}}`$. We calibrate finite pixelization effects by comparing the Minkowski functionals for the WMAP data to Monte Carlo simulations.

The WMAP data are a superposition of sky signal and instrument noise, each with a different morphology. The Minkowski functionals transform monotonically (although not linearly) between the limiting cases of a sky signal with no noise and a noise map with no sky signal. Unlike spatial analyses such as Fourier decomposition, different regions of the sky cannot be weighted by the signal-to-noise ratio, nor does the noise “average down” over many pixels. The choice of map pixelization becomes a tradeoff between resolution (favoring smaller pixels) versus signal-to-noise ratio (favoring larger pixels). We compute the Minkowski functionals at $`nside=16`$ through 256 (3072 to 786,432 pixels on the full sky). We use the WMAP Kp0 sky cut to reject pixels near the Galactic plane or contaminated by known sources. The cut sky has 1433 pixels at resolution $`nside=16`$ and 666,261 pixels at $`nside=256`$.

We compute the Minkowski functionals at 15 thresholds from $`3.5\sigma `$ to $`+3.5\sigma `$, and compare each functional to the simulations using a goodness-of-fit statistic,

$$\chi ^2=\underset{\nu _1\nu _2}{}[F_{\mathrm{WMAP}}^iF_{\mathrm{sim}}^i]_{\nu _1}\mathrm{\Sigma }_{\nu _1\nu _2}^1[F_{\mathrm{WMAP}}^iF_{\mathrm{sim}}^i]_{\nu _2},$$
(12)
where $`F_{\mathrm{WMAP}}^i`$ is a Minkowski functional from the WMAP data (the index $`i`$ denotes a kind of functional), $`F_{\mathrm{sim}}^i`$ is the mean from the Monte Carlo simulations, and $`\mathrm{\Sigma }_{\nu _1\nu _2}`$ is the bin-to-bin covariance matrix from the simulations.

### 2.3 Monte Carlo Simulations

Monte Carlo simulations are used to estimate the statistical significance of the non-Gaussian signals. One kind of simulation generates Gaussian random realizations of CMB sky maps for the angular power spectrum, window functions, and noise properties of the WMAP 1-year data. This simulation quantifies the uncertainty arising from Gaussian fields, or the uncertainty in the absence of non-Gaussian fluctuations. The other kind generates non-Gaussian CMB sky maps from primordial fluctuations of the form of equation (1) (see Appendix A for our method for simulating non-Gaussian maps). This simulation quantifies the uncertainty more accurately and consistently in the presence of non-Gaussian fluctuations.

In principle, one should always use the non-Gaussian simulations to characterize the uncertainty in $`f_{\mathrm{NL}}`$; however, the uncertainty estimated from the Gaussian realizations is good approximation to that from the non-Gaussian ones as long as $`\left|f_{\mathrm{NL}}\right|<500`$. Our non-Gaussian simulations verify that the distribution of $`f_{\mathrm{NL}}`$ and $`b_{\mathrm{src}}`$ around the mean is the same for Gaussian and non-Gaussian realizations (see Figure 8 in Appendix A for an example of $`f_{\mathrm{NL}}=100`$). The Gaussian simulations have the advantage of being much faster than the non-Gaussian ones. The former takes only a few seconds to simulate one map whereas the latter takes 3 hours on a single processor of an SGI Origin 300. Also, simulating non-Gaussian maps at $`nside=512`$ requires 17 GB of physical memory. We therefore use Gaussian simulations to estimate the uncertainty in measured $`f_{\mathrm{NL}}`$ and $`b_{\mathrm{src}}`$.

### 2.4 Limits to Primordial Non-Gaussianity

Figure 1 shows $`f_{\mathrm{NL}}`$ measured from the Q$`+`$V$`+`$W coadded map using the cubic statistic \[Eq (8)\], as a function of the maximum multipole $`l_{\mathrm{max}}`$. We find the best estimate of $`f_{\mathrm{NL}}=38\pm 48`$ (68%) for $`l_{\mathrm{max}}=265`$. The distribution of $`f_{\mathrm{NL}}`$ is close to a Gaussian, as suggested by Monte Carlo simulations (see Figure 8 in Appendix A). The 95% confidence interval is $`58<f_{\mathrm{NL}}<134`$. There is no significant detection of $`f_{\mathrm{NL}}`$ at any angular scale. The r.m.s. error, estimated from 500 Gaussian simulations, initially decreases as $`l_{\mathrm{max}}^1`$, although $`f_{\mathrm{NL}}`$ for $`l_{\mathrm{max}}=265`$ has a smaller error than that for $`l_{\mathrm{max}}=512`$ because the latter is dominated by the instrumental noise. Since all the pixels outside the cut region are uniformly weighted, the inhomogeneous noise in the map (pixels on the ecliptic equator are noisier than those on the north and south poles) is not accounted for. This leads to a noisier estimator than a minimum variance estimator. The constraint on $`f_{\mathrm{NL}}`$ for $`l_{\mathrm{max}}=512`$ will improve with more appropriate pixel-weighting schemes (Heavens, 1998; Santos et al., 2002). The simple inverse noise ($`N^1`$) weighting makes the constraints much worse than the uniform weighting, as it increases errors on large angular scales where the CMB signal dominates over the instrumental noise. The uniform weighting is thus closer to optimal. Note that for the power spectrum, one can simply use the uniform weighting to measure $`C_l`$ at small $`l`$ and the $`N^1`$ weighting at large $`l`$. For the bispectrum, however, this decomposition is not simple, as the bispectrum $`_{l_1l_2l_3}`$ measures the mode coupling from $`l_1`$ to $`l_2`$ and $`l_3`$ and vice versa. This property makes it difficult to use different weighting schemes on different angular scales. The first column of Table 1 shows $`f_{\mathrm{NL}}`$ measured in Q, V, and W bands separately. There is no a significant band-to-band variation, or a significant detection in any band.

Figure 2 shows the Minkowski functionals at $`nside=128`$ (147,594 high-latitude pixels, each $`28^{}`$ in diameter). The gray band shows the 68% confidence region derived from 1000 Gaussian simulations. Table 2 shows the $`\chi ^2`$ values \[Eq.(12)\]. The data are in excellent agreement with the Gaussian simulations at all resolutions. The individual Minkowski functionals are highly correlated with each other (e.g., Shandarin et al. (2002)). We account for this using a simultaneous analysis of all three Minkowski functionals, replacing the 15-element vectors $`F_{\mathrm{WMAP},\nu }^i`$ and $`F_{\mathrm{sim},\nu }^i`$ in equation (12) (the index $`i`$ denotes each Minkowski functional) with 45-element vectors $`F_\nu =[F^1,F^2,F^3]_\nu `$ = $`[`$Area, Contour, Genus$`]_\nu `$ and using the covariance of this larger vector as derived from the simulations. We compute $`\chi ^2`$ for values $`f_{\mathrm{NL}}`$ = 0 to 1000, comparing the results from WMAP to similar $`\chi ^2`$ values computed from non-Gaussian realizations. Figure 3 shows the result. We find a best-fit value $`f_{\mathrm{NL}}=22\pm 81`$ (68%), with 95% confidence upper limit $`f_{\mathrm{NL}}<139`$, in agreement with the cubic statistic.

### 2.5 Implications of the WMAP Limits on $`f_{\mathrm{NL}}`$

#### 2.5.1 Inflation

The limits on $`f_{\mathrm{NL}}`$ are consistent with simple inflation models: models based on a slowly rolling scalar field typically give $`\left|f_{\mathrm{NL}}\right|10^210^1`$ (Salopek & Bond, 1990, 1991; Falk et al., 1993; Gangui et al., 1994; Acquaviva et al., 2002; Maldacena, 2002), three to four orders of magnitude below our limits. Measuring $`f_{\mathrm{NL}}`$ at this level is difficult because of the cosmic variance. There are alternative models which allow larger amplitudes of non-Gaussiantiy in the primordial fluctuations, which we explore below.

A large $`f_{\mathrm{NL}}`$ may be produced when the following condition is met. Suppose that $`\mathrm{\Phi }`$ is given by $`\mathrm{\Phi }=ϵx`$ where $`ϵ`$ is a transfer function that converts $`x`$ to $`\mathrm{\Phi }`$ and $`x=x^{(1)}+x^{(2)}+𝒪(x^{(3)})`$ denotes a fluctuating field expanded into a series of $`x^{(i)}=f_ix^{(i1)}x^{(1)}`$ with $`f_1=1`$. Then, $`f_{\mathrm{NL}}=ϵ^1f_2`$. Inflation predicts the amplitude of $`x^{(i)}`$ and the form of $`f_i`$ which eventually depends upon the scalar field potential; thus, $`x^{(i)}`$ would be of order $`(H/m_{\mathrm{planck}})^i`$ ($`H`$ is the Hubble parameter during inflation) for $`H<m_{\mathrm{planck}}`$, and the leading order term is $`ϵH/m_{\mathrm{planck}}10^5ϵ`$. In this way $`ϵ`$ “suppresses” the amplitude of fluctuations, allowing a larger amplitude for $`H/m_{\mathrm{planck}}10^5ϵ^1`$. What does this mean? If $`H10^2m_{\mathrm{planck}}`$, then $`ϵ10^3`$ and $`f_{\mathrm{NL}}10^3f_2`$. The amplitude of $`f_{\mathrm{NL}}`$ is thus large enough to detect for $`f_20.1`$. This suppression factor, $`ϵ`$, seems necessary for one to obtain a large $`f_{\mathrm{NL}}`$ in the context of the slow-roll inflation. The suppression also helps us to avoid a “fine-tuning problem” of inflation models, as it allows $`H/m_{\mathrm{planck}}`$ to be of order slightly less than unity (which one might think natural) rather than forcing it to be of order $`10^5`$.

Curvatons proposed by Lyth & Wands (2002) provide an example of a supression mechanism. A curvaton is a scalar field, $`\sigma `$, having mass, $`m_\sigma `$, that develops fluctuations, $`\delta \sigma `$, during inflation with its energy density, $`\rho _\sigma V(\sigma )`$, tiny compared to that of the inflaton field that drives inflation. After inflation ends, radiation is produced as the inflaton decays, generating entropy perturbations between $`\sigma `$ and radiation, $`S_{\sigma \gamma }=\delta \rho _\sigma /\rho _\sigma \frac{3}{4}\delta \rho _\gamma /\rho _\gamma `$. When $`H`$ decreases to become comparable to $`m_\sigma `$, oscillations of $`\sigma `$ at the bottom of $`V(\sigma )`$ give $`\rho _\sigma m_\sigma ^2\sigma ^2`$. In the limit of “cold inflation” for which $`\delta \rho _\gamma /\rho _\gamma `$ is nearly zero, one finds $`S_{\sigma \gamma }\delta \rho _\sigma /\rho _\sigma 2\delta \sigma /\sigma +(\delta \sigma /\sigma )^2`$. As long as $`\sigma `$ survives after the production of $`S_{\sigma \gamma }`$, the curvature perturbation $`\mathrm{\Phi }`$ is generated as $`\mathrm{\Phi }=\frac{1}{2}ϵS_{\sigma \gamma }ϵ[x^{(1)}+\frac{1}{2}(x^{(1)})^2]`$ where $`x^{(1)}=\delta \sigma /\sigma `$ (i.e., $`f_2=\frac{1}{2}`$). The generation of $`\mathrm{\Phi }`$ continues until $`\sigma `$ decays, and $`\mathrm{\Phi }`$ is essentially determined by a ratio of $`\rho _\sigma `$ to the total energy density, $`\mathrm{\Omega }_\sigma `$, at the time of the decay. Lyth et al. (2002) numerically evolved perturbations to find $`ϵ\frac{2}{5}\mathrm{\Omega }_\sigma `$ at the time of the decay. The smaller the curvaton energy density is, the less efficient the $`S_{\sigma \gamma }`$ to $`\mathrm{\Phi }`$ conversion becomes (or the more efficient the supression becomes). The small $`\mathrm{\Omega }_\sigma `$ thus leads to the large $`f_{\mathrm{NL}}`$, as $`f_{\mathrm{NL}}=ϵ^1f_2\frac{5}{4}\mathrm{\Omega }_\sigma ^1`$ (i.e., $`f_{\mathrm{NL}}`$ is always positive in this model). Assuming the curvaton exists and is entirely responsible for the observed CMB anisotropy, our limits on $`f_{\mathrm{NL}}`$ imply $`\mathrm{\Omega }_\sigma >9\times 10^3`$ at the time of the curvaton decay. (However, the lower limit to $`\mathrm{\Omega }_\sigma `$ does not mean that we need the curvatons. This constraint makes sense only when the curvaton exists and is entierly producing the observed fluctuations.)

Features in an inflaton potential can generate significant non-Gaussian fluctuations (Kofman et al., 1991; Wang & Kamionkowski, 2000), and it is expected that measurements of non-Gaussianity can place constrains on a class of the feature models. In Appendix C, we calculate the angular bispectrum from a sudden step in a potential of the form in equation (C2). This step is motivated by a class of supergravity models yielding the steps as a consequence of successive spontaneous symmetry-breaking phase transitions of many fields coupled to the inflaton (Adams et al., 1997, 2001). One step generates two distinct regions in $`l`$ space where $`\left|f_{\mathrm{NL}}\right|`$ is very large: a positive $`f_{\mathrm{NL}}`$ is predicted at $`l<l_\mathrm{f}`$, while a negative $`f_{\mathrm{NL}}`$ at $`l>l_\mathrm{f}`$, where $`l_\mathrm{f}`$ is the projected location of the step. Our calculations suggest that the two regions are separated in $`l`$ by less than a factor of 2, and one cannot resolve them without knowing $`l_\mathrm{f}`$. The average of many $`l`$ modes further smeares out the signals. The averaged $`f_{\mathrm{NL}}`$ thus nearly cancels out to give only small signals, being hidden in our constraints in Figure 1. Peiris et al. (2003) argue that some sharp features in the WMAP angular power spectrum producing large $`\chi ^2`$ values may arise from features in the inflaton potential. If this is true, then one may be able to see non-Gaussian signals associated with the features by measuring the bispectrum at the scales of the sharp features of the power spectrum.

#### 2.5.2 Massive cluster abundance at high redshift

Massive halos, like clusters of galaxies at high redshift, are such rare objects in the universe that their abundance is sensitive to the presence of non-Gaussianity in the distribution function of primordial density fluctuations. Several authors have pointed out the power of the halo abundance as a tool for finding primordial non-Gaussianity (Lucchin & Matarrese, 1988; Robinson & Baker, 2000; Matarrese et al., 2000; Benson et al., 2002); however, the power of this method is extremely sensitive to the accuracy of the mass determinations of halos. It is necessary to go to redshifts of $`z1`$ to obtain tight constraints on primordial non-Gaussianity, as constraints from low and intermediate redshifts appear to be weak (Koyama et al., 1999; Robinson et al., 2000) (see also Figure 4 and 5). Due to the difficulty of measuring the mass of a high-redshift cluster the current constraints are not yet conclusive (Willick, 2000). The limited number of clusters observed at high redshift also limits the current sensitivity. In this section, we translate our constraints on $`f_{\mathrm{NL}}`$ from the WMAP 1-year CMB data into the effects on the massive halos in the high-redshift universe, showing the extent to which future cluster surveys would see signatures of non-Gaussian fluctuations.

We adopt the method of Matarrese et al. (2000) to calculate the dark-matter halo mass function $`dn/dM`$ for a given $`f_{\mathrm{NL}}`$, using the $`\mathrm{\Lambda }`$CDM with the running spectral index model best-fit to the WMAP data and the large-scale structure data. This set of parameters is best suited for the calculations of the cluster abundance. The parameters are in the rightmost column of Table 8 of Spergel et al. (2003). We calculate

$$\frac{dn}{dM}=2\frac{\rho _{\mathrm{m0}}}{M}\left|\frac{dP}{dM}\right|,$$
(13)
where $`\rho _{\mathrm{m0}}=2.775\times 10^{11}(\mathrm{\Omega }_\mathrm{m}h^2)M_{}\mathrm{Mpc}^3=3.7\times 10^{10}M_{}\mathrm{Mpc}^3`$ is the present-day mean mass density of the universe, $`P(M,z)`$ is the probability for halos of mass $`M`$ to collapse at redshift $`z`$, and $`dP/dM`$ is given by

$$\frac{dP}{dM}_0^{\mathrm{}}\frac{\lambda d\lambda }{2\pi }\left(\frac{d\sigma ^2}{dM}\mathrm{sin}\theta _\lambda \frac{\lambda }{3}\frac{d\mu _3}{dM}\mathrm{cos}\theta _\lambda \right)e^{\lambda ^2\sigma ^2/2},$$
(14)
where the angle $`\theta _\lambda `$ is given by $`\theta _\lambda \lambda \delta _\mathrm{c}+\lambda ^3\mu _3/6`$, and $`\delta _c(z)`$ is the threshold overdensity of spherical collapse (Lacey & Cole, 1993; Nakamura & Suto, 1997). The variance of the mass fluctuations as a function of $`z`$ is given by $`\sigma ^2(M,z)=D^2(z)\sigma ^2(M,0)`$, where $`D(z)`$ is the growth factor of linear density fluctuations, $`\sigma ^2(M,0)=_0^{\mathrm{}}𝑑kk^1F_M^2(k)\mathrm{\Delta }^2(k)`$, $`\mathrm{\Delta }^2(k)(2\pi ^2)^1k^3P(k)`$ is the dimensionless power spectrum of the Bardeen curvature perturbations, $`F_M(k)g(k)T(k)W(kR_M)`$ a filter function, $`g(k)\frac{2}{3}(k/H_0)^2\mathrm{\Omega }_{\mathrm{m0}}^1`$ a conversion factor from $`\mathrm{\Phi }`$ to density fluctuations, $`T(k)`$ the transfer function of linear density perturbations, $`W(x)3j_1(x)/x`$ the spherical top-hat window smoothing density fields, and $`R_M[3M/(4\pi \rho _{\mathrm{m0}})]^{1/3}`$ the spherical top-hat radius enclosing a mass $`M`$. The skewness $`\mu _3(M,z)=D^3(z)\mu _3(M,0)`$, where

$$\mu _3(M,0)=6f_{\mathrm{NL}}_0^{\mathrm{}}\frac{dk_1}{k_1}F_M(k_1)\mathrm{\Delta }^2(k_1)_0^{\mathrm{}}\frac{dk_2}{k_2}F_M(k_2)\mathrm{\Delta }^2(k_2)_0^1𝑑\mu F_M(\sqrt{k_1^2+k_2^2+2k_1k_2\mu }),$$
(15)
arises from the primordial non-Gaussianity. We use a Monte Carlo integration routine called vegas (Press et al., 1992) to evaluate the triple integral in equation (15). It follows from equation (15) that a positive $`f_{\mathrm{NL}}`$ gives a positive $`\mu _3`$, positively skewed density fluctuations. Also this $`dn/dM`$ reduces to the Press–Schechter form (Press & Schechter, 1974) in the limit of $`f_{\mathrm{NL}}0`$. Although the Press–Schechter form predicts significantly fewer massive halos than $`N`$-body simulations (Jenkins et al., 2001), we assume that a predicted ratio of the non-Gaussian $`dn/dM`$ to the Gaussian $`dn/dM`$ is still reasonably accurate, as the primordial non-Gaussianity does not affect the dynamics of halo formations which causes the difference between the Press–Schechter form of $`dn/dM`$ and the $`N`$-body simulations.

Figure 4 shows the WMAP constraints on the ratio of non-Gaussian $`dn/dM`$ to the Gaussian one, as a function of $`M`$ and $`z`$. We find that the WMAP constraint on $`f_{\mathrm{NL}}`$ strongly limits the amplitude of changes in $`dn/dM`$ due to the non-Gaussianity. At $`z=0`$, $`dn/dM`$ is changed by no more than 20% even for $`4\times 10^{15}M_{}`$ clusters. The number of clusters that would be newly found at $`z=1`$ for $`M<10^{15}M_{}`$ should be within $`{}_{10}{}^{+40}\%`$ of the value predicted from the Gaussian theory. At $`z=3`$, however, much larger effects are still allowed: $`dn/dM`$ can be increased by up to a factor of 2.5 for $`2\times 10^{14}M_{}`$.

Predictions for actual cluster surveys are made clearer by computing the source number counts as a function of $`z`$,

$$\frac{dN}{dz}\frac{dV}{dz}_{M_{\mathrm{lim}}}^{\mathrm{}}𝑑M\frac{dn}{dM},$$
(16)
where $`V(z)`$ is the comoving volume per steradian, and $`M_{\mathrm{lim}}`$ is the limiting mass that a survey can reach. In practice $`M_{\mathrm{lim}}`$ would depend on $`z`$ due to, for example, the redshift dimming of X-ray surface brightness; however, a constant $`M_{\mathrm{lim}}`$ turns out to be a good approximation for surveys of the Sunyaev–Zel’dovich (SZ) effect (Carlstrom et al., 2002). Figure 5 shows the ratio for $`dN/dz`$ as a function of $`z`$ and $`M_{\mathrm{lim}}`$. A source-detection sensitivity of $`S_{\mathrm{lim}}=0.5\mathrm{Jy}`$ roughly corresponds to $`M_{\mathrm{lim}}=1.4\times 10^{14}M_{}`$ (Carlstrom et al., 2002), for which $`dN/dz`$ should follow the prediction of the Gaussian theory out to $`z1`$ to within 10%, but $`dN/dz`$ at $`z=3`$ can be increased by up to a factor of 2. As $`M_{\mathrm{lim}}`$ increases, the impact on $`dN/dz`$ rapidly increases.

The SZ angular power spectrum $`C_l^{\mathrm{SZ}}`$ is so sensitive to $`\sigma _8`$ that we can use $`C_l^{\mathrm{SZ}}`$ to measure $`\sigma _8`$ (Komatsu & Kitayama, 1999). The sensitivity arises largely from massive ($`M>10^{14}M_{}`$) clusters at $`z1`$. From this fact one might argue that $`C_l^{\mathrm{SZ}}`$ is also sensitive to the primordial non-Gaussianity. We use a method of Komatsu & Seljak (2002) with $`dn/dM`$ replaced by equation (13) to compute $`C_l^{\mathrm{SZ}}`$ for the WMAP limits on $`f_{\mathrm{NL}}`$. We find that $`C_l^{\mathrm{SZ}}`$ should follow the prediction from the Gaussian theory to within 10% for $`100<l<10000`$. This is consistent with $`C_l^{\mathrm{SZ}}`$ being primarily sensitive to halos at $`z1`$, where the effect on $`dN/dz`$ is not too strong (see Figure 5). Since $`C_l^{\mathrm{SZ}}\sigma _8^7(\mathrm{\Omega }_\mathrm{b}h)^2`$ (Komatsu & Seljak, 2002), $`\sigma _8`$ can be determined from $`C_l^{\mathrm{SZ}}`$ to within 2% accuracy at a fixed $`\mathrm{\Omega }_\mathrm{b}h`$ using the Gaussian theory. The current theoretical uncertainty in the predictions of $`C_l^{\mathrm{SZ}}`$ is a factor of 2 in $`C_l^{\mathrm{SZ}}`$ (10% in $`\sigma _8`$), still much larger than the effect of the non-Gaussianity.

## 3 LIMITS TO RESIDUAL POINT SOURCES

### 3.1 Point-source Angular Power Spectrum and Bispectrum

Radio point sources distributed across the sky generate non-Gaussian signals, giving a positive bispectrum, $`b_{\mathrm{src}}`$ (Komatsu & Spergel, 2001). In addition, the point sources contribute significantly to the angular spectrum on small angular scales (Tegmark & Efstathiou, 1996), contaminating the cosmological angular power spectrum. It is thus important to understand how much of the measured angular power spectrum is due to sources. We constrain the source contribution to the angular power spectrum, $`c_{\mathrm{src}}`$, by measuring $`b_{\mathrm{src}}`$. Komatsu & Spergel (2001) have shown that WMAP can detect $`b_{\mathrm{src}}`$ even after subtracting all (bright) sources detected in the sky maps. Fortunately, there is no degeneracy between $`f_{\mathrm{NL}}`$ and $`b_{\mathrm{src}}`$, as shown later in Appendix A.

In this section we measure the amplitude of non-Gaussianity from “residual” point sources which are fainter than a certain flux threshold, $`S_\mathrm{c}`$, and left unmasked in the sky maps. The bispectrum $`b_{\mathrm{src}}`$ is related to the number of sources brighter than $`S_\mathrm{c}`$ per solid angle $`N(>S_\mathrm{c})`$:

$$b_{\mathrm{src}}(S_\mathrm{c})=_0^{S_\mathrm{c}}dS\frac{dN}{dS}[g(\nu )S]^3=N(>S_\mathrm{c})[g(\nu )S_\mathrm{c}]^3+3_0^{S_\mathrm{c}}\frac{dS}{S}N(>S)[g(\nu )S]^3,$$
(17)
where $`g(\nu )`$ is a conversion factor from $`\mathrm{Jy}\mathrm{sr}^1`$ to $`\mu `$K which depends on observing frequency $`\nu `$ as $`g(\nu )=(24.76\mathrm{Jy}\mu \mathrm{K}^1\mathrm{sr}^1)^1[(\mathrm{sinh}x/2)/x^2]^2`$, $`xh\nu /k_\mathrm{B}T_0\nu /(56.78\mathrm{GHz})`$ for $`T_0=2.725\mathrm{K}`$ (Mather et al., 1999), and $`dN/dS`$ is the differential source count per solid angle. The residual point sources also contribute to the point-source power spectrum $`c_{\mathrm{src}}`$ as

$$c_{\mathrm{src}}(S_\mathrm{c})=_0^{S_\mathrm{c}}dS\frac{dN}{dS}[g(\nu )S]^2=N(>S_\mathrm{c})[g(\nu )S_\mathrm{c}]^2+2_0^{S_\mathrm{c}}\frac{dS}{S}N(>S)[g(\nu )S]^2.$$
(18)
By combining equation (17) and (18) we find a relation between $`b_{\mathrm{src}}`$ and $`c_{\mathrm{src}}`$,

$$c_{\mathrm{src}}(S_\mathrm{c})=b_{\mathrm{src}}(S_\mathrm{c})[g(\nu )S_\mathrm{c}]^1+_0^{S_\mathrm{c}}\frac{dS}{S}b_{\mathrm{src}}(S)[g(\nu )S]^1.$$
(19)
We can use this equation combined with the measured $`b_{\mathrm{src}}`$ as a function of $`S_\mathrm{c}`$ to directly determine $`c_{\mathrm{src}}`$ as a function of $`S_\mathrm{c}`$, without relying on any extrapolations. When the source counts obey a power-law like $`dN/dSS^\beta `$, one finds $`b_{\mathrm{src}}(S)S^{4+\beta }`$; thus, brighter sources contribute more to the integral in equation (19) than fainter ones as long as $`\beta >3`$, which is the case for fluxes of interest. Bennett et al. (2003a) have found $`\beta =2.6\pm 0.2`$ for $`S=210`$ Jy in Q band. Below 1 Jy, $`\beta `$ becomes even flatter (Toffolatti et al., 1998), implying that one does not have to go down to the very faint end to obtain reasonable estimates of the integral. In practice, we use equation (17) with $`N(>S)`$ of the Toffolatti et al. (1998) model (hereafter T98) at 44 GHz to compute $`b_{\mathrm{src}}(S<0.5\mathrm{Jy})`$, inserting it into the integral to avoid missing faint sources and underestimating the integral.

### 3.2 Measurement of the Point-source Angular Bispectrum

The reduced point-source angular bispectrum, $`b_{\mathrm{src}}`$, is measured by a cubic statistic for point sources (Komatsu et al., 2003),

$$𝒮_{ps}=m_3^1\frac{d^2\widehat{𝐧}}{4\pi }D^3(\widehat{𝐧}),$$
(20)
where the filtered map $`D(\widehat{𝐧})`$ is given by

$$D(\widehat{𝐧})\underset{l=2}{\overset{l_{\mathrm{max}}}{}}\underset{m=l}{\overset{l}{}}\frac{b_l}{\stackrel{~}{C}_l}a_{lm}Y_{lm}(\widehat{𝐧}).$$
(21)
This statistic is even quicker ($`100`$ times) to compute than $`S_{\mathrm{prim}}`$ (eq.), as it involves only one integral over $`\widehat{𝐧}`$ and only one filtered map. This statistic also retains the same sensitivity to the point-source non-Gaussianity as the full bispectrum analysis. The cubic statistic $`S_{\mathrm{ps}}`$ gives $`b_{\mathrm{src}}`$ as

$$b_{\mathrm{src}}\left[\frac{3}{2\pi }\underset{l_1l_2l_3}{\overset{l_{\mathrm{max}}}{}}\frac{(_{l_1l_2l_3}^{\mathrm{ps}})^2}{𝒞_{l_1}𝒞_{l_2}𝒞_{l_3}}\right]^1𝒮_{\mathrm{ps}},$$
(22)
where $`_{l_1l_2l_3}^{\mathrm{ps}}`$ is the point-source bispectrum for $`b_{\mathrm{src}}=1`$ (Komatsu & Spergel, 2001) multiplied by $`b_{l_1}b_{l_2}b_{l_3}`$. While the uniform pixel-weighting outside the Galactic cut was used for $`f_{\mathrm{NL}}`$, we use here $`M(\widehat{𝐧})=\left[\sigma _{\mathrm{CMB}}^2+N(\widehat{𝐧})\right]^1`$ where $`\sigma _{\mathrm{CMB}}^2=(4\pi )^1_l(2l+1)C_lb_l^2`$ is the variance of CMB anisotropy and $`N(\widehat{𝐧})`$ is the variance of noise per pixel which varies across the sky. This weighting scheme is nearly optimal for measuring $`b_{\mathrm{src}}`$ as the signal comes from smaller angular scales where noise dominates. The factor of $`\sigma _{\mathrm{CMB}}^2`$ approximately takes into account the non-zero contribution to the variance from CMB anisotropy. This weight reduces uncertainties of $`b_{\mathrm{src}}`$ by 17%, 23%, and 31% in Q, V, and W bands, respectively, compared to the uniform weighting. We use the highest resolution level, $`nside=512`$, and integrate equation (22) up to $`l_{\mathrm{max}}=1024`$. In Appendix B, it is shown that this estimator is optimal and unbiased as long as very bright sources, which have contributions to $`\stackrel{~}{C}_l`$ too large to ignore, are masked. We cannot include $`c_{\mathrm{src}}`$ in the filter, as it is what we are trying to measure using $`b_{\mathrm{src}}`$.

The filled circles in the left panels of Figure 6 represent $`b_{\mathrm{src}}`$ measured in Q (top panel) and V (bottom panel) band. We have used source masks for various flux cuts, $`S_\mathrm{c}`$, defined at 4.85 GHz to make these measurements. (The masks are made from the GB6$`+`$PMN 5 GHz source catalogue.) We find that $`b_{\mathrm{src}}`$ increases as $`S_\mathrm{c}`$: the brighter sources unmasked, the more non-Gaussianity is detected. On the other hand one can make predictions for $`b_{\mathrm{src}}`$ using equation (17) for a given $`N(>S)`$. Comparing the measured values of $`b_{\mathrm{src}}`$ with the predicted values from $`N(>S)`$ of T98 (dashed lines) at 44 GHz, one finds that the measured values are smaller than the predicted values by a factor of 0.65. The solid lines show the predictions multiplied by 0.65. Both errors in the T98 predictions and a non-flat energy spectrum of sources easily cause this factor. (If sources have a non-flat spectrum like $`S\nu ^\alpha `$ where $`\alpha 0`$, then $`S_\mathrm{c}`$ at Q or V band is different from that at 4.85 GHz.) Bennett et al. (2003a) find that the majority of the radio sources detected in Q band have a flat spectrum, $`\alpha =0.0\pm 0.2`$. Our value for the correction factor matches well the one obtained from the WMAP source counts for $`210\mathrm{Jy}`$ in Q band (Bennett et al., 2003a).

Equation (18) combined with the measured $`b_{\mathrm{src}}`$ is used to estimate the point-source angular power spectrum $`c_{\mathrm{src}}`$. The right panels of Figure 6 show the estimated $`c_{\mathrm{src}}`$ as filled circles. These estimates agree well with predictions from equation (18) with $`N(>S)`$ of T98 multiplied by a factor of 0.65 (solid lines). For $`S_\mathrm{c}=1\mathrm{Jy}`$ at Q band, $`\widehat{c}_{\mathrm{src}}=(19\pm 5)\times 10^3\mu \mathrm{K}^2\mathrm{sr}`$, and matches well the value estimated from the WMAP source counts at the same flux threshold (Bennett et al., 2003a), which corresponds to the solid lines in the figure. At V band, $`\widehat{c}_{\mathrm{src}}=(5\pm 4)\times 10^3\mu \mathrm{K}^2\mathrm{sr}`$. Here, the hat denotes that these values do not represent $`c_{\mathrm{src}}`$ for the standard source mask used by Hinshaw et al. (2003b) for estimating the cosmological angular power spectrum. Since the standard source mask is made of several source catalogues with different selection thresholds, it is difficult to clearly identify a mask flux cut. We give the standard mask an “effective” flux cut threshold at 4.85 GHz by comparing $`b_{\mathrm{src}}`$ measured from the standard source mask (shaded areas in Figure 6; see the second column of Table 1 for actual values) with those from the GB6$`+`$PMN masks defined at 4.85 GHz. The measurements agree when $`S_\mathrm{c}0.75\mathrm{Jy}`$ in Q band. Using this effective threshold, one expects $`c_{\mathrm{src}}`$ for the standard source mask as $`c_{\mathrm{src}}=(15\pm 6)\times 10^3\mu \mathrm{K}^2\mathrm{sr}`$ in Q band. This value agrees with the excess power seen on small angular scales, $`(15.5\pm 1.7)\times 10^3\mu \mathrm{K}^2\mathrm{sr}`$ (Hinshaw et al., 2003b), as well as the value extrapolated from the WMAP source counts in Q band, $`(15.0\pm 1.4)\times 10^3\mu \mathrm{K}^2\mathrm{sr}`$ (Bennett et al., 2003a). In V band, $`c_{\mathrm{src}}=(4.5\pm 4)\times 10^3\mu \mathrm{K}^2\mathrm{sr}`$.

The source number counts, angular power spectrum, and bispectrum measure the first-, second-, and third-order moments of $`dN/dS`$, respectively. The good agreement among these three different estimates of $`c_{\mathrm{src}}`$ indicates the validity of the estimate of the effects of the residual point sources in Q band. There is no visible contribution to the angular power spectrum from the sources in V and W bands. We conclude that our understanding of the amplitude of the residual point sources is satisfactory for the analysis of the angular power spectrum not to be contaminated by the sources.

## 4 CONCLUSIONS

We use cubic (bispectrum) statistics and the Minkowski functionals to measure non-Gaussian fluctuations in the WMAP 1-year sky maps. The cubic statistic \[Eq.(7)\] and the Minkowski functionals place limits on the non-linear coupling parameter $`f_{\mathrm{NL}}`$, which characterizes the amplitude of a quadratic term in the Bardeen curvature perturbations \[Eq. (1)\]. It is important to remove the best-fit foreground templates from the WMAP maps in order to reduce the non-Gaussian Galactic foreground emission. The cubic statistic measures phase correlations of temperature fluctuations to find the best estimate of $`f_{\mathrm{NL}}`$ from the foreground-removed, weighted average of Q$`+`$V$`+`$W maps as $`f_{\mathrm{NL}}=38\pm 48`$ (68%) and $`58<f_{\mathrm{NL}}<134`$ (95%). The Minkowski functions measure morphological structures to find $`f_{\mathrm{NL}}=22\pm 81`$ (68%) and $`f_{\mathrm{NL}}<139`$ (95%), in good agreement with the cubic statistic. These two completely different statistics give consistent results, validating the robustness of our limits. Our limits are 20$``$30 times better than the previous ones (Komatsu et al., 2002; Santos et al., 2002; Cayón et al., 2002), and constrain the relative contribution from the non-linear term to the r.m.s. amplitude of $`\mathrm{\Phi }`$ to be smaller than $`2\times 10^5`$ (95%), much smaller than the limits on systematic errors in the WMAP observations. This validates that the angular power spectrum can fully characterize statistical properties of the WMAP CMB sky maps. We conclude that the WMAP 1-year data do not show evidence for significant primordial non-Gaussianity of the form in equation (1). Our limits are consistent with predictions from inflation models based upon a slowly rolling scalar field, $`\left|f_{\mathrm{NL}}\right|=10^210^1`$. The span of all non-Gaussian models, however, is large, and there are models which cannot be parametrized by equation (1) (e.g., Bernardeau & Uzan (2002b, a)). Other forms such as multi-field inflation models and topological defects will be tested in the future.

The non-Gaussianity also affects the dark-matter halo mass function $`dn/dM`$, since the massive halos at high redshift are sensitive to changes in the tail of the distribution function of density fluctuations. Our limits show that the number of clusters that would be newly found at $`z=1`$ for $`M<10^{15}M_{}`$ should be within $`{}_{10}{}^{+40}\%`$ of the value predicted from the Gaussian theory. At higher redshifts, however, much larger effects are still allowed. The number counts $`dN/dz`$ at $`z=3`$ with the limiting mass of $`3\times 10^{14}M_{}`$ can be reduced by a factor of 2, or increased by more than a factor of 3. Since the SZ angular power spectrum is primarily sensitive to massive halos at $`z1`$ where the impact of non-Gaussianity is constrained to be within 10%, a measurement of $`\sigma _8`$ from the SZ angular power spectrum is changed by no more than 2%. Our results on $`dn/dM`$ derived in this paper should be taken as the current observational limits to non-Gaussian effects on $`dn/dM`$. In other words, this is the uncertainty that we currently have in $`dn/dM`$ when the assumption of Gaussian fluctuations is relaxed.

The limits on $`f_{\mathrm{NL}}`$ will improve as the WMAP satellite acquires more data. Monte Carlo simulations show that the 4-year data will achieve 95% limit of 80. This value will further improve with a more proper pixel-weighting function that becomes the uniform weighting in the signal-dominated regime (large angular scales) and becomes the $`N^1`$ weighting in the noise-dominated regime (small angular scales). There is little hope of testing the expected levels of $`f_{\mathrm{NL}}=10^210^1`$ from simple inflation models, but some non-standard models can be excluded.

We have detected non-Gaussian signals arising from the residual radio point sources left unmasked at Q band, characterized by the reduced point-source angular bispectrum $`b_{\mathrm{src}}=(9.5\pm 4.4)\times 10^5\mu \mathrm{K}^3\mathrm{sr}^2`$, which, in turn, gives the point-source angular power spectrum $`c_{\mathrm{src}}=(15\pm 6)\times 10^3\mu \mathrm{K}^2\mathrm{sr}`$. This value agrees well with those from the source number counts (Bennett et al., 2003a) and the angular power spectrum analysis (Hinshaw et al., 2003b), giving us confidence on our understanding of the amplitude of the residual point sources. Since $`b_{\mathrm{src}}`$ directly measures $`c_{\mathrm{src}}`$ without relying on extrapolations, any CMB experiments which suffer from the point-source contamination should use $`b_{\mathrm{src}}`$ to quantify $`c_{\mathrm{src}}`$ to obtain an improved estimate of the CMB angular power spectrum for the cosmological-parameter determinations.

Hinshaw et al. (2003b) found that the best-fit power spectrum to the WMAP temperature data has a relatively large $`\chi ^2`$ value, corresponding to a chance probability of 3%. While still acceptable fit, there may be missing components in the error propagations over the Fisher matrix. Since the Fisher matrix is the four-point function of the temperature fluctuations, those missing components (e.g., gravitational lensing effects) may not be apparent in the bispectrum, the three-point function. The point-source non-Gaussianity contributes to the Fisher matrix by only a negligible amount, as it is dominated by the Gaussian instrumental noise. Non-Gaussianity in the instrumental noise due to the $`1/f`$ striping may have additional contributions to the Fisher matrix; however, since the Minkowski functionals, which are sensitive to higher-order moments of temperature fluctuations and instrumental noise, do not find significant non-Gaussian signals, non-Gaussianity in the instrumental noise is constrained to be very small.

The WMAP mission is made possible by the support of the Office of Space Sciences at NASA Headquarters and by the hard and capable work of scores of scientists, engineers, technicians, machinists, data analysts, budget analysts, managers, administrative staff, and reviewers. LV is supportd by NASA through Chandra Fellowship PF2-30022 issued by the Chandra X-ray Observatory center, which is operated by the Smithsonian Astrophysical Observatory for an on behalf of NASA under contract NAS8-39073.

## Appendix A SIMULATING CMB SKY MAPS FROM PRIMORDIAL FLUCTUATIONS

In this appendix, we describe how to simulate CMB sky maps from generic primordial fluctuations. As a specific example, we choose to use the primordial Bardeen curvature perturbations $`\mathrm{\Phi }(𝐱)`$, which generate CMB anisotropy at a given position of the sky $`\mathrm{\Delta }T(\widehat{𝐧})`$ as (Komatsu et al., 2003)

$$\mathrm{\Delta }T(\widehat{𝐧})=T_0\underset{lm}{}Y_{lm}(\widehat{𝐧})r^2𝑑r\mathrm{\Phi }_{lm}(r)\alpha _l(r),$$
(A1)
where $`\mathrm{\Phi }_{lm}(r)`$ is the harmonic transform of $`\mathrm{\Phi }(𝐱)`$ at a given comoving distance $`r\left|𝐱\right|`$, $`\mathrm{\Phi }_{lm}(r)d^2\widehat{𝐧}\mathrm{\Phi }(r,\widehat{𝐧})Y_{lm}^{}(\widehat{𝐧})`$, and $`\alpha _l(r)`$ was defined previously (Eq.). We can instead use isocurvature fluctuations or a mixture of the two. Equation (A1) suggests that $`\alpha _l(r)`$ is a transfer function projecting $`\mathrm{\Phi }(𝐱)`$ onto $`\mathrm{\Delta }T(\widehat{𝐧})`$ through the integral over the line of sight. Since $`\alpha _l(r)`$ is just a mathematical function, we pre-compute and store it for a given cosmology, reducing the computational time of a batch of simulations. We can thus use or extend equation (A1) to compute $`\mathrm{\Delta }T(\widehat{𝐧})`$ for generic primordial fluctuations.

We simulate CMB sky maps using a non-Gaussian model of the form in equation (1) as follows. (1) We generate $`\stackrel{~}{\mathrm{\Phi }}_\mathrm{L}(𝐤)`$ as a Gaussian random field in Fourier space for a given initial power spectrum $`P(k)`$, and transform it back to real space to obtain $`\mathrm{\Phi }_\mathrm{L}(𝐱)`$. (2) We transform from Cartesian to spherical coordinates to obtain $`\mathrm{\Phi }_\mathrm{L}(r,\widehat{𝐧})`$, compute its harmonic coefficients $`\mathrm{\Phi }_{lm}(r)`$, and obtain a temperature map of the Gaussian part $`\mathrm{\Delta }T_\mathrm{\Phi }(\widehat{𝐧})`$ by integrating equation (A1). (3) We repeat this procedure for $`\mathrm{\Phi }_\mathrm{L}^2(𝐱)V_x^1d^3𝐱\mathrm{\Phi }_\mathrm{L}^2(𝐱)`$ to obtain a temperature map of the non-Gaussian part $`\mathrm{\Delta }T_{\mathrm{\Phi }^2}(\widehat{𝐧})`$. (4) By combining these two temperature maps, we obtain non-Gaussian sky maps for any values of $`f_{\mathrm{NL}}`$,

$$\mathrm{\Delta }T(\widehat{𝐧})=\mathrm{\Delta }T_\mathrm{\Phi }(\widehat{𝐧})+f_{\mathrm{NL}}\mathrm{\Delta }T_{\mathrm{\Phi }^2}(\widehat{𝐧}).$$
(A2)
We do not need to run many simulations individually for different values of $`f_{\mathrm{NL}}`$, but run only twice to obtain $`\mathrm{\Delta }T_\mathrm{\Phi }(\widehat{𝐧})`$ and $`\mathrm{\Delta }T_{\mathrm{\Phi }^2}(\widehat{𝐧})`$ for a given initial random number seed. Also, we can combine $`\mathrm{\Delta }T_\mathrm{\Phi }(\widehat{𝐧})`$ for one seed with $`\mathrm{\Delta }T_{\mathrm{\Phi }^2}(\widehat{𝐧})`$ for the other to make realizations for a particular kind of two-field inflation models. We can apply the same procedure to isocurvature fluctuations with or without $`\mathrm{\Phi }(𝐱)`$ correlations.

We need the simulation box of the size of the present-day cosmic horizon size $`L_{\mathrm{box}}=2c\tau _0`$, where $`\tau _0`$ is the present-day conformal time. For example, $`L_{\mathrm{box}}20h^1\mathrm{Gpc}`$ is needed for a flat universe with $`\mathrm{\Omega }_\mathrm{m}=0.3`$, whereas we need spatial resolution of at least $`20h^1\mathrm{Mpc}`$ to resolve the last-scattering surface accurately. From this constraint the number of grid points is at least $`N_{\mathrm{grid}}=1024^3`$, and the required amount of physical memory to store $`\mathrm{\Phi }(𝐱)`$ is at least 4.3 GB. Moreover, when we simulate a sky map having 786 432 pixels at $`nside=256`$, we need 1.6 GB to store a field in spherical coordinate $`\mathrm{\Phi }(r,\widehat{𝐧})`$, where the number of $`r`$ evaluated for $`N_{\mathrm{grid}}=1024^3`$ is 512. Since our algorithm for transforming Cartesian into spherical coordinates requires another 1.6 GB, in total we need at least 7.5 GB of physical memory to simulate one sky map.

We have generated 300 realizations of non-Gaussian sky maps with $`N_{\mathrm{grid}}=1024^3`$ and $`nside=256`$. It takes 3 hours on 1 processor of SGI Origin 300 to simulate $`\mathrm{\Delta }T_\mathrm{\Phi }(\widehat{𝐧})`$ and $`\mathrm{\Delta }T_{\mathrm{\Phi }^2}(\widehat{𝐧})`$. We have used 6 processors to simulate 300 maps in one week. Figure 7 shows the one-point probability density function (PDF) of temperature fluctuations measured from simulated non-Gaussian maps (without noise and beam smearing) compared with the r.m.s. scatter of Gaussian realizations. We find it difficult for the PDF alone to distinguish non-Gaussian maps of $`\left|f_{\mathrm{NL}}\right|<500`$ from Gaussian maps, whereas the cubic statistic $`S_{\mathrm{prim}}`$ (eq.) can easily detect $`f_{\mathrm{NL}}=100`$ in the same data sets.

We measure $`f_{\mathrm{NL}}`$ on the simulated maps using $`S_{\mathrm{prim}}`$ to see if it can accurately recover $`f_{\mathrm{NL}}`$. Similar tests show the Minkowski functionals to be unbiased and able to discriminate different $`f_{\mathrm{NL}}`$ values at levels consistent with the quoted uncertainties. We also measure the point-source angular bispectrum $`b_{\mathrm{src}}`$ to see if it returns null values as the simulations do not contain point sources. We have included noise properties and window functions in the simulations. Figure 8 shows histograms of $`f_{\mathrm{NL}}`$ and $`b_{\mathrm{src}}`$ measured from 300 simulated maps of $`f_{\mathrm{NL}}=100`$ (solid lines) and $`f_{\mathrm{NL}}=0`$ (dashed lines). Our statistics find correct values for $`f_{\mathrm{NL}}`$ and find null values for $`b_{\mathrm{src}}`$; thus, our statistics are unbiased, and $`f_{\mathrm{NL}}`$ and $`b_{\mathrm{src}}`$ are orthogonal to each other as pointed out by Komatsu & Spergel (2001).

## Appendix B POWER OF THE POINT-SOURCE BISPECTRUM

In this appendix, we test our estimator for $`b_{\mathrm{src}}`$ and $`c_{\mathrm{src}}`$ using simulated Q-band maps of point sources, CMB, and detector noise. The 44 GHz source count model of T98 was used to generate the source populations. The total source count in each realization was fixed to 9043, the number predicted by T98 to lie between $`S_{\mathrm{min}}=0.1\mathrm{Jy}`$ and $`S_{\mathrm{max}}=10\mathrm{Jy}`$. By generating uniform deviates $`u(0,1)`$ and transforming to flux $`S`$ via

$$u=\frac{N(>S_{\mathrm{min}})N(>S)}{N(>S_{\mathrm{min}})N(>S_{\mathrm{max}})},$$
(B1)
we obtain the desired spectrum. The sources were distributed evenly over the sky and convolved with a Gaussian profile approximating the Q-band beam. Flux was converted to peak brightness using the values in Table 8 of Page et al. (2003b). The CMB and noise realizations were not varied between realizations. The goal in this appendix is to prove that our estimator for $`b_{\mathrm{src}}`$ works well and is very powerful in estimating $`c_{\mathrm{src}}`$.

The left panel of Figure 9 compares the measured $`b_{\mathrm{src}}`$ from simulated maps with the expectations of the simulations. Black, dark-gray, and light-gray indicate three different realizations of point sources. The measurements agree well with the expectations at $`S_\mathrm{c}<1.75\mathrm{Jy}`$. They however show significant scatter at $`S_\mathrm{c}>1.75\mathrm{Jy}`$, because our filter for computing $`b_{\mathrm{src}}`$ \[eq. (21)\] does not include contribution from $`c_{\mathrm{src}}`$ to $`\stackrel{~}{C}_l`$, making the filter less optimal in the limit of “too many” unmasked point sources. We can see from the figure that $`c_{\mathrm{src}}`$ at $`S_\mathrm{c}>2\mathrm{Jy}`$ is comparable to or larger than the noise power spectrum for Q band, $`54\times 10^3\mu \mathrm{K}^2\mathrm{sr}`$.

Fortunately this is not a problem in practice, as we can detect and mask those bright sources which contribute significantly to $`\stackrel{~}{C}_l`$. The residual point sources that we cannot detect (therefore we want to quantify using $`b_{\mathrm{src}}`$) should be hidden in the noise having only a small contribution to $`\stackrel{~}{C}_l`$. In this faint-source regime $`b_{\mathrm{src}}`$ works well in measuring the amplitude of residual point sources, offering a promising way for estimating $`c_{\mathrm{src}}`$. The right panel of Figure 9 compares $`c_{\mathrm{src}}`$ estimated from $`b_{\mathrm{src}}`$ \[Eq. (17)\] with the expectations. The agreement is good for $`S_\mathrm{c}<1.75\mathrm{Jy}`$, proving that estimates of $`c_{\mathrm{src}}`$ from $`b_{\mathrm{src}}`$ are unbiased and powerful. Since $`b_{\mathrm{src}}`$ measures $`c_{\mathrm{src}}`$ directly, we can use it for any CMB experiments which suffer from the effect of residual point sources. While we have considered the bispectrum only here, the forth order moment may also be used to increase our sensitivity to the point-source non-Gaussianity (Pierpaoli, 2003).

## Appendix C THE ANGULAR BISPECTRUM FROM A POTENTIAL STEP

A scalar-field potential $`V(\varphi )`$ with features can generate large non-Gaussian fluctuations in CMB by breaking the slow-roll conditions at the location of the features (Kofman et al., 1991; Wang & Kamionkowski, 2000). We estimate the impact of the features by using a scale-dependent $`f_{\mathrm{NL}}`$,

$$f_{\mathrm{NL}}(\varphi )=\frac{5}{24\pi G}\left(\frac{^2\mathrm{ln}H}{\varphi ^2}\right),$$
(C1)
which is calculated from a non-linear transformation between the curvature perturbations in the comoving gauge and the scalar-field fluctuations in the spatially flat gauge (Salopek & Bond, 1990, 1991). This expression does not assume the slow-roll conditions. Although this expression does not include all effects contributing to $`f_{\mathrm{NL}}`$ during inflation driven by a single field (Maldacena, 2002), we assume that an order-of-magnitude estimate can still be obtained.

A sharp feature in $`V(\varphi )`$ at $`\varphi _\mathrm{f}`$ produces a significantly scale-dependent $`f_{\mathrm{NL}}(\varphi )`$ near $`\varphi _\mathrm{f}`$ through the derivatives of $`H`$ in equation (C1). We illustrate the effects of the steps using the potential features proposed by Adams et al. (2001),

$$V(\varphi )=\frac{1}{2}m_\varphi ^2\varphi ^2\left[1+c\mathrm{tanh}\left(\frac{\varphi \varphi _\mathrm{f}}{d}\right)\right],$$
(C2)
which has a step in $`V(\varphi )`$ at $`\varphi _\mathrm{f}`$ with the height $`c`$ and the slope $`d^1`$. Adams et al. (1997) show that the steps are created by a class of supergravity models in which symmetry-breaking phase transitions of many fields in flat directions gravitationally coupled to $`\varphi `$ continuously generate steps in $`V(\varphi )`$ every $`1015`$ $`e`$-folds, giving a chance for a step to exist within the observable region of $`V(\varphi )`$.

It is instructive to evaluate equation (C1) combined with equation (C2) in the slow-roll limit, $`^2\mathrm{ln}H/\varphi ^2\frac{1}{2}^2\mathrm{ln}V/\varphi ^2`$. For $`\left|c\right|1`$, one obtains

$$f_{\mathrm{NL}}(\varphi )\frac{5}{24\pi G}\left(\frac{1}{\varphi ^2}+\frac{c}{d^2}\frac{\mathrm{tanh}x}{\mathrm{cosh}^2x}\right),$$
(C3)
where $`x(\varphi \varphi _\mathrm{f})/d`$. The first term corresponds to a standard, nearly scale-independent prediction giving $`7.4\times 10^3`$ at $`\varphi =3m_{\mathrm{planck}}`$, while the second term reveals a significant scale-dependence. The function $`\mathrm{tanh}x/\mathrm{cosh}^2x`$ is a symmetric odd function about $`x=0`$ with extrema of $`\pm 0.385`$ at $`x\pm 0.66`$. The picture is the following: as $`\varphi `$ rolls down $`V(\varphi )`$ from a positive $`x>0.66`$, $`\varphi `$ gets accelerated at $`x0.66`$, reaches constant velocity at $`x=0`$, decelerates at $`x0.66`$, and finally reaches slow roll at $`x<0.66`$. The ratio of the second term in equation (C3) to the first at the extrema is $`\pm 0.385c(\varphi /d)^2`$. For example, $`c=0.02`$ and $`\varphi /d=300`$ (i.e., $`d=0.01m_{\mathrm{planck}}`$) make the amplitude of the second term 700 times larger than the first, giving $`\left|f_{\mathrm{NL}}\right|5`$ at the extrema. Despite the slow-roll conditions having a tendency to underestimate $`f_{\mathrm{NL}}`$, it is possible to obtain $`\left|f_{\mathrm{NL}}\right|>1`$. Neglecting the first term in equation (C3) and converting $`\varphi `$ for $`k`$, one obtains

$$f_{\mathrm{NL}}(k)\frac{5c}{24\pi Gd^2}h_{\mathrm{step}}(k)\frac{5c}{24\pi Gd^2}\frac{\mathrm{tanh}x_k}{\mathrm{cosh}^2x_k},$$
(C4)
where $`x_kd^1(\varphi /\mathrm{ln}k)_\mathrm{f}(k/k_\mathrm{f}1)=d^1(\dot{\varphi }/H)_\mathrm{f}(k/k_\mathrm{f}1)`$ for $`kk_\mathrm{f}k_\mathrm{f}`$. The slow-roll approximation gives $`x_k(4\pi G\varphi _\mathrm{f}d)^1(k/k_\mathrm{f}1)`$. Finally, following the method of Komatsu & Spergel (2001), we obtain the reduced bispectrum of a potential step model, $`b_{l_1l_2l_3}^{\mathrm{step}}`$, as

$$b_{l_1l_2l_3}^{\mathrm{step}}=2\left(\frac{5c}{24\pi Gd^2}\right)_0^{\mathrm{}}r^2𝑑r\left[\beta _{l_1}(r)\beta _{l_2}(r)\alpha _{l_3}^{\mathrm{step}}(r)+\text{(2 permutations)}\right],$$
(C5)
where $`\beta _l(r)`$ is given by equation (6), and

$$\alpha _l^{\mathrm{step}}(r)\frac{2}{\pi }k^2𝑑kh_{\mathrm{step}}(k)g_{\mathrm{T}l}(k)j_l(kr).$$
(C6)
The amplitude is thus proportional to $`c/d^2`$: a bigger (larger $`c`$) and steeper (smaller $`d`$) step gives a larger bispectrum. The steepness affects the amplitude more, because the non-Gaussianity is generated by breaking the slow-roll conditions.

Since $`b_{l_1l_2l_3}^{\mathrm{step}}`$ linearly scales as $`c`$ for a fixed $`d`$, we can fit for $`c`$ by using exactly the same method as for the scale-independent $`f_{\mathrm{NL}}`$, but with $`\alpha _l(r)`$ in equation (3) replaced by $`\alpha _l^{\mathrm{step}}(r)`$. The exact form of the fitting parameter in the slow-roll limit is $`5c/(24\pi Gd^2)`$. A reason for the similarity between the two models in methods for the measurement is explained as follows. Komatsu et al. (2003) have shown that $`B(\widehat{𝐧},r)`$ \[Eq.(4)\] is a Wiener-filtered, reconstructed map of the primordial fluctuations $`\mathrm{\Phi }(\widehat{𝐧},r)`$. Our cubic statistic \[Eq.(7)\] effectively measures the skewness of the reconstructed $`\mathrm{\Phi }`$ field, maximizing the sensitivity to the primordial non-Gaussianity. One of the three maps comprising the cubic statistic is however not $`B(\widehat{𝐧},r)`$, but $`A(\widehat{𝐧},r)`$ given by equation (3). This map defines what kind of non-Gaussianity to look for, or more detailed form of the bispectrum. For the potential step case, $`A_{\mathrm{step}}(\widehat{𝐧},r)`$ made of $`\alpha _l^{\mathrm{step}}(r)`$ picks up the location of the step to measure $`5c/(24\pi Gd^2)`$ near $`k_\mathrm{f}`$, while for the form in equation (1), $`A(\widehat{𝐧},r)`$ explores all scales on equal footing to measure the scale-independent $`f_{\mathrm{NL}}`$.

The distinct features in $`k`$ space are often smeared out in $`l`$ space via the projection. This effect is estimated from equation (C6) as follows. The function $`h_{\mathrm{step}}(k)`$ near $`k_\mathrm{f}`$ is accurately approximated by $`h_{\mathrm{step}}(k)0.385\mathrm{sin}(2x_k)`$, which has a period of $`\mathrm{\Delta }k=4\pi ^2G\varphi _\mathrm{f}dk_\mathrm{f}`$. On the other hand, the radiation transfer function $`g_{\mathrm{T}l}(k)`$ behaves as $`j_l(kr_{})`$ where $`r_{}`$ is the comoving distance to the photon decoupling epoch, and $`g_{\mathrm{T}l}(k)j_l(kr)`$ behaves as $`j_l^2(kr_{})`$ (the integral is very small when $`rr_{}`$). The oscillation period of this part is thus $`\mathrm{\Delta }k=\pi /r_{}`$ for $`kr_{}>l`$. A ratio of the period of $`h_{\mathrm{step}}(k)`$ to that of $`g_{\mathrm{T}l}(k)j_l(kr)`$ is then estimated as $`4\pi G\varphi _\mathrm{f}dr_{}k_\mathrm{f}(l_\mathrm{f}/3)(d/0.01m_{\mathrm{planck}})(\varphi _\mathrm{f}/3m_{\mathrm{planck}})`$, where $`l_\mathrm{f}k_\mathrm{f}r_{}`$ is the angular wave number of the location of the step. We thus find that $`h_{\mathrm{step}}(k)`$ oscillates much more slowly than the rest of the integrand in equation (C6) for $`l_\mathrm{f}1`$.

What does it mean? It means that the results would look as if there were two distinct regions in $`l`$ space where $`f_{\mathrm{NL}}`$ is very large: a positive $`f_{\mathrm{NL}}`$ is found at $`l<l_\mathrm{f}`$ and a negative one at $`l>l_\mathrm{f}`$. The estimated location is $`l/l_\mathrm{f}1\pm 0.66(4\pi G\varphi _\mathrm{f}d)1\pm 0.2(d/0.01m_{\mathrm{planck}})(\varphi _\mathrm{f}/3m_{\mathrm{planck}})`$; thus, the positive and negative regions are separated in $`l`$ by only 40%, making the detection difficult when many $`l`$ modes are combined to improve the signal-to-noise ratio. The two extrema would cancel out to give only small signals. In other words, it is still possible that non-Gaussianity from a potential step is “hidden” in our measurements shown in Figure 1. Note that the cancellation occurs because of the point symmetry of $`h_{\mathrm{step}}(k)`$ about $`k=k_\mathrm{f}`$. If the function has a knee instead of a step, then the cancellation does not occur and there would be a single region in $`l`$ space where $`\left|f_{\mathrm{NL}}\right|`$ is large (Wang & Kamionkowski, 2000). Note that our estimate in this Appendix was based upon equation (C3), which uses the slow-roll approximations. While instructive, since the slow-roll approximations break down near the features, our estimate may not be very accrate. One needs to integrate the equation of motion of the scalar field to evaluate equation (C1) for more accurate estimations of the effect.


format: text

text:
# Relativistic many-body calculations of excitation energies and transition rates from core-excited states in copperlike ions

## I Introduction

This work further develops the application of the relativistic many-body perturbation theory (RMBPT) to the studies of atomic characteristics of particle-particle-hole excitations of closed-shell ions. Recently, RMBPT calculations of energies Safronova et al. (2002a) and transition rates Safronova et al. (2002b) in Na-like ions have been performed. The present paper focuses on the RMBPT calculations of energies and transition rates in Cu-like ions. These transitions form satellite lines to the brightest Ni-like ions and are of a great importance for M-shell diagnostics of heavy ions.

The second-order RMBPT calculations for Cu-like ions start from a $`1s^22s^22p^63s^23p^63d^{10}`$ \[Ni\] Dirac-Fock potential. All possible $`3l`$ holes and $`4l4l^{}`$ particles leading to 563 odd-parity and 571 even-parity $`3l^14l^{}4l\mathrm{"}(J)`$ states are considered. The energies of the 1134 states and transition probabilities of the 2294 electric-dipole lines are calculated for Cu-like ions with $`Z`$ = 30-100. Transitions from the $`(3s^23p^63d^94l4l^{})`$, $`(3s^23p^53d^{10}4l4l^{})`$, and $`(3s3p^63d^{10}4l4l^{})`$ states to the ground ($`3s^23p^63d^{10}4s`$) or singly excited ($`3s^23p^63d^{10}4l`$) states (with $`l=p,d,f`$) form satellite lines to the bright electric-dipole (E1) lines in Ni-like ions. These core-excited states (or often called doubly-excited states) in copperlike ions have been studied extensively both experimentally and theoretically in the past 15-20 years.

Experimentally, these configurations have been studied by photon and electron emission spectroscopy. To our knowledge, the first measurements of $`3d4p`$ transitions in Cu-like W and Tm were done by Klapisch et al. (1981) by classification of x-ray spectra from laser produced plasmas in the range 6 - 9 Å. It was shown that most of the Cu-like $`3d4p`$ line radiation came from $`(3s^23p^63d^{10}4s3s^23p^63d^94s4p)`$, $`(3s^23p^63d^{10}4p3s^23p^63d^94p^2)`$, $`(3s^23p^63d^{10}4d3s^23p^63d^94p4d)`$, and $`(3s^23p^63d^{10}4f3s^23p^63d^94p4f)`$ transitions. Wavelengths and transition probabilities were calculated in Klapisch et al. (1981) by the relativistic parametric potential method Klapisch et al. (1977). The same method was extended by Mandelbaum et al. (1983) to study x-ray spectra from laser produced plasmas of atoms from Tm ($`Z`$ = 69) up to Pt ($`Z`$ = 78). An extended analysis of x-ray spectra of laser produced gold plasmas has been performed by Busquet et al. (1985) including $`3d4f`$, $`3d4p`$, and $`3p4s`$ transitions in Ni-, Co-, Cu-, Zn-, and Ga-like Au. The Ni-, Co-, Cu-, Zn-, Ga-, Ge-, and As-like isoelectronic sequences have been considered in Ref. Tragin et al. (1988) to investigate the x-ray spectra of laser-irradiated elements from tantalum ($`Z`$ = 73) to lead ($`Z`$ = 82). In addition to the above mentioned isoelectronic sequences, the Fe- and Mn-like states have been included by Zigler et al. (1994) to analyze x-ray spectra emmited by laser-produced plasmas of lanthanum ($`Z`$ = 57) and praseodymium ($`Z`$ = 59). The wavelengths and transition probabilities have been calculated in Busquet et al. (1985); Tragin et al. (1988); Zigler et al. (1994) by the relativistic parametric potential method Klapisch et al. (1977). Investigation of the x-ray spectrum emitted by laser produced barium plasmas has been recently described by Doron et al. (1998a, b). The RELAC relativistic atomic code Klapisch et al. (1977) has been used to identify $`3dnl`$ ($`n`$ = 4 to 8), $`3p4s`$, and $`3p4d`$ transitions of Ni-like Ba28+ and corresponding satellite transitions in neighboring ionization states: Fe-, Co-, Cu-, Zn-, Ga-, and Ge-like ions. The X-ray spectrum emitted by a laser produced cerium plasma in the range 7.5 - 12 Å has been recently investigated in detail by Doron et al. (2000, 2001). The RELAC computer code Klapisch et al. (1977) has been used to study x-ray spectra from highly charged tungsten ions in tokamak plasmas in the range 7 - 10 Å Neu et al. (1997).

In the present paper, RMBPT is implemented to determine energies of $`(3s^23p^63d^94l4l^{})`$, $`(3s^23p^53d^{10}4l4l^{})`$, and $`(3s3p^63d^{10}4l4l^{})`$ states for Cu-like ions with nuclear charges in the range of $`Z`$ = 30 –100. The calculations are carried out to second order in perturbation theory and include second-order Coulomb and Breit interactions. Corrections for the frequency-dependent Breit interaction are taken into account in the lowest order. Screened self-energy and vacuum polarization data given by Kim et al. (1991) are used to determine the QED correction.

RMBPT is used to evaluate matrix elements, line strengths, and transition rates for 2294 allowed and forbidden electric-dipole transitions between the 1332 even-parity core-excited states \[Ni\]($`3d^1+3s^1)(4s^2+4s4d+4p^2+4d^2+4p4f+4f^2)`$ \+ \[Ni\]($`3p^1)(4s4p+4s4f+4p4d+4d4f)`$, and the singly excited \[Ni\]($`4p+4f`$) states and the 962 odd-parity core-excited states \[Ni\]($`3p^1)(4s^2+4s4d+4p^2+4d^2+4p4f+4f^2)`$ \+ \[Ni\]($`3d^1+3s^1)(4s4p+4s4f+4p4d+4d4f)`$ and the ground state \[Ni\]($`4s`$) together with the singly excited \[Ni\]($`4d`$) states in Cu-like ions with nuclear charges ranging from $`Z`$ = 30 to 100. Retarded E1 matrix elements are evaluated in both length and velocity forms. These calculations start from a \[Ni\] Dirac-Fock potential. First-order perturbation theory is used to obtain intermediate coupling coefficients and second-order RMBPT is used to determine transition matrix elements. The transition energies employed in the calculations of line strengths and transition rates are derived from second-order RMBPT.

## II Method

Details of the RMBPT method were presented in Ref. Safronova et al. (2002a) for calculation of energies of particle-particle-hole states and in Ref. Safronova et al. (2002b) for calculation of radiative transition rates from particle-particle-hole state to one-particle state. Differences between calculations for Na-like and Cu-like ions are due to the increased size of the model space ($`4l^{}4l^{\prime \prime }3l^1`$ instead of $`3l^{}3l^{\prime \prime }2l^1`$) and differences in the Dirac-Fock potential ($`1s^22s^22p^63s^23p^63d^{10}`$ instead of $`1s^22s^22p^6`$), leading to 1134 states instead of 106 and more laborious numerical calculations.

As a first step, we determine and store the single-particle contributions to the energies for five $`n`$=3 hole states ($`3s`$, $`3p_{1/2}`$, $`3p_{3/2}`$, $`3d_{3/2}`$, and $`3d_{5/2}`$) and the seven $`n`$=4 valence states ($`4s`$, $`4p_{1/2}`$, $`4p_{3/2}`$, $`4d_{3/2}`$, $`4d_{5/2}`$, $`4f_{5/2}`$, and $`4f_{7/2}`$) in lowest, first, and second orders. Next, we evaluate and store the 664 two-particle $`4l4l^{}J|H^{\mathrm{eff}}|4l^{\prime \prime }4l^{\prime \prime \prime }J`$ matrix elements and the 1127 hole-particle $`3l4l^{}J|H^{\mathrm{eff}}|3l^{\prime \prime }4l^{\prime \prime \prime }J`$ matrix elements of the effective Hamiltonian in first and second orders. It should be noted that these one-particle, two-particle, and hole-particle matrix elements were used previously to evaluate energies of the $`4l4l^{}`$ levels in Zn-like ions Safronova (2000) and energies of the $`3l^14l^{}`$ levels in Ni-like ions Safronova et al. (2000). Finally, second-order particle-particle-hole matrix elements are evaluated Safronova et al. (2002a). Combining these data using the method described below, we calculate one-, two-, and three-body contributions to the energies of Cu-like ions.

The calculations are carried out using sets of basis Dirac-Hartree-Fock (DHF) orbitals. The orbitals used in the present calculation are obtained as linear combinations of B-splines. These B-spline basis orbitals are determined using the method described in Ref. Johnson et al. (1988). Forty B-splines of order of eight for each single-particle angular momentum state are used and all orbitals with orbital angular momentum $`l7`$ are included in the basis set.

### II.1 Model space

The model space for core-excited $`4l4l^{}3l^1`$ states of copperlike ions includes 563 odd-parity states consisting of 78 $`J`$=1/2 states, 131 $`J`$=3/2 states, 143 $`J`$=5/2 states, 125 $`J`$=7/2 states, and 86 $`J`$=9/2 states. Additionally, there are 571 even-parity states consisting of 78 $`J`$=1/2 states, 131 $`J`$=3/2 states, 148 $`J`$=5/2 states, 125 $`J`$=7/2 states, and 89 $`J`$=9/2 states. The distribution of some of the 1134 states in the model space is summarized in Table I of the accompanying EPAPS document EPA .

### II.2 Energy-matrix elements

The evaluation of the second-order energies for the $`4l4l^{}(J_1)3l^{\prime \prime }(J)`$ states in Cu-like ions follows the pattern of the corresponding calculation for Zn-like and Ni-like ions given in Refs. Safronova (2000); Safronova et al. (2000). In particular, we use the second-order one- and two-particle matrix elements for Zn-like ions calculated in Safronova (2000) and hole-particle matrix elements for Ni-like ions calculated in Safronova et al. (2000), but recoupled as described below, to obtain the one- and two-particle contributions for Cu-like ions. We will discuss how these matrix elements are combined to obtain the one- and two-particle contributions to energies of Cu-like ions. We refer the reader to Ref. Safronova (2000); Safronova et al. (2000) for a discussion of the how the basic one- and two-particle matrix elements were evaluated. An intrinsic particle-particle-hole diagram also contributes to the second-order energy for Cu-like ions. It should be noted that the angular part of the expression for the particle-particle-hole diagram differs from the corresponding expression for the three-particle matrix elements given in Ref. Safronova et al. (1996). A detailed discussion of this difference is given in Ref. Safronova et al. (2002a).

Table 1 (see also Table II of EPA ) provides an illustration of various contributions to the second-order energies for the special case of Cu-like tungsten, $`Z`$ = 74. In this table, we show the one-body, two-body and three-body second-order Coulomb contributions to the energy matrix labelled as $`E_i^{(2)}`$, $`i`$ = 1, 2, 3, The one-body second-order energy, $`E_1^{(2)}`$ is obtained as the sum of the two valence and one hole $`E_v^{(2)}`$ energies. The values of $`E_1^{(2)}`$ are non-zero only for diagonal matrix elements. Even for odd-parity states with $`J`$=1/2 there are 78 diagonal and 6006 non-diagonal matrix elements for $`4l4l^{}(J_1)3l^{\prime \prime }(J)`$ core-excited states. We list data only for the first ten diagonal matrix elements of odd-parity states with $`J`$=1/2 in Table 1 (a more complete set of data is given in Table II of EPA ). It can be seen from the table that two-body and three-body second-order contributions are positive, when the one-body contributions are negative. The three-body contributions give about 20% in the total second-order contributions. The values of the $`E_3^{(2)}`$ and $`E_2^{(2)}`$ non-diagonal contributions are smaller than values of diagonal contributions by factor of 3-5.

After evaluating the energy matrices, eigenvalues and eigenvectors are calculated for states with given values of $`J`$ and parity. There are two possible methods to carry out the diagonalization: (a) diagonalize the sum of zeroth- and first-order matrices, then calculate the second-order contributions using the resulting eigenvectors; or (b) diagonalize the sum of the zeroth-, first- and second-order matrices together. Following Ref. Safronova et al. (1996), we choose the second method here.

Table 2 lists the energies of ten excited states of W45+ from Table 1 including the total energies $`E_{\mathrm{tot}}`$. The latter is the sum of the following contributions: $`E^{(0+1)}`$ = $`E^{(0)}`$ +$`E^{(1)}`$+$`B^{(1)}`$, the second-order Coulomb and Breit energy $`E^{(2)}`$, and the QED correction $`E_{\mathrm{LAMB}}`$. The QED correction is approximated as the sum of the one-electron self energy and the first-order vacuum-polarization energy. The screened self-energy and vacuum polarization data given by Kim et al. Kim et al. (1991) are used to determine the QED correction $`E_{\mathrm{LAMB}}`$ (see, for detail Ref. Safronova et al. (2002a)). As can be seen, the second order correction contributes to the total energy from 0.4 % for lowest levels up to 8 % for high-excited levels. The levels in this table (see also Table III of EPA ) could be divided into groups corresponding to excited $`4l4l^{}`$ states and $`3l`$ hole states.

### II.3 Dipole matrix element

We designate the first-order dipole matrix element by $`Z^{(1)}`$, the Coulomb-correction to the second-order matrix element $`Z^{(2)}`$, and the second-order Breit correction $`B^{(2)}`$. The evaluation of $`Z^{(1)}`$, $`Z^{(2)}`$, and $`B^{(2)}`$ for Cu-like ions follows the pattern of the corresponding calculation for Na-like ions in Ref. Safronova et al. (2002b). These matrix elements are calculated in both length and velocity gauges.

Table 3 lists values of uncoupled first-and second-order dipole matrix elements $`Z^{(1)}`$, $`Z^{(2)}`$, $`B^{(2)}`$, together with derivative terms $`P^{(\mathrm{derv})}`$ for Cu-like tungsten, $`Z`$ = 74. For simplicity, only the values for the selected dipole transitions between odd-parity states with $`J`$ = 1/2 and the ground $`4s`$ and excited $`4d_{3/2}`$ states are presented. The more comprehensive set of data is given in Table IV of EPA ). The derivative terms shown in Table 3 arise because transition amplitudes depend on energy, and the transition energy changes order-by-order in RMBPT calculations. Both length ($`L`$) and velocity ($`V`$) forms are given for the matrix elements. We find that the first-order matrix elements $`Z_L^{(1)}`$ and $`Z_V^{(1)}`$ differ by 10%; the $`L`$ \- $`V`$ differences between second-order matrix elements are much larger for some transitions. The term $`P^{(\mathrm{derv})}`$ in length form almost equals $`Z^{(1)}`$ in length form but in velocity form is smaller than $`Z^{(1)}`$ in length form by three to four orders of magnitude.

Values of line strengths calculated as a square of coupled reduced matrix element Safronova et al. (2002b) in length and velocity forms are given in Table 4 for the selected dipole transitions between odd-parity states with $`J`$ = 1/2 and the ground $`4s`$ and excited $`4d_{3/2}`$ states. A more complete set of data is given Table V of EPA . Although we use an intermediate-coupling scheme, it is nevertheless convenient to label the physical states using the $`jj`$ scheme. We see that $`L`$ and $`V`$ forms of the coupled matrix elements in Table 4 differ only in the third or fourth digits. These $`L`$$`V`$ differences arise because our RMBPT calculations start with a non-local Dirac-Fock (DF) potential. If we were to replace the DF potential by a local potential, the differences would disappear completely. The first two columns in Table 4 show $`L`$ and $`V`$ values of line strengths calculated in the first-order approximation which differ by factor of 10. Thr last two columns indicate that including of the second-order contribution almost removes the $`LV`$ differences.

## III Results and comparisons with other theory and experiment

We calculate energies of core-excited even-parity $`4lj4l^{}j^{}(J_1)3l^{\prime \prime }j^{\prime \prime }`$ as well as the singly excited $`4lj`$ states in Cu-like ions with $`Z`$ = 30 -100. Reduced matrix elements, line strengths, oscillator strengths, and transition rates are also determined for electric dipole transitions between the above mentioned core-excited and singly excited states in Cu-like ions for the same range of $`Z`$.

Table 5 lists theoretical data for selected transitions from core-excited levels to the ground $`4s`$ and a singly excited $`4p_{1/2}`$ levels with largest values of transition rates in W45+. Also, the comparison of theoretical data produced by the different methods and codes is inlcuded. Specifically, three values for wavelengths: $`\lambda _{\mathrm{MCDF}}`$, $`\lambda _{\mathrm{RMBPT}}`$, and $`\lambda _{\mathrm{cowan}}`$ are compared. The first values, $`\lambda _{\mathrm{MCDF}}`$, were obtained as the first-order RMBPT values. It can be seen from Table 5 that the results obtained by cowan code cow , $`\lambda _{\mathrm{cowan}}`$, better agree with the MCDF results, $`\lambda _{\mathrm{MCDF}}`$ than with the RMBPT results, $`\lambda _{\mathrm{RMBPT}}`$.

cowan code cow  gives results which are generally in good agreement with experimental energies by scaling the electrostatic Slater parameters to include the correlation effects (Refs. Sataka et al. (2002)) and here the scaling factor of 0.85 was used. Also, in cowan code the $`LS`$ coupling scheme is implemented and corresponding $`LS`$ designations were added to $`jj`$ designations in all tables which include comparisons with theoretical data produced by cowan code.

A comprehensive set of theoretical data and comparisons for transitions with largest values of transition rates in W45+ similar to presented in Table 5 is given in Table VI of the accompanying EPAPS document EPA . It includes transitions from core-excited levels not only to the ground $`4s`$ and to one singly excited level $`4p_1/2`$ as in Table 5 but to all singly excited $`4l_j`$ levels.

The similar comprehensive set of theoretical data and comparisons for transitions between core-excited states $`4lj4l^{}j^{}(J_1)3l^{\prime \prime }J`$ to the ground $`4s`$ and singly excited $`4l_j`$ states is presented for Ce29+ ion in Table VII of Ref. EPA  and for Au50+ ion in Table VIII of Ref. EPA . Comparisons with theoretical data obtained by using RMBPT, MCDF, and cowan codes show that the difference in results can be explained by the second-order corrections to energy and dipole matrix elements included in RMBPT.

The trends of the $`Z`$-dependence of transition rates for the transitions from core-excited odd-parity states with $`J`$ = 1/2 to the ground state $`4s^2S_{1/2}`$ states are presented in Fig. 1 . The trends for the transitions from core-excited even-parity states with $`J`$ = 1/2 to the singly excited state $`4p^2P_{3/2}`$ are shown in Fig. 2. More transition rate figures are given in EPA .

We see from the graphs that transitions with smooth $`Z`$-dependence are rarer than transitions with sharp features. Smooth $`Z`$-dependence occurs for transitions from doublet and quartet core-excited states. Usually, singularities happen in the intermediate interval of $`Z`$ = 40 - 60 when neither $`LS`$ nor $`jj`$ coupling schemes describe the states of these ions properly. One general conclusion that can be derived from these graphs is that the smooth $`Z`$-dependence occurs more frequently for transitions from low-lying core-excited states.

Singularities in the transition-rate curves have two distinct origins: avoided level crossings and zeros in dipole matrix elements. Avoided level crossings result in changes of the dominant configuration of a state at a particular value of $`Z`$ and lead to abrupt changes in the transition rate curves when the partial rates associated with the dominant configurations below and above the crossing point are significantly different. Zeros in transition matrix elements lead to cusp-like minima in the transition rate curves. Examples of each of these two singularity types can be seen in Figs. 1 and 2, as well as in Figs. 1 - 10 presented in EPAPS document EPA .

In Table 6, wavelengths and electric-dipole transition rates along with comparison with experimental data are presented for transitions in Cu-like Ce, W, and Au. The table is limited to identification of experimentally measured transitions given in Refs. Doron et al. (2001); Mandelbaum et al. (1983); Busquet et al. (1985). A more comprehensive set of theoretical data assigned to each experimental transition is presented in Table IX of EPAPS document EPA . We mentioned previously that all possible $`4l_1j_14l_2j_2(J_1)3l_3j_3(J)`$ \- $`4l_j`$ transitions produce 2294 spectrum lines. These lines in Ce29+, W45+, and Au50+ cover the spectral regions from 6.8 – 21.9 Å , 3.5 – 9.8 Å , and 2.9 – 8.7 Å , respectively. A number of spectral lines becomes smaller by a factor of ten when we consider transitions with largest values of transition rates, $`gA_r`$. The number of transitions with $`gA_r>`$ 1013 s-1 is about 200 for Ce29+, and a number of transitions with $`gA_r>`$ 1014 s-1 is about 130 for W45+ and about 190 for Au50+. Those transitions with the largest values of $`gA_r`$ cover the spectral regions smaller by factor of 3-4 than that mentioned for all lines. We can see that the number of predicted spectral lines even with the largest values of $`gA_r`$ still much larger than a number of experimental lines shown in Table 6. Also, the interval between these lines is about 0.01 Å or less which is comparable with the accuracy of experimental measurements. In this case, it could be reasonable to assign not a single transition as in Table 6 but 3-5 transitions to identify the experimental peak, as it is demonstrated in Table IX of EPAPS document EPA .

These atomic data are important in modeling of M-shell radiation spectra of heavy ions generated in electron beam ion trap experiments and in M-shell diagnostics of plasmas. For example, x-ray M-shell spectra of W ions collected at the different energies of the electron beam at LLNL electron beam ion trap indicate the existence of strong Cu-like spectral features Shlyaptseva et al. (2003). Also,the distinct features of x-ray M-shell spectra of W ions in a spectral region from 5 to 6 Å produced by laboratory plasmas are the brightest Ni-like lines and Co- and Cu-like spectral features. In particular, Cu-like satellite lines constitute most of the Cu-like spectral features, and Cu-like autoionization levels make an important contribution in ionization balance calculations Shlyaptseva et al. (2002). Synthetic x-ray spectra of Cu-like W computed with different resolution of 0.02 Å and 0.002 Å are presented in Fig. 3 and 4, respectively. It was assumed that spectral lines have the intensities proportional to the weighted transition probabilities and are fitted with the Gaussian profile. Specifically, the spectrum in Fig. 3 includes $`3l4l^{}`$ transitions and covers the spectral region from 4 Å to 6.4 Å. The most intense peaks at about 5.72 Å and 5.9 Å are formed by $`3d4f`$ transitions and are shown with better resolution in Fig. 4.

## IV Conclusion

In summary, a systematic second-order RMBPT study of the energies and transition rates for \[$`4l_1j_14l_2j_2(J_1)3l_3j_3(J)`$ \- $`4lj`$\] electric-dipole transitions in Cu-like ions with the nuclear charges $`Z`$ ranging from 30 to 100 has been presented. The retarded $`E1`$ matrix elements included correlation corrections from Coulomb and Breit interactions. Both length and velocity forms of the matrix elements were evaluated and small differences (0.4 - 1 %), caused by the non locality of the starting DF potential, were found between the two forms. Second-order RMBPT transition energies were used in the evaluation of transition rates. These calculations were compared with other calculations and with available experimental data. For $`Z40`$, we believe that the present theoretical data are more accurate than other theoretical or experimental data for transitions between the $`4l_1j_14l_2j_2(J_1)3l_3j_3(J)`$ core-excited states and the $`4lj`$ singly excited states in Cu-like ions. The results could be further improved by including third-order correlation corrections.These calculations are presented as a theoretical benchmark for comparison with experiment and theory. In addition, the application of generated atomic data to modeling and interpreting of x-ray M-shell spectra of heavy ions is discussed.

###### Acknowledgements.

The work of W.R.J. was supported in part by National Science Foundation Grant No. PHY-01-39928. U.I.S. acknowledges partial support by Grant No. B516165 from Lawrence Livermore National Laboratory. The work of A.S. was supported by the DOE-NNSA/NV Cooperative Agreement DE-FC08-01NV14050 and Sandia National Laboratories.