text:
# Imaging and Spatially Resolved Spectroscopy of AFGL 2688 in the Thermal Infrared Region 1footnote 11footnote 1Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.

## 1 Introduction

Low- to intermediate-mass stars (2–5 M) experience rapid transition from the asymptotic giant branch (AGB) to planetary nebula (PN) at the end of their life. An evolved star has high-mass loss when ascending the AGB, but at the tip of the AGB the star loses almost all the hydrogen envelope and cannot sustain high-mass loss. It then begins to evolve toward higher temperatures in H-R diagram while keeping almost constant luminosity as a proto-planetary nebula (PPN). The central star of the PPN eventually becomes hot enough ($`>`$30,000 K) to ionize the surrounding material ejected during the preceding mass-loss phase, and the PN is reached. The transition phase from PPN to PN is as short as about 5000 yr (Schönberner, 1983). PPN share several observational properties in common (Kwok, 1993, 2000), including a significant infrared excess at far-infrared wavelengths due to the warm circumstellar dust ($`T_{\mathrm{dust}}`$ = 150–250 K), a low-velocity molecular outflow (10–20 km s-1), a central star that is typically an F- to G- type supergiant, and extended or bipolar nebulosity in the visible.

AFGL 2688 is a prototypical carbon-rich PPN, one of the most extensively studied so far since the discovery of reflection nebulosity (Ney et al., 1975). AFGL 2688 is located 1.0–1.5 kpc away from us (Cohen & Kuhi, 1977) (hereafter we assume the distance to AFGL 2688 to be 1.2 kpc). From the visible spectroscopy at the reflection nebulosity, Crampton, Cowley, & Humphreys (1975) found that AFGL 2688 harbors an F5Ia star in the center of the nebula. It may have left the AGB only 200 years ago (Skinner et al., 1997).

In spite of many studies of AFGL 2688, the fundamental structure of the nebula is still actively debated. The old standard model requires equatorial density enhancement lying perpendicular to the polar axis to shape the optical bipolar nebula (Yusef-Zadeh, Morris, & White, 1984; Latter et al., 1993; Lopez & Perrin, 2000). The polar axis is defined by the optical bipolar nebula and lies at a position angle of about 15. However, recent millimeter and centimeter radio observations have revealed features hard to reconcile with the standard model. These structures including a high velocity $`{}_{}{}^{13}\mathrm{CO}`$ flow (Yamamura et al., 1995), a compact expanding shell (Cox et al., 2000), and elongated emission (Jura et al., 2000), and all have a PA $``$ 55. No present model can explain this apparent dual symmetry in AFGL 2688. Also surprising and unexplained is the quadrupole shape of the H2 emission (Latter et al., 1993; Sahai et al., 1998a).

The goal of our observations is to better understand these apparently contradictory morphological features. We infer that the features at PA $``$ 55 are young because they have relatively high velocity or are compact. To investigate the origin of these features, we must observe close to the central star. However, the inner 5″ of the nebula suffers from heavy obscuration, and is virtually invisible at wavelengths shorter than 2 µm. Imaging in the thermal infrared region with a 8 m class telescope is essential to access the innermost region of AFGL 2688, as well as to provide high spatial resolution comparable with the visible image.

This paper is organized as follows: In §2 we describe the observations and reduction of the data. In §3 we summarize current understanding of AFGL 2688 and discuss several new and recently found morphological features.

## 2 Observations and Data Reduction

### 2.1 $`L^{}`$ and $`M^{}`$ Imaging

The imaging observations were made on UT 2001 July 14 at the Subaru Telescope with the Infrared Camera and Spectrograph (IRCS) (Tokunaga et al., 1998; Kobayashi et al., 2000). IRCS is a cryogenic camera and spectrograph equipped with two 1024 $`\times `$ 1024 InSb arrays. We used special $`L^{}`$\- and $`M^{}`$-band filters optimized for the atmospheric transmission curve at Mauna Kea (Tokunaga, Simons, & Vacca, 2001). Images of $`L^{}`$ were recorded by 58 mas pixel-1 camera (60″$`\times `$ 60″ field of view) at five different locations with the array offset by 30″ to each time so that median sky frames could be made. The $`M^{}`$ images were recorded in the same manner, but with the 22 mas pixel-1 camera (23″$`\times `$ 23″ field of view) and smaller nodding offset (10″). The total on-source integration times were 270 and 360 s in $`L^{}`$\- and $`M^{}`$, respectively. Flux calibration was facilitated by observing a photometric reference star, HD 201941, selected from Elias et al. (1982). The seeing was good ($``$0$`\stackrel{}{\mathrm{.}}`$3 in FWHM in $`L^{}`$), and stable throughout the observations. The actual spatial resolution at $`M^{}`$ is somewhat inferior to the seeing FWHM because the high-resolution camera of IRCS was internally slightly out of focus.

The standard image reduction procedure was performed with IRAF.333IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation. After sky subtraction, flat-fielding by dividing by a sky flat, and bad pixel correction, image frames were registered and averaged. We show the final images in Figures Imaging and Spatially Resolved Spectroscopy of AFGL 2688 in the Thermal Infrared Region 111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. and Imaging and Spatially Resolved Spectroscopy of AFGL 2688 in the Thermal Infrared Region 111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. in linear and logarithmic scales. Image convolution was applied with a 2D Gaussian filter of 1.4 pixels FWHM to show the faint nebulosity.

The $`L^{}`$ and $`M^{}`$ images look quite similar to each other. In addition to the spindle features commonly observed at the shorter wavelengths, they reveal great details in the innermost region for the first time. The upper part of south lobe close to the center is very faint in the 2.15 µm image taken by HST/NICMOS (Sahai et al., 1998a), but it is the brightest region in the southern lobe at the thermal infrared wavelength. The north lobe appears by far brighter near the center and harbors several bright blobs. A clear nebulosity enclosing the blobs extends farther to the east of the north lobe somewhat independently of the spindle shape. The peak surface brightness is measured at 6.2 mag arcsec-2 and 4.4 mag arcsec-2 at the southernmost bright blob at $`L^{}`$ and $`M^{}`$, respectively. The total magnitudes are $`L^{}`$ = 5.5 mag and $`M^{}`$ = 3.5 mag. The central star is invisible in either wavelength, but a narrow dark lane lying 140 east of north in the center of the nebula was clearly detected in our image. The faint nebulosity extending about 6″ east of the north lobe roughly corresponds to the CO molecular outflow observed at 1.3 mm (Cox et al., 2000), suggesting shocked H2 emission falls on the $`L^{}`$ filter coverage.

We used HST/NICMOS archival data of AFGL 2688 obtained with H2 continuum (2.15 µm) and $`H`$-band (1.65 µm) filters and our $`L^{}`$ image to create a three-color composite image. The three images were registered at the brightest blob in the northern lobe. This could be inaccurate because the precise position of the surface brightness peak can be intrinsically different in each color, but we did not find any significant offset in the spindle pattern in the registered three images. The resultant composite image is shown in Figure Imaging and Spatially Resolved Spectroscopy of AFGL 2688 in the Thermal Infrared Region 111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. in different stretches.

### 2.2 3 µm Spectroscopy

The spectroscopic observation was made on UT 2000 September 24 with the same instrument and telescope using the medium-resolution grism in the IRCS camera section. The slit was put along north to south at three different locations on the several bright blobs in the northern lobe. The slit positions are illustrated in Figure Imaging and Spatially Resolved Spectroscopy of AFGL 2688 in the Thermal Infrared Region 111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.. We used the 0$`\stackrel{}{\mathrm{.}}`$3 slit to match the excellent seeing of the night ($``$0$`\stackrel{}{\mathrm{.}}`$25 in FWHM) to obtain low-resolution (R = 600–800) spectra from 2.84 to 4.18 µm. The sky background was subtracted by nodding the telescope by 3″ along the slit length. Special attention was paid to stay on one position for no longer than 60 s to ensure good sky subtraction. The total on-source integration time was 140 s at each slit position. We observed an F5V star (HR 1687, V = 5.0) as a spectroscopic standard at nearly the same airmass to cancel out the telluric absorption lines. The difference in the airmass between the object and the standard star was kept no larger than 0.1 so that the telluric transmission lines canceled well. Spectroscopic flat frames were obtained using a halogen lamp at the end of the night.

After we performed sky subtraction, flat-fielding, and correction of bad pixels in 2D spectrogram images, one-dimensional spectra were obtained using the aperture extraction package of IRAF. The width of the aperture was 0$`\stackrel{}{\mathrm{.}}`$233 (4 pixels), about the same size as the seeing. We defined the nine apertures within the slit offset by 0$`\stackrel{}{\mathrm{.}}`$233 from each other covering $`\pm `$1″ from the center of the slit (Fig. Imaging and Spatially Resolved Spectroscopy of AFGL 2688 in the Thermal Infrared Region 111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.).

Wavelength calibration was performed using over 60 telluric absorption lines. First, we calculated the model telluric transmission curve with matched spectral resolution using the ATRAN atmosphere modeling software (Lord, 1992). Then the local minimums of the model curve were picked up to make a list of telluric lines. The result of fitting the observed atmospheric absorption lines with a calculated line list shows differences no more than $`\pm `$0.001 µm in peak-to-peak all across the spectra. The intrinsic stellar lines of the spectral standard (mostly atomic hydrogen) were hard to remove because of the blending with the telluric absorption lines. Only Br$`\alpha `$ (4.05 µm) and Pf$`\gamma `$ (3.74 µm) were fit with a Lorentzian profile and subtracted before dividing. The one-dimensional spectra were registered to the standard star spectrum by subpixel shifting. The small discrepancy in the airmass between the standard and the object was corrected by rescaling the spectra according to Beer’s law. Since the 0$`\stackrel{}{\mathrm{.}}`$3 slit width corresponds to 5 pixels, the spectra shown in Figure Imaging and Spatially Resolved Spectroscopy of AFGL 2688 in the Thermal Infrared Region 111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. were binned by 2 pixels. The spectrum of the standard star was assumed to be a Planck function of $`T_{\mathrm{eff}}=`$ 6530 K, and the $`L^{}`$ magnitude of the standard star was estimated using the color of an F5V star ($`VL`$ = 1.35) in Tokunaga (2000).

## 3 Discussion

### 3.1 Current Understanding of the Nebula’s Structure

A summary of our present understanding of AFGL 2688 is sketched in Figure Imaging and Spatially Resolved Spectroscopy of AFGL 2688 in the Thermal Infrared Region 111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.. We focused our attention on recent observational studies of the morphology and kinematics of the nebula. A detailed discussion of models, such as the one by Skinner et al. (1997), are beyond the scope of this paper.

In the past decades, this nebula has been the one of the primary targets for new observing facilities and techniques. The central star has not been detected to date at any wavelength (but see Cox et al. 2000, who found a millimeter continuum peak). However, the location of the illumination source in the nebula is determined reasonably well by the center of the “searchlight beams” (Sahai et al., 1998b) and the symmetric center of the polarization map in the reflection nebula (Weintraub et al., 2000). In addition,the polarization center roughly coincides with the intersection of lines connecting blue- and redshifted components of multiple CO outflows (Cox et al., 2000). Low-velocity (10–20 km s-1) winds associated with the mass loss have been detected in the study of molecular line kinematics at radio wavelengths (Young et al., 1992; Yamamura et al., 1995). The visible nebulosity is thought to be a remnant of mass loss when the star was still in the AGB phase.

In the context of the current standard model, the horns extending out of the lobes have been interpreted as limb brightening at the inner wall of the conical opening. However, Sahai et al. (1998b) found the edges of these searchlight beams are too sharp, and the contrast between searchlights and regions between the two searchlight beams are too high, to be reproducible by limb brightening. They argue that the sharp edge is not a product of a material distribution. They propose an obscuring “dust cocoon” with a nonuniform annular opening through which the central source is illuminating the entire nebula to produce the observed intensity distribution of the searchlight beams.

The concentric arcs around the nebula have been taken to be evidence of material ejected in the periodic mass loss caused by either unstable pulsation of the latest AGB phase or recent interaction with a close binary. With faint and detailed structures revealed by HST/WFPC2, Sahai et al. (1998b) found the arcs are almost partial circles with no additional elongation. They argue that the perfect circularity is not consistent with latitudinal density variations because the arcs should proceed faster in the lower-density medium to form more elongated curves along polar axis, which is not apparent in their images. The absence of latitudinal density variation in Sahai’s model is not consistent with the current standard model.

In the $`H`$ and $`K`$ the nebula lobes appears as two spindles capped with shocked H2 emission blobs (Latter et al., 1993). The clear and bright spindle edge indicates bubbles invading the low-velocity wind sweeping up the circumstellar material. The axisymmetrical bubbles are probably infrared counterparts of the moderate-velocity wind ($``$40 km s-1) detected in the radio wavelength (Young et al., 1992). They argued against radiation pressure acceleration to drive such a fast outflow by the dragging of dust grains.

The blue and red components of the radio outflow found in 13CO observations are found along PA = 60 (Yamamura et al., 1995), rather than PA = 15 in the optical. Jura et al. (2000) have discovered the thermal emission tail of extremely large dust grains ($`a`$ 0.5 cm) at centimeter wavelengths. The emitting region is elliptical in shape, with the major axis along PA = 53. An expanding shell-like structure of 2″ diameter was found by Cox et al. (2000) in CO. The shell is slightly elongated to PA = 54. One of the puzzling features of the nebula is that the axis of symmetry of the optical and infrared lobes is at PA = 15, but the axis of symmetry at radio wavelengths is at PA = 53–60.

### 3.2 Sharp Edge in the Southern Lobe

What is recognized first in our three-color composite image (Fig. Imaging and Spatially Resolved Spectroscopy of AFGL 2688 in the Thermal Infrared Region 111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.) is a clear border in the southern lobe with a striking color change in the infrared. The crosscut at the south edge along the optical bipolar axis is shown in Figure Imaging and Spatially Resolved Spectroscopy of AFGL 2688 in the Thermal Infrared Region 111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.. The $`H`$ and 2.15 µm flux drops sharply at 1″ south from the center. The location of the border roughly corresponds to the edge of the dust cocoon proposed by Sahai et al. (1998a). The red color ($`KL^{}>`$ 3) suggests heavy reddening in the intervening dust cocoon. The sharpness of the limb implies an expanding cocoon into the circumstellar envelope pushing away the surrounding material. To allow bulletlike ejection of the bipolar bubbles, the cocoon would have to have formed after the lobes.

We noticed the spatial variation of the flux in Figure Imaging and Spatially Resolved Spectroscopy of AFGL 2688 in the Thermal Infrared Region 111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. in $`L^{}`$ and $`M^{}`$ is almost the same, and the similarity also holds for between $`H`$ and 2.15 µm, but it exhibits a huge difference between 2.15 µm and $`L^{}`$, indicating the specific grain size responsible for the reddening. To estimate the size of the grain and the amount of the extinction, we compare the flux density distribution observed at the cocoon and outside of the cocoon (Fig. Imaging and Spatially Resolved Spectroscopy of AFGL 2688 in the Thermal Infrared Region 111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.). It is evident the visible light is totally obscured in the cocoon, while there is no significant color change attributable to the cocoon in the thermal infrared region, except that the flux density increases in the cocoon by order of 2. This indicates that the light from the central region is first scattered by the dust grains in the lobe, and then subject to extinction in the foreground dust cocoon surrounding the lobes. The flux density enhancement at thermal wavelengths could be attributed to the higher particle density closer to the center.

It is obvious that the extinction curve is quite different from the average interstellar medium. We calculated the extinction efficiency of amorphous carbon grains of various sizes using the optical constants for amorphous carbon retrieved from Zubko et al. (1996). We found the reddening property is best reproduced by 0.2 µm single-size spherical amorphous carbon grains, and it has a flat extinction in the visible wavelength, a large drop between 2 to 3 µm region, and an equally negligible extinction at $`L^{}`$ and $`M^{}`$. We applied the reddening by the amorphous carbon to the flux density distribution observed outside the cocoon with column density of $`N_{\mathrm{dust}}=1.5\times 10^9`$ cm-2. We found the resulting density distribution matches that inside the cocoon well (Fig. Imaging and Spatially Resolved Spectroscopy of AFGL 2688 in the Thermal Infrared Region 111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.). The extinction applied in the cocoon corresponds to $`A_V`$ = 6.3 mag in visible and $`A_K`$ = 2.7 mag in $`K`$, where the optical depth was estimated from $`\tau _\lambda =N_{\mathrm{dust}}Q_{\mathrm{ext}}\pi a^2`$ and $`a`$ is the particle radius.

### 3.3 Dark Lane at the Center of the Nebula

The smaller extinction at $`L^{}`$ allows a view close of the central star through the heavily obscuring dust cocoon. At this wavelength a dark lane in the center of the nebula is evident (Fig. Imaging and Spatially Resolved Spectroscopy of AFGL 2688 in the Thermal Infrared Region 111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.). This dark lane, which is seen at $`L^{}`$ and $`M^{}`$, has a different orientation from the dust cocoon observed by Sahai et al. (1998b), and it appears to be perpendicular to the major axis of the centimeter continuum emission region observed by Jura et al. (2000).

The location of the central star proposed by the polarization map, searchlights, and the multiple outflow of CO molecules falls close to the dark lane, but there is no pointlike source at $`L^{}`$ or $`M^{}`$ corresponding to the central star. We infer the dark lane is a dust disk lying almost edge-on obscuring the central star. The similar appearance at $`L^{}`$ and $`M^{}`$ bands indicates it is absorbing light evenly at 3–5 µm. Dust grains in the disk should be cold enough compared to the surrounding nebulosity to produce no detectable emission at $`L^{}`$ or $`M^{}`$. Hence the size of the grain should be quite large ($`a>`$ 1 µm).

Although the nature of the disk is uncertain, there might be a physical association of the disk orientation to the morphological features observed in the radio wavelength oriented to PA $``$ 55 (Fig. Imaging and Spatially Resolved Spectroscopy of AFGL 2688 in the Thermal Infrared Region 111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.). The disk would collimate all these high-velocity compact components in a direction offset from the optical polar axis. The high velocity compact components should be the youngest structures in the nebula. We speculate that the central star is now experiencing a change in the symmetry axis. Consequently, the observed position angle of each component near the center may changes depending on the time it formed. It is tempting to note the series of blobs in the northern lobe seems to curve from the northwest to the northeast from the center (Fig. Imaging and Spatially Resolved Spectroscopy of AFGL 2688 in the Thermal Infrared Region 111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.). The sequence might be a tracer of explosive ejection from the central star in the order of formation. We might be witnessing the central star gradually changing the ejection axis from PA = 15 to PA $``$ 55.

### 3.4 The Nature of Peak A

An enigmatic object near the center of the nebula is peak A, a pointlike source at the northeast of the central star (Fig. Imaging and Spatially Resolved Spectroscopy of AFGL 2688 in the Thermal Infrared Region 111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.). Based on the distorted polarization pattern near the peak, Weintraub et al. (2000) argued the peak A is a self-luminous infrared source that is a binary companion to the invisible central star. At first glance the extremely red color of peak A (Fig. Imaging and Spatially Resolved Spectroscopy of AFGL 2688 in the Thermal Infrared Region 111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.) appears to be a clear signature of the thermal emission source, yet we find evidence to the contrary. We will discuss several observational features to suggest that the peak A is not self-luminous, but a dense dusty blob reflecting the infrared light from the central star.

First, although peak A is indeed the brightest spot, the overall look of the SED does not deviate from the surrounding reflection nebula. To see this, we compared the SED with those of adjacent regions (Fig. Imaging and Spatially Resolved Spectroscopy of AFGL 2688 in the Thermal Infrared Region 111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.). The infrared (1.65–4.70 µm) SED of peak A and those extracted from other locations are almost identical. If peak A has an illuminating source inside seen through the enshrouding dust in the thermal infrared, it would be hard to account for the identical SED in the infrared as far as 2$`\stackrel{}{\mathrm{.}}`$3 ($``$2800 AU) away from peak A. The direct light should be much redder than the scattered one. The only plausible interpretation of the identical SEDs is the identical radiation mechanism, namely, the relevant region is illuminated by a single central infrared source, and the reflected light becomes redder in the foreground dust cocoon. There is obviously no sign of a pointlike source in the southern lobe where a similar SED is observed, and thus we conclude peak A is not a self-luminous point source either, but a dense dusty blob reflecting the illumination from the central infrared source.

Second, if peak A were indeed self-luminous, what kind of source could it be? Peak A is located $``$600 AU away from the central star in the projected distance, so it cannot be thermal emission of dust heated up by the central star ($`L_{}`$ = 1.8 $`\times `$ 104 $`L_{\mathrm{}}`$) in equilibrium (thermal equilibrium temperature is $`<`$50 K). We may eliminate shock heating, because there is no shocked feature at the peak A in H2 2.122 µm narrowband imaging by Sahai et al. (1998a). Hence the remaining possibility is an enshrouded star. The apparent earlier phase of evolution relative to the central star (no nebulosity association, no hydrogen emission lines) suggests the possible stellar source should be a G- to M-type main-sequence or AGB star. However, to date we have no observational support for the presence of such a late type star there. Spatially resolved spectroscopy in the near-infrared (1.1–2.4 µm) made by Hora & Latter. (1994) covers the peak A location, but there is no significant CO absorption band-head at 2.3 µm in their spectra.

Third, our spectroscopy at peak A shows prominent emission features of hydrocarbon dust. Optically thick thermal emission cores are typically featureless, or at most exhibit molecular absorption lines in the 3.0 µm region due to C2H2 or HCN (Goto et al., 1997; Ridgway, Carbon , & Hall, 1978) that are not apparent in our spectrum. As no AGB stars have been reported to have IEF in the 3 µm region, the hydrocarbon features are a consequence of scattered light originating from the central star.

Fourth, the observed hydrocarbon features are all uniform around peak A. The spectra presented in Figure Imaging and Spatially Resolved Spectroscopy of AFGL 2688 in the Thermal Infrared Region 111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. are normalized to the power-law continuum to see if there is any spatial variation in the hydrocarbon feature at 3.3 and 3.4 µm, and this is shown in Figure Imaging and Spatially Resolved Spectroscopy of AFGL 2688 in the Thermal Infrared Region 111Based on data collected at Subaru Telescope, which is operated by the National Astronomical Observatory of Japan.. For comparison with the spectral feature extracted at each location along the slit, a template spectra is created by combining all the data available. We found all the normalized spectra are very similar, and there is no significant deviation from the template at any location within the observed area. This is again consistent with the view that the IEF emission source is remotely illuminating the dusty blobs, and we are observing the scattered light in the 3 µm region.

When scattered light is dominant at 3 µm, it is surely so at the shorter wavelengths. However, the polarization map at 2.2 µm presented by Weintraub et al. (2000) shows lower polarization close to peak A. We suspect that multiple scattering in a dense dust blob can reproduce the same polarization pattern, namely, less polarization at the brighter location. Polarization imaging in the thermal infrared band, where the multiple scattering effects are smaller, would be more conclusive. We conclude that peak A is probably not a self-luminous source, but is instead one of the dust blobs reflecting infrared illumination from the central source.

## 4 Summary

We presented high-resolution ($``$0$`\stackrel{}{\mathrm{.}}`$3) $`L^{}`$ and $`M^{}`$ images of AFGL 2688. Imaging in the thermal infrared wavelength at high angular resolution proved to be a powerful technique to penetrate deep into the central region of a bipolar proto-planetary nebula with heavy extinction and that is invisible at shorter wavelengths. We conclude the following:

1. A dust cocoon described by Sahai et al. (1998b) appears as a very red region in the south lobe with a dramatic color change between 2 to 3 µm. The amount of extinction in the dust cocoon was estimated to be $`A_V`$ = 6.3 mag, assuming 0.2 µm amorphous carbon grains are responsible for the extinction.

2. We found a dark lane at the core of the nebula at PA = 140. The normal position angle is more consistent with a high-velocity outflows rather than an optical reflection nebula. The central star is invisible at both $`L^{}`$ or $`M^{}`$.

3. Peak A, postulated to be a self-luminous object, has an almost identical SED to the nearby heavily obscured regions, and we conclude that it is not self-luminous.

4. From the 3 µm spectra we find that the IEFs at the bright blobs look uniform, and exhibit no significant spatial variation in the feature-to-continuum ratio. This supports the idea that the bright blobs are reflecting light from the invisible central star.

We are grateful to M. Weber and R. Potter for their assistance when these observations were made. We appreciate many useful discussions with H. Suto. We thank the staff and crew of the Subaru Telescope for their valuable assistance in obtaining these data. M. Goto is supported by Japan Society for the Promotion of Science fellowship.


format: text

text:
# 1 INTRODUCTION

## 1 INTRODUCTION

The origin of galaxies and of the observed structure of their inhomogeneous distribution in FRW cosmology is generally ascribed to the development of gravitational instability on the matter - dominant stage, leading to the slow growth of small initial density fluctuations on the nearly homogeneous expanding background. According to modern cosmology the spectrum of initial density fluctuations is generated on the inflationary stage. Particle theory offers a physical basis for both inflation and dark matter, dominating in the development of gravitational instability. However, the impact of particle physics in the problem of galaxy formation seems to be much stronger, adding new non-trivial elements to the theory of cosmological structure formation. One such new element, which is a rather general cosmological consequence of the particle symmetry breaking pattern, is the structure of strong primordial large scale inhomogeneities, tracing the large scale structure of cosmological defects, generated in second order phase transitions in the inflationary Universe.

The existence of topological defects reflects the non-trivial topological properties of broken particle symmetry. So, a magnetic monopole solution follows from compact Lie group breaking, leaving unbroken the U(1) symmetry of electromagnetism, which is the general case for GUT models. In some of these, the conditions for cosmic string and domain wall solutions are realized. The cosmological effect of these topologically stable solutions is widely recognized as an effective tool to probe the scenarios of very early Universe and the particle models underlying them. Recall that magnetic monopole overproduction  has put severe constraints on GUT models and stimulated the creation of inflationary cosmology , whereas the domain wall problem  constrained the models of spontaneous CP violation . On the other hand, topological defects represent a new element of galaxy formation, as it takes place in cosmic string cosmology .

However, the actual amount of primordial magnetic monopoles produced in the inflationary Universe is still a big question, whereas the conditions for stable wall and string solutions are not very general for particle models (see the review ). A much wider class of these models possesses the symmetry breaking pattern, which can be effectively described by pseudo-Nambu–Goldstone (PNG) field and which corresponds to the formation of an unstable topological defect structure in the early Universe.

The Nambu–Goldstone nature in such an effective description reflects the spontaneous breaking of global symmetry, resulting in continuous degeneracy of vacua. The explicit symmetry breaking at smaller energy scale changes this continuous degeneracy by discrete vacuum degeneracy. At high temperatures such a symmetry breaking pattern implies the succession of second order phase transitions. In the first transition, continuous degeneracy of vacua leads, at scales exceeding the correlation length, to the formation of topological defects in the form of a string network; in the second phase transition, continuous transitions in space between degenerated vacua form the surfaces: domain walls surrounded by strings. This last structure is unstable, but, as was shown in the example of the invisible axion , it is reflected in the large scale inhomogeneity of distribution of energy density of coherent PNG (axion) field oscillations. The role of inflation in the creation of this inhomogeneity is indirect, since it provides identical initial conditions of expansion in causally disconnected regions, which provides the simultaneous high temperature phase transition and the formation of structure of topological defects (axionic strings) that spreads far beyond the cosmological horizon in that period.

Second order phase transitions on the inflational stage result in more non-trivial influence of superluminal expansion on the forms of topological defects and properties of their structure. The result depends on the relationship between the symmetry breaking scales and Hubble constant on inflational stage; it leads to various observational consequences, related to various strong primordial inhomogeneities. Among these consequences, systematic study of which we approach, there appears a new interesting way of formation of primordial galactic nuclei, to be discussed in the present paper.

Now there is no doubt that the centres of almost all galaxies contain massive black holes . An original explanation of the formation of such supermassive black holes assumes the collapse of a large number of stars in the galaxy centres. However, the mechanism of the galactic nuclei formation is still unclear. According to , there are serious grounds to believe that the formation of stars and galaxies proceeded simultaneously. On the other hand, in the work  a model of galaxy formation around a massive black hole was considered and arguments in its favour were presented. Each of the two approaches has certain advantages, while being neither free of drawbacks.

We will consider below a new mechanism describing the very early formation of primordial black holes (PBHs), which serve as the nucleation centres in the subsequent formation of galaxies. This mechanism may prove to be free from disadvantages inherent in the models based on the concept of a single PBH being the nucleus of the future galaxy.

Previously  we proposed a mechanism of the PBHs formation that opens the possibility of massive black hole formation in the early Universe. The mechanism is based on the fact that black holes can be created as a result of a collapse of closed vacuum walls formed during a second order phase transition. The masses of such black holes may vary within broad limits, up to the level $`10^8M_{}`$.

To figure out the idea , let us assume that a potential of some scalar field possesses at least two different vacuum states. In such a situation there are two quantitatively different possibilities to distribute these states in the early Universe. The first possibility implies that the Universe contains approximately equal amounts of both vacuum states, which is typically the case at the usual thermal phase transition. The alternative possibility corresponds to the case when the two vacuum states are populated with different probabilities. In this case, islands of a less probable vacuum state surrounded by a sea of another, more probable vacuum state appear. As was shown in , an important condition for such an asymmetric distribution is the existence of effectively flat valleys in the scalar field potential during inflation. Under this condition the scalar field can be considered as a massless scalar field additional to inflaton, existing on the inflational de Sitter background. It is well known that the inflational fluctuations can essentially change the value of a massless scalar field, while the inflation itself blows up wavelengths of these fluctuations. These two factors define the space distribution of a scalar field, which is still massless until a certain moment. Thus, even though the phase transition itself takes place only after the end of inflation, deeply in the Friedman–Robertson–Walker (FRW) epoch, the initial distribution of the scalar field is already defined at the inflation period, in such a way that there are eventually, islands representing one vacuum in the sea of another vacuum.

It is supposed that after the phase transition, the two vacuum states are separated by a vacuum wall. The initial distribution of the scalar field formed during the inflation stage allows the formation of closed walls of sizes significantly exceeding the cosmological horizon at the moment when the phase transition just takes place. At some instant after crossing the horizon, such walls become causally connected as a whole and begin to contract because of the surface tension. As a result, provided that friction is small and the wall does not radiate a considerable part of its energy in the form of scalar waves, almost all the energy of a closed wall may be focused within a small volume inside the gravitational radius. This is the necessary condition for a black hole formation. The mass spectrum of black holes formed by this mechanism depends on parameters of the scalar field potential determining the direction and size of the potential valley during inflation and the post-inflation phase transition. Although we deal here with the so-called PNG field, the proposed mechanism is quite general. The presence of massive PBHs is a new factor in the development of gravitational instability in the surrounding matter and may serve as a basis for new scenarios of the formation and evolution of galaxies.

## 2 PRIMORDIAL BLACK HOLE FORMATION

Now we will describe a mechanism accounting for the appearance of massive walls of a size essentially greater than the horizon at the end of inflation. Let us consider a complex scalar field with the potential

$$V(\phi )=\lambda (\left|\phi \right|^2f^2/2)^2+\delta V(\theta ),$$
(1)
where $`\phi =re^{i\theta }`$. This field coexists with an inflaton field which drives the Hubble constant $`H`$ during the inflational stage. The term

$$\delta V(\theta )=\mathrm{\Lambda }^4\left(1\mathrm{cos}\theta \right),$$
(2)
reflecting the contribution of effects to the Lagrangian renormalization (see for example ), is negligible on the inflational stage and during some period in the FRW expansion. In other words, the parameter $`\mathrm{\Lambda }`$ vanishes with respect to $`H`$. The omitted term (2) begins to play a significant role only at the moment, after inflation, when the Hubble parameter sharply decreases with time ($`H=1/2t`$ during the radiation dominated epoch). Also, we assume the mass of the radial field component $`r`$ always to be sufficiently large with respect to $`H`$, which means that the complex field is in the ground state even before the end of inflation. Since the term (2) is negligible during inflation, the field has the form $`\phi f/\sqrt{2}e^{i\theta }`$, the quantity $`f\theta `$ acquiring the meaning of a massless field.

At the same time, the well established behaviour of quantum field fluctuations on the de Sitter background  implies that the wavelength of a vacuum fluctuation of every scalar field grows exponentially, having a fixed amplitude. Namely, when the wavelength of a particular fluctuation, in the inflating Universe, becomes greater than $`H^1`$, the average amplitude of this fluctuation freezes out at some non-zero value because of the large friction term in the equation of motion of the scalar field, whereas its wavelength grows exponentially. Such a frozen fluctuation is equivalent to the appearance of a classical field that does not vanish after averaging over macroscopic space intervals. Because the vacuum must contain fluctuations of every wavelength, inflation leads to the creation of more and more new regions containing a classical field of different amplitudes with scale greater than $`H^1`$. In the case of an effectively massless Nambu–Goldstone field considered here, the averaged amplitude of phase fluctuations generated during each e-fold (time interval $`H^1`$) is given by

$$\delta \theta =H/2\pi f.$$
(3)
Let us assume that the part of the Universe observed inside the contemporary horizon $`H_0^1=3000h^1`$Mpc was inflating, over $`N_U60`$ e-folds, out of a single causally connected domain of size $`H^1`$, which contains some average value of phase $`\theta _0`$ over it. When inflation begins in this region, after one e-fold, the volume of the Universe increases by a factor $`e^3`$ . The typical wavelength of the fluctuation $`\delta \theta `$ generated during every e-fold is equal to $`H^1`$. Thus, the whole domain $`H^1`$, containing $`\theta _0`$, after the first e-fold effectively becomes divided into $`e^3`$ separate, causally disconnected domains of size $`H^1`$. Each domain contains almost homogeneous phase value $`\theta _0\pm \delta \theta `$. Thereby, more and more domains appear with time, in which the phase differs significantly from the initial value $`\theta _0`$. A principally important point is the appearance of domains with phase $`\theta >\pi `$. Appearing only after a certain period of time during which the Universe exhibited exponential expansion, these domains turn out to be surrounded by a space with phase $`\theta <\pi `$. The coexistence of domains with phases $`\theta <\pi `$ and $`\theta >\pi `$ leads, in the following, to the formation of a large-scale structure of topological defects.

The potential (1) possesses a $`U(1)`$ symmetry, which is spontaneously broken, at least, after some period of inflation. Note that the phase fluctuations during the first e-folds may, generally speaking, transform eventually into fluctuations of the cosmic microwave radiation, which will lead to imposing restrictions on the scaling parameter $`f`$. This difficulty can be avoided by taking into account the interaction of the field $`\phi `$ with the inflaton field (i.e. by making parameter $`f`$ a variable ). This spontaneous breakdown is holding by the condition of smallness of the radial mass, $`m_r=\sqrt{\lambda _\varphi }>H`$. At the same time the condition

$$m_\theta =\frac{2f}{\mathrm{\Lambda }}^2H$$
(4)
on the angular mass provides the freezing out of the phase distribution until some moment of the FRW epoch. After the violation of condition (4) the term (2) contributes significantly to the potential (1) and explicitly breaks the continuous symmetry along the angular direction. Thus, potential (1) eventually has a number of discrete degenerate minima in the angular direction at the points $`\theta _{min}=0,\pm 2\pi ,\pm 4\pi ,\mathrm{}`$ .

As soon as the angular mass $`m_\theta `$ is of the order of the Hubble rate, the phase starts oscillating about the potential minimum, initial values being different in various space domains. Moreover, in the domains with the initial phase $`\pi <\theta <2\pi `$, the oscillations proceed around the potential minimum at $`\theta _{min}=2\pi `$, whereas the phase in the surrounding space tends to a minimum at the point $`\theta _{min}=0`$. Upon ceasing of the decaying phase oscillations, the system contains domains characterized by the phase $`\theta _{min}=2\pi `$ surrounded by space with $`\theta _{min}=0`$. Apparently, on moving in any direction from inside to outside of the domain, we will unavoidably pass through a point where $`\theta =\pi `$ because the phase varies continuously. This implies that a closed surface characterized by the phase $`\theta _{wall}=\pi `$ must exist. The size of this surface depends on the moment of domain formation in the inflation period, while the shape of the surface may be arbitrary. The principal point for the subsequent considerations is that the surface is closed. After reheating of the Universe, the evolution of domains with the phase $`\theta >\pi `$ proceeds on the background of the Friedman expansion and is described by the relativistic equation of state. When the temperature falls down to $`T_{}\mathrm{\Lambda }`$, an equilibrium state between the ”vacuum” phase $`\theta _{vac}=2\pi `$ inside the domain and the $`\theta _{vac}=0`$ phase outside it is established. Since the equation of motion corresponding to potential (2) admits a kink-like solution (see  and references therein), which interpolates between two adjacent vacua $`\theta _{vac}=0`$ and $`\theta _{vac}=2\pi `$, a closed wall corresponding to the transition region at $`\theta =\pi `$ is formed. The surface energy density of a wall of width $`1/mf/\mathrm{\Lambda }^2`$ is of the order of $`f\mathrm{\Lambda }^2`$ 111The existence of such domain walls in theory of the invisible axion was first pointed out in ..

Note that if the coherent phase oscillations do not decay for a long time, their energy density can play the role of CDM. This is the case, for example, in the cosmology of the invisible axion (see  and references therein).

It is clear that immediately after the end of inflation, the size of domains which contains a phase $`\theta _{vac}>2\pi `$ essentially exceeds the horizon size. This situation is replicated in the size distribution of vacuum walls, which appear at the temperature $`T_{}\mathrm{\Lambda }`$ whence the angular mass $`m_\theta `$ starts to build up. Those walls, which are larger than the cosmological horizon, still follow the general FRW expansion until the moment when they get causally connected as a whole; this happens as soon as the size of a wall becomes equal to the horizon size $`R_h`$. Evidently, internal stresses developed in the wall after crossing the horizon initiate processes tending to minimize the wall surface. This implies that the wall tends, first, to acquire a spherical shape and, second, to contract toward the centre. For simplicity, we will consider below the motion of closed spherical walls 222The motion of closed vacuum walls has been driven analytically in ..

The wall energy is proportional to its area at the instant of crossing the horizon. At the moment of maximum contraction, this energy is almost completely converted into kinetic energy. Should the wall at the same moment be localized within the gravitational radius, a PBH is formed.

Detailed consideration of BH formation was performed in . We proceed below to study the formation of a PBH cluster in the early Universe.

## 3 CORRELATIONS IN THE BLACK HOLE
DISTRIBUTION

Previously  we had studied a new process involving the formation of uncorrelated PBHs in the Universe. It was demonstrated that a model with reasonable parameters readily provides for the formation of $`10^{11}`$ massive (with the mass $`10^{30}`$-$`10^{40}`$g. each) black holes, which coincides with the number of galaxies in the observed part of the Universe. In that analysis, we did not take into account correlations (inherent in this mechanism) between the formation of a massive black hole and the appearance of smaller black holes surrounding it. This correlation is related primarily to certain features of the above-discussed process of the formation of domains with phases $`\theta >\pi `$. It seems that the appearance of such domains creates prerequisites for the formation of new, smaller domains inside. Let us estimate the mass distribution of these daughter domains. Consider a region with a size of the order of $`H^1`$ and a phase within $`\pi <\theta _0<\pi +\delta `$ (where $`\delta =H/2\pi f`$ is the average phase jump during the $`H^1`$ time period), formed during the inflation period as a result of fluctuation in a certain region of space of phase $`\theta <\pi `$. During the next e-fold, this space domain will be divided in $`e^3`$ subdomains with size $`H^1`$, and some of these will acquire a phase $`\theta _1`$ in the interval $`\pi \delta <\theta _1<\pi `$. Upon the subsequent phase transition, these domains will be separated by walls from the external region. Similar transitions, when crossing the phase $`\theta =\pi `$ in the reverse direction, will take place in each subdomain during the next e-fold. Thus, a structure of the fractal type appears 333A fractal structure of unclosed axionic walls, which is generated on the inflationary stage has been discussed in . which reproduces itself in each time step on a decreasing scale.

Let $`\zeta `$ denotes the number of subdomains formed in each step, around which a wall may appear with time. Apparently, this value obeys the inequality $`1<\zeta e^3`$. In the subsequent estimates, we will assume that $`\zeta 2`$$`3`$. Since each causally connected domain touches approximately six neighboring domains, we can hardly expect $`\zeta `$ to be greater for a total number of $`e^320`$. The mass of the future black hole (if this would actually form) is determined by the area of a closed surface with the phase $`\theta =\pi `$. The ratio of areas of the initial (mother) and daughter domains is readily estimated: the initial area, after a single e-fold, is $`S_0e^2H^2`$, and the daughter subdomain area is $`S_1H^2`$. Therefore, the ratio of masses of the black holes belonging to two sequential generations is

$$M_j/M_{j+1}S_j/S_{j+1}e^2,$$
(5)
for a relative number of them assumed to be

$$N_{j+1}/N_j=\zeta .$$
(6)
As is readily seen, the number and mass of black holes appearing upon the $`j`$-th e-fold after the initial domain formation are determined by parameters of the largest black hole genetically related to the primary domain in which the phase originally exceeded $`\pi `$. It is evident that

$$N_j\zeta ^j;M_jM_0/e^{2j}.$$
(7)
Excluding $`j`$ from these relationships, we obtain the desired black hole mass distribution in a cluster:

$$N_{cl}(M)\left(M_0/M\right)^{\frac{1}{2}\mathrm{ln}\zeta }.$$
(8)
The total mass of the cluster can be expressed through the mass $`M_0`$ of the largest initial black hole:

$$\begin{array}{c}M_{tot}=M_0+\zeta M_1+\zeta ^2M_2+\mathrm{}=\hfill \\ M_0+\zeta e^2M_0+\left(\zeta e^2\right)^2M_0+\mathrm{}=M_0[1\zeta /e^2]^1.\hfill \end{array}$$
(9)
As is seen, the total mass of the black hole cluster is only one and a half to two times greater than the largest initial black hole mass. The number of daughter black holes depends on the factors considered in the next section.

The inflationary mechanism described above leads to the occurrence of a fractal structure of the closed walls. After the end of inflation, as soon as the size of the horizon becomes larger than the characteristic size of closed walls, the walls begin to shrink. The energy of each wall, proportional to the area of their surface, concentrates in small spatial domains (in the following they are considered as point-like objects). These high density clots of energy could serve in the following for star and/or galaxy formation (see for example  and references there). Hence, according to the given models, the distribution of stars and galaxies should have fractal properties as well. It is important to note that the total surface of the walls in definite volume is proportional to the total energy of an object, while the number of walls is equal to the number of dense clots.

Suppose, for a characteristic time $`1/H`$, that several closed walls appear in a causally connected area of size $`R`$, and the phase is still at the maximum of the potential. Denote the number of walls by $`N`$ and their average size by $`\xi R`$, $`\xi >1/e`$ ($`\xi 1/e`$ because of a possible merging of causally disconnected subdomains with one common wall). In each of these subdomains, $`N`$ new, smaller closed walls of size $`\xi ^2R`$ arise during the next time step. Denote by $`a`$ the minimal size of such a wall that we are able to distinguish. This means that we may terminate the process after the $`n`$-th step such that $`a\xi ^nR`$. The total area of the closed walls in the initial volume is the sum of the areas with closed walls of size greater than $`a`$. This simple summation leads to the following result:

$$SR^2q(q^n1)/(q1),q\xi ^2N.$$
(10)
This expression can be written in the form

$$S(R/a)^D,$$
(11)
where $`D`$ is the fractal dimension. Equating these two expressions, one obtains

$$D=2+\frac{\mathrm{ln}\left(q\frac{q^{\frac{\mathrm{ln}\left(a/R\right)}{\mathrm{ln}\xi }}1}{q1}\right)}{\mathrm{ln}(R/a)}.$$
(12)
This quantity is constant only when the ratio $`R/a`$ is large; it is different for $`q<1`$ and $`q>1`$. It can be easily verified that $`D2`$ for $`q<1`$, while for $`q>1`$, $`D2+3\mathrm{ln}(q)/\mathrm{ln}(4N)`$. To get an estimate, suppose that the number of closed domains is $`N4`$, and $`\xi 1/e`$. The value of the parameter $`q`$ can be easily calculated, $`q0.5`$. Hence, the fractal dimension of the system of closed walls $`D2`$.

So, if quantum fluctuations lead to the formation of spatial areas, with the phase taking a value near a potential maximum, its further evolution results in a system of enclosing walls. The characteristic size of the next generations of walls differs from the previous one by a factor of approximately $`e`$. The fractal dimension of such a system is $`D2`$.

According to this scenario, it is interesting to find the number of walls inside a sphere of radius $`R`$, which is given by

$$N_{tot}=\underset{i=1}{\overset{n}{}}N^i=N\frac{N^n1}{N1}\frac{N^{n+1}}{N1}.$$
(13)
By analogy with the previous calculations and using Eq. (13), one obtains the distribution of point–like dense objects with fractal dimension $`D^{}\mathrm{ln}N/\mathrm{ln}(1/\xi )`$. For the realistic values $`N4and\xi 1/e`$, we find $`D^{}1.4`$, which differs somewhat from the value $`D2`$ previously obtained. This is not surprising since in the first case we measure the area of wall surfaces within a certain volume, while in the second case we measure the number of walls.

Let us compare our calculations with observational data of spatial distribution of galaxies and of stars in those galaxies. Recent data indicate that the distribution of stars and galaxies indeed carries a fractal character. So the number of galaxies inside a sphere of radius $`R`$ is $`N(R)(R)^{2.2\pm 0.2}`$ up to a sizes of 200 Mpc .

The distribution of stars inside galaxies also carries a fractal character. In Ref. this fractal dimension was determined by averaging observational data from ten galaxies; it was found to be equal to $`D2.3`$.

Evidently, the observable fractal dimension $`D`$ in distributions of stars and galaxies is in agreement with predictions of our model. Of course, other mechanisms, at a later stage, may contribute to the distribution and change the fractal dimension somewhat, but the scenario of PBHs formation, discussed here, could give a primordial reason of fractality in the galaxy and star distribution.

The mechanism of fractal structure production discussed in this paper is not unique. Another example is based on hybrid inflation, one of the most promising models of inflation . In the standard version of hybrid inflation the potential contains two fields

$$V=V_0+\frac{1}{2}m_\phi ^2\phi ^2+\frac{1}{2}\lambda _1\phi ^2\psi ^2\frac{1}{2}m_\psi ^2\psi ^2+\frac{1}{2}\lambda _2\psi ^4.$$
(14)
During inflation, the field $`\phi `$ rolls down along a valley $`\psi =0`$. In the meantime field fluctuations around the critical line $`\psi =0`$ lead to the formation of a fractal structure of domains. This critical line plays the same role as the critical point $`\pi `$ in the previous discussion. Just after passing the critical point $`\phi =m_\psi ^2/\lambda _1`$, the state $`\psi =0`$ becomes unstable and field $`\psi `$ moves (on average) to one of the new stable minima. These minima are separated by a potential maximum and we again inevitably come to the fractal structure of domain walls.

## 4 DISCUSSION

In the preceding sections, we considered only the principal possibility of the formation of domain walls connecting adjacent vacuum states. We have used the formulas derived above to estimate the efficiency of the proposed mechanism of the black hole cluster formation. The numerical calculations  were performed for the following values of the parameters (which are consistent with the observed anisotropy in the cosmic microwave radiation): the Hubble constant at the end of inflation, $`H=10^{13}`$Gev; Lagrangian parameters, $`f=1.77H`$ and $`\mathrm{\Lambda }=5`$ GeV. The initial phase, at which the visible part of the Universe is formed between the time $`t_U60H^1`$ to the end of inflation, controls the number of domains and, accordingly, the number of closed walls formed in the post-inflation stage. This random value, not related to the Lagrangian parameters, must be selected, taking into account the numerous restrictions on the abundance of PBHs in the Universe (see for example  and references there). We will use the numerical value $`\theta _U=0.05\pi `$, which ensures a sufficiently large number of massive black holes, while the presence of numerous smaller ones does not contradict experimental restrictions. As is seen in Fig.1, the PBH masses fall within the range from $`10^{25}`$ to $`10^{35}`$g. The initial phase $`\theta _U`$ was selected so as to provide for the number of massive PBHs (with mass $`10^{35}g`$) to be equal to the number of galaxies in the visible part of the Universe. The total mass of black holes is supposed to be $`1\%`$ of the contemporary baryonic contribution.

The results of these calculations are sensitive to changes in the parameter $`\mathrm{\Lambda }`$ and the initial phase $`\theta _U`$. As the $`\mathrm{\Lambda }`$ value decreases to $`1`$GeV, still greater PBHs appear with masses of up to $`10^{40}`$ g. A change in the initial phase leads to sharp variations in the total number of black holes. As was shown above, each domain generates a family of subdomains in the close vicinity. The total mass of such a cluster is only 1.5–2 times that of the largest initial black hole in this space region. Thus, our calculations confirm the possibility of formation of clusters of massive PBHs ( $`100M_{}`$ and above) in the earliest stages of the evolution of the Universe at a temperature of $`110`$GeV. These clusters represent stable energy density fluctuations around which increased baryonic (and cold dark matter) density may concentrate in the subsequent stages, followed by the evolution into galaxies.

It should be noted that additional energy density is supplied by closed walls of small sizes. Indeed, because the smallness of their gravitational radius, they do not collapse into BHs. After several oscillations such walls disappear, leaving coherent fluctuations of the PNG field. These fluctuations contribute to a local energy density excess, thus facilitating the formation of galaxies.

## 5 CONCLUSIONS

This paper proposes a new mechanism for the formation of protogalaxies, which is based on the cosmological inferences of the elementary particle models predicting non-equilibrium second order phase transition  in the inflation stage period and the formation of a domain wall upon the end of inflation. The presence of closed domain walls of a size markedly exceeding the cosmological horizon in their period of formation leads to the collapse of the wall in the post-inflation epoch (when the wall size becomes comparable with the cosmological horizon); this results of a formation of massive black hole clusters that can serve as nuclei for the future galaxies. The mass spectrum of PBHs do not contradict the available restrictions. The number of black holes with $`M100M_{}`$ and above is comparable with the number of galaxies in the visible Universe. Processes of deceleration of the wall motion, which have been considered in , may affect only the dynamics of collapse of supermassive walls, which are beyond of the scope of this paper.

A development of the proposed approach gives ground for a principally new scenario of the galaxy formation in the model of the hot Universe. Traditionally, the hot Universe model assumes a homogeneous distribution of matter on all scales, whereas the appearance of observed inhomogeneities is related to the growth of small initial density perturbations. However, the analysis of the cosmological consequences of the particle theory indicates the possible existence of strongly inhomogeneous primordial structures in the distribution of both the dark matter and baryons. These primordial structures represent a new factor in galaxy formation theory. Topological defects such as the cosmological walls and filaments, primordial black holes, archioles in the models of axionic CDM , and essentially inhomogeneous baryosynthesis (leading to the formation of antimatter domains in the baryon-asymmetric Universe)  offer by no means a complete list of possible primary inhomogeneities inferred from the existing elementary particle models.

The proposed approach discloses a number of interesting aspects in this direction. Indeed, our model offers a possibility of quantitative analysis of correlations in the formation of massive PBHs and the primary inhomogeneity of dark matter and baryons. Originally inherent in this mechanism is the inhomogeneous phase distribution, which eventually acquires (much as what takes place in the invisible axion cosmology) a dynamical sense of the initial amplitude of the coherent oscillations of a scalar field. Irrespective of the efficiency of dissipation of the energy of these oscillations, the regions of closed wall formation must be correlated with the regions of maximum energy density of the latent mass. If these oscillations are not decaying, their energy density may provide for the contemporary CDM density. Inhomogeneity in the initial amplitude of these oscillations would then imply the inhomogeneity in the initial energy density and, hence, the regions of black hole formation would become the regions of increased dark matter density.

A qualitatively similar effect (albeit not as pronounced) takes place in the dissipation of coherent oscillations at the expense of particle production. An increase in the oscillation energy density transforms into a local increase in the density of latent mass particles produced in this region.

The development of the proposed approach may thus lead to a number of interesting scenarios of the initial stages in the formation of protogalaxies, depending on the selection of particular elementary particle models and their parameters. This study presents the first step in this direction.

## 6 ACKNOWLEDGMENTS

The authors are grateful to E.D. Zhizhin and A.A. Starobinsky for their interest in this study. The work of M.Yu.Kh. and S.G.R. was partly performed within the framework of the ”Cosmoparticle Physics” Project of the State Scientific - Technological Program ”Astronomy. Fundamental Space Research”, and it was supported by the Cosmion-ETHZ, Epcos-AMS international cooperations and Scientific School Support Grant no. 00-15-96699. One of us (M.Yu.Kh.) expresses his gratitude to IHES and ICTP for hospitality, while finishing the work on this paper. We are also grateful to J. Ellis for useful comments.


text:
# New summing algorithm using ensemble computing

## I Introduction

In this paper, we propose an ensemble algorithm for evaluating and summing an arbitrary function, as an alternative to the quantum algorithms that are currently believed to be the most efficient algorithms available, in terms of query complexity.

The possibility of speeding up the evaluation and summing of a large number of function samples was first noted by Abrams and Williams Abrams99 , who suggested calculating numerical integrals and stochastic processes using quantum algorithms. Quantum algorithms exploit the inherent parallelism offered by entangled quantum states, to perform certain computational tasks much more efficiently than classical devices, using either pure or pseudopure quantum states DiVincenzo00 ; Vandersypen01  that are tensor products of multiple qubits. The numerical value of an integral is evaluated by employing either the mean estimation algorithm devised by Grover Grover97b  to calculate the mean of a discrete set of numbers, or by using the quantum counting algorithm proposed by Brassard, Hoyer, and Tapp Brassard98  to determine the number of elements that fulfill a specified condition. Both of these approaches rely on a generalization of Grover’s search algorithm, resulting in a quadratic speedup in comparison with classical randomized (Monte Carlo) algorithms, and an exponential speedup in comparison with classical deterministic algorithms for a single processor.

A systematic comparison of optimal summation of finite sequences and continuous-function integration for deterministic, randomized, and quantum algorithms has been done by Heinrich and Novak Heinrich2001a ; Heinrich2001b . They have examined the query complexity of quantum integration for different classes of integrand functions, assuming that the critical quantum speedup is obtained by using one of the two quantum summing algorithms mentioned above.

Recently an alternative paradigm for computing has been suggested by Madi, Bruschweiler, and Ernst, which operates on ensembles i.e., mixed states of identical spin sytems, using a spin Liouville space formalism Madi98 ; Bruschweiler00 . These types of ensemble algorithms are not quantum algorithms, insofar they do not involve entanglement of quantum states. Throughout this paper, ”mixed states” describe a statistical ensemble, not individual systems, so that each element of the ensemble performs part of the computation, in the same way as a classical parallel computer.

This new paradigm exploits the parallelism offered by simultaneously acting on linear combinations of many different input states in an ensemble of spins. Thus in general, ensemble computing requires an exponentially larger set of memory resources to encode the same number of distinct input states compared to quantum computing with pure states.

While ensemble computing requires more physical resources, it holds the important advantage that it is insensitive to the decoherence time of the spins, which is an outstanding limiting factor for quantum computations involving entangled states. Moreover, ensemble algorithms can be exponentially faster than the equivalent quantum algorithms, for adequate measurement sensitivities, so a trade-off exists between memory and speed capabilities.

We proceed to present a new approach to summing up function samples using an ensemble algorithm, and discuss its query complexity. We only consider the query complexity of the algorithm, i.e., the number of function evaluations that are performed, since the overall computational complexity will depend on the actual function that is being evaluated. At present, the most feasible physical implementation of this summing algorithm would rely on NMR technology, though any physical system of spins can be used in principle. The ideal physical system would allow us to have full control over a very large number of spins, in order to satisfy the large memory requirements. In the Discussion section, we comment on the application of the proposed summing algorithm to evaluating the mean of a continuous function, and as a corollary, on estimating the definite integral of a continuous multi-dimensional function.

## II Ensemble summing algorithm

### II.1 Statement of the Problem

Let $`f:\{1,2,\mathrm{},N\}[0,1]`$ be a real-valued function defined on a discrete set of samples comprised of $`N=2^n`$ points. The function $`f`$ may be known analytically or it may be the result of an explicit or hidden numerical computation. The latter case is known as an oracle. We want to evaluate efficiently the sum $`S_N`$,

$$S_N=\underset{i=1}{\overset{N}{}}f(i).$$
(1)
Here efficiency is understood in relation to the query complexity of the algorithm. Indeed, when $`N`$ is large and the function evaluation is costly in terms of computational complexity, reducing the number of function evaluations is critical.

We assume that the algorithm is to be implemented in a physically realizable system consisting of a finite number of two-valued spins. To accomodate the $`N`$ input values, we need $`n`$ spins in the input register. The finiteness of the system and the discreteness of the spin states implies that we have to approximate the set of function values, $`\{f(i)\}`$, with a set of finite-precision values $`\{f_i[0,1]\}`$ for $`i=1,2,\mathrm{},N`$. For the sake of simplicity, we shall consistently use the same notation for $`f(i)`$ and its $`k`$-digit approximation. Since the meaning of the formulas is explained in the text, there should be no ambiguity. The number of spins $`k`$ available in the output register will specify the minimal precision, $`\delta =2^k`$, for these values. Therefore we are actually evaluating the sum $`S_{N,k}`$,

$$S_{N,k}=\underset{i=1}{\overset{N}{}}f_i,$$
(2)
which converges exponentially fast to the sum $`S_N`$, as we increase the number of spins $`k`$ in the output register. Thus if we can evaluate $`S_{N,k}`$ efficiently, we can also evaluate $`S_N`$ efficiently.

### II.2 Outline of the Algorithm

The proposed algorithm has three main steps. The first step consists of preparing an ensemble mixture of input states representing the numbers $`i=1,2,\mathrm{},N`$. In the second step, the function $`f`$ is applied to the input states, using a single transformation $`U_f`$ to perform the function evaluation for every input state $`i`$ at once. This parallel application results in an ensemble mixture which contains all of the values $`f_i`$ in the output register. Finally, measurement of the output register automatically averages the contributions from the entire ensemble, yielding a signal proportional to the approximate sum, $`S_{N,k}`$.

### Step 1 - Initialization

We initialize the $`n`$-spin input register in an equally-weighted mixed state $`\rho _{in}^{(n)}`$,

$$\rho _{in}^{(n)}=\frac{1}{N}\underset{i=1}{\overset{N}{}}|i>_n<i|_n,$$
(3)
which accounts for all $`N=2^n`$ possible states. The mixed state $`\rho _{in}^{(n)}`$ is a density operator, which can be represented in spin Liouville space by a density matrix that has non-zero elements only on its diagonal. The off-diagonal elements are all zero, indicating the absence of quantum coherence between any of the states $`|i>_n`$.

We note again that this mixed state describes a statistical ensemble, thus the initialization is equivalent to assigning each of the input values $`i=1,2,\mathrm{},N`$, to one of $`N`$ classical processors in a parallel computer.

For example, in an NMR implementation, the states $`|i>_n`$ correspond to the eigenstates of the Zeeman Hamiltonian created by a strong external magnetic field Ernst87 . The ket states, $`|i>_n`$, can also be written in terms of individual spins,

$$|i>_n=|a_{i1}>|a_{i2}>\mathrm{}|a_{in}>$$
(4)
where $`(a_{i1},a_{i2},\mathrm{},a_{in})\{0,1\}`$ are the digits of the number $`(i1)`$ in binary format. The bra states, $`<i|_n`$, are the dual of the ket states. The state $`|0>`$ denotes a spin “up” and the state $`|1>`$ denotes a spin “down”. At room temperature, the thermal equilibrium state of an $`n`$-spin ensemble in an NMR experiment closely approximates the desired initial state $`\rho _{in}^{(n)}`$ Gershenfeld97 . The thermal state is equal to the sum of $`\rho _{in}^{(n)}`$ and a traceless deviation density matrix, with zero off-diagonal terms. The error introduced by using a thermal state instead of the equally-weighted mixed state $`\rho _{in}^{(n)}`$ is addressed in the Discussion.

We also assume that we have available an output register with $`k`$ spins, which is capable of encoding the real-numbered values of the series $`f_i[0,1]`$ with precision $`\delta =2^k`$. All of the states of the output register are initially set to zero, so the state of the entire ensemble (input and output registers) is given by $`\rho _{in}^{(n)}\rho _{out}^{(k)}`$,

$$\rho _{in}^{(n)}\rho _{out}^{(k)}=\frac{1}{N}\underset{i=1}{\overset{N}{}}|i>_n<i|_n|0>_k<0|_k.$$
(5)
### Step 2 - Function Evaluation

The function $`f`$, analogous to the oracle in Grover’s search algorithm, is evaluated by applying a reversible unitary transformation $`U_f`$. The transformation has no effect on the eigenstates $`|i>_n`$, but partitions the output register into a set of subensembles, $`|f_i>_k`$ , defined as the sets of $`k`$ spins which share the same state. Therefore we have,

$$U_f|i>_n|0>_k|i>_n|f_i>_k.$$
(6)
In a physical implementation, $`U_f`$ would be the product of a sequence of fundamental unitary transformations for each of the $`k`$ spins in the output register. The unitary extension of $`U_f`$ to arbitrary output states uses the bitwise XOR operator $``$:

$$U_f|i>_n|j>_k|i>_n|jf_i>_k.$$
(7)
The transformation $`U_f`$ is applied to the system, in order to evaluate the function $`f`$ simultaneously on the linear combination of all sample points $`i`$ given by the initial state in Eq.(5). This is equivalent to evaluating the function $`f`$ concurrently on $`N`$ classical processors in a parallel computer.

Since the initial mixed state is a density operator, the action of $`U_f`$ can be written as

$$U_f(\rho _{in}^{(n)}\rho _{out}^{(k)})U_f^{}=\frac{1}{N}\underset{i=1}{\overset{N}{}}|i>_n<i|_n|f_i>_k<f_i|_k$$

This operation transforms the state of the output register to a mixture that represents all of the approximate function values $`f_i`$, for $`i=1,2,\mathrm{},N`$. We remind the reader that the finite-precision values $`f_i`$ are represented by a binary string of $`k`$ spins. The states of the output register can be written in terms of individual spins,

$$|f_i>_k=|b_{f_i1}>|b_{f_i2}>\mathrm{}|b_{f_ik}>,$$
(8)
where $`(b_{f_i1},b_{f_i2},\mathrm{},b_{f_ik})\{0,1\}`$ are the digits of the approximate function value $`f_i`$ in binary format.

We use the following binary encoding scheme to approximate the set of function values $`f(i)`$ using the $`2^k`$ states available in the output register:

$`|0_1>|0_2>\mathrm{}|0_k>`$ $``$ $`f(i)[0,\delta )`$ (9)
$`|1_1>|0_2>\mathrm{}|0_k>`$ $``$ $`f(i)[\delta ,2\delta )`$ (10)
$`|0_1>|1_2>\mathrm{}|0_k>`$ $``$ $`f(i)[2\delta ,3\delta )`$
$`\mathrm{}`$
$`|1_1>|1_2>\mathrm{}|1_k>`$ $``$ $`f(i)[1\delta ,1].`$ (12)
Alternatively, the approximate function values $`f_i`$ can be defined directly in terms of the individual spin values, for example,

$$f_i=\delta \underset{j=1}{\overset{k}{}}2^{j1}b_{f_ij},$$
(13)
to set $`f_i`$ equal to the start of the range intervals given in the encoding scheme above.

Up to this point, the query complexity of the summing algorithm is one i.e., only one function invocation is required. An important, but separate issue that naturally arises is whether $`U_f`$ can be implemented efficiently. Despite a query complexity of $`O(1)`$, the computational complexity of the algorithm can be much higher, if the evaluation of the function $`f`$ is costly. We will not discuss the computational complexity in this paper, since we want to keep $`f`$ as general as possible, but we note that if $`f`$ is a classically efficiently computable function, then $`U_f`$ can be implemented with comparable complexity, as discussed by Nielsen and Chuang Nielsen . Moreover, functions that cannot be computed efficiently by classical devices may be rendered tractable in the future by other quantum algorithms.

### Step 3 - Measurement

In the last step of the algorithm, we measure the average value of the output register in the final ensemble given by Eq.(II). This measurement is an analog process, which performs a single concurrent evaluation, unlike the recursive evaluation on conventional digital computers. The result of the measurement is an ensemble average of all the approximate function values $`f_i`$,

$$\overline{f_i}=\frac{1}{N}\underset{i=1}{\overset{N}{}}f_i,$$
(14)
which depends on the distribution of $`f_i`$ and the precision $`\delta =2^k`$ of the encoding. Finally the desired sum $`S_{N,k}`$ can be obtained by multiplying the average value $`\overline{f_i}`$ and the total number of sample points $`N`$.

Note that in the statement of the problem we have assumed that the number of function samples is a power of two, $`N=2^n`$. This ensures that memory resources are employed optimally, by using every possible state of the input register to encode the sample points $`i`$. However, in general, the number of function samples can be arbitrary, in which case only a subset of the input register states is used to represent the sample points.

In the physical implementation of the algorithm, each of the spins in the output register generates an output signal $`\gamma _j`$, proportional to the number of spin subensembles that have the $`j`$-th spin in the state $`|1>`$. Each signal $`\gamma _j`$ can be transformed into a fraction $`\overline{\gamma }_j[0,1]`$,

$$\overline{\gamma }_j=\frac{\gamma _j}{\mathrm{\Gamma }_j},$$
(15)
by calibration against the maximum output signal $`\mathrm{\Gamma }_j`$, which is obtained when spin $`j`$ of the output register is set to $`|1>`$ for all subensembles. The normalized output signals $`\overline{\gamma }_j`$ are then multiplied by the corresponding binary weight $`2^{j1}`$, cf., Eq.(13), for $`j=1,2,\mathrm{}k`$, to give the ensemble average of the output register,

$$\overline{f_i}=\frac{1}{2^k}\underset{j=1}{\overset{k}{}}2^{j1}\overline{\gamma }_j,$$
(16)
and hence the sum $`S_{N,k}`$.

If the measurement sensitivity of the experiment is adequate to distinguish between distinct normalized output signals with a precision equal to or better than $`1/N`$, then the query complexity remains $`O(1)`$.

However, as the number of sample points $`N`$ increases, the sensitivity will eventually become inadequate. At that point, we cannot reliably measure the output signals for each spin in the output register, and significant differences between normalized output signals, differences larger than $`1/N`$, will not be detectable in a single experimental trial. We note that the error in the ensemble average value $`\overline{f_i}`$ is given by the weighted sum of the measurement errors for each of the spins in the output register, using the exponentially increasing weights in Eq.(16). To enhance the measurement sensitivity, the algorithm is repeated a number of times. For example, in an NMR implementation the proposed algorithm will have to be repeated $`N^2`$ times, taking into account the square-root scaling of the signal-to-noise ratio $`S\sqrt{N_e}`$ with the number of experimental trials, $`N_e`$.

## III Discussion

The ensemble summing algorithm proposed in this paper uses a radically different approach from previous summing algorithms, which rely either on the extrinsic parallelism of conventional digital computers, or on the intrinsic parallelism of entangled quantum states. Instead, the proposed algorithm uses the parallelism of mixed states in an ensemble of spins to evaluate a given function once, and then extracts the measurement result by ensemble averaging. We note that the lower bound derived by Nayak and Wu Nayak98  for the query complexity of quantum algorithms that calculate the mean of a function, does not apply to our ensemble algorithm, since the polynomial method Beals98  used in their derivation applies to pure states, not mixed states as in the ensemble case. The appropriate lower bound for the ensemble summing algorithm is obtained by considering the query complexity for a classical parallel computer with $`N`$ processors, which is $`O(1)`$.

Table 1 shows the query complexities of the present ensemble summing algorithm and the ensemble search algorithm proposed by Bruschweiler Bruschweiler00 . We compare both ensemble algorithms with Grover’s search algorithm Grover97a , which provides the critical speedup in existing quantum summing algorithms.

The ensemble summing algorithm has an exponential advantage in terms of query complexity, relative to the implementation of the quantum summing algorithm using Grover’s search algorithm with pseudopure states in NMR, as the two algorithms have the same scaling for the measurement sensitivity. Similarly, the query complexity of the ensemble search algorithm is exponentially smaller than that required for Grover’s search algorithm using pseudopure states in NMR.

However, a comparison with the (theoretical) implementation of Grover’s search algorithm using pure states shows that both the ensemble summing algorithm and the ensemble search algorithm will be more efficient only for a total number of samples below a threshold value determined by the measurement sensitivity. Bruschweiler Bruschweiler00  estimates that the ensemble search algorithm is more efficient for databases of size $`N`$ which fulfill the condition

$$N\sqrt{N}\mathrm{log}_2N<S^2,$$
(17)
where $`S`$ is the signal-to-noise ratio of measurements in an NMR implementation.

The same reasoning leads to the conclusion that the ensemble summing algorithm is more efficient than the quantum summing algorithm using Grover’s search algorithm with pure states, for a number of function samples $`N`$ given by

$$N\sqrt{N}<S^2,$$
(18)
with respect to the signal-to-noise ratio $`S`$. The best available signal-to-noise ratio in present NMR technology is $`S10^4`$, which results in an efficiency threshold value $`N_{max}2\times 10^5`$. However for values of $`N>N_{max}`$, the query complexity of the ensemble summing algorithm is $`O(N^2)`$, which is greater than the query complexity for quantum and classical summing algorithms.

A different type of error occurs in the initialization step of the algorithm, if a room-temperature thermal state is used instead of the initial state $`\rho _{in}^{(n)}`$ in Eq. (3). For an NMR implementation, Gershenfeld and Chuang Gershenfeld97  give the thermal state for $`n`$ spins in the form

$$\rho _{th}^{(n)}=\frac{1}{N}\underset{i=1}{\overset{N}{}}|i>_n<i|_n+\frac{\alpha }{N}\underset{i=1}{\overset{N}{}}\chi _i|i>_n<i|_n,$$
(19)
where $`\alpha =\frac{\mathrm{}\omega }{2k_BT}`$ ($`10^6`$ at room temperature) is the Boltzmann factor of the deviation density matrix, and the coefficients $`\chi _i[n,n]`$ represent the net integer sum of various spin-up and spin-down combinations of the $`n`$ spins.

The effect of using a thermal state instead of the equally-weighted mixed state $`\rho _{in}^{(n)}`$ becomes present in the measurement step of the algorithm, which results in an ensemble average

$$\overline{f_i}^{}=\frac{1}{N}\underset{i=1}{\overset{N}{}}f_i+\frac{\alpha }{N}\underset{i=1}{\overset{N}{}}\chi _if_i.$$
(20)
The error term in this average imposes an additional limit on the accuracy of the algorithm. In the worst case, the difference between the measured and desired values is bounded by

$$|\overline{f_i}^{}\overline{f_i}|\frac{n\alpha }{2},$$
(21)
since $`|\chi _i|n`$ and $`f_i1`$. In general, the error is much smaller than $`\frac{n\alpha }{2}`$, as the deviation matrix is traceless, i.e., $`\chi _i=0`$. The Boltzmann factor $`\alpha `$ can also be reduced by either raising the temperature $`T`$, or by choosing a lower average resonant frequency $`\omega `$ for the spins.

The ensemble summing algorithm presented in this paper can be applied to estimating the mean and/or the definite integral of a multi-dimensional function. The validity of the algorithm holds for rather general classes of functions, ranging from continuous to the Lebesgue measurable and integrable classes, $`L^p,1q<\mathrm{}`$. The former is based on the convergence of the Riemann sums to the Riemann integral of ”ordinary” functions, while the latter is based on the density of simple (Boolean) functions in $`L^p,1q<\mathrm{}`$ DS . Unfortunately, these results simply state that in the limit of infinite number of terms, the approximating sums coincide with the desired integrals. An evaluation of the error made when using a finite number of terms in the sum is impossible, in general. To estimate this error, one has to resort to the specific properties of the approximated (sub-classes of) functions. For instance, for any Lipschitz function $`f`$ with Lipschitz constant $`L`$, integrated over the finite interval $`[a,b]`$, the error, $`E_N`$, between the integral of $`f`$ and the Riemann sum with $`N`$ terms evaluated at equidistant points is bounded by $`(ba)L/N`$. This error can now be combined with the error made when estimating the discrete sum with $`N`$ terms (see Section II.1) and an efficient algorithm can then be devised for the estimation of the integral. For general functions, the expression of the error as a function of $`N`$ is unknown, although it is known that $`lim_N\mathrm{}E_N=0`$. However, if $`E_N`$ decreases with $`N`$ in a much slower fashion, say like $`E_N(lnN)^1`$, this would translate into a significant increase of the number of terms in the sum and therefore an increased complexity, to achieve a given overall precision for the integral. The relationship between various functional classes and their approximants by Boolean functions is an active research topic that addresses such notions as the complexity, capacity, and entropy of a function, which go beyond the scope of the present paper. The interested reader is referred to Refs. Heinrich2001a ; Heinrich2001b ; L ; TWW ; V .

We conclude by pointing out the two main advantages of the proposed ensemble summing algorithm. First, there is no need to maintain quantum coherences for ensemble algorithms, so they are easier to implement than their quantum counterparts. Indeed, the scaling of the resources required to maintain entanglement in pure state-based algorithms for the duration of the computation remains an open question, and could have a potentially large negative effect on the exact threshold value $`N_{max}`$ of the number of function samples.

Second, for a restricted number of function samples, $`NN_{max}`$, which is determined by the measurement sensitivity, ensemble algorithms may give an exponential speedup over all known quantum and classical summing algorithms. In this regime, the proposed ensemble summing algorithm requires only a single invocation of the function $`f`$.

## Acknowledgments

This research was partly supported by the U.S. Department of Energy, Office of Basic Energy Sciences. The Oak Ridge National Laboratory is managed for the U.S. DOE by UT-Battelle, LLC, under contract No. DE-AC05-00OR22725.