AKARI OBSERVATION OF THE FLUCTUATION OF THE NEAR-INFRARED BACKGROUND

Observations of the near-infrared background are believed to be very important for investigating the formation and evolution of the first stars (Santos et al. 2002). A search for the cosmic near-infrared background has been commenced using data from the Diffuse Infrared Background Experiment (DIRBE) on board Cosmic Background Explorer (COBE) and the Near-infrared Spectrometer (NIRS) on board Infrared Telescope in Space (IRTS; Hauser et al. 1998; Cambrésy et al. 2001; Arendt & Dwek 2003; Matsumoto et al. 2005; Levenson et al. 2007). The results have consistently shown that there remains an isotropic emission that can hardly be explained by known sources, as long as the same model is applied for zodiacal light. One possible origin of this excess emission is Population III (Pop III) stars, the first stars of the universe that caused its reionization (Salvaterra & Ferrara 2003; Dwek et al. 2005a). The observed excess emission, however, is much brighter than expected, and requires an extremely high rate of star formation during the Pop III era (Madau & Silk 2005). Furthermore, TeV-γ observations of distant blazars favor a low level of near-infrared extragalactic background light (Dwek et al. 2005b; Aharonian et al. 2006; Raue et al. 2009). It has been pointed out that uncertainties in the modeling of zodiacal light may explain the ambiguity in the result, since zodiacal light is the brightest foreground component in the near-infrared range (Dwek et al. 2005a). In order to obtain firm results on the cosmic near-infrared background, some researchers looked for fluctuations (Kashlinsky et al. 2002; Matsumoto et al. 2005), as zodiacal light is very smooth. Furthermore, fluctuation measurements are important for understanding the structure formation in the Pop III era. Recently, Kashlinsky et al. (2005) reported significant fluctuations at angular scales of 100''–300'' based on Spitzer data and concluded that they were caused by Pop III stars. However, a few studies have claimed that the origin of this fluctuation was the clustering of low-redshift galaxies (Cooray et al. 2007; Chary et al. 2008). Thompson et al. (2007a) detected a significant fluctuation at 1.1 μm and 1.6 μm with the Near-Infrared Camera and Multi-Object Spectrometer (NICMOS) on board Hubble Space Telescope (HST), although their observation was restricted to angular scales smaller than 100''. However, they claimed that the observed fluctuation could be due to faint galaxies at redshift z < 8 (Thompson et al. 2007b). In this paper, we present new observations of fluctuations in sky brightness using the Japanese infrared astronomical satellite AKARI. AKARI has a cold shutter to measure dark current precisely as well as a short wavelength band at 2.4 μm, which provided new observational evidence on the study of the sky fluctuation.

This paper is organized as follows. In Section 2, we provide detailed information on the observation and data reduction procedure. In Section 3, we describe the fluctuation analysis of the near-infrared sky, and we discuss the possible origins of the fluctuation in Section 4. Discussions of the implications of our data in view of the first stars are given in Section 5. The final section summarizes our results.

AKARI is the first Japanese infrared astronomy satellite with a cooled 68 cm aperture telescope, launched into a Sun-synchronous orbit on 2006 February 22 (Murakami et al. 2007). The Infrared Camera (IRC) is one of the instruments on board AKARI (Onaka et al. 2007), and has a frame of 10' × 10' with a pixel scale of 15. The IRC performed a large area survey around the north ecliptic pole (NEP; Lee et al. 2009), and monitored a fixed position in the sky near the NEP to evaluate its performance (Wada et al. 2007). The observed region of the sky was named NEP monitor field (α₂₀₀₀ = 17^h55^m24^s, δ₂₀₀₀ = 66:37:32). Pointed observations were carried out twice a month, and three sets of imaging observations were performed for all wavelength bands of the IRC during one pointed observation. Dithering of ∼10'' was applied between these three sets of imaging observations. Since maneuvering of the satellite sometimes commenced before completing the third set of imaging observations, the last image at 4.1 μm was usually useless. The integration time for one image was 44.41 s, and the data were retrieved with Fowler 4 samplings. Systematic errors due to uncertainty in the absolute calibration were 2.8%, 2.5%, and 3.4%, which were small compared to other errors. The details of observation mode and calibration are shown in the IRC user's manual (Lorente et al. 2008).

In order to detect the fluctuations in sky brightness, we used near-infrared images taken in 14 pointed observations between 2006 September and 2007 March, during which no radiation from the Earth entered the baffle and stray light was negligible. Incidence of earthshine into the telescope tube was detected from the temperature sensor on the baffle and confirmed by the IRC simultaneously as a spurious stray radiation. Furthermore, we did not use images when the attitude control was less stable than usual, as well as those images that suffered from a large number of cosmic-ray hits. After all, we chose 40, 39, and 24 images from 2.4, 3.2, and 4.1 μm bands, respectively, for the analysis (hereafter, we list data at these three bands in order of wavelength). Total integration times were 29.6, 28.9, and 17.8 minutes for each band, which are rather shorter than Spitzer observations (Kashlinsky et al. 2005). AKARI/IRC, however, has a larger field of view, 10' than Spitzer/IRAC, 52. Since pixel throughputs are almost the same, Spitzer/IRAC needs ∼4 times more observing time than AKARI/IRC to observe the same area of the sky with the same signal-to-noise ratio (S/N).

We analyzed the sky image following the IRC user's manual (Lorente et al. 2008) with some improvements. We estimated the dark current of the image field from that of the masked region on the array by applying the new method developed by Tsumura & Wada (2011). This method provides a more reliable estimate of the dark current for those frames suffering from an aftereffect due to strong cosmic-ray hits occurring in the South Atlantic Anomaly region. We also constructed better flat fields using all images obtained during the period as used for the observation of the NEP monitor field. In total, 774, 771, and 516 images were stacked, and significantly better flat fielding was achieved compared to the ordinary pipeline. For the correction of the aspect ratio, we used the "nearest" option as interpolation type in the MAGNIFY task of the Image Reduction and Analysis Facility (IRAF; Tody 1986) to avoid a stripe-like artifact.

Finally, we stacked the images and obtained the maps of the NEP monitor field. Due to constraints on attitude control, the position angles of the images were rotated by 1° per day, and by 200° during the observation period, resulting in circular images with 10' diameters. This improves the quality of the stacked images compared to usual stacking, since flat-field errors, fluctuations in sensitivity between pixels, residuals in dark frame subtraction, fluctuation of the zodiacal light (if any), etc., are reduced, except for the central part of the image and any axially symmetric structure. We also obtained dark images with the cold shutter closed. A dark image corresponding to each sky image was retrieved before or after the pointed observation. As a result, we obtained equal numbers of sky and dark images for each band. For the dark images, we applied the same procedure as for the sky images and stacked them to obtain dark maps for the three wavelength bands.

In order to obtain the fluctuations in sky brightness, we masked pixels whose signal exceeded 2σ level. We repeated this 10 times, after which additional clipping did not show significant changes. However, bright pixels remained around the masked region due to the extended point-spread function (PSF). In order to subtract this effect, we searched for foreground sources with a detection threshold of 2σ in the unmasked stacked images, and removed their contributions using the DAOPHOT package in IRAF. At first, we carefully modeled the PSF based on the beam profiles of the bright point sources in the unmasked stacked images. Using the clearly identified point sources (sources well fitted by the PSF), we estimated the contribution of sources in the outer region using the model PSF and subtracted them from the sky images. We were left with some point sources (sources rejected by the PSF fit) that cannot be fitted with the PSF due to faintness, overlapping of sources, or other reasons. For these sources, we simply masked eight neighboring pixels around their centers. We identified extended objects using ground-based data for the NEP region (Imai et al. 2007), which have a better spatial resolution. We created convolved images using AKARI's PSF and subtracted them from the sky images. For these PSF-subtracted images, we masked the same region in the 2σ clipping image. In order to minimize contributions from PSFs, we further masked a layer of 1 pixel around the masked regions in which the pixel of 2σ brighter level exists. For the dark maps, we masked the same region as with the sky maps and repeated 2σ clipping 10 times. Table 1 shows some parameters used and obtained in this analysis. Percentages of remaining pixels are ∼47% for all wavelength bands, and average sky brightness are 114, 73, and 105 nW m⁻² sr⁻¹. As for the sky brightness, systematic errors are not taken into account in order to show the statistical significance clearly.

Figure 1 shows the final fluctuation maps. The upper and lower panels indicate the sky and dark maps, respectively. The sky maps clearly show structure, while the dark maps look like random noise. As is well known, a major emission component in the near-infrared sky is zodiacal light. We have to keep in mind the fact that the contribution of Galactic stars is negligible, but there remains integrated light of galaxies fainter than the limiting magnitudes, and isotropic emission.

Zoom In Zoom Out Reset image size

Fluctuation analysis was first performed by Kashlinsky et al. (2002) and Odenwald et al. (2003) using Two Micron All Sky Survey (2MASS) data, and followed by Kashlinsky (2005), Kashlinsky et al. (2007a), and Arendt et al. (2010) for Spitzer data and by Thompson et al. (2007a, 2007b) for HST data. In order to examine the angular dependence of the fluctuation observed with AKARI, we performed a power spectrum analysis for both the sky and the dark maps (Figure 1), following the works of previous authors. A two-dimensional Fourier transformation was computed for the fluctuation field δF(x) as a function of the two-dimensional angular wavenumber, q, and the two-dimensional coordinate vector of the image, x:

$$F({\bf q}) = \int \delta F({\bf x}) \exp ^{(-i{\bf x} \cdot {\bf q})} d^{2}x.$$

The two-dimensional power spectrum and fluctuation spectrum are given by P₂(q) ≡ |F(q)|², and [q²P₂(q)/(2π)]^1/2. Figure 2 shows the amplitude maps of the two-dimensional fluctuation spectra in Fourier space. The maps are symmetric but faint annular patterns appear in 2.4 and 3.2 μm maps. These patterns correspond to small peaks in the fluctuation spectra (Figure 3) at ∼5'' and are probably artifacts caused by the masking pattern, since the same features appear in both sky and dark images. However, they do not affect the conclusion of this paper. The one-dimensional fluctuation spectrum, [q²P₂/(2π)]^1/2, is obtained by taking ensemble averages in the upper halves of the concentric rings of Figure 2 with equal spacing in the angular wavenumber, q.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 3 shows the fluctuation spectra. The upper panel shows the fluctuation spectra for both the sky (filled circles) and dark (open triangles) maps. Statistical errors, given by standard deviations divided by the square root of the number of the data points, are plotted. The fluctuation spectra of the sky maps are significantly larger than those of the dark maps.

We subtracted the fluctuation spectra of the dark maps from that of the sky maps in quadrature, and show the result as filled circles in the lower panel of Figure 3. The solid lines indicate the fluctuation spectra due to the shot noise of faint galaxies that were estimated by simulation. We first obtained galaxy counts for the stacked images before masking and extended them to the faint end assuming a slope of 0.23 for logarithmic counts (Maihara et al. 2001). We constructed fluctuation maps by randomly distributing the galaxies fainter than the limiting magnitudes, and performed the power spectrum analysis. The fluctuation spectra of the shot noise of faint galaxies fit well the observed power spectra with limiting magnitudes (AB) of 22.9, 23.2, and 23.8. The decline at small angle is caused by the PSF, while the absolute value depends on limiting magnitude and galaxy counts. Figure 3 shows that a significant residual fluctuation above the shot noise level remains at large angular scales. Table 2 shows the numerical values of fluctuation of the observed sky, fluctuation due to shot noise by faint galaxies, and excess fluctuation at angles larger than 100''.

We performed a subset analysis to confirm the celestial origin of the observed structure. We divided our set of images into two subsets corresponding to two cases. In case A, each subset is taken as a sequence alternating in time, while in case B the subsets correspond to the earlier and later halves of our set. For both cases, we obtained two stacked images, F1 and F2, by applying the same procedure as the one described previously in this section. The fluctuation spectra of the difference between these two stacked images are shown in Figure 4. Cases A and B are shown as squares and asterisks, respectively, while the fluctuation spectra of the dark maps are shown by triangles. The results for both subsets are consistent with those of the stacked dark maps. This indicates that the observed structure is indeed present in the original images and is of celestial origin.

Zoom In Zoom Out Reset image size

We also examined the impact of masking on the fluctuation spectra, since ∼ 53% of all pixels are masked. We constructed a common mask that includes all pixels masked in any of three wavelength bands. The fraction of remaining pixels in images with the common mask applied is ∼ 32%. We again performed the analysis for these images, and compared the obtained fluctuation spectra with the original ones shown in Figure 5. The fluctuation spectra of the images using the common mask have slightly lower values than the original ones, but are consistent with the original ones within the error bars. This result assures that the masks adopted do not cause significant artifacts and that the detected fluctuation is a real structure in the sky.

Zoom In Zoom Out Reset image size

The upper panel of Figure 6 shows smoothed sky maps that were obtained by taking averages within circles of 50'' in diameter centered at each pixel. The smoothed maps clearly show the existence of a significant extended spatial structure with an angular scale of a few hundred arcseconds for all wavelength bands, and the overall pattern is very similar. This is consistent with the excess fluctuation at large angular scales shown in Figure 3. The smoothed dark maps in the lower panel of Figure 6 do not show any structure besides random noise; thus, they reconfirm the celestial origin of the observed large-scale structure.

Zoom In Zoom Out Reset image size

In order to examine the fluctuation quantitatively, we took the averages of the excess fluctuation at angular scales between 100'' and 350''. Figure 7 shows the spectrum of the fluctuating component (filled circles) with 1σ errors. The characteristic feature of the fluctuating component is a blue-star-like spectrum, which is consistent with a Rayleigh–Jeans spectrum, ∼λ⁻³ (straight line in Figure 7). We note that the observed spectrum of the fluctuating component is similar to the spectrum of the isotropic excess emission observed by IRTS (Matsumoto et al. 2005). Open squares show the Spitzer results at 3.6 μm and 4.5 μm (Kashlinsky et al. 2005, 2007a), which are basically consistent with ours. We did not plot the Spitzer results at longer wavelength bands, since the data were scattered and appeared less reliable. At the wavelength bands at 1.1 μm and 1.6 μm (NICMOS/HST), Thompson et al. (2007b) detected fluctuation of 0.3 and 0.4 nW m⁻² sr⁻¹, respectively, at an angular scale of 85''. These values are lower than the values extrapolated from the AKARI data, assuming the Rayleigh–Jeans spectrum, ∼λ⁻³. However, the angular scale for NICMOS/HST is too small to have detected the large-scale structure shown in the smoothed sky maps (Figure 6). Therefore, we consider it is inappropriate to compare the NICMOS/HST results directly with those of AKARI and Spitzer.

Zoom In Zoom Out Reset image size

We performed a two-point correlation (autocorrelation) analysis, C_F(θ) = 〈δF(x + θ)δF(x)〉, both for the sky and dark maps. In order to save computation time, we reduced the angular resolution by adopting the average values groups of 2 × 2 pixels. The results, C_2.4(θ), C_3.2(θ), and C_4.1(θ) are shown in Figure 8 as a function of θ in units of arcseconds. The data were in groups of 4 pixels, corresponding to an angular resolution of 6''. Error bars indicate statistical errors, i.e., the standard deviation divided by the square root of the number of pairs. It must be noted that the sampling at angles larger than 300'' is not complete, that is, central parts of the images are not counted. Figure 8 clearly shows an excess low-frequency components, which are common for all three wavelength bands. The low-frequency components correspond to the excess fluctuations at large angles, as shown in Figure 3.

Zoom In Zoom Out Reset image size

We made a point-to-point correlation study between the wavelength bands using independent data in smoothed maps. The result is shown in Figure 9. The correlation between 3.2 and 2.4 μm is fairly strong (correlation coefficient ∼0.73), while that between 4.1 and 2.4 μm is somewhat weak (correlation coefficient ∼0.42). This is probably because shot noise due to faint galaxies becomes important at 4.1 μm. Colors of the correlating component (slope of Figure 9) normalized to be 1.0 at 2.4 μm (right ordinate) are also plotted in Figure 7 as open circles, which is similar to that of large-scale fluctuation (filled circles).

Zoom In Zoom Out Reset image size

We also obtained the cross-correlation function, C_F⊗G(θ) = 〈δF(x + θ)δG(x)〉 between the maps taken at different wavelengths. C_2.4⊗3.2(θ) and C_2.4⊗4.1(θ) are shown in Figure 10. Degradation of angular resolution and binning of the data are the same as for the two-point correlation analysis. Error bars indicate statistical errors, i.e., standard deviation divided by the square root of the number of pairs. The cross-correlation functions in Figure 10 are very similar to the 2.4 μm two-point correlation, C_2.4(θ). This is consistent with a point-to-point correlation analysis (Figure 9), and the ratios of cross-correlation functions to C_2.4(θ) are similar to the slopes in Figure 9. The result of our cross-correlation analysis implies that the fluctuation patterns are basically the same for all three wavelength bands.

Zoom In Zoom Out Reset image size

Being the main component of sky brightness, zodiacal light is a candidate for the origin of the fluctuation. The subset analysis provides evidence that zodiacal light is not the source of the observed fluctuation. Particularly in case B, the interplanetary dust observed during the first half period is almost independent from the one seen in the latter-half period, since the interplanetary dust in the column along the line of sight has different orbits and eccentricities. If the observed fluctuation were caused by zodiacal light, the fluctuation spectra of the differences between these two subsets should be higher than those of the dark images.

The detected fluctuation levels are also higher than those expected for zodiacal light. IRAS (Vrtilek & Hauser 1995), COBE (Kelsall et al. 1998), and Infrared Space Observatory (ISO; Abraham et al. 1997) attempted to detect the fluctuation of zodiacal emission on scales from arcminutes to degrees. The most stringent upper limit of 0.2% of the sky brightness was obtained by Abraham et al. (1997) who measured the fluctuation of zodiacal emission at 25 μm for a 45' × 45' field with a 3' pixel scale. Recently, Pyo et al. (2011, private communication) analyzed the smoothness of the mid-infrared sky using the same data set as the one we used in this paper. They applied the same power spectrum analysis for the mid-infrared bands, but did not detect any excess fluctuation over the random shot noise caused mainly by photon noise. We, therefore, assume the lowest fluctuation level at the 18 μm band to be an upper limit of fluctuations of the sky at an angle larger than 150'', that is, ∼0.02% of the sky brightness. Although this upper limit includes fluctuation due to zodiacal emission, diffuse Galactic light (DGL), integrated light of galaxies, and cosmic background, we simply assume that the fluctuation of the zodiacal emission is lower than ∼0.02% of the absolute sky brightness.

Fluctuation of zodiacal light is similar to that of zodiacal emission. Interplanetary dust clouds are thought to consist of three components (Kelsall et al. 1998). The first is a smooth cloud, which is broadly distributed in the solar system. The second is dust bands composed of asteroidal collisional debris (Nesvorny et al. 2003). The third is a circumsolar ring which is composed of dust resonantly trapped into orbit near 1 AU. We calculated the contributions of these components to the integrated brightness using the model by Kelsall et al. (1998). Figure 11 shows the volume emissivity of zodiacal light at 2.2 μm (scattering part, solid line) and that of zodiacal emission at 12 μm (thermal part, dashed line) along the line of sight toward the NEP. The volume emissivity is normalized such that the total emissivity at the origin is 1.0. Contributions of the dust bands (bottom), the circumsolar ring (middle), and the total emissivity (top) are shown separately. We used the model at the epoch 2006 September 15 when the contribution by the circumsolar ring reaches its maximum. Figure 11 shows that the smooth cloud is a main contributor for both zodiacal light and emission, while that of the dust bands is very low although they are conspicuous at low ecliptic latitudes. Figure 11 clearly shows that the spatial distribution of the volume emissivity is almost the same for scattered light (zodiacal light) and thermal emission (zodiacal emission). This is because scattering phase function and temperature of the dust in the column do not change much in this direction. Although there is a small deviation for the total emissivity at the distance of 1 AU and beyond, and a systematic difference between scattering and thermal part for the circumsolar ring, their contribution to the total surface brightness is only a few percent. We conclude that the difference of the fluctuation between zodiacal light and zodiacal emission is at most a few percent, even if a fluctuation of the dust density exists along the line of sight.

Zoom In Zoom Out Reset image size

More quantitatively, we estimated fluctuation levels to be expected due to zodiacal light by constraining the upper limit of the fluctuation of zodiacal light to 0.02%. The sky brightnesses corresponding to these upper limits are 0.023, 0.015, and 0.021 nW m⁻² sr⁻¹. The expected fluctuation of zodiacal light could be reduced further during the stacking procedure, since AKARI moved along the Earth's orbit and IRC observed different regions of the solar system during the observation period. The observed fluctuations at large angular scales, 0.25, 0.11, and 0.033 nW m⁻² sr⁻¹ (see Figure 7), are significantly larger than those expected for zodiacal light.

The spectrum of fluctuation observed with AKARI shows a very blue color up to 4.1 μm (Figure 7), while IRTS observations delineated that zodiacal emission becomes prominent at wavelengths longer than 3.5 μm (Ootsubo et al. 1998). This difference also indicates that the fluctuation observed by AKARI is not due to zodiacal light.

Based on these arguments, we conclude that the contribution of zodiacal light to the fluctuation is negligible.

The second possible origin of the fluctuation is DGL, i.e., scattered starlight and thermal emission of the interstellar dust. Since DGL is closely related to the dust column density, we examined its correlation with the far-infrared map. The Far-Infrared Surveyor (FIS) on board AKARI detected a clear cirrus feature toward NEP at 90 μm with ∼30'' angular resolution (Matsuura et al. 2011). Figure 12 shows the correlation between far-infrared and near-infrared maps with the angular resolution of the AKARI data degraded to 30''. No significant correlation was found, indicating that DGL is not the source of the observed near-infrared fluctuation.

Zoom In Zoom Out Reset image size

The third candidate is the clustering of faint galaxies. A few authors (Cooray et al. 2007; Chary et al. 2008) have claimed that the fluctuation observed with Spitzer is due to the clustering of low-redshift galaxies. The blue color observed by AKARI is not consistent with the proposal by Chary et al. (2008), which states that the fluctuating component originates in red dwarf galaxies. Kashlinsky et al. (2007b) carefully examined the correlation between the map of these galaxies and those observed with Spitzer and concluded that red dwarf galaxies are not responsible for the Spitzer result. Sullivan et al. (2007) examined the contribution of the clustering of faint galaxies to the sky fluctuation using the catalog of galaxies. Based on their work, we estimated that the fluctuation at 2.4 μm caused by the clustering of galaxies fainter than the limiting magnitude of this study at multipole, l ∼ 1000 (∼600''), to be ∼0.03 nW m⁻² sr⁻¹. This is significantly lower than the fluctuation detected at 2.4 μm. As for the L band, Figure 8 in Sullivan et al. (2007) shows that the fluctuation due to clustering is fairly weak compared with the one observed by Spitzer (Kashlinsky et al. 2005, 2007a). These results indicate that the clustering of low-redshift faint galaxies is not likely to be the source of the observed fluctuation.

The contribution of faint galaxies to the background is studied by several authors (Totani & Yoshii 2000; Totani et al. 2001; Maihara et al. 2001; Keenan et al. 2010). According to the recent observation by Keenan et al. (2010), the integrated light of galaxies in the K band amounts to 10.0 ± 0.8 nW m⁻² sr⁻¹. Galaxies of ∼18 mag are mostly responsible for the integrated brightness, and the galaxies fainter than that contribute even less. The integrated brightness at 2.4 μm due to galaxies fainter than the limiting magnitude of 22.9 mag is estimated to be at most 0.7 nW m⁻² sr⁻¹ (Keenan et al. 2010). Observed fluctuation at large angular scales, 0.25 nW m⁻² sr⁻¹, is a fraction of this brightness. Sullivan et al. (2007) estimated the fluctuation due to the clustering of galaxies with K magnitudes between 17 and 19 mag. The fluctuation obtained was 0.1 nW m⁻² sr⁻¹, while the integrated light of these galaxies is 4 nW m⁻² sr⁻¹ (Keenan et al. 2010), that is, the fluctuation is 1/40 of the integrated light of galaxies. Furthermore, the fainter galaxies contribute less to the clustering. Compared with this result, the fluctuation observed with AKARI is fairly large. The situation is more conspicuous for the case of NICMOS/HST observations (Thompson et al. 2007b), since their limiting magnitudes, 28.5 mag, are much fainter than the one of the AKARI observation. The integrated light in H band due to galaxies fainter than 28.5 mag is ∼0.04 nW m⁻² sr⁻¹ (Keenan et al. 2010), while the observed fluctuation is 0.4 nW m⁻² sr⁻¹. The observed fluctuation is 10 times brighter than the integrated light of galaxies, which indicates that the origin of the observed fluctuation is not the low-redshift galaxies. These considerations imply that a simple extrapolation of galaxy counts hardly explains fluctuation observed, and that there must exist a new component of infrared sources with magnitudes fainter than the current limit.

Accreting black holes could be another candidate if they are abundant in the early universe, since they show blue spectra. However, the bulk of the light will be in much shorter wavelengths like X-rays, since the accretion disks must be very hot.

To summarize, we found that the observed near-infrared fluctuation at large angular scales can hardly be explained by known sources. We now regard the most probable origin of the fluctuation to be Pop III stars, the first stars of the universe.

AKARI observations have shown that a large-scale structure existed even during the Pop III era. Since the fluctuation power is nearly constant at scales of a few hundred arcseconds, the large-scale structure in the Pop III era could be more extended. The co-moving distance for 10' at $\it z$ ∼10 is ∼30 Mpc, which is a scale similar to the local super cluster in the present universe. Theoretical studies of sky fluctuation due to the first stars have been predicted by several authors (Cooray et al. 2004; Kashlinsky et al. 2004; Fernandez et al. 2010).

Cooray et al. (2004) presented the first theoretical study on fluctuations caused by Pop III stars. They assumed biased star formation which traces the density distribution of the dark matter. Although the absolute level of the fluctuation spectrum has a wide range depending on the physical parameters, the overall shapes have almost the same convex feature with a flat peak at l ∼ 1000 (∼10') and a turnover toward large angles. We regard this flat peak as the most significant excess fluctuation observed by AKARI, which is mainly due to the Pop III stars. The fluctuating power observed with AKARI is marginally consistent with a pessimistic estimate (Figure 3 in Cooray et al. 2004). The wavelength dependence of fluctuation power of the model shows fairly blue color, which is basically consistent with AKARI observation. The fluctuation observed with AKARI could be interpreted within the framework of this model. However, detection of the turnover at an angular scale larger than 10' is essential.

Recently, Fernandez et al. (2010) made a numerical simulation of the cosmic near-infrared background caused by early populations (Pop II and III stars), assuming a halo mass of 2 × 10⁹ M_☉. They calculated both amplitude and fluctuation, taking into account stellar emission and emission from ionized gas surrounding stars. As for the fluctuation spectrum, they predicted a monotonic decrease of the fluctuation toward large angles, which is not consistent with our result. This discrepancy requires some kind of modification of the theory. However, their result provides valuable information on the physical conditions of early populations which are independent from halo mass.

A spectrum of the sky brightness at 2 ∼ 4 μm obtained by Fernandez et al. (2010) shows a blue color, and the main emission component is attributed to stars. Furthermore, the ratio of the fluctuating power to the absolute sky brightness, δI/I, does not depend on the wavelength. This justifies regarding the spectral shape of absolute brightness to be same as that of the fluctuation. In Figure 7, the spectrum of Figure 20 of Fernandez et al. (2010) is plotted as a dotted line. Their spectrum is a little redder than, but is qualitatively consistent with, the AKARI observation. A good correlation for the fluctuations between the wavelength bands observed with AKARI supports the idea that the emission source is stellar emission from early populations, and that the nearest stars provide the largest contribution to the sky brightness. δI/I estimated by Fernandez et al. (2010) is ∼0.02 at a few hundred arcsecond for the largest fluctuation. The observed fluctuation at 2.4 μm, ∼0.25 nW m⁻² sr⁻¹, renders the absolute sky brightness to be brighter than ∼13 nW m⁻² sr⁻¹, which is consistent with the excess brightness observed with IRTS and COBE (Matsumoto et al. 2005).

All theoretical models (Santos et al. 2002; Salvaterra & Ferrara 2003; Dwek et al. 2005a; Cooray et al. 2004; Fernandez et al. 2010) predict a clear peak due to the redshifted Lyα from the ionized gas surrounding first stars. The wavelength of this peak depends on the redshift of the nearest first stars. Since there is no signature of redshifted Lyα in AKARI data, it must appear below 2 μm, that is, the redshift of the nearest first stars is less than 15. Detection of the redshifted Lyα both in the amplitude and the fluctuation will be a key issue for delineating the epoch of formation of the first stars.

The fluctuation observed with AKARI provides evidence for light from the first stars. However, more observational results are required in order to investigate the physical conditions and the star formation history of the first stars. In a forthcoming paper, we will present a fluctuation analysis at larger angles than this paper based on the AKARI NEP wide field survey. The sounding rocket experiment, Cosmic Infrared Background Experiment (CIBER), will provide the spectrum of the sky brightness around 1 μm and the fluctuation at the 1° scale at the 0.8 μm and 1.6 μm band (Bock et al. 2006). The Korean infrared satellite, the Multi Purpose Infrared Imaging System (MIRIS; Han et al. 2010), will observe the fluctuations at an even larger angular scale (∼10°). These space observations will certainly improve our understanding of the first stars of the universe.

We performed a fluctuation analysis for a sky region with a diameter of 10' toward the NEP at 2.4, 3.2, and 4.1 μm based on observations with AKARI. A power spectrum analysis indicates that a significant excess fluctuation exists at an angular scale larger than 100'' that cannot be explained by zodiacal light, DGL, shot noise from faint galaxies, or clustering of low-redshift galaxies. The fluctuating component at large angles has a blue stellar spectrum, and shows a significant spatial correlation between adjacent wavelength bands.

The fluctuation observed with AKARI and NICMOS/HST (Thompson et al. 2007b) is considerably higher than the one expected from galaxy counts. This implies that a simple extrapolation of galaxy counts hardly explains the observed fluctuation, and that there must be a new component of infrared sources at the faint end.

The most probable origin of excess fluctuation is Pop III stars, the first stars of the universe. The detected fluctuation can be explained with biased star formation which traces dark matter (Cooray et al. 2004). A blue color of the observed fluctuation and a significant spatial correlation between wavelength bands imply that the source of emission is stellar emission from early populations (Fernandez et al. 2010).

The angular scale of 10' at $\it z$ ∼10 corresponds to a co-moving distance of 30 Mpc, a scale similar to the local super cluster in the present universe.

The detection of a turnover of fluctuation spectra at large angles and the spectral peak of redshifted Lyα at wavelengths shorter than 2 μm is a key issue for future observations in order to improve the investigation of the Pop III star formation.

We thank the AKARI/IRC team for their encouragement and helpful discussions. Thanks also to Professors A. Cooray, E. Komatsu, and A. Ferrara for valuable discussions and comments. T.M., S.M., and T.W. were supported by Grants-in-Aid from the Japan Society for the Promotion of Science (18204018 and 21111004). H.M.L. was supported by a grant from National Research Foundation of Korea (2006-341-C00018). W.-S.J. was supported by the STSAT-3 MIRIS project from the Ministry of Education, Science and Technology (MEST) of Korea.