text:
# Simulating the Sunyaev-Zel’dovich effect(s): including radiative cooling and energy injection by galactic winds

## 1. Introduction

On angular scales below $`10^{}`$, the microwave sky carries the imprint of large-scale structure in the low-$`z`$ universe. In particular, cosmic microwave background (CMB) photons propagating through the universe are inverse Compton or Doppler scattered by hot electrons along their path, either in dense structures such as clusters of galaxies or more generally in hot gas in the intergalactic medium. Inverse Compton scattering conserves the number of photons but preferentially increases their energy, leading to a spectral distortion whose amplitude is proportional to the product of the electron temperature and density (or pressure). Doppler scattering induces an intensity fluctuation with the same spectral shape as the CMB itself. These effects were first described by Sunyaev & Zel’dovich (SZ72 (1972, 1980)) and are known as the thermal and kinetic SZ effects, respectively (for recent reviews see Rephaeli Rep (1995) and Birkinshaw Bir (1999)).

The thermal effect is one of the primary sources of secondary anisotropies in the CMB on small angular scales. The change in the (thermodynamic) temperature of the CMB resulting from scattering off non-relativistic electrons is

$`{\displaystyle \frac{\mathrm{\Delta }T}{T}}`$ $`=`$ $`y\left(x{\displaystyle \frac{\mathrm{e}^x+1}{\mathrm{e}^x1}}4\right)`$ (1)
$``$ $`2y\text{for }x1,`$ (2)
where $`x=h\nu /kT_{\mathrm{CMB}}\nu /56.85`$GHz is the dimensionless frequency, and the second expression is valid in the Rayleigh-Jeans limit which we shall assume henceforth. The quantity $`y`$ is known as the Comptonization parameter and is given by

$$y\sigma _T𝑑\mathrm{}\frac{n_ek(T_eT_{\mathrm{CMB}})}{m_ec^2},$$
(3)
where the integral is performed along the photon path. Since $`T_eT_{\mathrm{CMB}}`$ the integrand is proportional to the integrated electron pressure along the line of sight.

The kinetic SZ effect arises from the motion of ionized gas with respect to the rest-frame of the CMB. The resulting temperature fluctuation is $`\mathrm{\Delta }T/T=b`$, where

$$b\sigma _T𝑑\mathrm{}n_e\frac{v_r}{c}$$
(4)
measures the magnitude of the effect along the line of sight if $`v_r`$ is the radial peculiar velocity of the gas (positive if the cluster is moving away from us, negative if it is moving towards us). The different dependence on frequency of the two effects can, in principle, be used to disentangle them observationally, though the primary CMB anisotropies provide “noise” for the kinetic SZ signal.

In this paper, we report results of a sequence of high resolution hydrodynamic simulations designed to study the SZ signal on small angular scales. The imprint of the SZ effects on the CMB has been studied by a number of authors (Ostriker & Vishniac OsVi (1986); Persi et al. PSCO (1995); da Silva et al. SBLT00 (2000, 2001); Refregier et al. RKSP (2000); Seljak et al. SelBurPen (2000); Kay, Liddle & Thomas KayLidTho (2001); Refregier & Teyssier RefTey (2001); da Silva et al. SKLTPB (2001); Zhang, Pen & Wang ZPW (2002)). Ours are the first such simulations to incorporate the effects of radiation and star formation, including gas cooling and feedback from supernovae and galactic winds. As such they extend our earlier work on the subject, reported in Springel, White & Hernquist (SprWhiHer (2001)).

The outline of this paper is as follows. In Section 2, we briefly describe our simulations and the techniques we use to compute SZ maps. In Section 3, we present our results for the power spectra of thermal and kinetic SZ effects and discuss the effects of non-Gaussianity on the sample variance. We show SZ source counts as a function of source strength and source size and how SZ source properties correlate with the intrinsic properties of the cluster. Finally, we summarize and discuss our results in Section 4.

## 2. Method

We chose as our cosmological model the currently favored $`\mathrm{\Lambda }`$CDM cosmology. In particular, we adopted $`\mathrm{\Omega }_\mathrm{m}=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`H_0=100h\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ with $`h=0.7`$, $`\mathrm{\Omega }_\mathrm{B}=0.04`$, $`n=1`$ and $`\sigma _8=0.9`$. This model yields a reasonable fit to the current suite of cosmological constraints and as such provides a good framework for making realistic predictions. The simulation domain was a $`100h^1`$Mpc cube with periodic boundary conditions. We ran smoothed particle hydrodynamic (SPH) simulations at two different mass resolutions corresponding to $`2\times 144^3`$ and $`2\times 216^3`$ particles, with initial conditions constructed such that the large-scale modes were equal for both resolutions, allowing us to perform a resolution study on an object by object basis without being affected by cosmic variance. For each mass resolution, we performed a simulation that computed only adiabatic111In what follows, we will refer to the runs that ignored radiative effects and star formation as ‘adiabatic’ but it should be kept in mind that these simulations did include shock heating and so are not strictly speaking adiabatic. gas physics and shock heating, and one that also followed the radiative cooling and heating processes of the gas in the presence of a UV background radiation field together with star formation and feedback processes in the dense collapsed gas. In Table 1, we list some basic numerical parameters of the four simulations that we carried out.

Radiative cooling and heating was modeled as in Katz, Weinberg & Hernquist (Ka96 (1996)), with a photo-ionizing flux similar to the one advocated by Haardt & Madau (Ha96 (1996)). Radiation from quasars reionizes the Universe at redshift $`z6`$ in this scenario (for details, see e.g. Davé, Hernquist, Katz & Weinberg Da99 (1999)). We have treated star formation in the framework of an effective multi-phase model for the interstellar medium (ISM), as detailed by Springel & Hernquist (2002a ). In this model, rapidly cooling gas of high overdensity is assumed to give rise to the formation of cold clouds, embedded in an ambient hot medium. The cloud material forms the reservoir of baryons available for star formation. Massive stars are taken to explode as supernovae on a short timescale, and to release their energy as heat to the ambient medium of the ISM. Supernovae are also assumed to evaporate cold clouds, thereby establishing a tight self-regulation cycle for the star-forming ISM.

In addition, we have included a phenomenological model for galactic outflows in our simulations, motivated by the large body of observational evidence for the ubiquitous existence of such galactic winds (e.g. Martin Mar99 (1999)). Galactic-scale outflows are thought to play a crucial role in the enrichment of the low-density intergalactic medium with heavy elements and for the global regulation of star formation, and they may have influenced the properties of intragroup and intracluster gas in a decisive way. Note that without the inclusion of such strong feedback processes, simulations tend to substantially overpredict the luminosity density of the universe (Balogh et al. Bal01 (2001)). However, the model for galactic winds that we use here is able to reduce the global efficiency of cooling and star formation sufficiently to match these constraints (Springel & Hernquist 2002c ).

In this model, each star-forming galaxy drives a wind with a mass-outflow rate equal to two times its star formation rate, and with a wind-velocity of $`v_w=484\mathrm{km}\mathrm{s}^1`$. These choices are deliberately ‘extreme’ in the sense that the total kinetic energy of the wind is of the order of the energy released by the supernovae. However, these parameters are quite typical for the observed properties of outflows from star-forming disks. Whether or not the wind can escape from a galaxy depends primarily on the depth of its dark matter potential well. Halos with virial temperatures below $`10^6`$K can lose some of their baryons to an outflowing wind, while more massive halos will contain the wind, such that the wind becomes progressively more unimportant for the regulation of star formation in more massive objects.

The simulations considered here are part of a larger simulation program which shows that the particular feedback model we employ allows a numerically converged prediction for the total star formation density of the universe as a function of epoch, provided all star-forming halos above virial temperatures of $`10^4`$K can be resolved. Of course, the simulations analyzed here cannot reach such a high mass resolution, because they need to model a volume large enough to contain rich clusters of galaxies. It is thus inevitable that part of the star formation is lost due to resolution limitations. This can be a sizable effect. Springel & Hernquist (2002c ) show that about 50% of all stars are born in halos with virial masses below $`10^{12}h^1\mathrm{M}_{}`$, which corresponds roughly to the mass of the smallest halos in our G4 simulation with reliable estimates for the star formation rate. Note that the ‘unresolved’ star formation in small halos is also the one that would predominantly be able to drive truly outflowing winds. It should therefore be kept in mind that the effects of galactic outflows on the intracluster medium are likely going to be underestimated by the present simulations.

The simulations were performed on an Athlon-MP cluster at the Center for Parallel Astrophysical Computing (CPAC) at the Harvard-Smithsonian Center for Astrophysics. We used a modified version of the GADGET-code (Springel, Yoshida & White SprYosWhi (2001)), and integrated the entropy as the independent thermodynamic variable (Springel & Hernquist 2002b ).

In all our runs, the full simulation box was output every $`100h^1`$Mpc along the line-of-sight between redshifts $`z19`$ and $`z=0`$. This corresponds to 77 dumps per model. We made maps222The missing factor of $`h`$ in our code which caused the maps from Springel et al. (SprWhiHer (2001)) to be in error has been corrected here. in a manner similar to da Silva et al. (SBLT00 (2000)) as described in our earlier paper (Springel et al. SprWhiHer (2001)). To recap briefly we stacked the various outputs from the boxes back along the line of sight, with each box randomly translated and oriented in the direction of either the $`x`$, $`y`$ or $`z`$ axis. Then, a grid of $`512^2`$ rays subtending a constant angle of $`1^\mathrm{o}`$ from the observer was traced through the boxes starting at $`z=19`$.

We produced maps of the $`y`$-parameter, the Doppler $`b`$-parameter, the projected gas density and the thermal bremsstrahlung emission. In the adiabatic simulations we do not track the ionization fraction, $`x_e`$, dynamically. Instead, we set it to 1.158, corresponding to full ionization, for temperatures above $`10^4`$K and take the gas at lower temperature to be neutral. The extra physics models self-consistently track $`x_e`$ and this is used in the computation.

It is also convenient to have a ‘cluster catalog’ from the simulations. We produced this by running a friends-of-friends group finder (e.g. Davis et al. DEFW (1985)) with a linking length of $`b=0.15`$ (in units of the mean inter-particle spacing). The FoF algorithm partitions the particles into equivalence classes by linking together all particle pairs separated by less than a distance $`b`$. We define the center of each halo as the position of the potential minimum, calculating the potential using only the particles in the FoF group. This proved to be more robust than using the center of mass, as the potential minimum coincided closely with the density maximum for all but the most disturbed clusters. Additionally, the position of the center was very insensitive to the presence or absence of the particles near the outskirts of the halo, and thus to the precise linking length used. For each of the halos we computed the mass, $`M_{200}`$, enclosed within a radius, $`r_{200}`$, interior to which the density contrast was 200 times the critical density. We shall sometimes refer to this as the virial radius, and denote its angular extent at the distance of the cluster by $`\theta _{200}`$.

## 3. Results

The evolution of the mean (mass weighted) temperature in the simulations is shown in Fig. 1. The formation of structure in the universe leads to a steady heating with the mean temperature today being approximately $`1/3`$keV. We can see the effects of cooling most clearly at early times in the suppression of $`T`$ in runs G3 and G4 with respect to G3a and G4a, but by $`z=0`$ most of the differences have disappeared. The slight jump in $`T`$ at $`z6`$ is due to reionization by quasars in this model, which raises the background radiation field and leads to increased photo-heating of the gas in the universe.

We show a typical example of our thermal SZ maps in Fig. 2. These maps were produced by coadding a large number of partial maps, each giving the contribution of one of the boxes that we stacked along the photons’ path. We made 15 maps for each of the 4 simulations, using the same set of random offsets and orientations for each of the 4 simulations. As we noted before (Springel et al. SprWhiHer (2001)) filaments, while prominent in the partial maps, are largely hidden in the high level of background arising from the summation over many of these structures (see also Croft et al. Cro (2001)).

The mean Comptonization we predict from our 15 maps is $`2.4\times 10^6`$ and $`2.5\times 10^6`$ for the lower and higher resolution adiabatic models respectively. This becomes $`2.2\times 10^6`$ and $`2.1\times 10^6`$ for the runs with cooling and feedback. Thus we see that the effect of cooling is more important than heating in aggregate for the gas contributing to the mean $`y`$ distortion and that our models easily pass the constraint from FIRAS (Fixsen et al. FIRAS (1996)).

### 3.1. Power spectra

For each map we computed a number of statistics, and we in general averaged the results for 15 random lines of sight to reduce field-to-field variance. Here we present the angular power spectrum of the thermal and kinetic SZ effects (Fig. 3).

A comparison of the low and high resolution simulations indicates that for the adiabatic case we have converged in the angular power spectrum out to $`\mathrm{}20,000`$. These results are also in excellent agreement with our earlier estimates (Springel, White & Hernquist SprWhiHer (2001)) once we take into account the slight differences in the cosmological model. For the runs including cooling we are more sensitive to the effects of finite mass and spatial resolution. The higher resolution simulation shows a lower power spectrum than the lower resolution spectrum, indicating that cooling has been able to proceed more efficiently when the simulation can correctly resolve denser structures. However, the final results are quite insensitive to extra physics included in these models, as might have been expected if the bulk of the signal is coming from gas outside the core of the cluster, and the gas content of the intracluster medium is not significantly depleted by gas cooling. In this respect our results differ from those of da Silva et al. (SKLTPB (2001)) who found a significant effect on the SZ signal if they include cooling in their simulations. The addition of star formation and energy injection in our simulations stabilizes the runaway cooling which otherwise occurs and significantly reduces the impact of the extra physics on the SZ effect. In the absence of such feedback the cooling is only mitigated by numerical resolution, leading to an excess of cold gas and stars in conflict with observations.

The power spectra we predict are in accord with recent measurements of arcminute angular scale anisotropy in the CMB at BIMA (Dawson et al. BIMA1 (2001); BIMA2 (2002)) though lower than the detection quoted by CBI (Mason et al. CBI1 (2002); Bond et al. Bond (2002)) from their deep field (the data from mosaic observations, Pearson et al. CBI2 (2002), do not significantly constrain the signal at high-$`\mathrm{}`$). To bring the adiabatic model into agreement with the CBI points would require us to e.g. increase $`\sigma _8`$ to 1–1.2 (or dramatically increase $`\mathrm{\Omega }_B`$; see e.g. Seljak, Burwell & Pen SelBurPen (2000) or Bond et al. Bond (2002)). The inclusion of cooling and feedback lowers our predictions by $`30\%`$ on CBI scales, which could be compensated for by an additional $`\mathrm{\hspace{0.17em}4}\%`$ change in $`\sigma _8`$.

We have also estimated the amount by which the evident non-Gaussian nature of the SZ signal increases the field-to-field variance in the power spectrum. For a Gaussian signal the fractional error per mode in the power spectrum is $`\sqrt{2}`$, so we expect the standard deviation in the binned power to be $`\sqrt{2/N}C_{\mathrm{}}`$ where $`N`$ is the number of modes in the bin. If we compare this expectation to the standard deviation computed from the 15 maps (Fig. 3), we find that the latter is higher by a factor of 2 (at $`\mathrm{}10^3`$) to 5 (at $`\mathrm{}10^4`$) due to the contributions from the 3- and 4-point functions (e.g. Meiksin & White MeiWhi (1999), §2). Since the non-Gaussianity leads to mode-coupling, the exact ratio depends slightly on the binning chosen – we show two choices. We caution the reader that all of the maps are created from the same underlying simulation volume, and so are not fully independent. This makes the estimates at lower $`\mathrm{}`$ somewhat uncertain. We do not have enough maps to accurately calculate the off-diagonal elements of the correlation matrix (which converge slowly compared to the diagonal elements), however the correlations are strong. As for the non-linear mass power spectrum, the primary fluctuation from field to field is in the amplitude, and not the shape, of the power spectrum. We show an example from 5 of our fields in Fig. 3.

Finally it is of interest to ask where the majority of the thermal SZ signal is coming from in the angular power spectrum – i.e. diffuse gas or gas within halos. To answer this we remade two of our maps excluding all gas with density smaller than $`100\overline{\rho }_{\mathrm{gas}}`$. Gas in halos should have higher density than this. The resulting maps are indistinguishable by eye, even ‘blinking’ between them. There is a small (25%) suppression of power for $`\mathrm{}<2000`$ in the maps with the density cut, but the suppression is $`<2\%`$ by $`\mathrm{}=6000`$ where the signal peaks. This suggests (da Silva et al. SKLTPB (2001)) that even the large-angle signal is coming primarily from relatively high density gas along the line-of-sight — the virialized halos – and not from the diffuse gas or filamentary structures (the ‘cosmic web’). If we increase the cut even further, to $`10^3\overline{\rho }_{\mathrm{gas}}`$, we start to significantly reduce the thermal SZ signal up to $`\mathrm{}10^4`$, indicating that the dense cores of the halos do not dominate the thermal SZ signal.

### 3.2. Source counts

A glance at Fig. 2 indicates that the SZ effect is dominated by discrete sources and that the pixel distribution is highly non-Gaussian. Thus a power spectral analysis contains only a small part of the information in the map. To obtain an understanding of the SZ source sizes and strengths we have used the source detection software SExtractor (Bertin & Arnouts Bertin96 (1996)) to count individual sources, and to measure their strength. We used default settings for the detection algorithm of the image processing software and let it estimate background and ‘noise’ levels automatically. Note that our maps are in principle noise-free, so what is interpreted as noise in this procedure effectively arises from limitations due to source confusion, which leads to non-detections of some of the faintest and smallest sources. In a typical analysis with SExtractor, we resolve slightly more than half of the total SZ signal into discrete sources.

We measure the strength of a source as the monochromatic brightness change

$$S_\nu =_\mathrm{\Omega }\mathrm{\Delta }B_\nu d\mathrm{\Omega }=f(x)B_\nu _\mathrm{\Omega }y(\theta )d\mathrm{\Omega }$$
(5)
of the CMB integrated over the solid angle of the source. Here $`B_\nu `$ is the Planck spectrum of the primary CMB, and $`f(x)`$ is the spectral function

$$f(x)=x\frac{\mathrm{e}^x}{\mathrm{e}^x1}\left(x\frac{\mathrm{e}^x+1}{\mathrm{e}^x1}4\right),$$
(6)
with $`x=h\nu /kT_{\mathrm{CMB}}`$ and we have taken $`\nu =30`$GHz.

In Fig. 4, we show the cumulative source counts per square degree found in the simulations. Again we find remarkably good agreement between all of the simulations, with the cooling runs predicting fewer sources at the bright end by 20-30% and at the faint end by 30-50% (at fixed flux). This is not terribly surprising at the bright end where the source counts are more massive clusters. In this case the bulk of the ICM is at too low a density to cool efficiently. However our calculations suggest that this holds true even for the relatively weak SZ sources. Cooling and feedback thus appear to be a small effect on the clustering and number counts of sources, with the gas dominating the SZ signal being more prone to cooling than heating in our model.

Each of these bright sources is associated with at least one massive halo in the simulation. It is instructive to pursue this association further. In Fig. 5 we show a scatter-plot of source strength vs. halo mass, both intrinsic and from the maps. For the ‘intrinsic’ effect we calculated the SZ signal coming from the most massive clusters in our box at $`z=0`$. Specifically we summed the SZ signal within half of $`r_{200}`$ to obtain the total decrement intrinsic to the cluster:

$$Y_{\mathrm{int}}=(\gamma 1)(1Y_p)\sigma _T\underset{j}{}\mu _jx_j\frac{m_j}{m_e}\frac{u_j}{c^2}$$
(7)
where $`\gamma =5/3`$, $`Y_p=24\%`$, $`\sigma _T`$ is the Thomson cross section, $`m_e`$ the electron mass, $`c`$ the speed of light, $`\mu _j`$ is the mean molecular weight of particle $`j`$, $`x_j`$ is the fractional ionization (relative to hydrogen), $`m_j`$ is the mass and $`u_j`$ the internal energy. The sum is over all particles within333The degree of correlation is almost independent of the exact aperture used. $`r_{200}/2`$ of the potential minimum in the group, which is the region dominating the SZ emission. Since $`y`$ is dimensionless we quote $`Y_{\mathrm{int}}`$ in units of $`(h^1\mathrm{Mpc})^2`$. As the total signal scales as $`M\times T`$ and for self-similar clusters $`TM^{2/3}`$ we expect $`Y_{\mathrm{int}}M^{5/3}`$ which is a good fit to our simulations. The upper panel of Fig. 5 shows that there is a very tight relation between $`Y_{\mathrm{int}}`$ and $`M_{200}`$, as has been noted before by e.g. Holder et al. (HMCEL (2000)). It is this relation which is behind the claims that a $`y`$-selected sample will be close to mass limited.

However, $`Y_{\mathrm{int}}`$ is not an observable. To mimic something which is closer to what an observation could probe we first found all of the clusters in the simulation which fell in the line-of-sight of a particular map. For each of these we integrated the Compton $`y`$ parameter values in pixels which lay within half $`\theta _{200}`$ of the cluster center to obtain $`Y_{\mathrm{map}}=𝑑\mathrm{\Omega }y`$. As the lower panel of Fig. 5 shows, there is a large amount of scatter in this $`Y_{\mathrm{map}}M`$ relation, even for our noise free, high resolution maps. This scatter is approximately constant as we vary the aperture used to define $`Y_{\mathrm{map}}`$ and can be traced to 3 prime sources. First, there is evolution in the $`MT`$ relation with redshift, and the clusters we have plotted are from a range of redshifts. Second, clusters are not spherical, leading to differences in signal when viewed from different orientations. Finally, the signal we see is the projection of SZ signal from the entire line-of-sight444This contribution was not included in the work of Holder et al. (HMCEL (2000)) who simulated only gas within the virial regions of their sources.. Since the signal is not dimmed by distance, objects all along the line-of-sight can contribute a significant amount to the signal in any pixel. A similar ‘confusion’ effect operates for optical and weak lensing cluster searches (White & Kochanek WhiKoc (2002); White, van Waerbeke & Mackey WvWM (2002)). Since clusters form in high density regions, with a lot of dense material nearby, there is a non-negligible chance that some lines of sight will be enhanced by non-cluster material. The extreme outliers in Fig. 5 are all clusters closer to the observer than about $`600h^1`$Mpc, whose virial radii subtend a large angle and are most likely to suffer from this projection effect. These clusters are marked with crosses in the figure.

To isolate this effect further we have computed the integrated SZ signal coming from spheres of gas, centered on rich clusters, viewed from many random directions. Specifically, we first extract all of the gas particles in the simulation within a distance $`r_{\mathrm{cut}}`$ of the minimum of the cluster potential. From this set we compute, using Eq. (7), the total $`Y`$ for all particles whose 2D projected radius $`R<r_{200}/2`$ for several lines-of-sight. We also compute $`Y`$ from the subset of the particles contained within our FoF group.

Our results are shown in Fig. 6 for $`r_{\mathrm{cut}}=5h^1`$Mpc and $`r_{\mathrm{cut}}=10h^1`$Mpc. For several clusters there is a large scatter between lines-of-sight depending on the presence or absence of material near the cluster. The scatter is larger for the larger radius, indicating that it is not just neighboring material which drives the scatter. If we include only those particles within the FoF group we still see a large scatter, though the SZ signal is noticeably depressed. Thus, some of the scatter is coming from asphericity of the cluster, but a non-negligible fraction of the signal is coming from non-cluster material along some lines-of-sight.

The large scatter that we see in the SZ signal at fixed mass has implications for SZ searches for galaxy clusters – any method which uses primarily SZ flux in target selection will not be mass limited. In fact, the steeply falling mass functions predicted by hierarchical theories would imply that some high-flux candidates are likely to be low mass groups with abnormally high SZ signal rather than the intrinsically rarer rich clusters. If this is true, a blind SZ search should turn up many groups in addition to the desired cluster sample.

## 4. Conclusions

We have performed a sequence of high resolution hydrodynamic simulations of structure formation in a $`\mathrm{\Lambda }`$CDM model to investigate the effect of extra physics on the thermal and kinetic Sunyaev-Zel’dovich effects. Including only adiabatic gas physics our simulations of the thermal effect are converged down to sub-arcminute scales and agree well with our earlier work. For the runs with extra physics we see an increased efficiency of cooling in the higher resolution simulation which leads to a lowering of the mean Compton $`y`$ parameter and a slight depression in the SZ angular power spectrum. Overall however the effects of including the extra physics on the 1- and 2-point functions are slight, indicating that the SZ effect is not overly sensitive to the (uncertain) thermal history of the gas.

The maps are noticeably non-Gaussian and dominated by discrete sources. This non-Gaussianity manifests itself in increased scatter, by a factor of 2-6, in $`C_{\mathrm{}}`$ from field-to-field compared to the Gaussian expectation.

The properties of the sources in the map are quite insensitive to the extra physics we have included also, and we are thus able to make robust statements about their detection. The brightest SZ sources have a density of $``$1 per square degree and their detection requires an angular resolution of $``$1, and the capability to separate brightness fluctuations of size $``$10$`\mathrm{mJy}`$ from the primary CMB fluctuations.

We have found that there is significant scatter in the SZ fluxes within apertures centered on clusters. This scatter has 3 sources: evolution of the mass temperature relation, asphericity in the matter distribution and line-of-sight projection. The existence of this scatter means that surveys using only SZ flux will select objects with a wide range of masses, including lower mass groups.

We would like to acknowledge useful conversations with Kyle Dawson, Bill Holzapfel, Adrian Lee and Saul Perlmutter. This work was supported by ACI 96-19019, AST 98-02568, AST 99-00877, AST 00-71019. M.W. was supported by a NASA Astrophysical Theory Grant and a Sloan Fellowship. The simulations were performed at the Center for Parallel Astrophysical Computing at the Harvard-Smithsonian Center for Astrophysics.


format: text

text:
# Imaging the Circumstellar Environments of Young Binaries in Southern Star-Forming Regions

## 1  Introduction

It has recently become clear that the great majority of the young low-mass stars in the nearest active star-forming regions (SFRs) are found in binary or multiple systems (Leinert et al. 1993; Ghez et al. 1993; Ghez et al. 1997). Although there is evidence that the mass in circumstellar material is smaller in binaries with separations $`\mathrm{¡}\mathrm{}100`$ AU than in wider binaries or single stars (Jensen, Mathieu, & Fuller 1996), a majority of young binaries emit excess radiation at infrared or submillimeter wavelengths which is typically presumed to arise in dusty disks. If the presence of one or more stellar companions does not severely hamper a disk’s ability to produce a planetary system, then binary and multiple systems probably contain the majority of the planets in the nearby SFRs.

At present, only two multiple-star systems are known to contain well-ordered, edge-on circumstellar disks visible because of their occlusion of the central star. The HK Tauri B and HV Tauri C disks exhibit a morphology much like that predicted by scattered-light models: A pair of bright lenslike reflection nebulae trace the upper and lower surfaces of the disk, while the dense region nearer the midplane is dark (Stapelfeldt et al. 1998; Koresko 1998; Monin & Bouvier 2000). There is no central point source seen in either the visible or near-infrared images. The near-infrared color of the HK Tauri B star+disk system is bluer than that of the primary star, presumably because the light we see is scattered by small particles.

A handful of visual binary companions to T Tauri stars are “infrared companions” (IRCs), objects unusually faint at visible wavelengths but bright in the infrared. The IRCs are mysterious because their spectral energy distributions show strong reddening and infrared excess reminiscent of stars at younger evolutionary phases than their primaries (e.g., Koresko, Herbst, & Leinert 1997). The structures seen in the two IRCs which have been resolved in published high-resolution images are not edge-on disk systems, but a diffuse and irregular nebula (Koresko et al. 1999) and a close double (Koresko 2000). Among the present sample, the companion to Elias 2-22 is known from previous observations (Chelli et al. 1988) to be an IRC.

Clearly, the circumstellar environments of young binary stars exhibit a range of phenomena which are accessible to modern imaging techniques. This paper describes the first results of an ongoing effort to characterize the circumstellar environments in a larger sample of young binary systems. The sources were chosen from the binaries detected in the visible-light imaging survey of Reipurth & Zinnecker (1993; hereafter RZ). They include all the binaries, with the exception of Sz 120, with separations between 1″.0 and 3″.7 and right ascensions between $`16^\mathrm{h}06^\mathrm{m}`$ and $`17^\mathrm{h}10^\mathrm{m}`$ (equinox 1950). This RA range includes all the RZ binaries in the Ophiuchus SFR, and several additional groups, for a total of 14 systems. The range of separations was chosen to avoid having the seeing-limited Point-Spread Functions (PSFs) overlap or fall off the edges of the detector. With one exception, the distances to the targets are given by RZ as 160 pc. Their physical separations are larger than 100 AU, so they are not expected to be strongly depleted of disk material via disk-companion interactions (Jensen, Mathieu, & Fuller 1996).

The sources are a somewhat mixed bag, consisting mainly of T Tauri stars but including at least one Herbig Ae/Be star (Elias 2-49) and a likely Post-T Tauri star (WSB 4). Many of the sources are well-studied, but a handful lack even published spectral types. Ten of the 14 systems are members of RZ’s “systematic survey sample”, which for the present purpose means that they were known to be members of dark clouds $`160`$ pc distant and are included in a list of young stars selected on the basis of a youth indicator such as near- or far-infrared excess or Ca II HK line emission. The remaining four systems are AS 205, both of whose components are Classical T Tauri stars (CTTS) (Prato et al. 1997); Elias 2-26, another CTTS (Martin et al. 1998); Elias 2-49, a Herbig Ae/Be star in Hillenbrand’s Group I (Hillenbrand 1992); and ESO H$`\alpha `$ 282, about which little seems to be known beyond the fact that it is an emission-line source.

The 14 pre-main sequence binary systems were observed using speckle holography at near-infrared wavelengths $`2\mu `$m to produce diffraction-limited images with dynamic range of $`2000`$ at 0″.3 from the star, which is much better than would have been possible with standard speckle imaging. The improved dynamic range is possible because the PSF is measured simultaneously with the science observations. Near-infrared photometry in the J, H, K, and L’ bands was obtained the following night to measure the infrared excess of each component as a tracer of warm dust.

Two new tertiary companions with separations $`0`$″.1 were detected, and an additional 7 objects were found to be marginally resolved. Of these, 5 are extended along only one direction, which is consistent with the presence of a very close tertiary companion or a narrow extended structure such as a compact edge-on disk. The remaining two marginally-resolved objects were extended in all directions, which would be consistent with extended nebulosity but not with tertiary companions.

No new edge-on disks were unambiguously detected. It is known from theoretical models (e.g., Wood et al. 1998) and from recent high-resolution images (e.g., Burrows et al. 1996) that pre-main sequence circumstellar disks with fiducial mass distributions are readily detectable in visible-light or near-infrared images only when they are viewed along a line of sight close to the disk midplane. In this orientation, the direct radiation from the central star is blocked, making it possible to see the much fainter scattered light which traces the disk surface. The range of viewing angles over which the disk hides the star, and therefore the likelihood of detecting it if its orientation is random, depends on the disk’s mass, radius, and flaring, and potentially on the degree to which the scattering dust has settled toward the midplane.

## 2  Observations

### 2.1  Holographic Imaging

The holography data were taken at the 10 m Keck 1 telescope on 1999 June 1 using the Near-Infrared Camera (NIRC; Matthews & Soifer 1994) with the NIRC Image Converter (Matthews et al. 1996) to produce a magnified pixel spacing of 20 mas, approximately Nyquist-sampling the $`50`$ mas diffraction limit at the $`2\mu `$m observing wavelength. The field of the $`256\times 256`$ pixel detector was 5″.1. The optical bandpass was set using either a standard K photometric filter, or, for stars which would saturate the detector in the K band, one of the narrower Ks or CH4 filters. The observations of each visual binary consisted of hundreds of frames with an integration time of 0.14 sec per frame. This short exposure time partially “froze” the atmospheric seeing, so that the raw PSF consisted of distinguishable speckles. The separations of these visual binaries were large enough compared to the 0″.5 typical seeing that only a small fraction of the frames were rejected due to poor separation of the stars. The data for each visual binary system were divided into $`7`$ series of frames, with $`140`$ frames in each series. Series taken on blank sky were used to estimate the background and bias level. A journal of the holographic observations is given in Table 1.

Individual frames were calibrated by subtracting mean sky frames, dividing by flatfield images, and “fixing” bad pixels. A model was computed for the detector “bleed” signal which extended along the readout direction, as described in the NIRC manual, and this was subtracted from the calibrated frame. For each frame, a subframe centered on each of the stars was extracted, with one subframe being treated as the science measurement and the other serving as a measurement of the instantaneous point-spread function (PSF). The subframes were either $`64\times 64`$ or $`48\times 48`$ pixels in size; the smaller subframes were necessary for doubles whose separations were less than 64 pixels. The Fourier power spectra of the target and PSF subframes, and the cross-spectrum (i.e., the Fourier transform of the cross-correlation) of the two subframes, were computed for each input frame. These quantities are meaningful within a roughly hexagonal region of the $`(u,v)`$ plane containing spatial frequencies below the diffraction limit of the hexagonal Keck primary mirror.

The cross-spectrum, the target power, and the PSF power were averaged over the whole series of good frames. An averaged Fourier power spectrum and Fourier phase for the target, corrected for the distortions due to the atmosphere and the telescope, were reconstructed from them as follows: The power spectrum was estimated as the ratio of the average target power to the average PSF power, after removing a noise bias term from each. The phase was estimated as the phase of the averaged cross-spectrum. Uncertainties in the power spectrum and phase were estimated from the scatter between the averages computed for the $`7`$ series for each visual binary system. The PSF star in each binary was chosen arbitrarily, and the choice was reversed if it showed signs of resolved structure. In no system was there clear evidence that both stars were resolved, although in principle the technique would be insensitive to this if the morphologies of the two stars were sufficiently similar.

The reconstructed Fourier transforms were compared with models of single stars, double stars, and elliptical Gaussian halos. These choices clearly cannot span the range of all possible astrophysical models, but for targets which are only marginally resolved, they are sufficient to fit the features in the data. There is a degree of degeneracy in these model fits, in the sense that an unresolved star can always be fitted by a double-star model (with either zero separation or zero brightness ratio) or by a halo model (with zero axis lengths). Further, an object which is marginally resolved along only one axis can be fitted by either a double-star model or a halo model with a zero semiminor axis, since either model produces a power spectrum having a single broad stripe, and a nearly featureless phase, in the measured frequency range.

The comparisons with models were done separately for each series of frames. Only frequencies below that for which the median uncertainty rose to four times its low-frequency value at were used. For the double-star and halo models, the fitting was done 11 times for each series, with the starting values of the parameters varied randomly, and the fit with the smallest $`\chi ^2`$ value was chosen to represent that series. This was necessary because for some starting parameters the fitting routine became trapped in local minima in the $`\chi ^2`$ surface. The final parameter values for each model were the averages of the results for the series, and the quoted uncertainties are the standard deviations of the means of those results. These statistical uncertainties are likely to be optimistic because they do not include the effects of systematic biases. The estimation of physical parameters is rather delicate and easily affected by those biases, especially for the marginally-resolved sources.

For the marginally-resolved sources, the final misfit ($`\chi ^2`$ per degree of freedom) for the power spectrum was close to unity. The power spectra are more useful than the phases for marginally-resolved systems because the power spectra begin to show strong parameter-dependent structure at lower spatial frequencies than the phases do. On the other hand, the position angles derived from power spectra suffer from an ambiguity of 180 deg (e.g., it is not possible to tell which star in a double is brighter). For the fully-resolved sources, which in the present sample are both primary-tertiary doubles, the fits to the power spectra tend to be poorer. This may be because the power spectra, being positive-definite, are susceptible to biases caused by zero-mean noise, while the phase tends to be unbiased by it. The phase fits also produce unambiguous position angles for resolved doubles. In summary, the phase data are more useful for resolved doubles, while the power spectra are more useful for marginally-resolved objects.

### 2.2  Near-Infrared Photometry

The photometric fluxes of the target visual binary systems were crudely measured in the J, H, K, and L’ filters at the Keck I telescope on 1999 June 2, by imaging them with NIRC. For the J, H, and K filters, subframes of $`256\times 128`$ pixels were used with the Image Converter, while for the L’ filter the Image Converter was removed and the subframe size was reduced to $`32\times 32`$ or $`64\times 32`$ pixels. These choices were driven largely by the need to read out the detector fast enough to avoid saturating it on the brighter stars or, at the longest wavelength, on the blank sky. The stars determined the peak signal levels in the J, H, and K filters, so the exposure times were adjusted between 68 msec and 4 sec, and no coadding was done. For the L’ filter the frames consisted of 1000 coadds of 5 msec each. A total of 10 frames were taken on each target in each filter. Periodically, frames were obtained on blank sky to estimate the background level. SR3 served as a photometric standard for the J, H, and K filters, and HD 161793 was used for the L’ filter.

The frames were calibrated in the usual way, with sky-subtraction, flatfielding, and pixel-fixing. To derive the fluxes in the frames, each binary was fitted using a model in which the two stars were represented by circular Gaussians of adjustable size, intensity, and position. Although the Gaussians were not very good models for the PSFs, the fact that the PSF in each image was essentially identical for the two stars permitted the ratio of their brightnesses and the distance between them to be estimated accurately using this technique. The brightness ratios were combined with total flux measurements to derive the fluxes for the individual stars. At J, H, and K, these total fluxes were computed by summing the values of all pixels within 32 pixels of either star, after subtracting the background level estimated as the median value of all pixels at least 64 pixels away from either star. In the L’ filter the total flux was simply the sum of all pixel-values, after subtraction of a background level. The photometric fluxes were calibrated in the usual way by comparison of these numbers with similar results on standard stars. No attempt was made to correct for the differences in airmass, as those corrections were found to be smaller than the statistical uncertainties. The final uncertainties in the magnitudes are approximately 0.2 mag.

Separations and position angles of the binaries were derived from the model fits in each filter. The platescale with IC out of the beam has been carefully measured by the Keck Observatory staff, and the pixel spacing is quoted in the instrument manual as 0″.1500. By contrast, the 0″.02 nominal platescale with the IC in is relatively uncertain. On the other hand, the measurements of the separation vectors taken with the smaller pixels are likely to be significantly more precise, in the sense of having smaller uncertainties if the platescale error is neglected. To take advantage of this, the assumed platescale in each of the J, H, and K bands was adjusted to make the mean separation for all the visual binaries (except for Sz 116, which was not measured in the L’ filter, and Elias 2-49, which is anomalous as discussed below) equal to the mean separation in the IC-out L’ images. The required IC-in pixel spacing of 0″.0207 was consistent between the J, H, and K images. Similarly, the assumed orientation of the IC-in images was adjusted to make the mean difference between the IC-in and IC-out position angle measurements vanish.

## 3  Results and Discussion

The results of the holographic imaging of the target binaries are summarized in Table 2, and the Fourier power spectra are displayed in Figure 1. The Table presents the values of $`\chi ^2`$ per degree of freedom (i.e., per frequency element in the $`(u,v)`$ plane) for unresolved-star, double-star, and halo models of the target in each of the 14 visual binary systems studied. For the two fully-resolved doubles, the fits were made to the Fourier phase, resulting in unambiguous position angle measurements, while for all other sources the fits were made to the power spectra. In the following discussion, a target is considered to be resolved if either the double-star or the halo model improves the $`\chi ^2`$ value by more than 0.5 compared to the unresolved-star model. This criterion provides a clean separation between the unresolved and extended targets.

The J,H,K,L’ photometric measurements for the 14 holographic imaging targets are given in Table 3. Measurements of the separation vectors of their components, which are the averages of those measured at J, H, and K with the corrected platescale, are included in the table. The brightness ratios for the visual binaries are plotted in Figure 2.

Out of the 14 visual binary systems (nominally 28 stars) observed, two are found to be well-resolved doubles, 5 to be marginally resolved along a single axis (consistent both with doubles with separations less than the diffraction limit, or with elongated halos), and two were marginally-resolved along all axes. The results for the individual resolved and marginally-resolved sources are discussed below.

### 3.1  Fully-Resolved
Primary-Tertiary Pairs

Two systems (WSB 18 and B59-1) have well-resolved tertiary companions separated from their primaries by $`0`$″.1. Models were fitted to their Fourier phases in order to derive the brightness ratio $`R`$, separation $`\rho `$, and position angle $`\theta `$ of the primary-tertiary pairs. The results are given in Table 4.

The presence of a tertiary companion to the primary in WSB 18 is noteworthy in light of the study by Brandner & Zinnecker (1997), in which the primary appeared somewhat more luminous than would have been predicted for a single star with its M2 spectral type and the $`5\times 10^6`$ yr age derived by comparing the secondary to the theoretical pre-main sequence tracks of D’Antona & Mazzitelli (1994). It is natural to suggest that the tertiary companion is responsible for the excess luminosity, and its detection strengthens the conclusion reached by Brandner & Zinnecker (1997) that the components of the binaries in their sample are coeval.

The projected separations of these primary-tertiary systems correspond to $`16`$ AU at 160 pc. Assuming that the true separations are not much larger, that the mass of each close pair is $`0.6\mathrm{M}_{}`$, and that the orbits are circular, results in an estimate of $`10^2`$ yr for their orbital periods. The systems should rotate by $`4`$ deg yr-1, making their orbital motion readily detectable from one year to the next.

### 3.2  Marginally-Resolved Sources

As discussed above, the nature of the 7
marginally resolved targets is somewhat ambiguous. Five of them appear to be extended along only one axis, making them consistent with both the double-star model (the double having a separation below the diffraction limit) and the halo model (the minor axis being much smaller than the major axis). In these cases, the double-star model is preferred for its astrophysical simplicity, although it may be possible to produce an elongated halo of the observed size by a suitable arrangement of circumstellar material such the inner regions of a nearly edge-on disk. The results of fitting the double-star model to these sources are included in Table 4.

The power spectrum of the Sz 116 secondary is consistent with a double whose separation is just below the diffraction limit. The best-fit model gives a projected separation of 40 mas, corresponding to 6 AU at 160 pc, making this the widest of the potential marginally-resolved doubles in the sample. The secondaries in AS 205, Elias 2-22, ESO H$`\alpha `$ 282, and the primary in B59-2 also have asymmetrical power spectra consistent with the presence of tertiary companions with separations ranging from only 6 to 13 mas, well below the diffraction limit. If the structures seen in these targets do reflect tertiary companions, their periods would be only $`120`$ yr. For these marginally-resolved objects, followup holographic observations to search for significant changes in the position angles of their extended structures would provide a good test of the tertiary-star interpretation.

The primary components in Elias 2-49 and WSB 71 are marginally-resolved but inconsistent with the double-star model because there is no position angle along which the power spectra do not decrease with frequency. In these objects the $`\chi ^2`$ is lower for the halo model than for the double-star model. The data are insufficient to determine the morphology of these sources in any detail, and it may be much more complex than a simple Gaussian halo. The marginally-resolved primaries in these systems are much redder than the unresolved secondaries (see Figure 2), which is not surprising if they are surrounded by dust envelopes, provided that the overall flux is not dominated by scattering. These sources merit followup observations at shorter wavelengths, e.g., visible-light imaging with the Hubble Space Telescope, to search for extended, low surface brightness scattered light to which the present holographic observations are insensitive.

### 3.3  Wavelength Dependence of Visual Binary Separation Vectors

If a star has a close companion with a significantly different infrared color, or if it is surrounded by an asymmetrical halo of non-gray dust, then its photocenter will be wavelength-dependent. A comparison between the separation vectors measured for a visual binary over a range of wavelengths can potentially detect this effect. The RMS scatter in the differences between the IC-in and IC-out separation values for the binaries, excluding Sz 116 (which was not measured at L’) and Elias 2-49, is 0″.014. This gives an indication of the precision of the individual separation measurements. For Elias 2-49 the separation decreases with wavelength from 1″.09 at J to 0″.97 at L’, a change which is much larger change than the estimated precision and suggestively similar to the estimated Gaussian size of the halo surrounding the primary star. Deeper high-resolution images of this system at visible or near-infrared wavelengths would be useful to reveal the structure of the circumstellar material. Interestingly, WSB 71, the other system which appears to have a halo around its primary, has a separation which is not dependent on wavelength between the J and L’ bands.

### 3.4  The Density of Tertiary Companions

If one assumes that all five marginally-resolved sources which are compatible with doubles do have physical tertiary companions at the separations derived from the model fits, then fully half of the nominal visual binaries in the sample are actually hierarchical triples. This perhaps surprising result is discussed below.

The angular surface density $`\mathrm{\Sigma }(\theta )`$ of companions to young stars in the Taurus, Ophiuchus, and Orion star-forming regions was computed by Simon (1997) for companion separations $`\theta `$ as small as 0″.01. In the “binary regime” in which $`\theta `$ is small enough that $`\mathrm{\Sigma }`$ is dominated by bound objects, the surface density in all three clusters was found to be $`\mathrm{\Sigma }(\theta )\theta ^b`$, where b is close to 2.0. A similar conclusion was reached for an X-ray selected sample of Taurus stars by Kohler & Leinert (1998; hereafter KL98) for separations between 0″.13 and 13″. The robustness of this result is impressive in light of range of physical conditions spanned by the three SFRs.

Suppose that this result applies within the present sample of visual binaries for tertiaries at separations small compared to the visual binary separation. This cannot be strictly true, both because the sample is statistically biased by the presence of the visual companions and because the companions’ gravitational influence restricts the range of separations for which a stable tertiary orbit may exist, but given the large ratio of the observed outer to inner separations and the steepness of the density function, it may not be terribly wrong. The orbital stability restriction has been estimated by Harrington (1972) and depends on the inclination of the inner orbit relative to the outer orbit, the eccentricity of the outer orbit, and the mass ratios of the components, but for an outer-orbit eccentricity of 0.5 and components of similar mass it requires the major axis of the outer orbit to be approximately a factor of 10 larger than that of the inner orbit. The average projected separation of the visual binaries in the present sample is 2″, so according to the Harrington criterion with the assumed parameters, and neglecting the effects of projection onto the plane of the sky, the systems on average should permit tertiary companions with separations no larger than about 200 mas. Normalizing $`\mathrm{\Sigma }`$ to produce the number of companions observed by KL98 within the range of radii to which their observations were sensitive, and integrating over an annulus extending from the approximate detection limit of 5 mas to the stability cutoff at 200 mas, leads to a prediction that 11 tertiaries should be found among the 28 primary and secondary stars in the sample. Given the small sample size, the crudeness of the above approximations, and the possibility that some additional tertiary companions were missed, this is consistent with the 7 detected here.

### 3.5  Infrared Brightness Ratios

The brightness ratios for the target visual binaries measured in all four filters, plotted against the filter wavelengths in Figure 2, show no clear tendency for the primary star in a visual binary to be redder than the secondary. Such a tendency would be manifested as a generally negative slope in the brightness-ratio curves. This contrasts with the result for the sample of Taurus visual binary TTS observed by Moneti & Zinnecker (1991), in which the primary stars were redder than the secondaries. An examination of the color-color plot in Figure 3 shows no clear correlations between the colors of the visual binary components, their structure in the holographic observations, or their relative brightness.

## 4  Conclusions

The present study is among the first to probe a sample of pre-main sequence visual binary systems with this combination of high resolution and dynamic range in the near-infrared. The results are not very consistent with the naive expectations that motivated it. More than half of the systems contain resolved or marginally-resolved components. Two tertiary companions were clearly resolved at $`100`$ mas separations from their primaries, and tertiaries at separations below the diffraction limit are likely to account for the structure seen in all but two of the marginally-resolved stars, including the infrared companion (IRC) in Elias 2-22. Followup observations to search for orbital motion would be useful. The presence of such a large number of tertiary companions would make it necessary to search binary systems carefully before attempting to use them to test pre-main sequence tracks, since the light of an undetected tertiary will make a star appear younger and more luminous that it is. Two stars showed extended structure which appears to be due to the presence of dust halos rather than tertiary companions, and one of these (Elias 2-49) has a separation which decreases with wavelength. Strikingly absent is a clear detection of an edge-on circumstellar disk in any of the binaries in the sample.

It is a pleasure to thank M. Creech-Eakman for her assistance with the observations, G. Blake, P. Goldreich, and J. Kwan for useful discussions, E. Gaidos for useful comments on a version of the manuscript, and to an anonymous referee for an exceptionally useful report. Data presented herein were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation. The author wishes to extend special thanks to those of Hawaiian ancestry on whose sacred mountain he was privileged to be a guest. Without their generous hospitality, and the hard work and dedication of the Keck Observatory staff, the observations presented herein would not have been possible. This research was supported by the National Aeronautics and Space Administration.


text:
# On the large-scale inhomogeneous Universe and the cosmological constant

We study the large-scale inhomogeneity of the Universe based on the averaging procedure of Buchert and Ehlers. The generalized Dyer-Roeder equation for the angular diameter distance of the inhomogeneous Universe is derived and solved for different cosmological models. We make a comparison of certain cosmic observables, such as the Hubble function, angular diameter distance, cosmological correction factor of homogeneous and inhomogeneous cosmological models, which are crucial ingredients in galaxy number counts and gravitational lenses.

It seems that present cosmological data consisting of supernova searches, structure formation, gravitational lenses, etc. suggest a rather unexpected equation of state for the simple homogeneous Friedmann-Lema$`\widehat{i}`$tre (FL) Universe with $`\mathrm{\Omega }_m0.3,\mathrm{\Omega }_\mathrm{\Lambda }0.7.`$ However, one should be cautious even when performing analyses of the cosmological data (for example, for supernova data see the reanalysis in ). Further concern is the low statistical confidence of all extracted observables. Nevertheless, this phenomenological result initiated a lot of speculative work on the origin of the positive cosmological constant.

Contrary to current investigations, there is only one derivation of the cosmological constant without any fine tuning or any obscure assumptions and this derivation is based on the Einstein-Cartan gravity . When the Universe is frozen at zero temperature $`T_\gamma =0`$ (spacelike infinity) and if all nonrelativistic fermionic matter is spinning matter, then owing to the same coupling between spacetime curvature and mass density on the one side and spacetime torsion and spin density of matter on the other side, the following relations emerge: $`\mathrm{\Omega }_m=2,\mathrm{\Omega }_\mathrm{\Lambda }=1`$ . It is important to strengthen that at $`T_\gamma =0`$ all spin densities of nonbaryonic and baryonic species act coherently to the total spin. The role of spin densities in the derivation of the primordial density fluctuation is described in Ref. . The existence of vorticity and acceleration  or nonvanishing shear  should be of fundamental importance for resolving certain cosmological problems, but the cosmological distance measures are not very much affected by relatively small deviations from the FL model . To conclude, let us notice that it is possible to construct a nonsingular and nonanomalous gauge theory with heavy and light neutrinos as cold and hot dark matter fermionic particles respectively, but without Higgs scalars and with no asymptotoic freedom in QCD .

In this paper we want to reconcile current measurements with large-scale inhomogeneous cosmic models. We assume that the Universe is inhomogeneous on large scales, thus clumpiness should not be characteristic only of small-scale structures, and that the Earth is placed in the region away from the centre of homogeneity.

Starting with the historical work of Lema$`\widehat{i}`$tre, Tolman and Bondi (LTB) on inhomogeneous models, one can account enormous activity in this field . However, a recent attempt  to exploit LTB models, following perturbative calculations in Ref. , was not successful at large redshifts.

The cellular structure of the Universe with a power-law distribution of matter of Ruffini et al , represents another approach to the problem of inhomogeneity.

The work of Zalaletdinov  is the first complete and consistent treatment of the averaged Einstein equations in an arbitrary Riemannian spacetime. However, even if the vorticity and acceleration do not vanish, they are small deviations (which we neglect in this paper) from the FL geometry, thus we choose the following line element and spatially averaged scalar quantites of Buchert and Ehlers :

$`ds^2=dt^2g_{ij}dX^idX^j,`$
$`\mathrm{\Psi }(t,X^i)_𝒟{\displaystyle \frac{1}{V_𝒟}}{\displaystyle _𝒟}d^3XJ\mathrm{\Psi }(t,X^i),`$
$`V_𝒟(t){\displaystyle _𝒟}Jd^3X,J(t,X^i)[det(g_{ij})]^{\frac{1}{2}}.`$
Introducing natural definitions

$`a_𝒟(t)({\displaystyle \frac{V_𝒟(t)}{V_{𝒟_0}}})^{\frac{1}{3}},`$
$`\theta _𝒟={\displaystyle \frac{\dot{V}_𝒟}{V_𝒟}}=3{\displaystyle \frac{\dot{a}_𝒟}{a_𝒟}},\dot{\psi }(t){\displaystyle \frac{d}{dt}}\psi (t),`$
$`M_𝒟={\displaystyle _𝒟}d^3XJ\rho (t,X^i),`$
$`\rho _𝒟={\displaystyle \frac{M_𝒟}{V_{𝒟_0}a_𝒟^3}}.`$
one can easily derive Raychaudhuri and constraint equations for averaged scalars from Einstein and conservation equations

$`3{\displaystyle \frac{\ddot{a}_𝒟}{a_𝒟}}+4\pi G_N{\displaystyle \frac{M_𝒟}{V_{𝒟_0}a_𝒟^3}}\mathrm{\Lambda }=Q_𝒟,`$ (1)
$`3({\displaystyle \frac{\dot{a}_𝒟}{a_𝒟}})^28\pi G_N{\displaystyle \frac{M_𝒟}{V_{𝒟_0}a_𝒟^3}}+{\displaystyle \frac{1}{2}}_𝒟\mathrm{\Lambda }={\displaystyle \frac{1}{2}}Q_𝒟,`$ (2)
$`Q_𝒟{\displaystyle \frac{2}{3}}(\theta \theta _𝒟)^2_𝒟2\sigma ^2_𝒟,`$
$`=_i^i,_{ij}spatialpartoftheRiccitensor.`$
In addition, we need a condition of integrability of these equations :

$`(a_𝒟^6Q_𝒟)^.+a_𝒟^4(a_𝒟^2_𝒟)^.=0.`$ (3)
Now we turn to the geometrical optics in order to derive the wanted observables of the inhomogeneous cosmic models.

Starting from the definition of the angular diameter distance and the focusing equation :

$`D({\displaystyle \frac{dA_S}{d\mathrm{\Omega }_O}})^{\frac{1}{2}},`$
$`{\displaystyle \frac{d^2}{dw^2}}A_S^{\frac{1}{2}}+{\displaystyle \frac{4\pi G_N}{H_0^2}}(1+z)^2\rho A_S^{\frac{1}{2}}=0,`$ (4)
$`wH_0\omega _0v,v=affineparameter.`$
We take into account that the shear of the congruence of photons vanishes for a spherically symmetric source . Then, performing the Buchert-Ehlers spatial averaging and neglecting the spacetime shear in the averaged equation for the change of redshift with respect to the affine parameter:

$`{\displaystyle \frac{dz}{dw}}=(1+z)^3{\displaystyle \frac{\dot{a}_𝒟(z)}{\dot{a}_𝒟(0)}},`$ (5)
we arrive at the following equation for the averaged angular diameter distance:

$`{\displaystyle \frac{d^2D_𝒟}{dz^2}}({\displaystyle \frac{dz}{dw}})^2+{\displaystyle \frac{dD_𝒟}{dz}}{\displaystyle \frac{d^2z}{dw^2}}+{\displaystyle \frac{4\pi G_N}{H_0^2}}(1+z)^2\rho D_𝒟=0.`$
Assuming the smooth space dependence of mass density and angular diameter distance, one can factorize the last term of the preceding equation as

$`\rho D_𝒟\rho (t,\overline{X^i})D_𝒟{\displaystyle \frac{M_𝒟}{V_𝒟}}D_𝒟,\overline{X^i}ϵ𝒟,`$
and we are left with the second-order generalized Dyer-Roeder equation with a suitable density normalization for the cosmological models:

$`\varphi (z){\displaystyle \frac{d^2\stackrel{~}{D}(z)}{dz^2}}+\chi (z){\displaystyle \frac{d\stackrel{~}{D}(z)}{dz}}+\eta (z)\stackrel{~}{D}(z)=0,`$ (6)
$`\stackrel{~}{D}D,\varphi (z)=({\displaystyle \frac{dz}{dw}})^2,\chi (z)={\displaystyle \frac{d^2z}{dw^2}},`$
$`{\displaystyle \frac{d^2z}{dw^2}}=({\displaystyle \frac{dz}{dw}})^2({\displaystyle \frac{3}{1+z}}+{\displaystyle \frac{1}{\dot{a}_𝒟(z)}}{\displaystyle \frac{d\dot{a}_𝒟(z)}{dz}}),`$
$`{\displaystyle \frac{d\dot{a}_𝒟}{dz}}={\displaystyle \frac{1}{H(z)(1+z)}}\ddot{a}_𝒟,`$
$`H(z){\displaystyle \frac{\dot{a}_𝒟(z)}{a_𝒟(z)}},a_𝒟(z)={\displaystyle \frac{1}{1+z}},`$
$`\eta (z)={\displaystyle \frac{3}{2}}(1+z)^5\mathrm{\Omega }_m,\mathrm{\Lambda }=3H_0^2\mathrm{\Omega }_\mathrm{\Lambda }.`$
Under adequate initial conditions  (Adams-Bashforth integration method used)

$`forintegrationfromz_1toz_2,z_2>z_1,`$
$`\stackrel{~}{D}(z_1,z_1)=0,{\displaystyle \frac{\stackrel{~}{D}(z_1,z)}{dz}}(z=z_1)={\displaystyle \frac{1}{H(z_1)}}{\displaystyle \frac{1}{1+z_1}},`$
we can calculate the averaged angular diameter distance for an arbitrary cosmological model $`\mathrm{\Omega }_m+\mathrm{\Omega }_\mathrm{\Lambda }=1`$. Needless to say, for the vanishing backreaction and curvature terms one recovers a homogeneous solution:

$`D(z)={\displaystyle \frac{1}{H_0(1+z)}}{\displaystyle _0^z}{\displaystyle \frac{d\zeta }{[\mathrm{\Omega }_m(1+\zeta )^3+\mathrm{\Omega }_\mathrm{\Lambda }]^{\frac{1}{2}}}}.`$
Numerical evaluations and comparisons are performed in the following manner: (1) choose some functional form for the backreaction Q(z), (2) evaluate curvature term R(z) from the integrability condition, (3) choose two inhomogeneous models with $`\mathrm{\Omega }_m=2`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=1`$ \[Einstein-Cartan (EC) model\] and $`\mathrm{\Omega }_m=1`$ , $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$ \[Einstein-deSitter (EdS) model\] and fit parameters of the backreaction and curvature to reach the input value for the Hubble constant $`H_0=\dot{a}(z=0)/a(z=0)`$ and to get the best fit of the $`\stackrel{~}{D}(z)`$ of the inhomogeneous models to the D(z) of the homogeneous model with
$`\mathrm{\Omega }_m=0.3`$ , $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ for $`z0.7`$.

One can easily check that for any $`0>\alpha >3`$ (we omit index $`𝒟`$)

$`Q(z)=Q_0a^\alpha (z),`$
$`_𝒟R(t)=a^2(t){\displaystyle _{t_0}^t}𝑑\tau a^4(\tau ){\displaystyle \frac{d}{d\tau }}[a^6(\tau )Q(\tau )],`$
$`R(z)={\displaystyle \frac{\alpha +6}{\alpha +2}}[Q(z)({\displaystyle \frac{a(z_0)}{a(z)}})^2Q(z_0)],\alpha 2,`$
$`R(z)=4Q(z)\mathrm{ln}{\displaystyle \frac{1+z}{1+z_0}},\alpha =2,`$
the curvature term power-law coefficient is within the same interval as that for the backreaction. Thus, this simple functional form garantees the homogeneity at very large scales $`Q(z),R(z)H_0^2\mathrm{\Omega }_ma_𝒟(z)^3,z>>1`$(it could also be worth including the Gaussian damping, if necessary).

For definite numerical comparisons of the three cosmological models we use this Ansatz:

$`Q(z)=Q_0[1+c_0a^1(z)],`$
$`R(z)=Q_0(1+z)^2[3({\displaystyle \frac{1}{(1+z)^2}}{\displaystyle \frac{1}{(1+z_0)^2}})+5c_0({\displaystyle \frac{1}{1+z}}{\displaystyle \frac{1}{1+z_0}})].`$
In order to match the input Hubble parameter at present time $`H_0`$ , $`z_0`$ has to fulfil the following equation:

$`Q(0)+R(0)=0`$
$`2(1+3c_0)(1+z_0)^25c_0(1+z_0)+3=0.`$
Searching for the best fit values of $`c_0`$ and $`Q_0`$, we find

$`H_0h_0u,u100kms^1Mpc^1,`$
$`(EC):h_0=0.6,Q_0=3.789u^2,z_0=\sqrt{6}/21,c_0=0,`$
$`(EdS):h_0=0.6,Q_0=1.2632u^2,z_0=\sqrt{6}/21,c_0=0,`$
$`(Hom):h_0=0.6,\mathrm{\Omega }_m=0.3,\mathrm{\Omega }_\mathrm{\Lambda }=0.7,Q_0=0.`$
For $`z0.5`$, the angular diameter distances differ less than $`5\%`$. Figs. 1 and 2 show the Hubble flow and the angular diameter distance for the three models, while Table 1 presents the cosmological correction terms for the difference in light travel time between gravitational lens images :

$`H_0t=Tf(\theta _{1,Obs},\theta _{2,Obs},z_d,z_s),`$
$`T=H_0{\displaystyle \frac{\stackrel{~}{D}^2(0,z_d)}{\stackrel{~}{D}}}(1+z_d){\displaystyle \frac{z_sz_d}{z_sz_d}},`$
$`\stackrel{~}{D}={\displaystyle \frac{\stackrel{~}{D}(z_d,z_s)\stackrel{~}{D}(0,z_d)}{\stackrel{~}{D}(0,z_s)}}.`$
|  | $`z_d=0.5,z_s=1`$ | $`z_d=0.5,z_s=1.5`$ | $`z_d=1,z_s=1.5`$ | $`z_d=1,z_s=2`$ |
| --- | --- | --- | --- | --- |
| T(EC) | 0.7550 | 0.6686 | 0.4977 | 0.4226 |
| T(EdS) | 0.9559 | 0.9350 | 0.8775 | 0.8471 |
| T(Hom) | 1.8543 | 2.1098 | 4.7948 | 8.9155 |

The galaxy number count depends essentially on the cosmological model through the comoving volume:

$`{\displaystyle \frac{dV(z)}{dz}}\stackrel{~}{D}(z)^2.`$
We conclude with two observations: (1) For the large-scale inhomogeneity of the Universe to be established one needs clear indications that none of homogeneous models can simultaneously fit and explain all data in cosmography, structure evolution, gravitational lenses, etc. for all redshifts. At present it is difficult to make any conclusion what is the source of some recently found disparities in gravitational lenses  or galaxy evolution . (2) If the inhomogeneity is established, then one could intend to model a backreaction function for the global fit of data, and ultimately attempt to derive it from the evolution of multicomponent imperfect fluid with photons, baryonic and nonbaryonic matter, with a necessary knowledge of all relevant masses and couplings.

* * *
This work was supported by the Ministry of Science and Technology of the Republic of Croatia under Contract No. 00980103.


text:
# Electroweak properties of the nucleon in a chiral constituent quark model

## Abstract

Results for all elastic electroweak nucleon form factors are presented for the chiral constituent quark model based on Goldstone-boson-exchange dynamics. The calculations are performed in a covariant framework using the point-form approach to relativistic quantum mechanics. The direct predictions of the model yield a remarkably consistent picture of the electroweak nucleon structure.

We present results for the elastic electroweak nucleon observables as a progress report of a more comprehensive programme aiming at a consistent description of the electroweak properties of baryons at low energy. The theoretical context is represented by the chiral Constituent Quark Model (CQM) based on the Goldstone-boson exchange (GBE) quark-quark interaction, that is induced by the spontaneous chiral symmetry breaking in QCD and that accurately reproduces the baryon spectrum of light and strange flavors . The dynamics of quarks inside the nucleon is essentially relativistic. Therefore, we have adopted the point-form realization of relativistic quantum mechanics, where the boost generators are interaction-free and make the theory manifestly covariant . The electromagnetic photon-quark interaction is assumed point-like, but in point-form the momentum delivered to the nucleon is different from the one delivered to the struck quark; hence, we will name this approach the Point-form Spectator Approximation (PFSA) . The quark wave functions deduced from fitting the baryon spectrum are used as input and no further parameter is introduced, since quarks are considered point-like and the point-form allows for an exact calculation of all boosts required by a covariant description. Results have recently been published for electromagnetic , axial  and pseudoscalar  nucleon form factors. They are summarized here in Figs. 1-2, and in Tab. 1.

The agreement with experimental data is remarkable and it indicates that by a proper choice of low-energy degrees of freedom a quark model is capable of describing the spectroscopy and the low-energy dynamics of baryons at the same time. However, a more detailed comparison with data shows that there is still room for quantitative improvements, e.g. by considering two-body electromagnetic current operators and constituent quark sizes.


text:
# Manifestation of photonic band structure in small clusters of spherical particles

## Abstract

We study the formation of the photonic band structure in small clusters of dielectric spheres. The first signs of the band structure, an attribute of an infinite crystal, can appear for clusters of 5 particles. Density of resonant states of a cluster of 32 spheres may exhibit a well defined structure similar to the density of electromagnetic states of the infinite photonic crystal. The resonant mode structure of finite-size aggregates is shown to be insensitive to random displacements of particles off the perfect lattice positions as large as half-radius of the particle. The results were obtained by an efficient numerical method, which relates the density of resonant states to the the scattering coefficients of the electromagnetic scattering problem. Generalized multisphere Mie (GMM) solution was used to obtain scattering matrix elements. These results are important to miniature photonic crystal design as well as understanding of light localization in dense random media.

The idea of employing photonic crystals, structures with periodically modulated dielectric constant, to manipulate the density of states (DOS) of electromagnetic waves (EM) Yablonovitch (1987); John (1987) led to an explosion of interest, both academic and practical Johnson and Joannopoulos (2002); vla . In a photonic crystal, in a given crystallographic direction, light can propagate only for the frequencies within photonic bands, described by a dispersion relation $`\omega =\omega (\text{k})`$, similar to electron de Broglie waves in crystals. If the gaps between the bands overlap for all directions in the crystal, photonic crystal posses a complete photonic band gap Ho et al. (1990). For frequencies inside the complete gap, density of EM states turns to zero, that leads to new phenomena such as suppression of spontaneous emission, light localization, zero-threshold lasers, all-optical transistors and circuitry, and anomalous nonlinear propertiesYablonovitch (1987); John (1987); Johnson and Joannopoulos (2002); vla ; Soukoulis (2001). For theoretical considerations, the size of the crystal is usually assumed infinite. In practice, however, one always deals with finite structures. Naturally, a question arises: How large a photonic crystal should be, in order to exhibit photonic band structure? In this letter, we show that the modification of the density of EM states may begin for a cluster of 5 spheres. This suggests that scattering of light by ordered aggregates of spheres can be important in the problem of light propagation in dense random media. Indeed, an ensemble of mono-disperse, resonant scatterers (e.g. spheres) is considered the most favorable Lagendijk and van Tiggelen (1996) for satisfying Ioffe-Regel criterion for Anderson localization $`kl_{scat}1`$ Sheng (1995); Lagendijk and van Tiggelen (1996); Soukoulis (2001); Chabanov and Genack (2001); Chabanov et al. (2000); Wiersma et al. (1997), where $`k`$ is the wavenumber in the medium and $`l_{scat}`$ is the scattering length. Closely-packed face centered cubic (FCC) arrangement often appear in self-assembled structures and it is shown to have lowest free energy. This fact points to a high probability of formation of ordered clusters in a collection of spheres. Even though the complete photonic gap can never be realized in an FCC structure of dielectric spheres in air, for some parameters, the total density of EM states can be significantly suppressed in certain frequency regions. Our results suggest that such suppression can occur due to the presence of the clusters even without long range order.

Treatment of photonic structures, with all 3 dimensions of the order of the wavelength has been a challenge. When the distance between scatterers is comparable to their size, analytical methods, such as single scattering approximation, dipole approximation become invalidLagendijk and van Tiggelen (1996). Numerical methods tailored for calculation of transmission through a slab Stefanou et al. (2000); Pendry and Mackinnon  cannot be applied because it requires the lateral size of the slab to be larger than all characteristic lengths. Finite difference time domain (FDTD) method is used extensively in photonics design A. Taflove (2000); Painter et al. (1999). Obtaining the density of states with FDTD-based ONYX method Pendry and Mackinnon  requires time-dependent solution. Combined with dense $`\lambda /20`$ spacial grids, it makes the problem computationally demanding in 3D, even for the smallest structure considered in this paper. However, the tendency to miniaturization of photonic devices calls for an efficient tool for dealing with such structures. An efficient alternative approach is needed.

For small truncated photonic crystal structures, scattering language becomes appropriate. We use the GMM solution l. Xu (1997) to obtain the elements of the scattering matrix of the collection of dielectric spheres. The density of resonant EM states of the cluster can be easily expressed in terms of the scattering matrix. Within this approach, one only needs to evaluate vector-spherical-function expansion coefficients, there is no need, as in FDTD method, to find EM fields at every spacial point. The method requires considerably less memory and computation time resources, that makes it possible to use on a personal computer for up to 300 particles. An important class of photonic crystals, opals, is treated in this framework as an example.

The main results of this paper are: (i) the signatures of the photonic band structure can appear for aggregates as small as 5 particles; (ii) in a cluster of 32 particles the density of (resonant) states can show strong resemblance to that of an infinite structure, with a pronounced depletion in the region of pseudo gap; (iii) photonic band structure of the finite clusters is tolerable to a certain degree of random displacements of the particles off their lattice sites.

We consider clusters of 5 and 32 spheres, with diameter $`0.95`$ cm and refractive index $`n=3.14`$. For both clusters we assumed FCC ordering with the distance between nearest neighbors of 1.9 cm. The system was chosen to match the setup used in the light localization experiment – Ref. Chabanov and Genack, 2001. Smaller cluster is comprised of three planar layers perpendicular to the 111 direction, with 1, 3, and 1 sphere in each layer. In the larger cluster, particles were arranged in 5 layers with 3, 7, 12, 7, and 3 spheres in each layer, with the same stacking as in the first cluster.

The cluster is an open system that does not support stable modes, quasi-states of EM field are leaky, and density of the states cannot be specified. For electrons Fyodorov and Sommers (1997) in open stochastic systems, Wigner delay time, $`\widehat{\tau }_w(E)=\mathrm{}/iln(\widehat{\text{S}})/E`$, expressed in terms of scattering matrix $`\widehat{\text{S}}`$, can be shown to be proportional to the density of resonant modes

$$\tau _w(\omega )\underset{j}{}\frac{\gamma _j/2}{(\omega \omega _j)^2+(\gamma _j/2)^2},$$
(1)
where $`\omega _j`$ and $`\gamma _j`$ are the energy and the linewidth of the $`j`$th resonance. For an infinite system, the poles of the scattering matrix $`\widehat{\text{S}}`$ approach the real energy axis and the stationary modes are formed, $`\gamma _j`$ become zero. For EM waves, one can define dwell time, $`\tau _d`$, that closely follows EM version of the Wigner time, giving asymptotically the same result on the resonances Lagendijk and van Tiggelen (1996). Dwell time can be readily found in terms of scattering coefficients as follows:

$$\tau _d(\omega )=\frac{1/N_i^N_{S_i}W(\text{r},\omega )𝑑\text{r}}{c\sigma _{scat}(\omega )},$$
(2)
where $`W(\text{r},\omega )`$ is electromagnetic energy normalized to 1 in vacuum, $`\sigma _{scat}`$ is the scattering cross-section, normalized by the total geometrical cross-section $`_i^N\pi r_i^2`$, and $`c`$ is speed of light. The integral in Eq. (2) is taken over the volume of each scatterer, $`S_i`$. It gives $`\tau _d`$ in units of $`r/c`$ and numerically coincides with cavity quality factor $`Q`$. The above expression can be understood physically as a ratio between the energy stored in a scatterer divided by outgoing energy current. $`_{S_i}W(\text{r},\omega )𝑑\text{r}`$ can be found by introducing an infinitesimally small absorption in the refractive index of the spheres $`n+i\kappa `$Lagendijk and van Tiggelen (1996)

$$_{S_i}W(\text{r},\omega )\frac{3}{8}n\underset{\kappa 0}{lim}\frac{Q_{absi}}{x_i\kappa },$$
(3)
where $`Q_{absi}`$ is the absorption cross-section, and $`x_i=2\pi r_i/\lambda `$ is size parameter of the $`i`$th sphere. The problem of finding quasi-DOS of the finite cluster is now reduced to finding the scattering coefficients $`Q_{absi},\sigma _{scat}`$. This step can be done by evaluating the expansion series of GMM solution l. Xu (1997). To obtain a measure of the total density of resonant states $`\tau _d^{tot}`$, expression Eq.(2) was averaged over 9 symmetrically non-equivalent orientations (altogether 108 directions, accounting symmetry). We would like to stress, that within this framework, the source of error are numerical accuracy of the evaluation of expansion coefficients and limited number of orientations used in angular averaging, both of them can be controlled.

One characteristic property of a photonic crystal is spacial dispersion, $`\omega (\text{k})`$. On Fig. 1 we compare dispersion curves of the infinite photonic crystal, obtained using MIT photonic bands code Johnson and Joannopoulos (2001), to the dwell times for two clusters, described above, calculated for 2 different incident angles corresponding to $`X`$ and $`L`$ crystallographic directions. It can be easily seen, that even for smaller cluster of 5 particles $`\tau _d`$ shows the formation of the resonant modes at the positions of photonic modes in $`L`$-direction. For the cluster of 32 particles, well defined mode structure appears in both directions. However it is more pronounced for $`L`$ direction, where 6 resonant states between 12 and 16 GHz can be traced to $`L`$-modes of the infinite photonic crystal. The shift of the modes can ascribed to the coupling between different modes owing to the finite size of the system. Indeed, in the infinite structure the conservation of momentum of EM wave, due to periodicity, would decouple $`L`$-modes, whereas the finite coupling in truncated crystal leads to the shift and widening of the resonances. Based on symmetry considerations, we can understand the stronger $`L`$-modes as compared to $`X`$-modes. The former are formed due to Bragg reflections from 3 or 5 (for the smaller and larger clusters respectively) planes formed by adjacent (nearest) neighbors in 111 direction, while in $`X`$ direction, Bragg planes are formed only by next-nearest neighbors, which would make them more susceptible to the truncation.

$`\tau _d^{tot}`$ and angularly (as for $`\tau _d^{tot}`$) averaged $`\sigma _{scat}`$ in the region of the first two single sphere Mie resonances are shown on Figs. 2(b,c). For comparison we included the total density of EM states, Fig. 2(a), calculated for the infinite crystal using the method described in Ref. Busch and John, 1998. Judging from scattering cross-section only, one cannot compare a finite cluster to the infinite structure, while $`\tau _d^{tot}`$ allows straightforward comparison with DOS of the photonic crystal

Fig. 2(b) shows twofold decrease of $`\sigma _{scat}`$ at Mie resonances in the cluster. This suppression of the scattering efficiency can be attributed to the hybridization of the single particle resonances Moroz (1999) due to multiple scattering. Dwell time is closely related to the near-fields Lagendijk and van Tiggelen (1996) inside the clusters, where formation of the Bragg standing modes leads to substantial modification of the spectrum – Fig. 2(c). For 32-sphere cluster one can see strong resemblance between the density of the resonant modes expressed by $`\tau _d^{tot}`$ and DOS in the photonic crystal, Fig. 2(a). It is interesting to compare our results to 1D case considered in Ref. Deych and Lisyansky, 2000. For a stack of periodically arranged quantum wells it was argued that with increase of the system size, so-called, subradiant EM modes formed stable modes of the photonic passband with large life-times, similar to our result for dwell times in the 3D clusters.

To assess the sensitivity of our results to disorder in the cluster, we performed the calculations of $`\sigma _{scat}`$ and $`\tau _d^{tot}`$, Fig 2(d,e), for: the cluster of 32 spheres displaced by the half of the particle radius in random directions (dashed lines) off their initial positions and the cluster of 32 spheres where 5 spheres replaced by the defect spheres of a different diameter (0.85 cm) - solid lines. One can see, that two types of disorder had different effect on the spectra. While positional (topological) disorder led only to smearing of only sharp resonances, the defects also introduced a number of new peaks. The later effect can be related to the new resonances introduced by the defect spheres. The effect of the topological disorder on photonic bandgaps was studied in (finite but comprised of large number of particles) 2D Sigalas et al. (1996); Asatryan et al. (2000) and 3D photonic crystals Li and Zhang (2000). Our conclusion on stability of the band structure to the topological disorder seems to be in line with Ref. Sigalas et al., 1996. However, besides obvious difference in the dimensionality of the space (we work in 3D), another important difference exists. In large systems, to explain the existence of the band structure one can resort to the periodic-on-average argument, while in such a small system as ours one cannot.

It is worth noting, that FCC opal structures are not usually considered as a good candidate for potential applications connected to complete photonic bandgap Leung and Liu (1990). However, as it can be seen from our example, for some parameters, DOS can be significantly suppressed in a wide spectral region. Fig. 2 also suggests that such structures can be tolerable to strong topological disorder.

As it was already mentioned in the introduction, our results are relevant for the problem of light propagation in a dense random media made of particles with a narrow size dispersion. Topological disorder does not prevent the occurrence of ordered clusters. In Ref. Chabanov and Genack, 2001, random displacement of scatterers off their lattice positions did not have a significant effect on the density of EM states. Our calculations performed for the clusters with similar random displacements (Fig. 2(d,e)), confirm small effect on the formation of the photonic band structure. This also suggests that the light localization observed in Ref. Chabanov and Genack, 2001 may have been facilitated by depletion of the density of electromagnetic states due to the presence of the pseudogap.

It can be further argued that when the scattering length $`l_{scat}`$ becomes of the same order as the mean distance between the particles, the occurrence of clusters cannot be ignored. Indeed, the scattering inside the clusters is dominated by Bragg reflection mechanism, while $`l_{Bragg}`$ can be as small as a few lattice constants for high index contrast structures. Furthermore, in Ref. Huang et al. (2001) it was shown that periodicity of a photonic crystal may also exhibit itself in coherent back-scattering effect. The modification of coherent back-scattering due to the presence of ordered clusters is an interesting problem that deserves a detailed consideration.

In conclusion, we calculated scattering cross-section and density of resonant states in small 3D ordered aggregates of dielectric spheres within framework of the rigorous GMM solution. In contrast with previous studies we made no assumptions about the strength, size, or separation between scatterers. All multiple scattering effects were automatically retained in the solution. The results suggest that the photonic band structure of the infinite crystal can show up for clusters as small as 5 particles. Density of resonant states is significantly perturbed compared to single particle case. For a cluster of 32 particle the density of resonant states can show a close resemblance to DOS of the infinite system.

###### Acknowledgements.

AY gratefully acknowledges numerous insightful discussions and comments by A. A. Chabanov, A. A. Lisyansky, and L. I. Deych. This work is supported by the National Science Foundation under Grant No DMR-0093949. HC acknowledges the support from the David and Lucile Packard Foundation and Alfred P. Sloan Foundation.