text:
# ISO deep far-infrared survey in the “Lockman Hole”

## 1 Introduction

In order to understand the history of galaxy formation and evolution, surveys at far-infrared(FIR) and submillimeter wavelengths are essentially important, since a large portion of star formation activity in the universe may be hidden by dust, prohibiting optical and near-infrared studies due to the enormous extinction. The recent detection of the far-infrared Cosmic Infrared Background(CIB) radiation by COBE, which appears to have comparable brightness to the total intensity in deep optical counts from the Hubble deep field (Fixen et al. Fix98 (1998); Hauser et al. Hau98 (1998); Lagache et al. Lag00b (2000)), indicates that the infrared-bright galaxies are responsible for roughly half of the energy released by nucleosynthesis. Various discrete source surveys are also being pursued from mid-infrared to submillimeter wavelengths. For example, source counts at 15 $`\mu `$m (Aussel et al. Aus98 (1998); Altieri et al. Alt98 (1998); Elbaz et al. Elb99 (1999)), at 170 $`\mu `$m (Kawara et al. Kaw98 (1998)(paper I) ; Puget et al. Pug99 (1999)), and at 850 $`\mu `$m (Blain et al. Bla99 (1999)) have been reported. Using both CIB and these new number count data, modeling of the cosmic star formation history has been attempted by many authors(Guiderdoni et al. Gui98 (1998); Dwek et al. Dwe98 (1998); Rowan-Robinson RR99 (1999); Tan et al. Tan99 (1999); Ishii et al. IT99 (1999)). All of them require a strong evolution in the star formation rate as we look back to high redshift(more than 10 times larger at $`z1`$).

Although the FIR deep survey is now considered a key observing method for the exploration of the “optically dark side” of the star formation history of galaxies(Guiderdoni et al. Gui97 (1997)), the detectivity of 1m-class space FIR telescopes is likely to be limited by the noise due to the fluctuation of the IR cirrus(Low et al. Low84 (1984)), emission from the interstellar dust, even at high Galactic latitude(Helou & Beichman HeBe90 (1990)). Hence a detailed study of the IR cirrus fluctuations is highly important, especially for the planning of the deep surveys intended with forthcoming IR space telescopes such as IRIS(Murakami Mur98 (1998)) and SIRTF(Fanson et al. Fan98 (1998)).

The spatial structure of the IR cirrus at 100 $`\mu `$m as measured by IRAS was extensively studied by Gautier et al. (Gau92 (1992)) and Abergel et al.(Abe96 (1996)). From a Fourier analysis of the brightness distribution, Gautier et al. found that the Power Spectral Density(PSD) of the brightness fluctuation at 100 $`\mu `$m follows a power-law function of the spatial frequency with an index of about 3 below the spatial frequency corresponding to the IRAS beam size(about $`0.25\mathrm{arcmin}^1`$). They also found that the PSD is proportional to $`B_0^3`$, where $`B_0`$ is the mean brightness of the IR cirrus. The ISOPHOT(Lemke et al. Lem96 (1996)) onboard the ISO(Kessler et al. Kes96 (1996)) is capable of observing the IR cirrus at wavelengths longer than 100 $`\mu `$m, with a better spatial resolution than that of IRAS. Herbstmeier et al. (Her98 (1998)) analyzed the spatial characteristics of four fields measured by ISOPHOT, and obtained similar power-law spectra for relatively bright cirrus regions. Recently Lagache & Puget (Lag00a (2000)) have made a power spectrum analysis of the 170 $`\mu `$m image of the Marano 1 field, and found a significant excess at $`f=0.250.6\mathrm{arcmin}^1`$, which they attributed to the fluctuations due to unresolved extra-galactic sources.

In this paper we investigate the characteristics of FIR brightness fluctuations in the Lockman Hole, a region with a uniquely low HI column density(Lockman et al. Loc86 (1986)), and thus minimal IR cirrus contribution to the fluctuations. Therefore the fluctuations in the Lockman Hole are likely to be dominated by faint, distant galaxies(Helou & Beichman HeBe90 (1990); Herbstmeier et al. Her98 (1998)) and the fluctuation analysis of FIR images will provide unique information on the number counts of infrared galaxies even below the source-confusion limit, thus constraining the parameters characterizing the number count models.

This paper is organized as follows: section 2 briefly describes the observations and data processing. Section 3 explains the power spectrum analysis and examines the contribution by the IR cirrus fluctuations and also describes a simulation of the images and the PSDs. Section 4 describes the constraints on the galaxy number counts. The nature of the sources responsible for the fluctuations is discussed in section 5, and section 6 gives the conclusions.

## 2 Observation and data processing

The FIR survey of the Lockman Hole, which was executed as a part of Japan/UH cosmology program using the ISAS guaranteed time, is described in Taniguchi et al. (Tan94 (1994)) and paper I. Two $`44\mathrm{}\times 44\mathrm{}`$ fields named LHEX and LHNW were mapped with two filters: C\_90(centered at 90 $`\mu `$m) and C\_160(170 $`\mu `$m). Each of the two fields is made up of 4 sub-fields.

Each sub-field map was produced from the edited raw data by the PHT Interactive Analysis(PIA; Gabriel et al. Gab97 (1997)) version 7.1 or 7.2, and is hereafter referred to as an AAP(Astronomical Analysis Processing of PIA) map. Each AAP map is either $`58\times 58`$ pixels($`23\mathrm{}/\mathrm{pixel}`$) for a 90 $`\mu `$m sub-field, or $`31\times 31`$ pixels($`46\mathrm{}/\mathrm{pixel}`$) for a 170 $`\mu `$m sub-field. To correct the drift in the responsivity of the detectors, we applied the median filter smoothing (see Paper I for details) to the AAP data.

Together with the final AAP maps we also produced the uncertainty maps, and found that the typical 3$`\sigma `$ noise is as low as 0.012MJy/sr( 0.60mJy/pixel for a $`46\mathrm{}\times 46\mathrm{}`$ pixel), indicating that the instrumental noise is negligible in the following results.

As was done in Paper I, the observed fluxes as well as the brightness of the images are scaled based on the fluxes of the brightest source F10507+5723(UGC 06009) measured with IRAS. We found that the flux calibration based on the FCS1 measurements underestimates the 90 $`\mu `$m flux of the IRAS source by a factor of 2.6, although it gives the mean brightness of the images consistent with the COBE/DIRBE brightness within 25 per cent. The 170 $`\mu `$m flux of the IRAS source is assumed to be 1133 mJy as described in Paper I. The flux calibration based on the FCS1 measurements again underestimates the 170 $`\mu `$m flux by a factor of 1.5. For the 170 $`\mu `$m flux, we found that the discrepancy is mostly due to an underestimate of the effective solid angle of the ISOPHOT detector. Puget et al. (Pug99 (1999)) and Lagache & Puget (Lag00a (2000)) used an effective solid angle at 170 $`\mu `$m derived from Saturn footprint measurements, which is significantly larger than that used in the PIA. Here we assume theoretical PSFs of a telescope with a 60 cm primary mirror and a 20 cm secondary mirror which agree well with those measured during the ISOPHOT calibration observations(Müller Mul99 (2000)). With this assumption, the factors of 2.6(90 $`\mu `$m) and 1.5(170 $`\mu `$m) discrepancies reduce to factors of 2.0(90 $`\mu `$m) and 1.2(170 $`\mu `$m). The origin of the residual discrepancies is still unknown. However, this is not problematic for the main results and conclusions of this paper as described in section 4 and 5 as long as the fluxes of F10507+5723 determined by the IRAS faint source survey are correct.

Fig. 1 shows the images of LHEX and LHNW used for the fluctuation analysis, each of which is the largest square area extracted from the mosaiced map(see Fig. 3 of Paper I) made up from four sub-field AAP maps. Each image is rebinned into $`103\times 103`$ pixels2 (C\_90), $`53\times 53`$ pixels2(C\_160 LHEX), or $`51\times 51`$ pixels2(C\_160 LHNW). The plate scale is $`23\mathrm{}/\mathrm{pixel}`$ for 90 $`\mu `$m maps, and $`46\mathrm{}/\mathrm{pixel}`$ for 170 $`\mu `$m maps.

## 3 The power spectra and the simulation

For each image shown in Fig. 1 the 2-dimensional angular correlation function $`C(x,y)`$ was calculated from the brightness distribution $`B(x_0,y_0)`$:

$$C(x,y)=<(B(x_0,y_0)\overline{B})(B(x_0+x,y_0+y)\overline{B})>$$
(1)
where bracket $`<>`$ represents the average over the whole area in each image shown in Fig. 1, and $`\overline{B}=<B>`$. Then, the coordinates were expressed in polar coordinates as $`C(r,\theta )`$, and the Fourier transforms in various radial directions were calculated (i.e. $`\theta `$ is fixed for each transform). One-dimensional power spectral density(PSD) $`P(f)`$ is calculated from

$$P(f)=\frac{2}{N}\underset{k=0}{\overset{N1}{}}C_kcos(2\pi fr_k),$$
(2)
where $`f`$ is spatial frequency($`1/(Nr_{\mathrm{pix}})f1/(2r_{\mathrm{pix}})`$), $`r_kr_{\mathrm{pix}}k`$ ($`k=0,1,\mathrm{},N1`$), $`r_{\mathrm{pix}}`$ is the pixel size of the map, $`Nr_{\mathrm{pix}}`$ is the largest angular size for which the angular correlation function is evaluated, and $`C_kC(r=r_k,\theta )`$ . Note that the units of the PSDs are the same as that of the angular correlation function: $`\mathrm{Jy}^2/\mathrm{sr}^2`$. Finally, the PSDs were averaged with respect to $`\theta `$ and are shown in Fig. 2. It is noteworthy that the fluctuations at high frequencies($`f0.8`$ arcmin-1 for 90 $`\mu `$m, and $`f0.4`$ arcmin-1 for 170 $`\mu `$m) are smoothed by the instrumental beam and therefore the PSDs decrease appreciably. The small error bars of PSDs represent the standard deviation among a set of the PSDs with different $`\theta `$, showing that the PSDs are almost independent of $`\theta `$.

In order to check the contributions from bright sources, the PSDs are derived by masking circular regions with a 4$`\times `$FWHM diameter around bright sources with fluxes above $`S_{\mathrm{max}}`$: $`S_{\mathrm{max}}=250`$ mJy at 170 $`\mu `$m111 As for 170 $`\mu `$m images, smaller values of $`S_{\mathrm{max}}`$ significantly reduce the amount of residual images., $`S_{\mathrm{max}}=150`$ mJy at 90 $`\mu `$m. In the following the resultant PSDs are called “residual PSDs”, and are also shown in Fig. 2. Interestingly, the residual PSD for each image is more than half of the PSD of corresponding original image with almost the same spectral shape, indicating that there remains significant contribution from randomly distributed point sources with fluxes below $`S_{\mathrm{max}}`$.

In Fig. 2 typical IR cirrus PSDs are also compared in order to check the contribution of the IR cirrus to the PSDs of the Lockman Hole222The formula for the IR cirrus PSD by Gautier et al. (Gau92 (1992)) cannot be directly compared with the present work, because the power spectrum analysis presented here is one-dimensional one while the formula of Gautier et al. is based on the two-dimensional analysis.. We examined several IRAS 100 $`\mu `$m maps of high-latitude clouds in Ursa Major($`l=145^{}`$, $`b=40^{}`$), which are reproduced from the IRAS Sky Survey Atlas(ISSA) by reducing the brightness by a factor of 0.72, following the COBE/DIRBE calibration(Wheelock et al. Whe94 (1994)). The average brightness of the cirrus is 2–3 MJy/sr, and each map is $`150\mathrm{}\times 150\mathrm{}`$ wide with $`1\mathrm{}.5/\mathrm{pixel}`$. The derived PSDs show a power-law distribution with an index of -1.5. Gautier et al.(Gau92 (1992)) noted that one-dimensional analysis of the IR cirrus yielded spectral indices near -2. These cirrus fluctuation spectra are much different from those obtained for the Lockman Hole images at 90 $`\mu `$m, which show rather flat spectra at lower spatial frequencies. The 170 $`\mu `$m spectra present a slope similar to the IR cirrus one, but this can be explained by the shape of the footprint power spectrum of ISOPHOT detectors. Moreover, the fluctuation powers are much larger than those estimated for the IR cirrus, which are also shown in Fig. 2. We assume that the cirrus PSD is proportional to $`B_0^3`$(Gautier et al. Gau92 (1992)), and taking the mean brightness $`B_0`$ of the IR cirrus in the Lockman Hole as 0.33 MJy/sr at 90 $`\mu `$m and 1.0 MJy/sr at 170 $`\mu `$m. These values are estimated from the atomic hydrogen column density of $`6\times 10^{19}\mathrm{atoms}\mathrm{cm}^2`$ in LHEX(Jahoda et al. Jah90 (1990)) and the COBE/DIRBE data analysis in the Lockman Hole by Lagache et al. (Lag99a (1999), in their Table 4). We can conclude that the IR cirrus contribution to the PSDs in the Lockman Hole is negligible over all spatial frequencies $`f0.05`$.

In the following we interpret the residual PSDs in terms of unresolved point sources, probably galaxies. We neglect the spatial correlation between galaxies, thus assuming that galaxies are randomly distributed in the images. Then the residual PSD will be the product of the footprint power spectrum of the ISOPHOT detector and the fluctuation power due to the point sources, which is a white power spectrum given by:

$$P_{\mathrm{source}}=_{S_{\mathrm{min}}}^{S_{\mathrm{max}}}S^2\frac{dN}{dS}𝑑S$$
(3)
where $`dN/dS`$ is differential source counts. From the residual PSDs observed, we derive $`P_{\mathrm{source}}=13000\pm 3000\mathrm{Jy}^2/\mathrm{sr}`$ at 90 $`\mu `$m ($`S_{\mathrm{max}}=150\mathrm{mJy}`$) and $`P_{\mathrm{source}}=12000\pm 2000\mathrm{Jy}^2/\mathrm{sr}`$ at 170 $`\mu `$m ($`S_{\mathrm{max}}=250\mathrm{mJy}`$), where the errors do not include systematic ones due to uncertainties in the flux calibration. Lagache & Puget (Lag00a (2000)) reported the detection of $`P_{\mathrm{source}}=7400\mathrm{Jy}^2/\mathrm{sr}`$ at 170 $`\mu `$m for the Marano 1 field, after sources brighter than 100 mJy are removed. Contribution from detected sources with fluxes between 100 mJy and 250 mJy is estimated to be about 7000 $`\mathrm{Jy}^2/\mathrm{sr}`$. Thus the fluctuation power in the Lockman Hole due to the sources fainter than 100 mJy is $`5000\mathrm{Jy}^2/\mathrm{sr}`$, which is comparable to that observed in the Marano 1 field.

A simulation was performed by making 90 $`\mu `$m and 170 $`\mu `$m images made up only by galaxies with fluxes between $`S_{\mathrm{min}}`$ and $`S_{\mathrm{max}}`$ and calculating their PSDs. Here galaxies are treated as point sources with a PSF specific to the respective wavelength band of ISOPHOT. We used the image of the bright IRAS source (F10507+5723) seen in LHEX images as the PSF. The number of sources and their flux densities are controlled by the source counts. We examined the non-evolution count model by Takeuchi et al. (Tak99 (1999)) and the scenario E by Guiderdoni et al.(Gui98 (1998)). We assume $`S_{\mathrm{min}}=10\mathrm{mJy}`$ because the fluctuations due to galaxies fainter than 10 mJy are negligible in case of these models. The resulted PSDs are compared with the residual PSDs in Fig. 3. These simulated PSDs are not sufficient to explain the observed PSDs although the spectral shapes are quite similar.

## 4 Constraints on the Number Counts

In this section we will investigate the source counts below $`S_{\mathrm{max}}`$, by fitting the simulated PSDs to the residual PSDs. In Paper I we presented the number density of sources brighter than 150 mJy at both filter bands. After Paper I was published, the data processing and the source extraction technique have been improved, and now we obtained source number counts down to $`S=70`$ mJy at 90 $`\mu `$m and $`S=100`$ mJy at 170 $`\mu `$m, although the counts at such low fluxes may be overestimated due to the source confusion. We found that the cumulative counts above 150 mJy show quite a steep increase as the source flux decreases: $`\beta 3`$ where $`N(>S)S^\beta `$. These results will be described in detail in Kawara et al. (Kaw00 (2000), paper III). If $`\beta `$ remains equal to -3 down to $`S_{\mathrm{min}}=1\mathrm{mJy}`$, then the fluctuation power($`P_{\mathrm{source}}`$) calculated by equation 3 exceeds the observed ones by an order of magnitude. Hence at a certain flux below $`S_{\mathrm{max}}`$, we expect that the slope of the counts must flatten to $`\beta <2`$ so that the predicted fluctuations do not exceed the observed ones. We thus considered simple double power-law count models and evaluated the simulated PSDs: in the flux range $`S_{\mathrm{max}}SS_{\mathrm{crit}}`$

$$N(>S)=N_{\mathrm{max}}\left(\frac{S}{S_{\mathrm{max}}}\right)^{\beta _0},$$
(4)
and in flux range $`S_{\mathrm{crit}}SS_{\mathrm{min}}`$

$$N(>S)=N_{\mathrm{crit}}\left(\frac{S}{S_{\mathrm{crit}}}\right)^{\beta _1},$$
(5)
where $`N_{\mathrm{crit}}=N_{\mathrm{max}}(S_{\mathrm{crit}}/S_{\mathrm{max}})^{\beta _0}`$, and $`\beta _0`$, $`\beta _1`$ are assumed to obey inequalities $`0\beta _1<\beta _04`$. $`N_{\mathrm{max}}`$ and its uncertainty are listed in Table 1. The uncertainty in $`N_{\mathrm{max}}`$ includes Poisson uncertainties based on the total number of sources with flux above $`S_{\mathrm{max}}`$ in all images, and systematic ones originating from the incompleteness due to the source confusion, which will be discussed in paper III. If $`\beta _1<2<\beta _0`$, the simulated PSD is not sensitive to $`S_{\mathrm{min}}`$ and is dominated by sources with fluxes around $`S_{\mathrm{crit}}`$, as discussed in Lagache and Puget (Lag00b (2000)). Various sets of parameters can be chosen so that the simulated PSDs fit to the residual PSDs. Examples of the simulated PSDs which fit well to the residual PSDs are also shown in Fig. 3, for which $`\beta _0=3.0`$, $`\beta _1=1.0`$ for both 90 $`\mu `$m and 170 $`\mu `$m counts, and $`S_{\mathrm{crit}}=80`$ mJy for 90 $`\mu `$m counts and $`S_{\mathrm{crit}}=150`$ mJy for 170 $`\mu `$m one. We found that the simulated PSDs can fit the residual PSDs over all spatial frequencies within $`\pm 22`$ per cent for the 90 $`\mu `$m PSDs and $`\pm 17`$ per cent for the 170 $`\mu `$m ones. These deviations of the residual PSDs from the simulated PSDs are also included as uncertainties in the residual PSDs.

We now derive the allowed regions in plots of the cumulative counts at fluxes below $`S_{\mathrm{max}}`$. The bottom border of the allowed regions is derived by examining the single power-law counts($`S_{\mathrm{crit}}=S_{\mathrm{min}}=1\mathrm{mJy}`$)333Lower value of $`S_{\mathrm{min}}`$ does not appreciably change the results: $`\beta _0`$=1.7 even if $`S_{\mathrm{min}}=0`$. and we obtain $`\beta _0`$=1.8 for 90 $`\mu `$m counts and $`\beta _0`$=1.7 for 170 $`\mu `$m ones. In case of double power-law counts we may choose a larger value of $`\beta _0`$ with a set of parameter pairs($`S_{\mathrm{crit}}`$, $`\beta _1`$), and among them $`S_{\mathrm{crit}}`$ for $`\beta _1=0`$ is the largest. Thus we determine the upper borders of the allowed regions by connecting ($`S_{\mathrm{crit}}`$, $`N_{\mathrm{crit}}`$) points of numerous models with $`\beta _1=0`$ and $`\beta _0=24`$. The results are shown by the shaded area in Fig. 4. As for the 170 $`\mu `$m counts, the bottom border of the allowed region is close to the Scenario E model by Guiderdoni et al. (Gui98 (1998)), and is consistent with the counts at 120 mJy and at 200 mJy obtained from the FIRBACK Marano 1 survey (Puget et al. Pug99 (1999)). On the other hand, the allowed region for 90 $`\mu `$m counts are significantly above those of any currently existing models.

The bottom border derived here cannot be applied below a certain flux $`S_\mathrm{L}`$, which is the flux at the intersection between an upper-border count model with the steepest slope ($`\beta _0=4`$, $`\beta _1=0`$), and the bottom border: $`S_\mathrm{L}=30mJy`$(90 $`\mu `$m) or 50 mJy(175 $`\mu `$m). By definition, the upper border model do not require sources below $`S_{\mathrm{crit}}`$($`>S_\mathrm{L}`$). On the other hand, in case of the bottom border model, contribution to the fluctuations from sources below $`S_\mathrm{L}`$ is still appreciable, about half of the total. Considering the uncertainties in $`N_{\mathrm{max}}`$ and the residual PSDs as discussed above, we finally give $`S_\mathrm{L}=35mJy`$ for 90 $`\mu `$m counts and $`S_\mathrm{L}=60mJy`$ for 170 $`\mu `$m counts.

## 5 Discussion

We attribute the origin of the residual PSDs to the random distributions of faint sources which can no longer be identified as individual sources due to the source confusion. Another piece of evidence which supports this interpretation is the 170 $`\mu `$m/90 $`\mu `$m brightness ratio of the FIR background emission. In Fig. 5 we show the point-to-point comparison between the 90 $`\mu `$m and the 170 $`\mu `$m sky brightness. Here we examine the residual images after masking pixels around the sources with $`S>S_{\mathrm{max}}`$. In order to match the spatial resolution between 90 $`\mu `$m($`76\mathrm{}`$ FWHM) and 170 $`\mu `$m($`144\mathrm{}`$ FWHM) images, one wavelength band image is convoluted with the other wavelength band beam profile. Although the scatter is large, the plot shows linear correlation with a slope of unity. This background color is quite different from that expected for the IR cirrus(about 3.1), and can be interpreted as a typical FIR color of galaxies contributing to the fluctuation. The fluctuation color due to faint galaxies depends on their redshift, and their SEDs. The relations between the FIR flux ratio and the redshift for SEDs of the cirrus dominated galaxy, the pure starburst galaxy(Efstathiou et al. Efs00 (2000); Efstathiou & Siebenmorgen ES00 (2000)), and the mixture of these two SEDs which represents a star-forming galaxy’s SED like IRAS F10507+5723, are shown in Fig. 6. The flux ratio of $`S_\nu (170\mu \mathrm{m})/S_\nu (90\mu \mathrm{m})=1`$ is attained at $`z0.7`$ for a pure starburst galaxy, while a cirrus dominated galaxy or a star-forming galaxy with a small contribution from starburst component must be local($`z0.1`$).

From Fig. 6 we can also justify the flux calibration based on the IRAS fluxes of F10507+5723. If we adopt the ISO 90 $`\mu `$m flux determined by PIA and the theoretical PSFs, then $`S_\nu (90\mu \mathrm{m})/S_\nu (60\mu \mathrm{m})=1.0`$ which is the color of a pure-starburst galaxy, while the ISO flux ratio $`S_\nu (170\mu \mathrm{m})/S_\nu (90\mu \mathrm{m})=1.6`$ is much larger than that expected for any nearby($`z0`$) galaxies.

The local cirrus dominated galaxies are not likely to be the dominant sources responsible for the fluctuations, since such galaxies are relatively bright at optical wavelength and hence the source number density in the Lockman Hole at optical wavelength would exceed the observed one by more than an order of magnitude: as shown in Table 1 we obtained source density of about $`2\times 10^5`$ sources/sr($`60\mathrm{sources}/\mathrm{deg}^2`$) at $`S100`$ mJy. Using the averaged SED of normal galaxies given by Schmidt et al. (Sch97 (1997)), the B magnitude of a 100 mJy source at 90 $`\mu `$m is estimated to be about 16 mag. Reported B-band source counts above B = 16 mag is only a few sources per square degrees(Kirschner et al. KOS79 (1979); Ellis Ell83 (1983)). Hence the galaxies responsible for the fluctuations must be heavily extincted at optical wavelengths, which is a well-known feature of starburst galaxies. In conclusion, the major source of the brightness fluctuations is most probably the star-forming galaxies located at $`z<1`$ extincted at optical wavelengths.

In the plot of the 90 $`\mu `$m counts in Fig. 4, the IRAS counts are also plotted. These were originally the 60 $`\mu `$m counts, in which the flux is scaled by the factor $`S_\nu `$(90 $`\mu `$m)/$`S_\nu `$(60 $`\mu `$m)=2.1. This scaling factor is derived by the IRAS source counts at relatively high fluxes($`S_\nu (100\mu \mathrm{m})>1\mathrm{Jy}`$) (Rowan-Robinson et al. Row86 (1986)), and is consistent with the color of normal galaxies with a partial contribution from starburst. Thus, if such galaxies dominate the 60 $`\mu `$m counts down to 110 mJy, this scaling is valid. We used the complete differential source counts down to 110mJy by Lonsdale et al. (Lon90 (1990)) and Bertin et al.(Ber97 (1997)), and as shown in Fig. 4 the scaled counts are in good agreement with the model counts given by Guiderdoni et al.(Gui98 (1998)). Since the integrated counts down to 150 mJy is much larger than the model counts, the integrated counts must show a steep rise between 150 mJy and 240 mJy, which, however, could not be observed in the IRAS 60 $`\mu `$m counts due to the source confusion(Hacking and Houck Hac87 (1987)). Further studies from space like IRIS and SIRTF will judge if this steep rise is real.

Finally, we mention on the impact of the present work to the cosmic infrared background(CIB). Hauser et al.(Hau98 (1998)) reported an upper limit of 1.1 MJy/sr at 100 $`\mu `$m and a positive detection of the CIB of 1.2 MJy/sr at 140 $`\mu `$m and 1.1 MJy/sr at 240 $`\mu `$m. Lagache et al. (Lag99a (1999)) reported lower values of the CIB: $`0.72\pm 0.30`$ MJy/sr at 140 $`\mu `$m and $`0.91\pm 0.15`$ MJy/sr at 240 $`\mu `$m by further subtracting the dust emission associated with the diffuse ionized gas. Lagache et al.(Lag00b (2000)) also reported the CIB at 100 $`\mu `$m: $`0.78\pm 0.21`$ MJy/sr. By summing the fluxes of all detected sources above 150 mJy, we obtained an integrated brightness of the CIB of 0.031 MJy/sr at 90 $`\mu `$m and 0.050 MJy/sr at 170 $`\mu `$m. If we consider the constraints on the number counts for $`SS_\mathrm{L}`$ shown in Fig. 4, then the integrated brightness would be 0.09 – 0.30 MJy/sr at 90 $`\mu `$m for sources above 35 mJy and 0.053 – 0.15 MJy/sr at 170 $`\mu `$m for source above 60 mJy. Hence, 5 – 40 per cent of CIB can now be attributed to the integrated light of discrete sources above $`S_\mathrm{L}`$ which are responsible for the fluctuations.

## 6 Conclusion

We have studied the characteristics of FIR brightness fluctuations at 90 $`\mu `$m and 170 $`\mu `$m in the Lockman Hole, using the power spectrum analysis over the spatial frequency range of $`f=0.051`$ arcmin-1. The spectra of PSDs are found to be rather flat at low frequencies, which differs from that expected for IR cirrus fluctuations. We interpret the spectral characteristics of PSDs as those of randomly distributed point sources. The fluctuations in the Lockman Hole are not dominated by the IR cirrus, and are instead most likely due to star-forming galaxies. We next showed the constraints on the galaxy number counts down to 35 mJy at 90 $`\mu `$m and 60 mJy at 170 $`\mu `$m, assuming that only galaxies contribute to the observed fluctuation powers. This analysis indicates the existence of steep rise in the integrated counts down to these fluxes and especially at 90 $`\mu `$m, the source density is much larger than that expected from the currently available number count models.

We also found a linear correlation between sky brightnesses at 90 $`\mu `$m and 170 $`\mu `$m, even if pixels around the bright sources are masked from the images. The slope of the linear fit is much lower than the color of the IR cirrus, indicating that the point sources are star-forming galaxies at low or moderate redshift($`z<1`$). These galaxies with fluxes above 35 mJy or 60 mJy also significantly contribute to the CIB recently reported by Hauser et al.(Hau98 (1998)) and Lagache et al.(Lag00b (2000)).

###### Acknowledgements.

The authors would like to thank Hiroshi Shibai, Tsutomu Takeuchi, Hiroyuki Hirashita, and Chris P. Pearson for their extremely useful comments. HM thanks especially Andreas Efstathiou, who kindly provided us with the cirrus galaxy SED before its publication.


format: text

text:
# USA and RXTE Observations of a Variable Low-Frequency QPO in XTE J1118+480

## 1. Background on XTE J1118+480

The transient source XTE J1118+480 (Remillard et al. (2000)) in Ursa Major appeared as an X-ray source in the first few days of 2000 and has continued to be detectable through the subsequent six months. Its X-ray activity to date has consisted of two distinct outbursts, the first of which rose to maximum rapidly and declined in 2000 January. The rise to the second maximum followed about 25 days after the decline from the first and it was during this second maximum that the All Sky Monitor (ASM) on the Rossi X-ray Timing Explorer (RXTE) discovered the source. The discovery announcement noted the earlier maximum based on pre-discovery observations.

XTE J1118+480 is tentatively classified as a black hole transient based on the detection of low-frequency quasi-periodic oscillations (QPOs), lack of high-frequency noise (Revnivtsev, Sunyaev, & Borozdin (2000)), and power-law photon spectrum extending to at least 120 keV (Wilson & McCollough (2000)). However, no dynamical mass measurement has yet been made for the compact object, preventing definitive identification as a black hole. The X-ray flux is low, not having exceeded 75 mCrab at any time to date. Tentative distance estimates to the putative optical counterpart based on companion spectral types in other X-ray transients (Uemura et al. (2000)) imply an unusually low X-ray luminosity ($`L_X10^{34}10^{35}`$ erg/s). Its Galactic latitude of 62 is extremely high for this source class, and the nominal distance would place it in the Galactic halo (Uemura et al. (2000)). The second outburst has prompted extensive multi-wavelength monitoring, including observations in radio, optical, EUV, soft X-rays and hard X-rays. A continuous spectrum is inferred from optical to hard X-rays with the photon spectral index determined in X-rays as $`1.8\pm 0.1`$ (Hynes et al. (2000)).

X-ray variability of the source is characterized in the time domain by repeated flaring on timescales of a few seconds. The power spectrum shows a variable, low-frequency QPO near 0.1 Hz, first reported from observations with the Proportional Counter Array (PCA) on RXTE (Revnivtsev, Sunyaev, & Borozdin (2000)). An outstanding feature of the source is that the QPO is seen in X-rays, EUV and optical (Haswell et al. (2000)). The optical source shows a photometric period, presumably orbital, of 4.1 hours (Uemura et al. (2000)). During the outburst the source has brightened from $`V\stackrel{>}{}18.8`$ in quiescence to $`V13`$ (Uemura et al. (2000)).

This Letter reports X-ray observations made with the Unconventional Stellar Aspect (USA) experiment on the US Air Force Advanced Research and Global Observation Satellite (ARGOS) and the PCA on RXTE. The $`0.1`$Hz QPO is detected in most USA and RXTE observations, providing a picture of the evolution of this QPO over more than 60 days with unprecedented coverage. The USA observing sequence was initiated as soon as possible after the first report (Remillard et al. (2000)) of the transient and has subsequently been maintained through the remainder of the second maximum of the outburst. The high revisit frequency of USA observations is meant to provide, among other things, a framework for linking up and correlating multi-wavelength observations.

## 2. Description of the USA experiment

USA is an X-ray timing experiment built jointly by the Naval Research Laboratory and the Stanford Linear Accelerator Center for the dual purposes of conducting studies of variability in X-ray sources and exploring applications of X-ray sensor technology (see Ray et al. 1999 for a more detailed description). USA was launched on 1999 February 23 on ARGOS into a nearly circular 830 km orbit at 98.8 inclination. It is a reflight of two proportional counter X-ray detectors flown previously on the NASA Spartan-1 mission (Kowalski et al. (1993)), after which they were recovered and refurbished. The primary observing targets are bright Galactic X-ray binaries, with a goal of obtaining large exposures on a modest number of sources.

The detector consists of two multiwire constant-flow proportional counters equipped with a 5.0 $`\mu `$m Mylar window and an additional 1.8 $`\mu `$m thick aluminized Mylar heat shield. The detector gas is a mixture of 90% argon and 10% methane (P-10) at a pressure of 16.1 psia (at room temperature). The detectors are sensitive in the range 1–15 keV with an effective area of about 1000 cm2 per detector at 3 keV. The collimators serve to support the window as well as to define the field of view, which is approximately 1.2 FWHM circular. The Crab Nebula gives about 4000 cts/s in one detector at the center of the field of view. The ARGOS spacecraft on which USA is mounted is three-axis-stabilized and nadir-pointed. Consequently, the X-ray detectors are mounted on a 2-axis gimballed platform to permit inertial pointing at celestial objects. The pitch and yaw drive capability is $`3.5`$/min (track) and $`20`$ /min (slew) and the common pitch/yaw pivot design allows 180 travel in each axis.

USA has five telemetry modes, four event and one spectral. Events are time-tagged to an onboard GPS receiver which provides an absolute time reference. Events are recorded with either 32 $`\mu `$s time resolution and 16 pulse height channels or 2 $`\mu `$s resolution with 8 channels. Event modes can be used up to count rates of $``$1000 cts/s at 40 kbps or $``$6000 cts/s at 128 kbps before events are discarded. Spectral mode records a 48 channel spectrum every 10 ms.

The USA instrument has performed well since activation began on 1999 April 30, but the mission has not been without difficulties. Approximately two weeks after launch the detector heat shields suffered from degradation which has imposed additional constraints on USA pointing with respect to the Sun. On 1999 June 8, Detector 2 developed a rapid gas leak, possibly caused by a micrometeor impact, and exhausted its gas supply leaving only Detector 1 to complete the mission, halving the effective area.

## 3. Observations and Data Analysis

USA has obtained a total of 425 ks of data on XTE J1118+480 between 2000 April 10 and June 13. These observations are continuing through the duration of the outburst. Using the 32 $`\mu `$s event mode, the source was observed between 5 and 11 times per day, providing the most consistent sampling of its X-ray behavior in existence. We have also utilized the 24 public target-of-opportunity (TOO) observations of XTE J1118+480 made by RXTE between 2000 April 13 and June 11. Analysis of the USA and RXTE data, including energy selection, binning, construction of the power spectra, and energy spectral fitting, was performed using FTOOLS v5.0.1. Typically the observed QPO has a fractional RMS amplitude of $`5`$% and a FWHM of $``$0.01 Hz. Additionally, in some cases there are suggestions of multiple peaks in the power spectra. Detailed results of our QPO analysis are summarized in Table 1.

### 3.1. USA

To search for QPOs in the USA data, we created 30 groups of observations, each of which contain 4–8 individual observations from within a continuous interval of 8–20 hours. As a result of the nearly polar orbit of USA, the usable time from each individual observation varies in length between 300 and 1100 seconds; hence the total on-source integration time contained in each group ranges from 1.8 to 4.9 ks. To construct the power spectrum for each group of USA observations, we selected channels 1–8, corresponding to an energy range of 1–10 keV, and binned the data by either 0.1 or 0.125 seconds. Power spectra of length 4096 were computed and averaged for each observation, the expected Poisson level subtracted, and normalized to fractional RMS2/Hz (see Nowak et al. 1999 and references therein for a detailed description of power spectral analysis of X-ray data). The QPO parameters were determined by fitting the resultant power spectra between 0.03 and 0.8 Hz to a power law plus a Gaussian QPO feature. The centroid and FWHM of the feature were determined from the fit parameters and the fractional RMS attributed to the QPO was determined from the integral of the Gaussian profile. USA average fluxes for each observation are determined by subtracting a background model and correcting for obscuration of the detector by the instrument support structure. The USA spectral calibration is still being determined and spectral results from USA will be presented in a future paper.

### 3.2. RXTE

The power spectra from the RXTE observations were generated from the Standard1 data which bin all of the good counts in each detector at 0.125 s (Jahoda et al. (1996)). These time series were summed for all active PCUs and the power spectra were computed and fitted as described above for the USA data. Spectral fits were also made using the PCA data, employing the background models and response matrices provided by the RXTE team. We fit a pure power law model between 2.5 and 20 keV to each of the RXTE observations listed in Table 1. This provided a reasonably good fit to the data, but the residuals typically show excesses below 4 keV and near 6 keV. Fitting a power law between 8 and 25 keV to avoid these features, we find that the photon spectral index is consistent with a constant value of 1.73$`\pm `$0.01.

## 4. Discussion

### 4.1. Source Models

Examining the options for the nature of this unusual X-ray transient, we note that its emission characteristics and timescales require a compact stellar object, but the hard X-ray spectrum excludes a white dwarf. A neutron star model has not been supported by detection of signatures such as pulsations, Type I X-ray bursts, or high frequency QPO activity. Therefore the observed X-ray characteristics of XTE J1118+480 seem to be best explained by identification of this source as a binary black hole transient in the canonical low hard state. However, this source is far from a typical example of that class. The luminosity is relatively low, and it is also unusual for a black hole transient to make its first appearance in the low hard state. (In this regard Uemura et al. (2000) have suggested that the source is viewed at high inclination, but there are no observed eclipses to support this suggestion.) The possible Galactic halo location of XTE J1118+480 and the X-ray/EUV/optical QPO are unique, although previous transients of this type have shown X-ray and optical correlated variability. The ratio of the X-ray to optical flux is also extremely low, leading to speculation that this outburst of XTE J1118+480 may be a “mini-outburst” of the type seen in the black hole candidate system GRO J0422+32 (Hynes et al. (2000)). We also note that short-timescale X-ray and optical variability similar to that seen in XTE J1118+480 has been observed in the hard state of the black hole candidate GX 339$``$4 (Motch et al. (1983)). Merloni et al. (2000) suggest that both the optical and hard X-ray variability in XTE J1118+480 are due to magnetic flares in a Comptonizing corona and predict X-ray/optical correlations which can be tested with ongoing multi-wavelength observations. Alternatively, Titarchuk and Osherovich (2000) argue that the low frequency QPO is due to a global oscillation mode in the accretion disk and show that these frequencies carry information about the system as a whole such as the size of the disk and the mass of the central object.

### 4.2. QPO Evolution

Figure 1 shows the evolution of the source flux as monitored by the ASM on RXTE (Levine et al. (1996)) and the QPO centroid as measured by USA and RXTE. The QPO moves upward in frequency over the 62 days of observations, apparently monotonically, but with significant changes in rate. During the same period, the 2–12 keV X-ray flux slowly rises, then begins to decrease. Figure 2 shows the QPO centroid frequency versus measured USA (1–10 keV) fluxes. Clearly, the QPO frequency fails to track the source intensity according to a simple proportionality.

This effect could be correlated with change in the source energy spectrum. However, the power-law component between 8 and 25 keV seen with the PCA stays essentially constant while the QPO frequency changes by a factor of two. It is possible that a soft component, which peaks below the PCA energy range, tracks the state of the accretion disk and thus the QPO frequency. The variable flux in the power-law component would then be due to inverse-Compton scattered seed photons from the variable disk component.

Low-frequency QPOs have been seen in several black hole binaries in the low hard state, including Cyg X-1, GRO J0422+32, and GX339$``$4. Unlike the QPO observed in XTE J1118+480 , the $``$0.3 Hz QPO observed in the low hard state of GX339$``$4 appears to be relatively stable, but like XTE J1118+480 no simple correlation between the QPO frequency and source flux was found (Nowak, Wilms, & Dove (1999)). Variable low-frequency QPOs have been observed during the high and intermediate states of the black hole transient XTE J1550$``$564. In this source, the QPO frequency generally increases as the disk flux increases (Sobczak et al. (1999)). However, we note that the time evolution of the variable QPO in XTE J1550$``$564 is more erratic than in XTE J1118+480 with the QPO frequency increasing and decreasing multiple times within a similar interval ($``$60 days).

## 5. Summary

The QPO in XTE J1118+480 shows an upward drift from 0.07 Hz to 0.15 Hz over a 62 day time interval while the source intensity slowly rises and then decreases. The effect is clear but the explanation seems to call for changes in the disk state that are not completely specified by measurement of the PCA X-ray spectrum, possibly by a variable soft disk component which is not well determined by these observations. We are further investigating the variability by studying flaring in the time domain. The flare events, which have a large signal-to-noise ratio and occur on the same timescale as the QPO, may shed light on the physical processes behind the QPO. Continued X-ray observations will be important for covering a larger dynamic range in the flux and continuing to track the QPO centroid.

This paper is dedicated to the memory of Terry E. Crandall who contributed many ingenious and creative solutions to the USA Experiment software. We would like to acknowledge Ganwise Fewtrell for his contributions to the USA data analysis software. Basic research in X-ray Astronomy at the Naval Research Laboratory is supported by the Office of Naval Research. Work at SLAC was supported by Department of Energy contract DE-AC03-76SF00515. This work was performed while RMB held a National Research Council Research Associateship Award at NRL. This paper made use of quick-look results provided by the ASM/RXTE team (see http://xte.mit.edu).


text:
# References

SOCIAL PERCOLATION ON

INHOMOGENEOUS SPANNING NETWORK

Abhijit Kar Gupta111Postdoctoral fellow at: Raman Research Institute, Bangalore-560 080, India and Dietrich Stauffer

Institute for Theoretical Physics, Cologne University,
50923 Köln, Germany
E-mail: abhi@rri.ernet.in, stauffer@thp.uni-koeln.de

## Abstract

The Social Percolation model recently proposed by Solomon et al.  is studied on the Ising correlated inhomogeneous network. The dynamics in this is studied so as to understand the role of correlations in the social structure. Thus the possible role of the structural social connectivity is examined.

Keywords: dynamics of social systems, percolation, statistical mechanics of model systems, self-organization.

In recent times, there has been a great interest in applications of Statistical Physics in Social Science . These works are wide ranged; from the stock market data analysis to the ‘microscopic’ models to understand the social dynamics. The study in social dynamics is concerned with the concept of information propagation and thus the formation of opinion in a complex social network. Various approaches have been developed in this respect. A class of social impact models based on social impact theory has been dealt with in ref.. One of the main premises in such models are the interactions among the social agents and groups; the agents are influenced by one another while taking decision. Another approach is the assumption that the agents take decisions of their own (not influenced by each other), only the spreading of information is necessary through some kind of contact process, the philosophy of epidemiology (see a brief review ). The primary interest is to study through simple models how a complex process of opinion formation (and regulation) takes place among the agents in a social space when the information is shared locally by ‘word of mouth’ (‘nearest neighbour’ connection may be). This is prior to a domination by powerful media (e.g., TV, Newspapers, Internet) which may influence the agents (people) quite homogeneously. Every social system has a kind of feedback mechanism. Thus it is reasonable to expect that the quality of a product, initially released in the market, is also adjusted according to the ever upgraded opinions (preferences) of the agents and the economic constraints.

In this paper we examine a model which has been recently proposed by Solomon et al.  (also see ref.) in a context of how an opinion is formed, when a commercial product is released in the market. The product is sometimes chosen by a spanning population of a region or sometimes fails to hit the target. Naturally, a concept of percolation  comes in and it has been termed as ‘Social Percolation’ in ref.).

The above ‘social percolation’ model is based on a social network which, for example, may be a two dimensional regular network. The agents are assumed to be situated on the sites of the network which provides the frame of fixed connectivity. They decide to buy some product (or watch a movie) if that qualifies according to their personal preferences. There is no other external influence in decision making. Thus the ‘quality’ and the personal ‘preference’ are two variables in the model, the former is the global and the later being a local variable.

The steps are the following:

$``$ To start with it is assumed that the agents along one boundary line are informed about the product. Each agent has a certain independent personal preference ($`p_i`$, for the $`i`$th agent) to begin with. A uniform random distribution is chosen for this purpose so that $`1>p_i>0`$, for all $`i`$. The quality $`q`$ is assigned a certain value ($`q<1`$) at the beginning.

$``$ Now if the quality is greater than the preference of an agent, $`p_i<q`$, the agent $`i`$ responds and the information is passed on to the nearest neighbour uninformed agents. Otherwise the agent remains inactive and the information is not propagated to any of its nearest neighbours. After it responds, the agent $`i`$ enhances its personal preference from $`p_i`$ to $`p_i+\delta p`$. Else the inactive agent reduces its preference from $`p_i`$ to $`p_i\delta p`$ at the next time step. Those who remain uninformed in the process keep their preferences unchanged.

$``$ At the end of one sweep (one time unit) of the whole system, it is checked if any of the agents on the other end has got the information and has responded. If so that means the information is propagated from one end to the other of the community. That is, in the language of percolation, it is checked if the responded agents span or percolate the system from one end to the other. This is a case of commercial success. So, in such a case, for financial constraints, the quality is reduced from $`q`$ to $`q\delta q`$ for the next time step. Otherwise (in case of failure) the quality is increased from $`q`$ to $`q+\delta q`$.

The above three steps constitute the execution of the model of Solomon et al. which shows an interesting ‘self-organized’ social dynamics . A Leath-type algorithm is implemented to upgrade the sites sequentially from a set of agents (sites) who are informed at the start.

In the above work  it is assumed that the population is distributed homogeneously on a social space and so the dynamics is tested on a regular square lattice. Here in the present work we ask the question of how the social structure plays a role in the social dynamics. It may happen that the population is not homogeneously widespread in a big locality. Because of geographical and other reasons, there may be certain clusterings of houses which may be connected by pathways. Thus there can be an entire inhomogeneous but correlated network of such population formed in a big region. Our interest is to examine how this correlated inhomogeneous structure plays role in the dynamics of the simple model system in contrast to the earlier studies  done on homogeneous random percolating network.

For the purpose of this work, in order to generate correlated inhomogeneous network, we create spanning clusters of up spins at the critical temperature $`J/K_BT_c=0.44`$ of a 2D Ising model at zero field in the sense of Fortuin-Kasteleyn-Coniglio-Klein-Swendsen-Wang cluster definitions . We removed bonds between two parallel (up) spins with a probability $`P=\mathrm{exp}(2J/K_BT)`$. For a typical $`101\times 101`$ and a $`1001\times 1001`$ systems we iterate up to $`10^4`$ Monte Carlo steps per spin.

Now, precisely, we check for a spanning cluster of neighbouring up spins connected by unremoved bonds. The agents are now situated on a tenuous and inhomogeneously connected network. Thus instead of random site percolation as in Solomon et al. we use Ising-correlated site-bond percolation to take into account the herding of human beings. In the rest of the work we refer the homogeneous network of ref. as ’random case’.

Below we present our observations on various aspects of the dynamics and we discuss the related scenario. For our entire work we chose $`\delta p=0.001`$ and $`\delta q=0.001`$ and all the quantities are calculated on typical networks of size $`101\times 101`$ and $`1001\times 1001`$. The initial value of $`q`$ was chosen to be 0.5. The late time dynamics does not depend on this choice.

We enumerate the time evolution of the quality $`q`$ and the average preference $`<p>=(1/N)\underset{i=1}{\overset{N}{}}p_i`$, where $`N`$ is number of agents on the spanning network. Both the quantities increase with time while $`q`$ shows a plateau type behaviour (not shown here) for an intermediate regime and $`<p>`$ increases initially in a power law. In Fig.1, we have plotted the difference of the above two, $`q<p>`$, with time $`t`$. To demonstrate a size effect we also plot the same for a network of size $`500\times 500`$ (indicated in the figure). For a comparison we plot the above quantity for the case of homogeneous random case. The difference in the two cases is apparent: in the Ising correlated case the difference $`(q<p>)`$ ever increases with time after a period of nonmonoticity, while in the random case the same goes down or remains constant in the same time domain. This reason may be attributed to the more frequent failures (no spanning) in this case, the quality of the product had to be increased while at the same time because of less effective communication through the tenuous network, many of the agents remain uninformed and thus average preference, although it increases with time, remains low compared to $`q`$.

To check the time evolution of the distribution of the agents personal preferences, we plot the distribution of $`p_i`$’s for three different times, $`t=200`$ (early time), $`t=1000`$ (intermediate time) and at $`t=2500`$ (a relatively late time) in Fig.2, Fig.3 and in Fig.4 respectively. The value of $`q`$ is indicated in the figures in each case for a guideline. In all cases the average of $`p_i`$’s ($`<p>`$) is always less than $`q`$. In the early time the distribution is more or less uniform and soon the preferences start clustering around a value, a sharp peak appears and that shifts to the right. This reflects the fact that the agents are again and again exposed to the information of the product so that they thrive for the high quality product and thus enhance their personal taste. However, apart from the sharp peak in the intermediate and late times, there appears to be a number of agents whose preferences remain far below. These are counted in the average $`<p>`$ which brings down the value significantly from the most probable value and thus keeps it far below the quality ($`q`$) of the product provided at that time. We have checked that this is even more prominent at very late times when the distribution is clearly divided into two parts, one around the sharp peak on the right and a number of agents’ preferences cluster around 0. Thus this further demonstrates the drifting of average preference away from the quality value, as time passes in the Ising correlated case. For comparison with the random homogeneous case Fig.5 is plotted corresponding to that in Fig.4. The left side of the peak in Fig.5 decays roughly exponentially and thus the contribution of the agents having such preference values to the average is much less in this case as compared to the case of Fig.4 for Ising correlated network.

We enumerate some additional quantities in order to understand the additional influences of the correlated inhomogeneous network on the dynamics. We look at a magnetization-like quantity and compare our present case with that of the random percolation model. Let us define a quantity $`d_i`$, which we may call as decision of the agent $`i`$: $`d_i=+1`$, if the agent $`i`$ chooses to buy the product or go for a movie ($`p_i<q`$) and $`d_i=1`$ otherwise ($`p_i>q`$). The average decision, $`<d>=(1/N)\underset{i=1}{\overset{N}{}}d_i`$, in the model, interestingly, oscillates in time between two well separated branches. Again for a comparison we plot the above for Ising correlated and for the random case, shown in Fig.6 and Fig.7 respectively. Although the qualitative behaviour is similar in two cases, there is a noticeable quantitative difference which may be important in determining some other quantities of interest. The average decision $`<d>`$ is always positive (high or low) as can be seen in Fig.6 whereas the same attains negative values in the lower branch in the random case (Fig.7). This presumably tells us that a comparatively larger fraction of agents (due to the many smaller non-spanning clusters) abstain from buying products at each time, while in the Ising correlated case the fraction of agents who abstain from buying the product at each time step is comparatively less. However, as may be expected, there are a certain fraction of agents which are always difficult to reach due to ill connectivity (the structure of the tenuous network) in Ising correlated case. Therefore, their preferences remain low compared to the most probable preferences in the social network. Now we see that this fact is reflected further in other quantities we define below for brevity.

Within the framework of the social percolation model one may define a concept of ‘profit’ (or ‘loss’) and the accumulation of ‘wealth’. Money spent by the producer/ manufacturer is proportional to the quality $`q`$ of the product. Therefore, the price of the product could be proportional to $`q/N`$, where $`N`$ is the number of agents targeted. But the product is sold among the agents whose preferences $`p_i`$ are less than the quality $`q`$ at that time. So we can define a quantity which we may call ‘profit’ (the negative of this may be called ’loss’): Profit = $`(q/N)\underset{p_i<q}{}p_iq`$. The cumulative sum of the so-called profit (or loss) with time may be thought of as ‘wealth’.

We numerically evaluate the above quantities in the random case (homogeneous structure) and later test those quantities for the inhomogeneous Ising correlated network. The distinct two branches shown in Fig.8 may be a company’s profit (and loss) with time. In Fig.9 we draw a distribution of these profit (or loss) values with time. This shows a bimodal distribution. The ’wealth’ is plotted in Fig.10 where we see that the wealth is negative in the entire time domain shown but it tries to recover the total loss with time, the curves goes through a dip minimum and then a sharp rise.

In Fig.11, Fig.12 and in Fig.13 we demonstrate the above quantities (as indicated in the figure) for a Ising correlated network corresponding to Fig.8, Fig.9 and in Fig.10 respectively. Now the quantitative as well as a qualitative differences are apparent. In this type of social network (under the scope of the present simple model) a company does not probably suffer much loss. Of course in this type of analysis the ‘price’ is not fixed, neither the quality or the preference values are bounded. This, however, indicates that the social structure (this means the connectivity among agents here) seems to play a definite role in the process of dynamics of any financial interest.

In conclusion, we may say that the indication of self-organized like behaviour (as that found in the same model but for a random homogeneous lattice) is probably absent here in the Ising correlated network. However, this correlated, inhomogeneous and tenuous network brings into the possibilities of many other interesting differences as compared to the random case. The spread of information in all directions is what dictates the quantitative and in some cases qualitative change in some relevant quantities.

Acknowledgment:

This work is supported by SFB 341 at the Institute for Theoretical Physics, Cologne University. Discussions with A. Novak and Z-F Huang are thankfully acknowledged.


text:
# Time in Quantum Geometrodynamics

## 1 Introduction.

A proper introduction of time is one of the central issues of any viable programme of gravity quantization. Its resolution is important conceptually as it determines in a profound way the meaning of the quantization procedure or the meaning of the basic structure that contains quantum gravity as a particular case or an approximation. It is also likely that a proper understanding of this issue will provide the clues to answering the questions motivating gravity quantization in the first place, such as avoidance of singularities, the issues of the final and the initial state of the Universe, etc. It might contribute to a better understanding of the aspects of modern quantum field theories that depend crucially on the causal structure of spacetime. One can be referred to the recent paper  by J. Butterfield and C. J. Isham for an exhaustive exposition of motivational and conceptual problems of quantum gravity, as well as for an overview of modern programmes of gravity quantization.

Quite understandably, the issue of time is formulated differently in different approaches to gravity quantization. In what follows we restrict ourselves to what can be thought of as a slight extension of the canonical gravity quantization programme.

The standard canonical quantum gravity approach  is based on the classical dynamic picture of the evolving 3–geometry of a slicing of a spacetime manifold. The slicing is essentially a reference foliation of the spacetime manifold (endowed with a 4–geometry) with respect to which the canonical variables are assigned. It is usually parametrized by the lapse function $`N`$ and the shift functions $`N^i`$. The canonical variables are the 3–metric $`g_{ik}`$ on a spatial slice $`\mathrm{\Sigma }`$ of the foliation induced by the spacetime 4–metric, and their canonical conjugate matrix $`\pi ^{ik}`$. The latter is related to the extrinsic curvature of $`\mathrm{\Sigma }`$ when it is considered as embedded in the spacetime.

The customary variational procedure applied to the Hilbert action expressed in terms of the canonical variables produces Hamilton equations describing the time evolution of the canonical variables, with the Hamiltonian given as $`N+N^i_i`$ where $``$ and $`_i`$ are functions of the canonical variables and their spatial derivatives. The procedure is not extended to the derivation of the Hamilton–Jacobi equation in the usual manner as such an equation is rendered to be meaningless with the chosen set of canonical variables (cf. ) when the general covariance of the theory is taken into account.

As a way out, the general covariance is introduced in the variational principle from the very onset as the requirement of the action to be invariant with respect to variations of the lapse and shift which leads to the constraint equations (to simplify notations, we omit indices on components of $`g`$ and $`\pi `$ in all equations below)

$$(g,\pi ;x)=0,$$
(1)
and

$$_i(g,\pi ;x)=0,$$
(2)
imposed on the canonical variables on each slice. An important feature of general relativity is that its dynamics is fully constrained. It can be shown that if the geometry of spacetime is such that the constraints are satisfied on all the slices of all spatial foliations of spacetime, then the canonical variables necessarily satisfy the Hamilton evolution equations. This feature is often referred to as a key property of general relativity  and is interpreted as an argument that the entire theory is coded in the constraints, with the conclusion that the Hamilton equations are redundant and can be ignored in dynamic considerations. Substitution of $`\delta S/\delta g`$ in the place of $`p`$ in constraint equations leads to the new set of equations

$$(g,\frac{\delta S}{\delta g};x)=0,$$
(3)
and

$$_i(g,\frac{\delta S}{\delta g};x)=0,$$
(4)
the first of which is considered to be the Hamilton–Jacobi equation. This assertion is supported by arguments appealing to the variational principle on superspace of 3–geometries (for detailed arguments and the interpretation of other equations cf. ).

Dirac’s procedure of canonical gravity quantization is based directly on this Hamilton–Jacobi equation and produces the quantum theory that consists of commutation relations imposed on all canonical variables and the Wheeler – DeWitt equation.

The ADM square root quantization procedure is also based entirely on constraints, but in this procedure the set of canonical variables is split in two subsets, embedding variables (four of them altogether; one slicing parameter $`\mathrm{\Omega }`$ and three coordinatization parameters $`\alpha `$) and true dynamic variables $`\beta `$ (two of them) , , . The set of constraints is solved with respect to the momenta conjugate to the embedding variables. After substituting $`\delta S/\delta \mathrm{\Omega }`$, $`\delta S/\delta \alpha `$ for $`p_\mathrm{\Omega }`$, $`p_\alpha `$, (where $`S`$ is the principal Hamilton function) one of the resulting equations (the equation for the momentum conjugate to the slicing parameter) is identified with the Hamilton–Jacobi equation, and its right hand side yields an expression for a new (square root) Hamiltonian. The quantization is based on this equation and produces the quantum theory that consists of the Schrödinger equation and commutation relations imposed on true dynamic variables and their conjugate momenta.

In both approaches, the description of time evolution of quantized gravitational fields or systems including such fields becomes extremely troublesome. Any attempt to introduce time that can be used in a way similar to that of time in quantum mechanics or in quantum field theory on a flat background invariably leads to the notorious problems of time , some of which are of a conceptual nature while others are technical. Attempts to introduce time in such systems in a universal way from outside, as a reading of a specially designed clock have been unsuccessful and there are all reasons to believe that it is impossible , whether the clock is believed to be gravitational (the readings depend only on the variables describing gravitational field) in its nature or it is a matter clock, for as long as it interacts with gravity.

The difficulties of the conceptual nature (the problem of functional evolution, and the multiple choice problem, in Kuchař’s terminology) emerge due to the dual nature of time parametrization in general relativity. If spacetime is considered as a manifold it can be coordinatized and sliced in any arbitrary manner. However, this is not sufficient for the description of geometrodynamic evolution. Both slicing and coordinatization need to be tied to the metric on the manifold. The standard way of doing it in classical geometrodynamics involves lapse and shift. These parameters express slicing condition in terms of readings of the clocks of resting test particles, and coordinatization conditions as the metric shift of the coordinate grid on the slice with time. In such a description, the reduction of the dynamic picture to the constraints is based entirely on existence of the unique spacetime metric (although the metric might be not known until the geometrodynamic problem is resolved). While this is not a problem in classical general relativity, there is, in general, no possibility to assign such a unique metric on spacetime in canonical quantum gravity.

With this in mind, keeping the constraints equations as a foundation of geometrodynamics becomes not very meaningful, to say the least. Even in the classical theory, such a procedure results in the treatment of the Hamilton–Jacobi equation, that looks quite artificial and totally different from the standard mechanical considerations. From the point of view of Hamilton dynamics, shift and lapse invariance is just a symmetry of the system and should not be used in deriving dynamic equations.

In fact, the derivation of Hamilton equations does not depend on assumptions of lapse and shift invariance of action (general covariance). It is only in the derivation of the Hamilton–Jacobi equation this invariance becomes involved in a fundamental way and essentially replaces standard dynamic considerations.

As a result, the Hamiltonian of Hamilton equations does not coincide with the Hamiltonian of the Hamilton–Jacobi equation, only the first being related to the Lagrangian in a standard way. Also, the Hamilton–Jacobi equation does not contain the reference to the time evolution at all.

Quantization of the dynamic picture based on the constraints is essentially equivalent to restricting the states of the resulting quantum systems to a “shell” determined by the constraints that are classical in their origin. An attempt to undertake a similar action in quantum mechanics or quantum field theory would be quite disastrous under all but very carefully selected conditions.

One way to resolve this dilemma would be to weaken the requirement of covariance, essentially discarding it in dynamic considerations and recovering it by imposing symmetries on solutions only to the extent and in the sense that is allowed by dynamics. The general covariance in its traditional meaning should be recovered in the classical limit. In a sense, this requirement should determine, at least partially, what constitutes the classical limit in quantum geometrodynamics.

In order to achieve this goal a formulation of geometrodynamic Hamilton–Jacobi theory independent of the symmetry assumptions is needed, in the same spirit as in a standard setting of mechanics. Is there a possibility of writing the Hamilton–Jacobi equation in a way that is closer to that encountered in mechanics? As we have mentioned, the answer is no, if one considers as the object of geometrodynamics the whole 3–metric (or 3–geometry) of the slice .

The situation changes dramatically if York’s analysis of gravitational degrees of freedom  is taken into account and actively utilized. It becomes possible to reformulate classical geometrodynamics in a standard way from the very beginning to the very end and to treat general covariance as a symmetry of gravitational systems. Although the Hamilton–Jacobi equation looks different and makes sense in its traditional form, the resulting description is equivalent to the commonly accepted in classical geometrodynamics. However, quantization based on this new Hamilton–Jacobi equation provides an appropriate interpretation of the conceptual problems of time making them quite natural statements concerning the properties of gravity quantization. It also seems to avoid the technical problems of time, such as the Hilbert space problem, and the spectral analysis problem, as it produces the Schrödinger equation for the state evolution, and the Hamiltonian does not include the square root operation. The procedure has been described previously elsewhere , , but we do not believe that it is widely known and provide a brief overview of it below.

In this setting, time can be introduced as a slicing parametrization on the spacetime manifold and tied to the metric structure without any contradictions. The metric interpretation of time is coupled with geometrodynamic evolution. The true meaning of time becomes completely determined only after the geometrodynamic evolution problem has been solved. In a sense, quantum geometrodynamic configuration and time emerge together and the meaning of the clock readings is influenced by quantum gravitational system.

We illustrate the emerging meaning of time by considering different examples and introducing times as counted by matter clocks and by gravitational clocks. These clocks are akin to the clocks measuring timelike intervals along the world lines of test particles (matter clock) and the York’s extrinsic time parameter. In the end, we discuss the general features of time parametrization and its metric interpretation in quantum geometrodynamics.

## 2 Geometrodynamic Quantization in General <br>Setting.

According to York’s analysis of gravitational degrees of freedom, the set of six parameters describing the slice 3–metric should be split in two subsets, $`\{\beta _1,\beta _2\}`$ (two functions) and $`\{\alpha _1,\alpha _2,\alpha _3,\mathrm{\Omega }\}`$. The first of these is treated as the set of true gravitational degrees of freedom (the initial values for them can be given freely), while the second is considered to be the set of embedding variables. The $`\alpha `$ parameters are often referred to as coordinatization parameters, while $`\mathrm{\Omega }`$ is called, depending on the context, the slicing parameter, the scale factor, or the many–fingered time parameter. Information relevant to the gravity field is carried by $`\beta `$ parameters, while $`\alpha `$ and $`\mathrm{\Omega }`$ essentially describe time. In the original York’s analysis the $`\beta `$ variables describe the conformal part of a slice 3–geometry, $`\mathrm{\Omega }`$ represents the scale factor and the $`\alpha `$ variables are determined by the choice of coordinatization of 3–slices. The true dynamic variables form what we call a dynamic superspace while the embedding variables are treated as functional parameters.

The idea is to develop geometrodynamics from the very beginning on the dynamic superspace instead of the superspace of 3–metrics or 3–geometries. The variational principle on the dynamic superspace or its phase space (formed by true dynamic variables $`\{\beta _1,\beta _2\}`$ and their conjugate momenta $`\{\pi _{\beta _1},\pi _{\beta _2}\}`$ yields the dynamic equations describing evolution of true dynamic variables. All of these equations depend on lapse and shift and contain embedding variables as functional parameters. These are treated as an external field and is determined by additional equations that do not follow from the variational principle on dynamic superspace. The quantization procedure is performed on the dynamic superspace (only $`\beta `$-s are quantized, i. e. generate commutation relations, while embedding variables form a classical field). The Schrödinger equation is obtained by a quantization procedure from the Hamilton–Jacobi equation on the dynamic superspace and describes the time evolution of the state functional on true dynamic superspace coupled with the external classical field determined by embedding variables. In this paper we introduce such a coupling via a procedure similar to that of Hartree–Fock.

In a more detailed and precise description that follows, we omit indices on variables $`\beta `$ and $`\alpha `$ for the sake of notational simplicity. They can be recovered easily whenever it becomes necessary.

We start from the standard Lagrangian $``$ (written in terms of the 3–metric, shift and lapse) and the associated action (with appropriate boundary terms, as needed, to remove the terms containing second time derivatives) and we introduce the momenta conjugate to the true dynamic variables

$$\pi _\beta =\frac{}{\dot{\beta }}.$$
(5)
We then use these $`\pi _\beta `$’s to form the geometrodynamic Hamiltonian $`_{dyn}`$,

$$_{dyn}=\pi _\beta \dot{\beta }.$$
(6)
The arguments of the Hamiltonian $`_{dyn}`$ are described by the expression

$$_{dyn}=_{dyn}(\beta ,\pi _\beta ;\mathrm{\Omega },\alpha ).$$
(7)
The variables following the semicolon are treated as describing an external field, while the ones preceding the semicolon are the coordinates and momenta of the gravitational true degrees of freedom, i.e. of the true geometrodynamics. The variation of $`\beta `$ and $`\pi _\beta `$ leads to the equations of geometrodynamics, i.e. to two pairs of Hamilton equations,

$$\dot{\beta }=\frac{_{dyn}}{\pi _\beta },$$
(8)
$$\dot{\pi }_\beta =\frac{_{dyn}}{\beta },$$
(9)
and, subsequently, to the Hamilton–Jacobi equation

$$\frac{\delta S}{\delta t}=_{DYN}(\beta ,\frac{\delta S}{\delta \beta };\mathrm{\Omega },\alpha ).$$
(10)
Here $`S`$ is a functional of $`\beta `$ and, in addition, a function of $`t`$,

$$S=S[\beta ;t).$$
(11)
and $`\frac{\delta }{\delta t}`$ is defined by

$$\frac{S}{t}=\frac{\delta S}{\delta t}d^3x.$$
(12)
Neither the Hamilton equations (8), (9) nor the Hamilton–Jacobi equation (10) are capable of providing any predictions as their solutions depend on the functional parameters $`\mathrm{\Omega }`$ and $`\alpha `$ which are not yet known. One can complete the system of equations by adding to the Hamilton equations, or to the Hamilton–Jacobi equation, the standard constraint equations of general relativity. They should be satisfied when the solution for $`\beta `$, $`\pi _\beta `$ of equations of true geometrodynamics (with appropriate initial data) is substituted in them (we use symbols $`[\beta ]_s`$, $`[\pi _\beta ]_s`$ for such a solution)

$$\begin{array}{c}^i([\beta ]_s,[\pi _\beta ]_s,\mathrm{\Omega },\alpha )=0\\ \\ ([\beta ]_s,[\pi _\beta ]_s,\mathrm{\Omega },\alpha )=0\end{array}$$
(13)
These constraint equations cannot be derived from variational principles on dynamic superspace. Rather, they should be treated as additional symmetries, or the equations for an external field. They do follow from the shift and lapse invariance of the action but their derivation in this new setting depends on the structure of the whole action integral (cf. section 4). As a result, they cannot replace the full set of equations for geometrodynamic evolution. However, the resulting complete system of equations (dynamic equations on conformal superspace and constraint equations) is equivalent to this of the standard geometrodynamics on the superspace of 3–geometries .

For the purpose of quantization, we make a transition to the corresponding Schrödinger equation based entirely on dynamics and ignoring the system symmetries

$$i\mathrm{}\frac{\delta \mathrm{\Psi }}{\delta t}=\widehat{}_{dyn}(\beta ,\widehat{\pi }_\beta ;\mathrm{\Omega },\alpha )\mathrm{\Psi }$$
(14)
where $`\widehat{\pi }_\beta =\frac{\mathrm{}}{i}\frac{\delta }{\delta \beta }`$. The Schrödinger equation (14) implies that commutation relations are imposed only on true dynamic variables and treats embedding variables as external classical fields. The state functional $`\mathrm{\Psi }`$ in this equation is a functional of $`\beta `$ and a function of $`t`$.

$$\mathrm{\Psi }=\mathrm{\Psi }[\beta ,t)$$
(15)
This Schrödinger equation (with specific initial data) can be solved (cf., for instance the example of the Bianchi 1A cosmological model below). The resulting solution $`\mathrm{\Psi }_s`$ of this Schrödinger equation is not capable of providing any definite predictions as it depends on four functional parameters $`\mathrm{\Omega }`$, $`\alpha `$ which remain at this stage undetermined. All expectations, such as the expectation values of $`\beta `$

$$<\beta >_s=\mathrm{\Psi }_s|\beta |\mathrm{\Psi }_s=\mathrm{\Psi }_s^{}\beta \mathrm{\Psi }_s𝒟\beta $$
(16)
or of $`\widehat{\pi }_\beta `$

$$<\pi _\beta >_s=\mathrm{\Psi }_s|\widehat{\pi }_\beta |\mathrm{\Psi }_s=\mathrm{\Psi }_s^{}\widehat{\pi }_\beta \mathrm{\Psi }_s𝒟\beta $$
(17)
also depend on these functional parameters. To specify these functions we resort to the constraint equations. The treatment of the constraints has nothing to do with quantization of geometrodynamics. It is merely introducing the coupling between already quantized geometrodynamics and the classical field determined by embedding variables. In other words, they take care of the symmetries, which are classical in their nature, to the extent they are capable of doing that.

As in case of classical geometrodynamics, we impose the constraints on the solution of the dynamic equations (Schrödinger equation) with appropriate initial data and in this way, determine the unique values of $`\mathrm{\Omega }`$ and $`\alpha `$. It is possible that there are several ways to couple the constraints to the quantization of the true dynamic variables, $`\beta `$. As per our previous proposal we impose the four constraints only on the expectation values of the conformal dynamics

$$\begin{array}{c}^i(<\beta >_s,<\pi _\beta >_s,\mathrm{\Omega },\alpha )=0\\ \\ (<\beta >_s,<\pi _\beta >_s,\mathrm{\Omega },\alpha )=0\end{array}$$
(18)
The way evolution occurs can be described as follows. Initial data at $`t=t_0`$ consist of the initial state functional $`\mathrm{\Psi }=\mathrm{\Psi }_0`$ and the initial values (functions) of embedding variables. In addition, lapse and shift are supposed to be given. Equations (16), (17) yield the expectation values (functions) of true dynamic variables and their conjugate momenta. The result are substituted in the constraints (18). After this, the constraints are solved with respect to the time derivatives of embedding variables. A step forward in time (say, with the increment $`\mathrm{\Delta }t`$) is performed by integration of obtained expressions to evolve embedding variables and and by integration of the Schrödinger equation (14) to evolve the state functional. This concludes one step forward in time. The next step is performed by repeating the same operations in the same order.

One can be referred to ,  for two particular examples illustrating such geometrodynamic evolution in cases of Bianchi 1A cosmology and Taub cosmology. The first one can and has been solved analytically, while the latter one has been solved numerically.

We provide an abbreviated description of only the first example (Bianchi 1A) as its analytical solution is more useful in presenting the issue of time in quantum geometrodynamics.

## 3 Geometrodynamic Quantization: Bianchi 1A Model.

The Bianchi 1A cosmological model is commonly referred to as the axisymmetric Kasner model . Its metric is determined by two parameters, the scale factor $`\mathrm{\Omega }`$ and the anisotropy parameter $`\beta `$

$$ds^2=dt^2+\mathrm{e}^{2\mathrm{\Omega }}\left(\mathrm{e}^{2\beta }dx^2+\mathrm{e}^{2\beta }dy^2+\mathrm{e}^{4\beta }dz^2\right).$$
(19)
The choice of this expression for the metric implies that we have chosen $`N^i=0`$ and $`N=1`$ values of shift and lapse for this example. As this cosmology is homogeneous the two functions $`\mathrm{\Omega }`$ and $`\beta `$ are the functions of the time parameter $`t`$ only. The scalar 4–curvature can be expressed in terms of these two functions to yield the Hilbert action and, after subtracting the boundary term, the cosmological action,

$$I_C=I_H+\frac{3V}{8\pi }\dot{\mathrm{\Omega }}\mathrm{e}^{3\mathrm{\Omega }}|_{t_0}^{t_f}=\frac{3V}{8\pi }\underset{t_0}{\overset{t_f}{}}\left(\dot{\beta }^2\dot{\mathrm{\Omega }}^2\right)\mathrm{e}^{3\mathrm{\Omega }}𝑑t,$$
(20)
where $`V=𝑑x𝑑y𝑑z`$ is the spatial volume element. As it is usually done in case of homogeneous cosmologies, we integrate appropriate quantities over the spatial volumes and work with integrated Lagrangian $`L`$, Hamiltonian $`H`$ and momenta $`p_\beta `$ rather than with their densities $`,,\pi _\beta `$.

We treat the scale factor $`\mathrm{\Omega }(t)`$ as the embedding variable and the anisotropy $`\beta (t)`$ as the dynamic degree of freedom. The momentum conjugate to $`\beta `$ is

$$p_\beta =\frac{L}{\dot{\beta }}=\frac{3V}{4\pi }\mathrm{e}^{3\mathrm{\Omega }}\dot{\beta }.$$
(21)
The Hamiltonian of the system in our approach can be expressed in terms of the momentum conjugate to $`\beta `$ and the Lagrangian.

$$H_{dyn}=p_\beta \dot{\beta }L=\frac{2\pi }{3V}\mathrm{e}^{3\mathrm{\Omega }}p_\beta ^2+\frac{3V}{8\pi }\dot{\mathrm{\Omega }}^2\mathrm{e}^{3\mathrm{\Omega }}.$$

In the classical theory this Hamiltonian can be used to produce either one pair of Hamilton equations or the equivalent Hamilton–Jacobi equation. In either case, the dynamics picture derived in this way is incomplete. To complete it, we impose the super-Hamiltonian constraint.

$$p_\beta ^2=\left(\frac{3V}{4\pi }\right)^2\mathrm{e}^{6\mathrm{\Omega }}\dot{\mathrm{\Omega }}^2.$$
(22)
Using the Hamilton–Jacobi equation,

$$\frac{S}{t}=H_{dyn}(\frac{S}{\beta },\mathrm{\Omega }(t),\dot{\mathrm{\Omega }}(t)),$$
(23)
together with the expression (3) for the Hamiltonian $`H_{dyn}`$, we obtain the Schrödinger equation for the axisymmetric Kasner model.

$$i\mathrm{}\frac{\mathrm{\Psi }}{t}=\frac{2\pi \mathrm{}^2}{3V}\mathrm{e}^{3\mathrm{\Omega }}\frac{^2\mathrm{\Psi }}{\beta ^2}+\frac{3V}{8\pi }\dot{\mathrm{\Omega }}^2\mathrm{e}^{3\mathrm{\Omega }}\mathrm{\Psi }.$$
(24)
The constant $`\mathrm{}`$ in this equation should be understood as the square of Planck’s length scale, rather than the standard Planck constant. The quantum picture based on the Schrödinger equation (24) is incomplete as the scale factor $`\mathrm{\Omega }`$ is so far an unknown function of time. To complete the dynamics picture we follow our prescription and impose, in addition to equation (24), the super-Hamiltonian constraint.

$$<p_\beta >_s^2=\left(\frac{4\pi }{3V}\right)^2\mathrm{e}^{6\mathrm{\Omega }}\dot{\mathrm{\Omega }}^2.$$
(25)
Here $`<p_\beta >_s`$ is the expectation value of the momentum $`\widehat{p}_\beta =\frac{\mathrm{}}{i}\frac{}{\beta }`$

$$<p_\beta >_s=\mathrm{\Psi }_s|\widehat{p}_\beta |\mathrm{\Psi }_s=\underset{\mathrm{}}{\overset{\mathrm{}}{}}\mathrm{\Psi }_s^{}(\beta ,t)\widehat{p}_\beta \mathrm{\Psi }_s(\beta ,t)𝑑\beta $$
(26)
where $`\mathrm{\Psi }_s`$ is the solution of the Schrödinger equation with specified initial data. The system of equations (24), (25) provide us with a complete quantum dynamic picture of the axisymmetric Kasner model evolution and, when augmented by appropriate initial and boundary conditions, can be solved analytically.

For instance, we can specify the initial data for the Schrödinger equation in the form

$$\mathrm{\Psi }(\beta ,t)|_{t_0}=\mathrm{\Psi }_s(\beta ,t_0)=\underset{\mathrm{}}{\overset{\mathrm{}}{}}A_k\mathrm{e}^{\frac{i}{\mathrm{}}k\beta }𝑑k$$
(27)
with

$$A_k=C\mathrm{e}^{a(kk_0)^2}.$$
(28)
Such a choice introduces the Gaussian wave packet centered initially at the value $`k_0`$ of $`k`$ (we will describe the meaning of $`k_0`$ later), and of the initial width determined by the constant $`a`$, with $`C`$ being merely the normalization constant, picked to satisfy the condition $`\mathrm{\Psi }(\beta ,t_0)|\mathrm{\Psi }(\beta ,t_0)=1`$ .

The solution of the Schrödinger equation with such initial data can be written as

$$\mathrm{\Psi }_s(\beta ,t)=C\sqrt{\pi }\left(a^2+\frac{f^2}{\mathrm{}^2}\right)^{\frac{1}{4}}\mathrm{exp}\left\{\frac{a}{4\left(a^2+\frac{f^2}{\mathrm{}^2}\right)}\frac{(\beta 2k_0f)^2}{\mathrm{}^2}\right\}\mathrm{e}^{iF},$$
(29)
where

$$f=f(t)=\frac{2\pi }{3V}\underset{t_0}{\overset{t}{}}\mathrm{e}^{3\mathrm{\Omega }}𝑑t,$$
(30)
and $`F=F(\beta ,t)`$ is rather involved real valued expression  that depends on $`\beta `$, $`\mathrm{\Omega }`$, and $`\dot{\mathrm{\Omega }}`$. However, it has a structure that makes computations of relevant expectation values easy. Such a computation of the expectation $`p_\beta _s`$ of the momentum $`\widehat{p}_\beta =\frac{\mathrm{}}{i}\frac{}{\beta }`$ yields

$$<p_\beta >_s=\mathrm{\Psi }_s|\widehat{p}_\beta |\mathrm{\Psi }_s=k_0.$$
(31)
In other words the expectation value of the momentum $`<p_\beta >_s`$ does not change with time. It is determined by the $`k`$–center of the packet at $`t=t_0`$.

It is clear that this solution of the Schrödinger equation describing the wave packet time evolution cannot provide any definite predictions as it contains as yet undetermined scale factor $`\mathrm{\Omega }`$. To find $`\mathrm{\Omega }(t)`$ we need to substitute this expectation value into the constraint (25) and to solve the resulting equation with respect to $`\mathrm{\Omega }`$. Substitution of $`p_\beta _s=k_0`$ in (25) yields

$$k_0^2=\left(\frac{3V}{4\pi }\right)^2\mathrm{e}^{6\mathrm{\Omega }}\dot{\mathrm{\Omega }}^2.$$
(32)
Once the solution of this equation is substituted in (29) the geometrodynamic problem (24), (25) for the wave packet (28) is solved completely.

The solution can be used to compute the expectation value for $`\beta `$:

$$<\beta >_s=\mathrm{\Psi }_s|\beta |\mathrm{\Psi }_s=2k_0f(t),$$
(33)
and the variance in $`\beta `$

$$<\left(\beta <\beta >\right)^2>_s=\mathrm{\Psi }_s|\left(\beta <\beta >\right)^2|\mathrm{\Psi }_s=\frac{\mathrm{}^2a^2+f^2}{a}$$
(34)
Thus “the center” of the wave packet evolves as the classical Kasner universe determined by the momentum value equal to $`k_0`$ would evolve, while the spread of the packet increases with time. The result is similar to that of the quantum mechanics of a free particle; after all the Bianchi I cosmology is the free–particle analogue of quantum cosmology.

## 4 Time in Quantum Geometrodynamics.

In sections 2, 3 we have demonstrated that a proper understanding of classical geometrodynamics as the dynamics on the superspace of true dynamic variables amended by the constraints attributed to the universal symmetries of gravitational systems (lapse and shift invariance) opens a way to circumvent the problems of time (for a full discussion, cf. , ).

The important point is that these problems disappear as soon as the proper object of quantization (true geometrodynamic variables) is chosen. The constraints themselves are of no primary significance in this process. Their presence in the theory reflects the fact that for gravity the true dynamics cannot be decoupled from the evolution of “embedding” variables and that within the classical theory the true dynamic variables picture is limited to a “shell” determined by the constraints. There is absolutely no reason to expect that the last feature will survive after quantization, except, perhaps, for some particular carefully chosen systems. We do not require it and thus avoid the problems of time.

The particular choice of lapse ($`N=1`$) and shift ($`N^i=0`$) in the previous section is quite sufficient to make this point. However, such a choice becomes an obstacle for understanding the subtle differences between our treatment of constraints and the standard one. Also it fixes a particular choice of time and its interpretation, thus precluding the study of alternatives in the important issue of time in quantum geometrodynamics. Essentially, we have become tied to the classical matter clocks at rest as determined by spacelike slices of the Kasner universe. We cannot even switch to another classical clock, such as the one that produces the trace of the extrinsic curvature as its readings, while such a clock might be relevant in resolving the issue of the final singularity for some cosmological models. Such a transition demands releasing of at least lapse $`N`$ and making it a function of time $`N=N(t)`$ ($`N=N(K)`$, if we want to assign $`t=K=\mathrm{Tr}(K)`$).

This amounts to the choice of metric in the form

$$ds^2=\left[N(t)dt\right]^2+\mathrm{e}^{2\mathrm{\Omega }(t)}\left[\mathrm{e}^{2\beta (t)}dx^2+\mathrm{e}^{2\beta (t)}dy^2+\mathrm{e}^{4\beta (t)}dz^2\right].$$
(35)
instead of (19) (we still retain zero shift $`N^i=0`$), which results in the expression for the Lagrangian

$$L=\frac{3v}{8\pi }\frac{1}{N}\left(\dot{\beta }^2\dot{\mathrm{\Omega }}^2\right)\mathrm{e}^{3\mathrm{\Omega }}$$
(36)
The momentum conjugate to $`\beta `$ becomes

$$p_\beta =\frac{L}{\dot{\beta }}=\frac{3V}{4\pi }\frac{1}{N}\mathrm{e}^{3\mathrm{\Omega }}\dot{\beta }.$$
(37)
The Hamiltonian of the system can be expressed now in terms of the momentum conjugate to $`\beta `$ and the Lagrangian.

$$H_{dyn}=p_\beta \dot{\beta }L=\frac{2\pi }{3V}N\mathrm{e}^{3\mathrm{\Omega }}p_\beta ^2+\frac{3V}{8\pi }\frac{1}{N}\dot{\mathrm{\Omega }}^2\mathrm{e}^{3\mathrm{\Omega }}.$$

The expression for the action integral

$$I=p_\beta 𝑑\beta H_{dyn}dt=p_\beta 𝑑\beta \left[\frac{2\pi }{3V}N\mathrm{e}^{3\mathrm{\Omega }}p_\beta ^2+\frac{3V}{8\pi }\frac{1}{N}\dot{\mathrm{\Omega }}^2\mathrm{e}^{3\mathrm{\Omega }}\right]dt$$
(38)
differs in its appearance from the standard one. In particular, lapse $`N`$ is not a mere factor in front of some expression anymore. However, it is easy to see that variation of the action with respect to $`N`$ produces the constraint

$$p_\beta ^2=\left(\frac{3V}{4\pi }\right)^2\frac{1}{N^2}\mathrm{e}^{6\mathrm{\Omega }}\dot{\mathrm{\Omega }}^2.$$
(39)
that coincides with the Hamiltonian constraint of the standard approach when the variables of two different approaches are properly identified.

Just as in previous section, we use the Hamilton–Jacobi equation (23) together with the expression (4) for the Hamiltonian $`H_{dyn}`$, to obtain the Schrödinger equation for the axisymmetric Kasner model.

$$i\mathrm{}\frac{\mathrm{\Psi }}{t}=\frac{2\pi \mathrm{}^2}{3V}N\mathrm{e}^{3\mathrm{\Omega }}\frac{^2\mathrm{\Psi }}{\beta ^2}+\frac{3V}{8\pi }\dot{\mathrm{\Omega }}^2\frac{1}{N}\mathrm{e}^{3\mathrm{\Omega }}\mathrm{\Psi }.$$
(40)
To complete the dynamic picture, this equation should be amended by the equation

$$<p_\beta >_s^2=\left(\frac{4\pi }{3V}\right)^2\frac{1}{N^2}\mathrm{e}^{6\mathrm{\Omega }}\dot{\mathrm{\Omega }}^2.$$
(41)
obtained from the constraint equation (39) via the same procedure as equation (25) of the previous section.

The solution of these equations describing “propagation” of a wave packet of Kasner universes (cf. the previous section) can be written down right away. A simple inspection of equations (40), (41) reveals that they are obtained from the similar equations (24), (25) by replacing everywhere the factor $`\mathrm{e}^{3\mathrm{\Omega }}`$ for the factor $`N\mathrm{e}^{3\mathrm{\Omega }}`$. The solution of the Schrödinger equation (39) is expressed by (29) with $`f(t)`$ given by

$$f(t)=\frac{2\pi }{3V}\underset{t_0}{\overset{t}{}}N\mathrm{e}^{3\mathrm{\Omega }}𝑑t,$$
(42)
instead of that given by (30). Correction of the equation (32) is just as trivial.

Up until now the lapse function $`N(t)`$ and, together with it, the meaning of the time parameter has not been specified. A variety of approaches can be used to fix both of them. For a case when shift is fixed (as in our example), one more condition including $`N(t)`$ should be imposed. One extreme case has been considered in the previous section. The lapse has been fixed and, after this, the meaning of the time parameter has been interpreted. Another possibility is to fix the interpretation of the time parameter, in which case the lapse function will be fixed implicitly and will be determined by an additional condition considered together with the dynamic evolution equations.

As an example we consider the choice of $`K=\mathrm{Tr}K`$ as the time parameter. This choice might provide advantages in particular problems (such as avoidance of singularities, etc.) but we are not concerned with this now. In this case equations (40) – (42) take the form

$$i\mathrm{}\frac{\mathrm{\Psi }}{K}=\frac{2\pi \mathrm{}^2}{3V}N\mathrm{e}^{3\mathrm{\Omega }}\frac{^2\mathrm{\Psi }}{\beta ^2}+\frac{3V}{8\pi }\left(\frac{d\mathrm{\Omega }}{dK}\right)^2\frac{1}{N}\mathrm{e}^{3\mathrm{\Omega }}\mathrm{\Psi }$$
(43)
$$<p_\beta >_s^2=\left(\frac{4\pi }{3V}\right)^2\frac{1}{N^2}\mathrm{e}^{6\mathrm{\Omega }}\left(\frac{d\mathrm{\Omega }}{dK}\right)^2.$$
(44)
and

$$f(K)=\frac{2\pi }{3V}\underset{K_0}{\overset{K}{}}N(K)\mathrm{e}^{3\mathrm{\Omega }(K)}𝑑K,$$
(45)
After solving the Schrödinger equation and computing the expectation value $`<p_\beta >_s`$, (44) yields

$$k_0^2=\left(\frac{3V}{4\pi }\right)^2\frac{1}{N^2}\mathrm{e}^{6\mathrm{\Omega }}\left(\frac{d\mathrm{\Omega }}{dK}\right)^2.$$
(46)
which, together with the standard expression (this expression holds for any slicing parametrization, including parametrization by $`K`$)

$$K=\frac{3}{N}\frac{d\mathrm{\Omega }}{dK}$$
(47)
provides the basis for computing both $`\mathrm{\Omega }(K)`$ and $`N(K)`$. Indeed, (47) implies

$$\frac{1}{N}\frac{d\mathrm{\Omega }}{dK}=\frac{K}{3}$$
(48)
substitution of which in (46) yields the equation for $`\mathrm{\Omega }(K)`$

$$k_0^2=\frac{V^2}{16\pi ^2}K^2\mathrm{e}^{6\mathrm{\Omega }}$$
(49)
It is remarkable that for this particular time parametrization the equation (49) for $`\mathrm{\Omega }(K)`$ is an algebraic equation (generically it should be differential). This equation implies

$$\mathrm{\Omega }(K)=\frac{1}{6}\mathrm{ln}\left(\frac{V^2}{16\pi ^2}\frac{K^2}{k_0^2}\right)$$
(50)
which, together with (47), yields the expression for $`N(K)`$

$$N(K)=\frac{3}{K}\frac{d\mathrm{\Omega }}{dK}=\frac{1}{K^2}$$
(51)
## 5 Discussion.

The issue of time in the standard approach to canonical quantum gravity does not seem to have a satisfactory resolution no matter which of the two quantization procedures is used. Both Dirac and ADM quantization do not seem to be able to handle the emerging difficulties of conceptual and technical nature. The difficulties can be traced to identification of the entire 3–geometry of a spacelike slicing as the principal dynamic object. The resulting constrained dynamics consists of proper dynamics that follows from variations of dynamic variables, and of the constraints that enforce the fundamental symmetries of gravitational systems (general covariance) and are obtained by varying lapse and shift. The constraints essentially restrict the evolution of gravitational systems to a “shell”. The peculiarity of gravitation is that on this shell all of the dynamics is essentially determined by the constraints, thus rendering the proper dynamic relations obsolete and unnecessary. This does not present any problem in the classical theory. However, quantization of such a classical theory is reduced to quantization of the constraint equations. This essentially amounts to restricting the system states to the shell and neglecting “off shell” contributions, which leads to seemingly intractable problems in describing the time evolution of gravitational systems.

Our suggestion is to separate the true dynamics of gravity fields from enforcing the symmetries. Practically, it is done via putting the dynamic object in the geometrodynamic superspace of true gravitational degrees of freedom and formulating of variational principles on this superspace rather than on the superspace of total 3–geometries, while leaving the embedding variables as free parameters reserved for enforcing of the symmetries on the solutions of proper dynamic equations at a suitable moment. The equations enforcing the symmetries following from the general covariance are not obtained by variations on the geometrodynamic superspace. They come from variations of the action with respect to lapse and shift, as in standard approach. Within the realm of the classical theory, the final outcome of the theory predictions does not change (compared to the standard approach). Only the order of operations is changed. First, we consider all the solutions of the proper dynamic equations (on and off shell), and then force them on the shell via adjusting the embedding parameters. Since the latter enter the solutions for true dynamic variables, the final expressions for these variables, too, depend on the outcome of forcing the system on the shell. In the classical domain, the resulting theory is indistinguishable from the standard one.

The quantum theory, however, becomes quite different. We consider the true geometrodynamics (on the superspace of true dynamic variables or its phase space) as the object of quantization. The resulting equation of the theory is a Schrödinger equation with the Hamiltonian that is not a square root Hamiltonian and thus avoids most of the conceptual as well as technical problems of time. The initial values as well as the solutions of these equations contain embedding variables as functional parameters. They do not know anything about the shell (one can say that they include both on and off shell contributions). By themselves, they cannot either violate or enforce constraints. This should be achieved via adjusting embedding parameters. In addition, the solutions contain lapse and shift functions. These are responsible for the interpretation of time. When imposing constraints, it is important to realize that they might not be enforceable exactly (by simply treating them as additional operator equations obtained by the same procedure as the Schrödinger equation and applied to the state functional) if their operator versions do not commute with the Hamiltonian of the Schrödinger equation and cannot be made to commute via adjusting embedding parameters. It only means that, in this particular problem, off shell contributions cannot be neglected. A weaker version of the constraints should be introduced. The version chosen by us in examples was to impose the constraints on expectations. In this way the covariance requirement is satisfied as much as it can be satisfied. To a reasonable extent, it is exactly satisfied in the sense that predictions on the final hypersurface remain the same if the surface is not changing. The question of the change of the prediction when the surface is changing can be also answered positively by simple count of equations and adjustable parameters. The same is true, under reasonable conditions, about the multiple choice problem of time.

As about interpretation of time (metric interpretation of the slicing parameter) our consideration indicates that, within our approach, it is equivalent to imposing one (or more, if shift is involved) additional relation that involves lapse (or lapse and shift), and will not create any problems if these additional relations do not involve true dynamic variables. We do not believe that even relations involving true dynamic variables can cause considerable difficulties if they are introduced in the way similar to the one described for constraints. In our examples the relations included only embedding variables, thus commuting with the Hamiltonian and not causing any troubles. It should be stressed, however, that all three components of the evolution description for quantum geometrodynamic system — quantum dynamics itself, constraints enforcing the symmetries, and the interpretation of time — emerge together as the solution of the total problem of geometrodynamic evolution.

## Acknowledgements

For discussion, advice, or judgment on one or another issue taken up in this manuscript, we are indebted to S. Chakrabarti, R. Fulp, D. Holz, L. Norris, M. Maclean, and J. A. Wheeler. One of us (AK) gratefully acknowledges the support of this research by the Theoretical Division of the LANL during the Summer, 1999.