text:
# Limits on isocurvature fluctuations from Boomerang and MAXIMA

## I Introduction

Two recent balloon-borne experiments, Boomerang and Maxima-1, have lent considerable support to the inflationary paradigm. The Boomerang experiment, which was flown over Antarctica for 10 days in 1999 and which measured the temperature fluctuations over 1% of the microwave sky, was the first to provide firm evidence for the first (and possibly second) acoustic peak. The angular power spectrum, obtained from a preliminary analysis of the Boomerang data, is compatible with approximate scale invariance predicted by cosmic inflation. The first acoustic peak is at the multipole $`l200`$, implying a flat universe with $`\mathrm{\Omega }\mathrm{\Omega }_m+\mathrm{\Omega }_\mathrm{\Lambda }1`$ in the case of adiabatic fluctuations . The Boomerang results have been confirmed by Maxima-1, which observed a 124 square degree patch of the sky. Boomerang and Maxima-1 are the first precision measurements of the cosmic microwave background (CMB) temperature fluctuations at multipoles $`l30`$. They will be followed by the second generation satellite experiments Microwave Anisotropy Probe (MAP)  and Planck.

The simplest one-field inflation model typically predicts a near scale invariant power spectrum with Gaussian fluctuations. The perturbations generated by quantum fluctuations are adiabatic, with the number density proportional to entropy density so that $`\delta (n/s)=0`$. This is so because the quantum fluctuations of the inflaton field are directly reflected in perturbations of the inflaton energy density.

More generally, in addition to adiabatic fluctuations, there is also another fundamental fluctuation mode. It does not perturb the total energy density on comoving hypersurfaces (orthogonal to fluid flow), so that $`\delta \rho =0`$. Since there is then no curvature perturbation of the comoving hypersurfaces, it is called an isocurvature fluctuation. Such fluctuations can arise in particle physics models where, during inflation, there appear also other low mass fields in addition to the inflaton. For instance, for a complex scalar field with a U(1) symmetry the fluctuations of the phase do not affect the energy density but would show up as perturbations in the conserved charge associated with the symmetry. This is effectively the case in the minimally supersymmetric model, and its extensions, which have several flat directions in the potential by virtue of gauge invariance and supersymmetry. The complex fields along these flat directions, called Affleck-Dine fields, will be subject to quantum fluctuations similar to the inflaton. As the Affleck-Dine fields typically carry global charges such as baryon number, in addition to adiabatic perturbations, there will also be isocurvature perturbations. The Affleck-Dine condensate itself is not stable but fragments into non-topological solitons, or Q-balls. These may be absolutely stable, in which case there will be isocurvature fluctuations in the baryon number, or unstable but long-lived, in which case isocurvature fluctuation can be imprinted also on the cold dark matter (CDM) particles. Other particle physics motivated sources for isocurvature fluctuations are axion models  or models of inflation with more than one field.

For isocurvature perturbations the overdensities in a given particle species are balanced by perturbations in other particle species, such as radiation. At the last scattering surface (LSS) the compensation for the isocurvature perturbations can be maintained only for scales larger than the horizon, effectively generating extra power to photon perturbations at small multipoles $`l`$. As a consequence, the CMB angular power spectrum, $`C_l`$, of isocurvature perturbations differs a great deal from adiabatic perturbations, and a purely isocurvature CDM perturbation with a scale free spectrum is clearly ruled out  on the basis of the Cosmic Background Explorer (COBE) normalization and $`\sigma _8`$, the observed amplitude of the rms mass fluctuations in an $`8h^1`$ Mpc sphere.

To match $`\sigma _8`$ with the COBE normalization requires a large “blue” tilt $`n_{\mathrm{iso}}2.2`$. Such a tilted isocurvature model can be obtained from a reasonable inflation model and agrees with most observational data. (The main difficulty before Boomerang was that it has more tilt at low multipoles than what COBE saw.) Thus it could be argued that a pure isocurvature model was still acceptable.

A mixed perturbation, with a small isocurvature component is certainly allowed and indeed motivated by particle physics, although the forthcoming satellite experiments are expected to be able to constrain the isocurvature amplitude with a high precision.

The purpose of the present paper is to find out the constraints on the isocurvature component as implied by the Boomerang and Maxima-1 data. The isocurvature perturbation is highly degenerate with the tensor perturbation, and for their separation one probably has to wait for the polarization data from Planck. Here we will make the simplifying assumption that the tensor contribution is negligible. This will be the case in so-called small-field inflation models. For large-field inflation, the tensor/scalar ratio is $`rC_2^T/C_2^{\mathrm{ad}}7(1n_{\mathrm{ad}})`$, so that the tensor contribution disappears for a scale-free spectrum, $`n_{\mathrm{ad}}1`$. As we shall see, such an assumption is consistent with the Boomerang and Maxima-1 data.

An isocurvature contribution is also somewhat degenerate with reionization. As we are looking for the upper limit to an isocurvature contribution, we make the simplifying assumption that reionization has a negligible effect.

## II Adiabatic vs isocurvature perturbations

Adiabatic perturbations are characterized by the quantity $``$, which is related to the curvature perturbation of the comoving hypersurface,

$$_k=\frac{1}{4}\left(\frac{a}{k}\right)^2R_k^{(3)},$$
(1)
where $`a`$ is the scale factor, $`k`$ is the comoving wave number, and $`R^{(3)}`$ is the scalar curvature of the hypersurface (for $`\mathrm{\Omega }=1`$, the unperturbed hypersurface is flat). For adiabatic perturbations $``$ is independent of time while outside the horizon.

Isocurvature perturbations are characterized by the entropy perturbation $`S\delta (n_c/s_\gamma )/(n_c/s_\gamma )=\delta _c(3/4)\delta _\gamma `$, where $`n_c`$ is the number density of CDM, $`s_\gamma `$ is the entropy density associated with photons, and $`\delta _c`$ and $`\delta _\gamma `$ are the relative overdensities in the CDM and photon energy densities. During radiation domination $`S\delta _c\delta n_c/n_c`$. Outside the horizon, $`S`$ does not change with time.

The terms “adiabatic” and “isocurvature” refer to “initial conditions” specified at a time when the universe is radiation-dominated and all relevant scales are well beyond the horizon. During this era $`S`$ and $``$ are time-independent. An adiabatic perturbation is then defined as one for which $`S=0`$, and an isocurvature perturbation as one for which $`=0`$. A general perturbation is a superposition of an adiabatic and an isocurvature perturbation (in fact, there can be several kinds of isocurvature perturbations as there are several kinds of matter; we consider here a CDM isocurvature perturbation).

If the perturbations are Gaussian, their statistical properties are fully described by their power spectra $`P_{}(k)|_k|^2`$ and $`P_S(k)|S_k|^2`$. The spectral index is defined by $`n(k)d\mathrm{ln}P(k)/dk+4`$. If the indices are scale independent, we can write

$`P_{}(k)=Ak^{n_{\mathrm{ad}}4},`$ (2)
$`P_S(k)=Bk^{n_{\mathrm{iso}}4},`$ (3)
where $`n_{\mathrm{ad}}`$ and $`n_{\mathrm{iso}}`$ are the spectral indices with $`n_{\mathrm{ad}}=1`$ and $`n_{\mathrm{iso}}=1`$ corresponding to a completely scale free spectrum. (For $`n_{\mathrm{ad}}`$, this definition is standard, for $`n_{\mathrm{iso}}`$, other definitions in use set the index for a scale-free spectrum to $`n_{\mathrm{iso}}=0`$ or $`n_{\mathrm{iso}}=3`$.)

While the perturbations are outside the horizon, $`S`$ remains time independent, but an isocurvature perturbation develops a nonzero $`=\frac{1}{3}S`$ when the universe becomes matter dominated. Thus, during last scattering a curvature perturbation is already present. During matter domination, the density perturbation is related to the curvature perturbation by

$$\delta _k=\frac{2}{5}\left(\frac{k}{aH}\right)^2_k,$$
(4)
while outside the horizon ($`kaH`$).

For low multipoles, $`l200`$, the corresponding scales are well outside the horizon during last scattering and the CMB anisotropy has a simple relation to the initial perturbations. For adiabatic perturbations, the anisotropy comes from the Sachs-Wolfe-effect,

$$\frac{\delta T}{T}=\frac{1}{5},$$
(5)
whereas for isocurvature perturbations there is also a direct contribution from the entropy perturbation, which is 5 times larger,

$$\frac{\delta T}{T}=\frac{1}{5}\frac{1}{3}S=\frac{6}{15}S=\frac{6}{5}.$$
(6)
(This factor 6 in $`\delta T/T`$ becomes a factor 36 in the $`C_l`$.)

For higher multipoles the perturbations have evolved by the time of last scattering and the perturbations contribute to the anisotropy by other mechanisms also. The full angular power spectra from given initial perturbations need to be calculated by computer codes such as CMBFAST.

## III Cosmological Models

We allow for different spectral indices (“tilts”) $`n_{\mathrm{ad}}`$ and $`n_{\mathrm{iso}}`$ for adiabatic and isocurvature perturbations, and define

$$\alpha \frac{C_2^{\mathrm{iso}}}{C_2^{\mathrm{iso}}+C_2^{\mathrm{ad}}}$$
(7)
to indicate the relative contribution of the isocurvature fluctuations. Here the $`C_2^{\mathrm{iso}}`$ and $`C_2^{\mathrm{ad}}`$ are the expectation values of the isocurvature and adiabatic contributions to the quadrupole ($`C_2=C_2^{\mathrm{iso}}+C_2^{\mathrm{ad}}`$).

Note that different authors use different definitions for this parameter. Since our definition fixes the isocurvature contribution at the low end of the multipole range, the relative isocurvature contribution at high multipoles can be quite different depending on $`n_{\mathrm{iso}}`$ and $`n_{\mathrm{ad}}`$. Since isocurvature perturbations produce much more power at low multipoles, the isocurvature contribution to high multipoles is much smaller than $`\alpha `$ for $`n_{\mathrm{iso}}n_{\mathrm{ad}}`$.

Motivated by inflation, we restrict ourselves to models with $`\mathrm{\Omega }=1`$. We also take $`r=0`$, assume $`\tau =0`$ for the optical depth due to reionization, and consider a seven-parameter space of cosmological models. These parameters are $`\mathrm{\Omega }_\mathrm{\Lambda }1\mathrm{\Omega }_m`$, $`\omega _b\mathrm{\Omega }_bh^2`$, $`\omega _c\mathrm{\Omega }_ch^2`$, the spectral indices $`n_{\mathrm{ad}}`$ and $`n_{\mathrm{iso}}`$, and the amplitudes $`A`$ and $`B`$ of the power spectra \[Eq. (3)\]. The last two translate into the isocurvature contribution parameter $`\alpha `$ and an overall normalization. The Hubble constant is given by $`h=\sqrt{(\omega _c+\omega _b)/(1\mathrm{\Omega }_\mathrm{\Lambda })}`$.

In the present paper we assume that the adiabatic and isocurvature perturbations are uncorrelated and that their spectral indices are constant over the cosmological range of scales. This is not necessarily the situation in all particle physics models, but the nature of correlation between isocurvature and adiabatic amplitudes is very much model dependent. In general, such correlations would result in bounds more stringent than presented here.

We use the COBE, Boomerang, and Maxima-1 data and find the best-fit models which satisfy $`h=0.45`$–0.90 and $`\sigma _8\mathrm{\Omega }_m^{0.56}`$ = 0.43–0.70 (top-hat priors). We do not use any prior for the baryon density $`\omega _b`$.

As noted in , especially the Boomerang data favors $`\omega _b0.03`$, which is above the range allowed by standard big-bang nucleosynthesis, $`\omega _b`$ = 0.006–0.023. If this situation persists in the future, when more accurate CMB data becomes available, one may need to seriously consider non-standard nucleosynthesis scenarios.

We first used a coarse grid in the range (following) $`\mathrm{\Omega }_\mathrm{\Lambda }`$ = 0–0.8, $`\omega _b`$ = 0.003–0.130, $`\omega _c`$ = 0.02–0.80. We then refined it to focus on the relevant region; for adiabatic and mixed models: $`\mathrm{\Omega }_\mathrm{\Lambda }`$ = 0.48–0.78, $`\omega _b`$ = 0.018–0.042, $`\omega _c`$ = 0.13–0.21, with about a dozen values for each parameter; and for pure isocurvature models: $`\mathrm{\Omega }_\mathrm{\Lambda }`$ = 0.25–0.75, $`\omega _b`$ = 0.004–0.028, $`\omega _c`$ = 0.10–0.40. For the spectral indices we used a denser grid with step 0.02; in the range $`n_{\mathrm{ad}}`$ = 0.5–2.0, $`n_{\mathrm{iso}}`$ = 0.5–2.5 for adiabatic and mixed models; and $`n_{\mathrm{iso}}`$ = 1.5–2.7 for pure isocurvature models.

## IV Results

The best-fit adiabatic model (model 1, see Table I) is close to scale free, $`n_{\mathrm{ad}}=0.96`$, and the fit has $`\chi ^2=23.0`$ for 30 data points and 5 parameters. (We did not attempt to refine the parameter values between our grid points, as this is not the focus of this paper. See for cosmological parameter estimation in adiabatic models for the Boomerang and Maxima-1 data.)

The best-fit isocurvature model (model 2) has $`n_{\mathrm{iso}}=2.10`$ and $`\chi ^2=126.4`$, for 5 parameters. Clearly a pure isocurvature model is ruled out.

The best fit mixed (adiabatic+isocurvature) model (model 3) has $`n_{\mathrm{ad}}`$ = 0.98, $`n_{\mathrm{iso}}`$ = 2.26, and $`\alpha `$ = 0.008. The fit has $`\chi ^2=22.5`$, for 7 parameters. The small improvement in the fit from allowing an isocurvature contribution indicates that the data does not suggest the presence of such a contribution (but neither rules it out). The big blue tilt $`n_{\mathrm{iso}}=2.26`$ in the best-fit isocurvature contribution is not significant, since the fit was almost as good for any $`n_{\mathrm{iso}}`$ between 0.5 and 2.5. The isocurvature contribution is so small that it does not affect the angular power spectrum very much, and therefore also its tilt is not constrained. The small shifts from changing $`n_{\mathrm{iso}}`$ can be approximated by changing the other parameters, within the accuracy of the present data.

For large $`n_{\mathrm{iso}}`$ the smallness of the parameter $`\alpha `$ is somewhat misleading, as we defined it as the contribution to the lowest multipole. To represent the isocurvature contribution in the region of the first peak, let us define

$$\alpha _{200}\frac{C_{200}^{\mathrm{iso}}}{C_{200}^{\mathrm{iso}}+C_{200}^{\mathrm{ad}}},$$
(8)
where for the purpose of this paper $`C_{200}`$ is the average $`l(l+1)C_l`$ in the 4th Boomerang bin $`l`$ = 176–225. For the above best-fit model, $`\alpha _{200}=0.04`$. Similarly, the isocurvature contribution to $`\sigma _8^2`$ is 8% in this model. The models with smaller $`n_{\mathrm{iso}}`$ that are almost as good fits naturally have larger $`\alpha `$ (see Fig. 1) but smaller $`\alpha _{200}`$. Accordingly, the best-fit scale-free ($`n_{\mathrm{iso}}`$,$`n_{\mathrm{ad}}`$) = (1,1) mixed model (model 4) is almost as good a fit, $`\chi ^2=23.0`$ for only 5 parameters, but it has a larger $`\alpha `$ = 0.09 and a smaller $`\alpha _{200}`$ = 0.002.

For the 2-$`\sigma `$ ($`\mathrm{\Delta }\chi ^2=4`$) upper limit to $`\alpha `$ we find $`\alpha 0.63`$ and $`\alpha _{200}0.13`$. Within scale-free models the upper limits are $`\alpha 0.34`$ and $`\alpha _{200}0.010`$.

In mixed models with a large $`n_{\mathrm{iso}}`$, the isocurvature contribution to $`\sigma _8`$ can be significant. The 2-$`\sigma `$ upper limit to the isocurvature contribution to $`\sigma _8^2`$ is 35%, but within scale-free models only 0.3%.

In Fig. 1 we show the allowed region in the ($`n_{\mathrm{iso}}`$,$`n_{\mathrm{ad}}`$)-plane, and the best-fit values of $`\alpha `$ for each ($`n_{\mathrm{iso}}`$,$`n_{\mathrm{ad}}`$). In Fig. 2 we show the allowed region ($`\mathrm{\Delta }\chi ^23.5`$) in the ($`n_{\mathrm{iso}}`$,$`n_{\mathrm{ad}}`$,$`\alpha `$)-space as contours of maximum $`\alpha `$ on the ($`n_{\mathrm{iso}}`$,$`n_{\mathrm{ad}}`$) plane. This indicates the magnitude of the isocurvature contribution which remains allowed by the present data for different spectral indices.

We see that a blue tilt, $`n_{\mathrm{ad}}>1`$, in the adiabatic spectrum favors a larger isocurvature contribution, since the low-multipole power from the isocurvature perturbations compensates for this tilt. In fact, the presence of an isocurvature contribution makes a larger $`n_{\mathrm{ad}}`$ acceptable than in the pure adiabatic case. A large blue tilt, $`n_{\mathrm{iso}}1`$, in the isocurvature spectrum makes it less helpful for this purpose, and the upper limit to $`\alpha `$ goes down again.

Since $`\alpha `$ measures the contribution to the lowest multipoles, the contribution to high multipoles for a given $`\alpha `$ is larger for large $`n_{\mathrm{iso}}`$. Thus, although the limits to $`\alpha `$ are tighter for large $`n_{\mathrm{iso}}`$, this does not mean that the contribution to high multipoles would have to be less. In fact, the largest allowed $`\alpha _{200}`$, above 0.1, are for spectral indices $`n_{\mathrm{iso}}>1.8`$ (see Fig. 3).

In Fig. 4 we show the angular power spectra for some of the models discussed above, together with the COBE, Boomerang, and Maxima-1 data. Figures 5 and 6 show separately the adiabatic and isocurvature contributions in models 5 and 6, which have large isocurvature contributions (large $`\alpha `$ and large $`\alpha _{200}`$, respectively) but are allowed by the present CMB data.

We see how the isocurvature contribution slightly modifies the shape of the angular power spectrum. In model 6 it fills the minima between the peaks in the adiabatic contribution. An isocurvature contribution with a large blue tilt could be the reason for the lack of prominence of the 2nd peak in the Boomerang data.

The Boomerang and Maxima teams estimate a calibration uncertainty of 10% and 4%, respectively, in their measurements. This uncertainty is not included in the above results. To study its effect, we multiplied the Boomerang measurements by $`1.1^2`$ and divided the Maxima-1 measurements by $`1.04^2`$. This brings the measurements closer to each other, since the Maxima-1 points tend to lie above the Boomerang points, with the exception of the $`ł_{\mathrm{eff}}=147`$ point. The best-fit mixed model becomes $`n_{\mathrm{ad}}=1.00`$, $`n_{\mathrm{iso}}=1.78`$, and $`\alpha =0.05`$ (model 9), with an improved fit, $`\chi ^2=20.9`$. The upper limit to the isocurvature contribution is tightened to $`\alpha 0.56`$ ($`\alpha _{200}0.13`$ stays the same), but the results do not change qualitatively.

The data is consistent with a purely adiabatic perturbation with $`n_{\mathrm{ad}}1`$, with the best fit $`n_{\mathrm{ad}}`$ = 0.96 (allowing for reionization would shift this even closer to 1). In large-field inflation models this translates into $`r0.3`$, so that the tensor contribution could be of the same order as the isocurvature contribution. Because of the degeneracy between isocurvature and tensor perturbation such a situation is not altogether surprising.

The COBE data together with $`\sigma _8`$ have an important role in constraining the spectral indices of the models. The Boomerang and Maxima-1 data are important in constraining the amplitude of the isocurvature contribution, as they outline the pattern of acoustic peaks and thus are able to distinguish between the different patterns of peaks in the isocurvature vs adiabatic models.

## V Conclusions

The main conclusion of the present work is that the Boomerang and Maxima-1 data definitely rule out all purely isocurvature models, including those with a large tilt . This is because Boomerang and Maxima-1 define the shape of the angular power spectrum in the region of the first acoustic peaks much more accurately than the previous data.

However, a significant isocurvature contribution remains allowed. Indeed, as much as about half of the power at low multipoles could come from the isocurvature contribution, although the adiabatic contribution must dominate the first acoustic peak. This certainly leaves much room for various particle physics models. We should however point out that in many models, such as two-field inflation models or Affleck-Dine models, isocurvature and adiabatic fluctuations are not completely independent, and much more stringent constraints on $`\alpha `$ (and possibly $`n_{\mathrm{iso}}`$) could be obtained in specific cases.

Note also that we cannot obtain any clear conclusion about the spectral index $`n_{\mathrm{iso}}`$; except for large $`n_{\mathrm{iso}}`$, the correlation between $`\alpha `$ and $`n_{\mathrm{iso}}`$ is small. MAP and Planck will provide new accurate data on high multipoles and therefore could further constrain $`n_{\mathrm{iso}}`$, but at moderate tilts we do not expect MAP to improve much on the Boomerang and Maxima-1 limit on $`\alpha `$. The situation with Planck is different, provided that high quality polarization data is achieved as expected. Polarization will resolve the degeneracy between isocurvature and tensor perturbations and push the limit on $`\alpha `$ lower by an order of magnitude . Nevertheless, as discussed in the present paper, already the balloon experiments on the temperature fluctuations of the microwave sky can yield interesting constraints on particle physics models.

## Acknowledgments

This work has been supported by the Academy of Finland under the contracts 101-35224 and 101-47213. We thank the Center for Scientific Computing (Finland) for computational resources. We acknowledge the use of the CMBFAST Boltzmann code developed by Uroš Seljak and Matias Zaldarriaga.


format: text

text:
# On the Spectrum of Ultrahigh Energy Cosmic Rays and the Gamma Ray Burst Origin Hypothesis

## 1 Introduction

Ultra-high energy cosmic rays (UHECRs) are among the most energetic form of radiation ever detected, with energies apparently reaching beyond 1020eV. Cosmic rays with ultrahigh energies were first detected over 35 years ago by the Volcano Ranch experiment . With typical event rates of one per km2 per century, subsequent experiments including the Fly’s Eye  and AGASA  experiments have observed only a double digit number of events. However, in spite of the paucity of data, a surprising characteristic in the UHECR spectrum has been observed, viz., its apparent smooth continuation beyond $`10^{20}`$ eV. UHECRs appear to be distributed isotropically, indicating an extragalactic origin. Protons traveling from intergalactic distances should experience energy losses owing to photopion production interactions with the 2.7K cosmic background radiation. For a uniform distribution of protons in the universe, even neglecting redshift effects and cosmological source evolution, these interactions would be expected to produce a cutoff in the observed spectrum of UHECRs at $`6\times 10^{19}`$ eV, commonly referred to the Greisen-Zatsepin-Kuzmin (GZK) cutoff. This is because of the strong energy dependence in the energy loss rate with energy from photomeson production interactions. As a result, above $`10^{20}`$ eV, only protons from within the relatively nearby distance of $`100`$ Mpc will survive at these energies to reach the Earth . The present data indicate, however, a flattening of the spectrum above $`10^{19}`$ eV and no obvious cutoff out to an energy of $`3\times 10^{20}`$ eV.

Two classes of theories have emerged to explain this phenomenon. In one class of theories, “top down”, a small branching ratio into high energy nucleons would result from the decays of unstable or meta-stable supermassive particles originating at the grand unified scale or in topological defects produced at that scale. (Most of the energy from these decays goes into the production of pions which, in turn, decay to produce much larger numbers of neutrinos, photons and electrons.)

In the other class of theories, “bottom up”, protons achieve their high energies through conventional electromagnetic acceleration processes. Only a few astrophysical sources are presently considered to have the possibility to accelerate particles to ultrahigh energies of the order of $`10^{20}`$ eV. Suggested possible acceleration sites include AGN jets and hot spots at the ends of radio galaxy jets. For reviews, see Ref. .

Among other origin hypotheses, it has been suggested that $`\gamma `$-ray bursts (GRBs) can produce the observed flux of UHECRs. In such scenarios, it has been speculated that UHECRs would be emitted by GRBs fireballs with roughly the same amount of total emitted energy as the observed $`\gamma `$-rays in the keV to MeV range. Based on this assumption, and also assuming that the GRB rate is independent of redshift, it has been argued that the flux and spectrum of UHECRs with energies above $`10^{19}`$ eV (a putative extragalactic component which could account for the flattening and possible composition change above $`10^{19}`$ eV) is consistent with their production by GRBs .

In recent years, X-ray, optical, and radio afterglows of about a dozen GRBs have been detected leading to the subsequent identification of the host galaxies of these objects and consequently, their redshifts. The regions in the host galaxies where the bursts are located have also been identified as regions of very active ongoing star formation. One of us  has therefore argued that a more appropriate redshift distribution to take for GRBs is that of the average star formation rate. One consequence of most GRBs being at moderate to high redshifts, is that the GZK effect in the UHECR spectrum produced by GRBs would materialize at lower energies and be more pronounced, in stark contradiction to the data. In this Letter, we present a detailed calculation of the expected UHECR spectrum which would result from GRBs with a more realistic redshift distribution, one following the star formation rate.

## 2 Gamma Ray Burst Redshift Evolution

To date, some 15 GRBs afterglows of long duration bursts have been detected with a subsequent identification of their host galaxies. At least 14 of the 15 are at moderate to high redshifts with the highest one (GRB000131) lying at a redshift of 4.50,.111The origin of GRB980425 is somewhat controversial; a possible X-ray source and an unusual nearby Type Ic supernova have both been put forward as candidates. Taking the supernova identification gives an energy release which is orders of magnitude smaller than that for a typical cosmological GRB. Redshift data on GRB host galaxies are presently available only for long duration bursts; the short duration bursts have not as yet been identified with any astronomical sources, however, their source-count distribution is also consistent with a cosmological origin .

The host galaxies of GRBs appear to be sites of active star formation. The colors and morphological types of these galaxies are also consistent with active star formation as is the detection of Ly$`\alpha `$ and \[OII\] in several of the host galaxies. Further evidence suggests that bursts themselves are directly associated with star forming regions within their host galaxies; their positions correspond to regions having significant hydrogen column densities with evidence of dust extinction. A good argument in favor of strong redshift evolution for GRBs is made in Ref. based on the nature of the host galaxies. Other recent analyses have also favored a GRB redshift distribution which follows the strong redshift evolution of the star formation rate ,.

One argument for strong redshift evolution for GRBs which is independent of the SFR assumption has been put forward very recently. It is based on correlating the time variability of GRB events with the absolute luminosity of the events. The relationship is based on seven events and has been extended to to a sample of 224 long, bright GRBs detected by BATSE. The GRB rate was found to scale as $`(1+z)^{3.6\pm 0.3}`$ which is roughly consistent with the SFR for $`z<1.5`$. For higher $`z`$, the SFR can not be accurately determined because it is strongly affected by dust. As we shall show in the next section, this uncertainty at such high redshifts does not strongly affect our result because contributions to the observed UHECR flux from such redshifts are negligible.

## 3 Calculations

Energy losses of high energy protons from their creation to the present epoch $`(z=0)`$ are a result of cosmological redshifting and by pair production ($`p+\gamma p+e^++e^{}`$) and pion production ($`p+\gamma \pi +n`$) through interactions with cosmic background radiation (CBR) photons. We shall assume for this calculation a flat ($`\mathrm{\Omega }=1`$) Einstein-de Sitter universe with a Hubble constant of H0 = 70 km sec-1 Mpc-1 ($`h=0.7`$) and taking $`\mathrm{\Lambda }=0`$.222Taking a $`\mathrm{\Lambda }`$ = 0.7 model does not significantly affect the results. The relationship between time and redshift is then given by

$`t={\displaystyle \frac{2}{3}}H_0^1(1+z)^{3/2},`$ (1)
and the energy loss owing to redshifting is

$`{\displaystyle \frac{1}{E}}{\displaystyle \frac{dE}{dt}}=H_0(1+z)^{3/2}.`$ (2)
Photopion production is the dominant interaction between UHECRs and the CMB for energies $`>`$ 1020 eV. Below energies of 3 $`\times `$ 1019 eV, protons will lose energy mainly through pair production. The combined energy loss rate for both pion and pair production for protons in collisions with photons of the 2.7K background at the present epoch ($`z=0`$) is defined as

$`{\displaystyle \frac{1}{E}}{\displaystyle \frac{dE}{dt}}\tau (E)^1.`$ (3)
The photomeson loss time is given by Stecker and that for pair production is given by Blumenthal. In order to compute the energy losses due to pion and pair production occurring at a redshift $`z`$, one must account for the increased photon densities and energies at larger $`z`$. The energy of photons increases as $`(1+z)`$ and the photon density increases as $`(1+z)^3`$. The energy loss for photons for any epoch $`z`$ can then be expressed as

$`{\displaystyle \frac{1}{E}}{\displaystyle \frac{dE}{dt}}=(1+z)^3\tau [(1+z)E]^1.`$ (4)
We can then calculate the energy E0(z) at which a proton is created at a redshift $`z`$ whose observed energy today is E.

In order to calculate the (number) flux of protons observed on earth of energy $`E`$, we follow the formalism derived in Berenzinsky and Grigor’eva . Because of their high energies, we will assume that these protons will propagate along roughly straight lines from their source to the observer, unaffected by small intergalactic magnetic fields. (Dropping this assumption increases their propagation time and energy losses so that it would only strengthen our argument.)

The observed flux from a volume element of the sphere $`dV=R^3(z)r^2drd\mathrm{\Omega }`$ is

$`{\displaystyle \frac{dj}{dE}}dE={\displaystyle \frac{F(E_0,z)dE_0n(z)dV}{(1+z)4\pi R_0^2r^2}}.`$ (5)
where F(E0,z) is the flux of particles emitted at $`z`$ of energy E0 and $`n(z)`$ is the density of GRBs as a function of redshift. We will take this function to be of the form

$`n(z)=(1+z)^{(3+q)}n(0).`$ (6)
Furthermore with $`R_0=R(z)(1+z)`$ and $`R(z)dr=cdt`$, changing time to redshift dependence utilizing (1), and by integrating (5) to the edge of the sphere (i.e. out to a maximum redshift $`z_{max}`$) one obtains

$`j(E)={\displaystyle \frac{3}{8\pi }}R_0n(0)F(E){\displaystyle _0^{z_{max}}}(1+z)^{(q5/2)}\left({\displaystyle \frac{E_0}{E}}\right)^\alpha {\displaystyle \frac{dE_0}{dE}}𝑑z.`$ (7)
assuming a power-law spectrum of index $`\alpha `$ for F(E) to be generated by the GRBs, where F(E) is the total number of particles per second with energy E observed.

Above $`10^{19}`$ eV, the Fly’s Eye and AGASA data are well fit by a spectrum with a power law $``$ E-2.75. The “GZK cutoff energy” for UHECRs versus redshift of origin is plotted in Figure 1. We have defined this as the energy at which the spectrum drops to $`1/e`$ of its original value owing to photomeson production. It can be seen from Figure 1 that the protons of energy above ($`>10^{20}`$ eV) can have originated no further than a redshift of $`z=.034`$ or a distance of 130 Mpc.

For the redshift distribution of GRBs, we have have considered three scenarios and corresponding values for the evolution index, $`q`$. In the first two cases, we assume a value for $`\alpha `$ of 2.75 which fits the data above 1019eV. For the first case, we take an isotropic comoving source density distribution corresponding to $`q=0`$ out to a maximum redshift of 2.5. In a second case, we take a redshift dependence proportional to the star formation rate corresponding to $`q=3`$ out to a redshift of 2.5.

For our third case, we take the redshift distribution for GRBs derived from the time variability-redshift relationship for GRBs. This corresponds to taking $`q=3.6`$ out to a redshift of 3.6. (It should be noted that there is no observable contribution to the cosmic ray flux from redshifts beyond 2.5 in any case.) For this case, we assume a flatter source spectrum with a smaller value for $`\alpha `$ of 2.35. This gives a better fit to the data between 1018 and 1019 eV, but requires a larger energy input. (If we take $`\alpha `$ = 2.75 for this case, we get a very similar result to the second case, but slightly lower fluxes.)

We have normalized our resulting spectra to the Fly’s Eye data at $`E=10^{18}`$ eV and calculated the corresponding UHECR spectra observed at Earth under the assumptions for the three scenarios described above. The results of these scenarios are plotted against the Fly’s Eye and AGASA data in figure 2. The failure of all three scenarios to fit the observed observed spectrum, particularly for energies above 1020 eV, is evident. In the scenarios where the source rate follows the star formation rate as a function of redshift, the spectrum is significantly steepened from the source spectrum above 1018 eV, a situation which does not occur in the case of a uniform distribution independent of redshift.

For the $`\alpha `$ = 2.75 cases, the required energy input in cosmic rays above 1019 eV is $`2\times 10^{45}`$ erg Mpc-3 yr-1. This is two orders of magnitude greater than the energy release rate from GRBs in the 0.01-1 MeV range . For the $`\alpha =2.35`$ case, the discrepency increases to three orders of magnitude, given the required normalization and the flatter source spectrum.

One may speculate that there are numerous nearby low-flux GRBs below the BATSE threshold. However, BATSE has only detected one burst which may have had an energy release significantly lower than $`10^{50}`$ erg, even this event being in doubt -. Thus, there is no reason to believe that such a scenario could overcome the energy discrepency shown above unless much more GRB energy is released in UHECRs than in low-energy $`\gamma `$-rays.

Our calculated spectrum above $`10^{20}`$ eV may be underestimated owing to the fact that fluctuations in the interaction rate from very nearby sources can allow more protons to reach the Earth than our continuous-energy-loss approximation indicates . However, even with the inclusion of this effect, the predicted energy fluxes from GRBs would be much lower than the observations indicate. With regard to fluctuations, it has been argued that a few nearby sources could produce the UHECRs  but statistical studies rule this out .

## 4 Conclusions

In this paper, we have presented the predicted spectra of UHECRs assuming that the UHECR sources are GRBs and that their source density distribution follows the star formation rate. Our motivation is the recent redshift studies of GRBs which suggest a strong redshift evolution for these objects. The results indicate that GRBs can not be responsible for the UHECRs observed at Earth because their redshift evolution leads to a predicted spectrum which is well below the data, particularly for energies above $`10^{20}`$ eV. In fact, the “effective GZK cutoff” for UHECRs of cosmological GRB origin is found to be $`3\times 10^{19}`$eV (see Figure 2).

It should be noted that, independent of spectral considerations, if the GRB redshift distribution follows the star formation rate, GRBs can be ruled out as candidates for the UHECRs based on energetics considerations. The present star formation rate coupled with the mean energy release per burst leads to an energy release rate per volume which is more than an order of magnitude below the rate need to account for the UHECRs with energies above $`10^{20}`$ eV, and, as shown in the previous section, two to three orders of magnitude below the rate needed to account for UHECRs above $`10^{19}`$ eV, contradicting the hypothesis advocated in Ref.  and stengthening the conclusion of Ref. .

The recently launched High Energy Transient Explorer (HETE-2) and the Swift satellite mission, projected to be launched in 2003, should eventually provide redshift information for thousands of bursts, including the short duration bursts. These missions therefore have the potential to conclusively exclude gamma ray bursts as the source of the UHECRs, should the density distribution of the bursts conclusively indicate redshift evolution.

Future and planned UHECR detectors such as the Pierre Auger Observatory, and the space-based orbiting wide-angle lens collector (OWL) satellite system should provide the coverage area necessary to determine the true particle spectrum above 1020eV. Confirmation of the lack of a pronounced GZK cutoff and an event rate similar to that indicated by present data would also spell the demise for a GRB origin of the UHECRs.


text:
# On the genera of 𝑋₀(𝑁)

## 1. Introduction

For each positive integer $`N`$, the modular curve $`X_0(N)`$ parametrizes elliptic curves together with a cyclic subgroup of order $`N`$ (for more details, see  and ). The genus $`g_0(N)`$ of $`X_0(N)`$ tends to infinity as the level $`N`$ increases. We will examine the sequence $`\{g_0(N)\}`$ more closely. Among other things, we will show that the average size of $`g_0(N)`$ is $`1.25N/\pi ^2`$; however, the values of $`g_0(N)`$ form a density-zero subset of the integers. These two properties imply that there is much collapsing under the map $`Ng_0(N)`$: for instance, there are integers whose preimage under this map is arbitrarily large.

The results of this note are:

1. Upper and lower bounds, including the asymptotic results
$`0<lim\; infg_0(N)/N<lim\; supg_0(N)/(N\mathrm{log}\mathrm{log}N)<\mathrm{}`$.
2. Average behavior: $`lim_B\mathrm{}(1/B)_{NB}g_0(N)/N=1.25/\pi ^2`$.
3. Natural density: $`\{g_0(N)\}`$ is a density zero subset of the integers.
4. Non-uniformity of $`g_0(N)`$ modulo a fixed prime $`p`$: for instance, $`g_0(N)`$ is odd with probability 1, and (for a fixed odd prime $`p`$) the probability that $`g_0(N)1(modp)`$ is much less than $`1/p`$.

This note arose from a question about curves over finite fields. Typically one estimates the number of points on such a curve by means of the Weil (upper) bound. However, in various situations this bound can be improved, which leads to the converse question of producing curves with many points (of a given genus, over a given field) in order to see how far the Weil bound can be improved. It is known that, if the prime $`p`$ does not divide $`N`$, then $`X_0(N)`$ has many rational points over the finite field $`𝔽_{p^2}`$. Covers of these curves have been used to show that there are curves with many points (over $`𝔽_{p^2}`$) in every genus . It is natural to ask whether covers are needed to prove the latter result, or whether the genera of $`X_0(N)`$ already achieve all sufficiently large integer values. The present note grew from our proof that infinitely many positive integers do not occur in the sequence $`\{g_0(N)\}`$.

The authors suspect that some of the results of this note have been discovered previously, but apparently they have not been published.111Indeed, we recently saw an unpublished manuscript of S. Wong which contains some of the results of this note.

## 2. The genus of $`X_0(N)`$

For a positive integer $`N`$, let $`g_0(N)`$ denote the genus of $`X_0(N)`$. By \[15, Prop. 1.43\], we have

$$g_0(N)=1+\frac{\mu }{12}\frac{\nu _2}{4}\frac{\nu _3}{3}\frac{\nu _{\mathrm{}}}{2}$$

where $`\nu _2`$ and $`\nu _3`$ are the numbers of solutions in $`/N`$ of the equations $`x^2+1=0`$ and $`x^2+x+1=0`$, respectively; $`\nu _{\mathrm{}}=_d\phi (\text{gcd}(d,N/d))`$; and $`\mu =N_a1/a`$. Here $`\phi `$ denotes Euler’s totient function, $`d`$ ranges over the positive divisors of $`N`$, and $`a`$ ranges over the squarefree positive divisors of $`N`$. Likewise, let $`p`$ range over the prime divisors of $`N`$, and write $`N=_pp^{r_p}`$. Then $`\mu =_p(p+1)p^{r_p1}`$. Also $`\nu _2`$ is 0 if $`4N`$ or if some $`p3(mod4)`$, and otherwise $`\nu _2=2^s`$ where $`s`$ is the number of $`p1(mod4)`$; similarly, $`\nu _3`$ is 0 if $`9N`$ or if some $`p2(mod3)`$, and otherwise $`\nu _3=2^t`$ where $`t`$ is the number of $`p1(mod3)`$. Finally, $`\nu _{\mathrm{}}=_p\theta (p,r_p)`$, where we define $`\theta (p,2R+1)=2p^R`$ and $`\theta (p,2R)=(p+1)p^{R1}`$.

## 3. Bounds on $`g_0(N)`$

We now give general upper and lower bounds on $`g_0(N)`$. For the lower bound we show that

$$g_0(N)(N5\sqrt{N}8)/12,$$

with equality if and only if $`N=p^2`$ where $`p`$ is a prime congruent to 1 (mod 12). For the upper bound we show the asymptotic result

$$\underset{N\mathrm{}}{lim\; sup}\frac{g_0(N)}{N\mathrm{log}\mathrm{log}N}=\frac{e^\gamma }{2\pi ^2}$$

and the explicit bound

$$g_0(N)<N\frac{e^\gamma }{2\pi ^2}(\mathrm{log}\mathrm{log}N+2/\mathrm{log}\mathrm{log}N)\text{for }N>2,$$

where $`\gamma =0.5772\mathrm{}`$ is Euler’s constant.

We start by proving the lower bound $`g_0(N)(N5\sqrt{N}8)/12`$. If $`N`$ is a prime power (or $`N=1`$), it is easy to prove this bound and to show that equality occurs precisely when $`N`$ is the square of a prime congruent to 1 (mod 12). It is also easy to check the lower bound for $`N`$ less than 2000, so we may assume in the sequel that $`N2000`$.

Suppose that $`N`$ has exactly two distinct prime factors $`p`$ and $`q`$. Then $`\mu N>N/p+N/q>2\sqrt{N}`$. Since $`\nu _{\mathrm{}}\mu /\sqrt{N}`$ and $`\nu _2,\nu _34`$, we have $`g_0(N)1+\mu /124/44/3\mu /(2\sqrt{N}).`$ Thus,

$`12g_0(N)`$ $`\mu (16/\sqrt{N})16`$
$`>(N+2\sqrt{N})(16/\sqrt{N})16`$
$`>N5\sqrt{N}8.`$
Finally, consider the case where $`N`$ has at least 3 distinct prime factors. In this case we observe that $`\mu N>3N^{2/3}`$ and $`\nu _2,\nu _3\nu _{\mathrm{}}`$. It follows that $`g_0(N)1+\mu /12(1/4+1/3+1/2)\nu _{\mathrm{}}`$; hence,

$`12g_0(N)`$ $`12+\mu 13\mu /\sqrt{N}`$
$`>(N5\sqrt{N}8)+(\mu N)(113/\sqrt{N})8\sqrt{N}+20`$
$`>(N5\sqrt{N}8)+3N^{4/6}39N^{1/6}8N^{3/6}+20.`$
The largest real root of $`3x^48x^339x+20`$ is $`3.5495\mathrm{}<2000^{1/6}`$. The lower bound on $`g_0(N)`$ follows.

We now prove the limsup result. We have already noted that $`\nu _2`$, $`\nu _3`$, and $`\nu _{\mathrm{}}`$ are each bounded by $`\mu /\sqrt{N}`$; thus, it suffices to show that

$$\underset{N\mathrm{}}{lim\; sup}\frac{\mu }{N\mathrm{log}\mathrm{log}N}=\frac{6}{\pi ^2}e^\gamma .$$

For any $`x2`$ let $`N_x=_{px}p`$. Among all $`N`$ with $`N_xN<N_{x+1}`$ (for fixed $`x`$), the case $`N=N_x`$ maximizes $`\mu /N`$ and minimizes $`\mathrm{log}\mathrm{log}N`$. Thus, we need only consider values $`N=N_x`$. For these $`N`$, Mertens showed in 1874 that $`\mu /N`$ is asymptotic to $`(\mathrm{log}x)e^\gamma 6/\pi ^2`$ \[7, p. 429\]; the limsup result follows, since $`\mathrm{log}\mathrm{log}N_x`$ is asymptotic to $`\mathrm{log}x`$ \[7, p. 341\].

The upper bound on $`g_0(N)`$ is easily verified for $`N<210`$, so we now assume $`N210`$. Note that $`g_0(N)\mu /12`$, so our upper bound follows from the inequality $`\mu /N<(\mathrm{log}\mathrm{log}N+2/\mathrm{log}\mathrm{log}N)e^\gamma /(2\pi ^2)`$. As above, it suffices to prove this inequality when $`N=N_x`$ (since $`N210=N_7`$). But in this case the inequality follows easily from results in .

## 4. Data for small $`N`$

We now determine the first few positive integers $`n`$ which do not occur as $`g_0(N)`$ for any $`N`$. The lower bound of the previous section implies that, if $`n=g_0(N)`$, then $`N<12n+18\sqrt{n}+40`$. So, to determine whether $`n`$ occurs as $`g_0(N)`$, we just check all levels $`N`$ up to this bound. Doing this by computer, we find the first few missed values (i.e., positive integers $`n`$ not of the form $`g_0(N)`$): 150, 180, 210, 286, 304, 312, … Note that all of these are even; in fact, the first several thousand missed values are even. It is easy to describe all positive integers $`N`$ for which $`g_0(N)`$ is even. They are given in the following list, where $`p`$ denotes a prime and $`r`$ denotes a positive integer:

1. $`N=1,2,3,4,8`$ or $`16`$;
2. $`N=p^r`$ where $`p5(mod8)`$;
3. $`N=p^r`$ where $`p7(mod8)`$ and $`r`$ odd;
4. $`N=p^r`$ where $`p3(mod8)`$ and $`r`$ even;
5. $`N=2p^r`$ where $`p\pm 3(mod8)`$;
6. $`N=4p^r`$ where $`p3(mod4)`$ and $`r`$ odd.

It follows that, if $`N`$ is a randomly chosen positive integer, then $`g_0(N)`$ is odd with probability 1. This shows that a randomly chosen even positive integer has probability zero of occurring as a value of $`g_0(N)`$.

Looking further in the list of integers not of the form $`g_0(N)`$, we do eventually find some odd values, the first one occuring at the $`3885`$th position. There are four such up to $`10^5`$ (out of 9035 total missed values), namely 49267, 74135, 94091, 96463. In Section 6 we will show that the paucity of odd missed values is not a general phenomenon, but instead an accident caused by the fact that ‘small’ numbers do not have enough prime factors. In particular, we will see that there are infinitely many positive odd integers not of the form $`g_0(N)`$, and in fact the set of integers of the form $`g_0(N)`$ has density zero in the set of all nonnegative integers.

## 5. Average size of $`g_0(N)`$

We have shown that $`g_0(N)`$ is sometimes as small as $`N/12`$, and sometimes as large as $`cN\mathrm{log}\mathrm{log}N`$. We now determine the usual behavior. More precisely, we show that

$$\frac{1}{B}\underset{N=1}{\overset{B}{}}g_0(N)=\frac{5}{8\pi ^2}B+o(B).$$

By Abel’s lemma \[14, p. 65\], this result is equivalent to the following:

(1)
$$\underset{B\mathrm{}}{lim}\frac{1}{B}\underset{N=1}{\overset{B}{}}\frac{g_0(N)}{N}=\frac{5}{4\pi ^2}=0.12665\mathrm{}$$

In Section 3, we showed that $`g_0(N)/N=(1/12)_{aN}(1/a)+o(1)`$, and we can ignore the error term since it contributes nothing to the limit. Here we use the convention that $`p`$ is always prime, and $`a`$ is always a squarefree positive integer. The following computation proves (1):

$`\underset{B\mathrm{}}{lim}(1/B){\displaystyle \underset{N=1}{\overset{B}{}}}{\displaystyle \underset{aN}{}}1/a=\underset{B\mathrm{}}{lim}(1/B){\displaystyle \underset{aB}{}}(B/a)(1/a)`$
$`={\displaystyle \underset{a<\mathrm{}}{}}1/a^2=\zeta (2)/\zeta (4)=15/\pi ^2.`$
A. M. Odlyzko showed us the following alternative proof of (1): for any complex $`s`$ with $`\mathrm{Re}(s)>1`$, the series $`F(s):=_{N=1}^{\mathrm{}}(\mu /N)N^s`$ converges, and we have

$`F(s)`$ $`={\displaystyle \underset{N=1}{\overset{\mathrm{}}{}}}N^s{\displaystyle \underset{pN}{}}(1+1/p)={\displaystyle \underset{p}{}}\left(1+(1+1/p)(p^s+p^{2s}+\mathrm{})\right)`$
$`=\zeta (s){\displaystyle \underset{p}{}}(1+1/p^{s+1})=\zeta (s)\zeta (s+1)/\zeta (2s+2).`$
It follows that $`F(s)`$ has a simple pole at $`s=1`$ with residue $`\zeta (2)/\zeta (4)`$, and that $`F(s)`$ is regular on $`\mathrm{Re}(s)=1`$ except at $`s=1`$. Now (1) follows at once from Ikehara’s Tauberian theorem \[8, p. 311\].

## 6. Density of $`\{g_0(N)\}`$

In this section we consider the (natural) density of the set $`\{g_0(N)\}`$ as a subset of the non-negative integers. We have already seen that $`g_0(N)`$ is almost never even, so this density (if it exists) is at most $`1/2`$. In fact, we will show that the set $`\{g_0(N)\}`$ has density zero—in other words, we will show that $`lim_x\mathrm{}(\mathrm{\#}S(x))/x=0`$, where

$$S(x)=\{n:nx\text{ and }n=g_0(N)\text{ for some }N\}.$$

Let $`N`$ be a positive integer and suppose that $`N`$ is divisible by at least $`s>2`$ distinct odd primes. The formulas in Section 2 imply that $`2^{s1}`$ divides $`\nu _3`$ and $`2^s`$ divides each of $`\mu `$, $`\nu _2`$, and $`\nu _{\mathrm{}}`$. It follows that $`g_0(N)1(mod2^{s2})`$.

We now show that $`\mathrm{\#}S(x)x/2^d+o_d(x)`$ for each positive integer $`d`$; this implies that $`\mathrm{\#}S(x)=o(x)`$, as desired. Fix a positive integer $`d`$. Clearly the number of $`nS(x)`$ with $`n1(mod2^d)`$ is less than $`x/2^d+1`$. It remains to show that the number of $`nS(x)`$ with $`n1(mod2^d)`$ is $`o(x)`$. Each such $`n`$ has the form $`g_0(N)`$ where $`N`$ has at most $`d+1`$ distinct odd prime divisors and $`N<12x+18\sqrt{x}+40`$. And the number of such $`N`$ (hence also the number of such $`n`$) is well-known to be $`o(x)`$ \[7, p. 356\].

###### Remark.

We have shown that $`\mathrm{\#}S(x)=o(x)`$. It would be interesting to obtain a more precise estimate for $`\mathrm{\#}S(x)`$. We make some preliminary observations in this vein. First, the number of even integers in $`S(x)`$ is asymptotic to $`9x/\mathrm{log}x`$. (This follows from the classification of $`N`$ for which $`g_0(N)`$ is even, together with the fact that, for any fixed nonzero integers $`a,b`$, there are only $`o(x/\mathrm{log}x)`$ primes $`p<x`$ for which $`|ap+b|`$ is prime \[5, Cor. 2.4.1\].) We can prove the following upper and lower bounds on $`\mathrm{\#}S(x)`$ for $`x3`$:

$$\frac{x}{\mathrm{log}x}e^{a(\mathrm{log}\mathrm{log}\mathrm{log}x)^2}\mathrm{\#}S(x)\frac{x}{(\mathrm{log}x)^b(\mathrm{log}\mathrm{log}x)^c}.$$

(For positive $`f(x),g(x)`$, the notation $`f(x)g(x)`$ is equivalent to $`f(x)O(g(x))`$.) Here $`a`$ is any constant less than $`a_0`$, where the constants $`a_0=0.8168146\mathrm{}`$ and $`b=0.2587966\mathrm{}`$ and $`c=0.2064969\mathrm{}`$ are defined as follows. Let $`B`$ be the unique root of $`1/B+\mathrm{log}B=1+\mathrm{log}2`$ in the interval $`(0,1)`$, let $`A`$ be the unique root of $`_{n=1}^{\mathrm{}}A^n((n+1)\mathrm{log}(n+1)n\mathrm{log}n1)=1`$ in the interval $`(0,1)`$, and put $`a_0=1/(2\mathrm{log}A)`$ and $`b=B\mathrm{log}2`$ and $`c=(B\mathrm{log}2)/(22B)`$. The lower bound is proved by considering numbers $`N`$ which are products of a fixed number $`k`$ of distinct primes, the least of which is $`11`$; for such $`N`$ we have $`g_0(N)=12^{k1}+(1/12)_{p|N}(p+1)`$, and the lower bound follows from . The upper bound is proved by optimizing the choice of $`d`$ in our density-zero proof, after one has modified that proof by replacing the $`o(x)`$ bound from \[7, p. 356\] by the more precise bound from . Similar bounds have been proved for the number $`\mathrm{\#}V(x)`$ of distinct values of Euler’s $`\phi `$-function not exceeding $`x`$. In this setting, our upper bound is essentially an optimization of the paper . In , upper and lower bounds are proved for $`\mathrm{\#}V(x)`$, both of which have the same shape as the above lower bound (but the upper bound is for any $`b>b_0`$). The precise order of magnitude of $`\mathrm{\#}V(x)`$ is determined in the excellent paper , which (together with ) also contains several other interesting results whose analogues would be interesting to study in our situation. However, the multiplicativity of $`\phi (n)`$ plays a crucial role in all proofs giving better upper bounds for $`\mathrm{\#}V(x)`$ than our upper bound for $`\mathrm{\#}S(x)`$ above. Since $`g_0(N)`$ (and $`g_0(N)1`$) is not multiplicative, it seems that these methods do not apply to $`\mathrm{\#}S(x)`$. We do not know which of our upper and lower bounds is closer to the truth.

## 7. Distribution of $`g_0(N)`$ modulo primes

We now study the distribution of $`g_0(N)`$ modulo a fixed prime $`\mathrm{}`$ as $`N`$ varies. We will see that some residue classes occur more often than others. Let $`\mathrm{}`$ be a fixed prime. We may restrict to levels $`N`$ having a prime factor $`p`$ congruent to $`1`$ modulo $`12\mathrm{}`$, since the set of such $`N`$’s has density 1 in the set of positive integers (as shown by Dirichlet, cf. \[14, p. 73\]). Our assumption on $`N`$ forces $`\nu _2=\nu _3=0`$ and $`12\mathrm{}\mu `$, so

$$g_0(N)1\nu _{\mathrm{}}/2(mod\mathrm{}).$$

First consider the case $`\mathrm{}=2`$: for our restricted class of $`N`$’s, we have $`g_0(N)1(mod2)`$ unless $`N`$ is either $`p^{2r}`$ or $`4p^{2r}`$, so certainly $`g_0(N)`$ is odd with probability 1 (as was observed in Section 4).

We have seen that the sequence of residues mod 2 of $`g_0(N)`$ is biased. We now show a similar result modulo other primes $`\mathrm{}`$. We compute the density of positive integers $`N`$ for which $`g_0(N)1(mod\mathrm{})`$. As above, we may assume $`N`$ has a prime factor $`p`$ with $`12\mathrm{}(p+1)`$; then $`\mathrm{}(g_0(N)1)`$ is equivalent to $`\mathrm{}\nu _{\mathrm{}}`$, and this occurs precisely when either $`\mathrm{}^3N`$ or some prime congruent to $`1`$ mod $`\mathrm{}`$ divides $`N`$ with (positive) even multiplicity. By a standard argument (similar to \[10, Thm. 2.18\]), the probability that $`N`$ does not satisfy either of these conditions is

$$(11/\mathrm{}^3)\underset{\begin{array}{c}s1\text{ (mod }\mathrm{}\text{)}\\ s\text{ prime}\end{array}}{}(11/s^2+1/s^31/s^4+1/s^5\mathrm{}),$$

which equals $`(11/\mathrm{}^3)(11/(s^2+s))`$. The following table gives upper bounds for the probability $`P(\mathrm{})`$ that $`g_0(N)1(mod\mathrm{})`$:

| $`\mathrm{}`$ | 3 | 5 | 7 | 11 | 13 | 17 | 19 | 23 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $`P(\mathrm{})<`$ | 1/4 | 1/78 | 1/105 | 1/653 | 1/1542 | 1/1793 | 1/978 | 1/5821 |

In these first few cases, we see that the probability is much less than $`1/\mathrm{}`$. More generally, for every $`\mathrm{}`$ we have $`P(\mathrm{})<3/\mathrm{}^2`$ (where the constant ‘3’ is not optimal):

$`P(\mathrm{})`$ $`=1(11/\mathrm{}^3){\displaystyle \underset{\begin{array}{c}s1\text{ (mod }\mathrm{}\text{)}\\ s\text{ prime}\end{array}}{}}(11/(s^2+s))`$
$`<1(11/\mathrm{}^3){\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(11/((n\mathrm{})^2n\mathrm{}))`$
$`<1(11/\mathrm{}^3){\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(11.5/(n\mathrm{})^2)`$
$`<1(11/\mathrm{}^3)(1{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}(1.5/\mathrm{}^2)/n^2)`$
$`=1(11/\mathrm{}^3)(1\pi ^2/(4\mathrm{}^2))`$
$`<3\mathrm{}^2.`$
Also note that other cosets mod $`\mathrm{}`$ will have special behavior as well. If $`N`$ is squarefree, $`\nu _{\mathrm{}}`$ is a power of 2, so that $`g_0(N)`$ is congruent modulo $`\mathrm{}`$ to a number of the form $`12^k`$. Recall that $`N`$ is squarefree with probability $`6/\pi ^2=0.6079\mathrm{}`$ \[7, p. 269\]. Thus, if 2 is not a primitive root mod $`\mathrm{}`$, then the residue classes (mod $`\mathrm{}`$) of the integers $`12^k`$ will occur more frequently than the other residue classes. For example, for $`\mathrm{}=7`$, the classes $`0,4,6`$ occur much more frequently than the classes $`2,3,5`$.

Closing remark: readers wishing to see further details of various arguments in this note are invited to look at the comments in the source code.


text:
# Morphologies and ages of star cluster pairs and multiplets in the Small Magellanic Cloud

## 1 Introduction

Star cluster pairs are common objects in the Magellanic Clouds and it is important to understand their formation and evolution processes. A list of 30 cluster pairs in the SMC was first presented by Hatzidimitriou & Bhatia (1990). Bica & Schmitt (1995, hereafter BS95) revised previous data on SMC extended objects (star clusters, associations and emission nebulae) and identified new ones using Sky Survey ESO/SERC R and J films. They presented a list of 40 pairs and 2 triple star clusters. Pietrzyński & Udalski (1999a) reported 23 pairs and 4 triplets derived from Pietrzyński et al.’s (1998) star cluster catalogue in the OGLE survey area. Bica & Dutra (2000) updated BS95’s catalogue considering the entries in Pietrzyński et al’s catalogue. Bica & Dutra (2000) indicated 75 pairs and multiplets comprising 176 individual objects.

In recent years there has been growing evidence of interacting star clusters in the Magellanic Clouds, especially in the LMC. Bhatia & Hatzidimitriou (1988) concluded that more than 50% of LMC pairs must be physical systems. Bhatia & McGillivray (1988) found evidence that NGC2214 is a merging binary star cluster, based on the presence of a flattened core and an extended halo. Indeed Sagar et al. (1991) detected two turnoff points revealing the presence of two interacting clusters. Bica et al. (1992) studied cluster pairs and multiplets in the LMC bar by means of integrated colours and found systems which resulted coeval and some with age differences. Vallenari et al (1998) confirmed such scenarios by means of colour-magnitude diagrams. They also presented isophotal contours indicating physical interaction. Several other studies have found binarity evidence in LMC cluster pairs (e.g. Kontizas et al. 1993, Dieball & Grebel 2000a, 2000b).

Bhatia et al. (1991) presented a photographic atlas of binary star cluster candidates in the LMC. For the SMC no morphological atlas is available and isophotal maps are required to test possible physical interactions. Indeed comparisons of isophotes with isopleth maps of N-body simulations proved to be a useful tool (Rodrigues et al. 1994, de Oliveira et al. 1998, hereafter ODB98), since the simulations produce features such as bridges, common envelopes and extensions. The observational importance of tidal tails as interaction signatures was also indicated by Leon et al. (1999).

In this work we provide isophotal maps for SMC pairs and multiplets to study the following properties of these candidate physical systems: (i) angular distribution; (ii) projected centre-to-centre separation of members; (iii) isophotal structures using the Digitized Sky Survey (hereafter DSS); (iv) ages derived by means of isochrone fitting, when possible. We discuss candidate physical systems and infer a scenario for their formation and evolution.

In Sect. 2 we gather the objects providing coordinates, sizes, centre-to-centre separations and other details for the SMC pairs and multiplets. In Sect. 3 we provide the isophotal contour atlas together with classifications of interaction features whenever present. A preliminary version of the present isophotal SMC atlas together with one for the LMC was given in de Oliveira (1996). In Sect. 4 we derive cluster ages by means of Colour-Magnitude diagrams extracted from the OGLE-II BVI photometric survey (Udalski et al. 1998). In Sect. 5 we discuss the properties of the systems and the possible scenarios for their origin and evolution. Finally, concluding remarks are given in Sect. 6.

## 2 SMC pairs and Multiplets

Table 1 shows data for the 176 SMC objects which form 75 star cluster pairs and multiplets (Bica & Dutra 2000), considering that the pair BS63/B67 in the latter study is possibly a triplet with NGC294 (ODB98). By columns: (1) object cross-identification in the different catalogues: N (Henize 1956), K (Kron 1956), L (Lindsay 1958), H (Hodge 1960), SL (Shapley & Lindsay 1963), B (Brück 1976), DEM (Davies et al. 1976), ESO (Lauberts 1982), H86- (Hodge 1986), MA (Meyssonier & Azzopardi 1993), BS (Bica & Schmitt 1995), OGLE (Pietrzyński et al. 1998). Note that some are embedded objects, named after the corresponding HII region; (2) and (3) right ascension and declination for the epoch J2000, respectively; (4) object type following BS95: C for star cluster, NC for small HII regions with embedded star clusters, CN for clusters which show some traces of emission, A for associations, CA and/or AC for star clusters of low density and objects with intermediate properties, AN for associations which show some traces of emission and NA for HII regions with embedded associations; (5) and (6) major and minor sizes, respectively; (7) position angle of major axis ($`0^{}`$=N, $`90^{}`$=E); (8) centre-to-centre angular separation (1 arcsec = 0.28 pc assuming an absolute distance modulus (m-M) = 18.9 for the SMC, Westerlund 1990). For triplets and multiplets we measured the separation between the two main members, as indicated in the corresponding object line in the table; (9) remarks: $`m6,m5,m4,mT`$ and $`mP`$ indicate groups with six, five and four members, triplets and pairs, respectively. A running number identifies each candidate system. Abbreviations ‘br\*’ indicates that a bright star is present; ‘att’ means attached to. A hierarchical indication is given for objects embedded in or superimposed on larger ones: ‘in’ suggests a possible physical connection while ‘sup’ suggests a projection. There are 56 pairs, 15 triplets, 2 quadruplets, 1 quintuplet and 1 sextuplet. The quintuplet is located in the star forming region NGC395 with dimensions 25 pc x 18 pc which in turn is embedded in the HII complex SMC-DEM126 with 66 pc x 35 pc according to angular sizes in BS95. The sextuplet is mostly contained in the OB association H-A35 (Hodge 1985) with dimensions 57 pc x 47 pc. One quadruplet is located in H-A60 with 44 pc x 24 pc. Note that dimensions (columns 5 and 6) and position angle (column 7) are from Bica & Dutra (2000) and were measured on plates. The dimensions are in general larger than the isophotal sizes in the subsequent analyses.

## 3 The Isophotal Atlas

For isophotal analysis purposes we extracted digitized images of pairs and multiplets from the DSS. The plates are from the SERC Southern Sky Survey and include IIIa-J long (3600s), V band medium (1200s) and V band short (300s) exposures. The PDS pixel values correspond to photographic density measures from the original plates, and are not calibrated. The images were processed with the IRAF package at the Instituto de Física - UFRGS, applying a 2-d Gaussian filter to smooth out individual stars, producing isodensity maps.

We show the isophotal atlas of SMC star cluster pairs and multiplets in Figs. 1 throughout 7, where morphological evidence of interactions can be searched for. In each panel we point members and provide designations.

### 3.1 Isophotal Distortions

Cluster pairs and multiplets with evidence of physical interaction are marked in column 6 of Table 2. Their isophotal maps show features such as isophotal distortions, common envelope, isophotal twistings, etc. We classified these isophotal maps according to the following criteria: (i) isophotes of the members are detached (d), but showing isophotal distortions; (ii) isophotes of the members are connected by a “bridge” (b); (iii) the members are embedded in the same isophotal envelope (e). These classification criteria are an important tool for a selection of interacting cluster pair candidates, since such isophotal morphologies are predominant in N-body model encounters (Rodrigues et al. 1994, ODB98). The classifications above were made for objects not disturbed by gas emission. It is worth noting that both bridges and envelopes imply a common outer isophotal, however the former includes a dimension which is considerably smaller than the smaller cluster diameter. The pair B59/SMC-N46 would be an envelope limiting case (Fig. 4).

As examples, the pair H86-186/IC1612 (Fig. 5) is embedded in the same isophotal envelope and it is classified as e. The pair H86-140/H86-139 (Fig. 4) shows a bridge linking its members and is classificated as b. In Fig. 6 the pair NGC376/BS114 shows isophotal distortions, however the outer isophotes of the components are not connected and it is thus classified as d.

We detected relevant isophotal features for about 25% of the sample distributed as follows: 6 candidate cluster systems with common envelope, 7 with bridge and 5 detached cases.

High isophotal densities such as those observed for B78/L51 (Fig. 5) are related to the high surface brightness which often occurs in blue clusters.

It is important to note that non relevant isophotal features can appear in the maps caused by bright stars or background fluctuations. By means of previous inspections in the original images we verified the non relevant objects. For example, the compact object in the lower right corner of the B39/BS30 isophotal map (Fig. 2) is a bright star, however the compact one to the right of the main cluster can be a star or a compact cluster, as seen in the original form. Higher resolution images would be necessary to check the presence of a third cluster.

## 4 Ages of cluster members

Udalski et al. (1998) provided a BVI photometric survey of the SMC central $``$ 2.4 square degrees with results for about 2.2 million stars. The data were collected during the OGLE-II microlensing search project (Udalski et al. 1997).

Pietrzyński & Udalski (1999b) studied colour-magnitude diagrams (CMDs) extracted from the BVI database for 93 SMC star clusters, taking into account neighbouring field extractions for comparisons. They estimated reddening values for each cluster using: (i) the mean I band magnitude of red clump stars in the cluster neighbourhood; (ii) the assumed extinction-free magnitude of the SMC’s red clump stars, I = 18.34 mag (Udalski 1998). They derived ages using Bertelli et al.’s (1994) isochrones, adopting as a rule for the SMC a metallicity Z = 0.004, and an absolute distance modulus (m-M) = 18.65.

Considering the OGLE-II angular coverage we conclude that out of the present 176 objects 133 are therein included. We obtained from the database V and I band CMDs for these objects, using their coordinates and diameters to define a box extraction. Guided by DSS images of each system we selected representative field regions to extract stars and construct CMDs for comparison purposes. Ages and reddenings for the clusters were estimated by fitting the Padova isochrones. We adopted Bertelli et al.’s (1994) isochrones rather than the new isochrones by Girardi et al. (2000), since the previous set includes younger ages. The isochrone age grid suitable for 4, 5, 6.3, 7.9, 10, 12, 16, 20, 25, 32, 40, 50, 63, 79, 100, 126, 160, 199, 251, 316, 400, 500, 630, 790 and 1000 Myr. The isochrone age range and resolution allow one to estimate errors taking into account the stellar statistics in the cluster and field CMDs. We also assumed a SMC metallicity of Z = 0.004, a foreground galactic reddening of E(B-V)f = 0.03 in the SMC direction and E(V-I)/E(B-V) = 1.31. Table 2 shows the ages, in column 2 and reddening values (foreground plus SMC internal) in column 3 determined in the present work via isochrone fitting. Ages and reddening values derived by Pietrzyński & Udalski (1999b) are in column 4 and 5 respectively, when available for comparisons. In our sample we could estimate CMD ages for 91 objects, 40 of them in common with Pietrzyński & Udalski (1999b). We conclude that despite the reddening method and distance modulus differences between the two studies, there is good overall agreement for the age determinations. We also included ages for 21 objects embedded in HII regions, classified as NA and NC (see column 4 of Table 1, and also BS95). We assumed for them an age of 3 Myr.

## 5 Discussion

We overlap two age histograms in Fig. 8. The dotted line histogram shows objects from Table 2 with ages via CMDs. There are three peaks, at about 80, 220 and 450 Myr. The first two peaks appear to be present in Grebel et al.’s (2000) sample for 200 SMC clusters based on OGLE-II data, reported at t $``$ 100 and 200 Myr respectively. They suggested them as enhanced star formation epochs. Considering Pietrzyński & Udalski’s (1999b) sample with CMD ages for 93 SMC clusters, no peak is seen at t = 80-100 Myr, having their histogram a maximum at the youngest bin. In their histogram there occurs a local maximum at t $``$ 250-300 Myr, but the statistics is low. The third peak in the present study has no counterpart in both previous analyses, and it is probably an artifact from the low isochrone age step resolution in that range.

The solid line histogram in Fig. 8 shows additionally the embedded objects in HII regions for which we assumed an age of 3 Myr, and the objects from Table 2 with lower or upper limits by assuming them as the ages themselves. We can observe again a three peak distribution but the first peak is now shifted to the youngest bin similarly to Pietrzyński & Udalski’s (1999b) histogram. This suggests that only the 200 Myr peak is relevant, being related to the SMC/LMC last encounter (Gardiner et al. 1994, Grebel et al.’s 2000). The maximum seen in the youngest bin is possibly related to the cluster formation/destruction rates. Since the present sample deals with pairs and multiplets a fast destruction rate might be caused by the internal dynamical evolution in each cluster complex, caused by merger and/or other effects.

Approximately 55% of cluster pairs and multiplets in Table 2 present similar ages between their members indicating that they are coeval. This suggests that most of the pairs and multiplets had a common origin, possibly from the same molecular complex. Note that about 60% of the pairs and multiplets embedded in OB associations (H-A, Hodge 1985), as indicated in column 9 of Table 1, have comparable ages between their members. This could be an explanation for the origin of cluster systems.

Considering triplet and multiplet members, we found that about 70% of them are younger than 100 Myr. These results suggest a possible binary cluster formation scenario: clusters can be born in multiplet systems and coalesce by mergers and tidal disruptions forming binary clusters in a timescale of $``$ 100 Myr. This time is in agreement with dynamical times required for an interacting pair to merge into a single cluster (ODB98, de Oliveira et al. 2000).

The pairs with a bridge in the isophotal maps have comparable ages for their components (Table 2). As examples, the pairs NGC241/NGC242 and B39/BS30 in Fig. 2 show a bridge linking their members which could be interpreted as a sign of interaction (see the similarity with the N-body simulation model in Fig. 11 of ODB98). A typical timescale for the bridge phenomenon is $``$ 30 Myr, as deduced from a series of N-body simulations related to bridge formation and evolution (ODB98).

Another interesting isophotal feature is related to the cluster triplet NGC220/NGC222/B23 which shows distortions for the small cluster in a direction almost perpendicular to the line connecting itself to the large components NGC220 and NGC222. This configuration and morphologies are compatible with a fast hyperbolic encounter with small impact parameter (e.g. Fig. 12 of ODB98).

Fig. 9 shows the angular distribution of pairs and multiplets together with SMC HI column density isophotes from Mathewson & Ford (1984). It can be seen that most of the objects are concentrated in the SMC main body, close to the higher concentration of HI, so it is not unexpected that in general they result young (Sect. 4). The objects appear to be distributed along an axis. Such distribution is present in the overall SMC cluster sample and there is growing evidence that it is related to a nearly edge-on disk containing the bulk of the young stellar population in the SMC (Bica et al. 1999, Westerlund 1990).

A nearly edge-on disk in a low internal reddening galaxy like the SMC would imply an increase of projection effects as compared to a simulation such as that carried out by Bhatia & Hatzidimitriou (1988) for the nearly face-on LMC disk where the physical pairs would be about 50%. Consequently the fraction of physical pairs in the SMC would be lower. The present approach including morphological evidence of interaction is an attempt to constrain this aspect. Indeed the fraction with isophotal distortions is only 25% (Sect. 3.1). Projection effects can be responsible for the age spread in some multiplets. For example the sextuplet (Table 2) has component ages 80 Myr (B72), 200 Myr (H86-143), 50 Myr (BS257), 3 Myr (SMC-N52A and SMC-N52B) and 400 Myr (H86-148). Possibly only the 3 or 4 younger components could be related to OB-Association H-A35, the remaining objects would be captures or projection effects. This age spread is also present in IC1611’s quadruplet and in some triplets. On the other hand the quintuplet in the star forming complex NGC395 has all its members with the same age (3 Myr) thus forming a physical system.

In Fig. 10 we show the distribution of the diameter ratio between members for all pairs in the sample. The diameter ratio is mostly in the range $`12`$, with a peak at $`1`$ indicating that the majority of pair members have a comparable size. This effect was also observed in the LMC (Bhatia et al. 1991).

Fig. 11 shows the distribution of the centre-to-centre angular separation between pair members. The separation range is $``$ $`1080`$ arcsec ($``$3 - 22 pc) with a pronounced peak at $``$ $`45`$ arcsec (13 pc). A similar peak was also observed by Bhatia et al. (1991) and de Oliveira (1996) who found a bimodal distribution for the LMC pairs with peaks at $``$ 5 and 13 pc. The observed upper limit of the projected centre-to-centre linear separation $``$ 23 pc (80 arcsec) is comparable to Bhatia & Hatzidimitriou’s (1988) separation criterion for pairs in the LMC (18.7 pc). The frequent separation value around 13 pc may reflect a preferred survival distance for the systems, combined to projection effects.

## 6 Conclusions

We presented an isophotal atlas for 75 star cluster pairs and multiplets in the Small Magellanic Cloud, comprising 176 objects.

It was possible to derive ages from Colour-Magnitude Diagrams using the OGLE-II photometric survey for 91 objects. In addition we included in the analysis ages for 21 embedded objects in HII regions. The age distribution has a maximum in the youngest bin with a profile related to the cluster formation/destruction rates, in particular cluster multiple systems can have a fast destruction rate caused by their internal dynamical evolution. There is a second peak around 220 Myr which is probably related to the SMC/LMC last encounter.

We find that 55% of the pairs and multiplets in the sample are in general coeval, indicating that captures are a rare phenomenon. Most of the cluster multiplets occur in OB stellar associations and/or HII region complexes which indicates a common origin and suggests that multiplets coalesce into pairs or single clusters in a short time scale ($``$ 100 Myr).

The majority of the cluster members have comparable sizes, with a diameter ratio ranging mostly between $`12`$. The projected separation distribution between the members of a pair has a pronounced peak at $``$ 13 pc. These observational results are important constraints to theoretical models of star cluster pair encounters and could be related with the formation process and subsequent dynamical evolution of cluster systems.

The angular distribution of cluster pairs and multiplets shows that most of the objects are located in the SMC main body. The overall SMC cluster sample presents a similar distribution and there is evidence that it is related to a nearly edge-on disk in the SMC. Considering this, it is expected an increase of projection effects as compared to estimates for the LMC disk where physical pairs would be about 50% (Bhatia & Hatzidimitriou 1988).

The atlas shows that around 25% of the isophote maps present relevant structures like bridges, common envelopes or detached distorted isophotes. N-body simulations have indicated that these structures arise from interactions between the members of the cluster systems. Indeed cluster pairs as NGC241/NGC242 and B39/BS30 show in their isophotal maps bridges linking their members and have comparable ages.

We conclude that multiplicity may have an important role in the early dynamical evolution of star clusters in general, and signatures of that may survive in the long term structure of large single clusters (de Oliveira et al. 2000).

###### Acknowledgements.

We thank the referee Dr. B. Westerlund for interesting remarks. We acknowledge support from the Brazilian institutions CNPq, CAPES and FINEP.


text:
# Clustering and fluidization in a one-dimensional granular system: molecular dynamics and direct-simulation Monte Carlo method

## I Introduction

Granular systems have been extensively studied due both to the theoretical challenges they present (for a recent review see ) and to the applications of industrial importance that spring from the rich phenomena they exhibit (see  and references therein). These systems are characterized by an energy loss in collisions. This loss is at the base of many interesting phenomena, such as inelastic collapse, where the particles collide infinitely often in finite time, and clustering (for a sample of theoretical, simulational, and experimental approaches see ). Different methods for keeping the system from collapsing have been devised, such as subjecting the particles to Brownian forces, and forcing through the boundaries by putting the system in a box with one or more thermal-like walls (see for example ). This work focuses on the latter method.

Being one of the simplest types of forcing, several authors have studied a one-dimensional system in a box with one or two heated (stochastic) walls. Of these,  study cluster formation, although  are not strictly one-dimensional.

This article studies a quasielastic one-dimensional system not subject to gravity between two thermalizing walls. We focus on two control parameters: the total inelasticity parameter $`qNN(1r)/2`$, where $`N`$ is the number of particles and $`r`$ is the restitution coefficient, and the externally imposed temperature gradient. The parameter $`qN`$ has been shown to be relevant for the quasielastic system. By varying these parameters we determine the region in parameter space where clustering is fully inhibited, obtaining a fluidized state. We present a singular feature of the distribution function for the clustering regime, and then study how this feature is modified for the fluidized state.

In  the authors study a one-dimensional system of point-like particles between an elastic and a heated wall. They emphasize that a cluster inevitably forms away from the heated wall, regardless of how elastic the system is (as long as it is not perfectly elastic). They also study the same system, but with both walls expelling the particles with a fixed velocity. In this case they find that the cluster forms away from the walls and roams slowly about the system, with two groups of fast particles connecting the cluster with the “heated” walls.

In  the same system is studied for different types of boundary conditions at the heated wall. The stochastic boundary condition studied has the form of a power of the velocity times the “thermal” condition (the one that produces a Maxwell-Boltzmann distribution in the elastic case). The authors show that when the power that multiplies the thermal condition is positive the test-particle equation (derived from the Boltzmann equation) has a steady-state solution. Thus the thermal case does not have a steady state and develops a cluster away from the heated wall. The mechanism for the growth of the cluster is explained and verified numerically.

In  a similar system is studied: a long thin pipe of inelastic hard disks with heated walls (at the same temperature) at the ends of the pipe and periodic side walls. The pipe is thin enough for the particle order to be preserved. The probability distribution for the distance between the central particles is studied. This distribution gives a markedly denser system near the center than in the elastic case, although the limit to the elastic case is smooth, unlike the strictly one-dimensional case of . In  the same author studies the velocity correlations that this system develops as inelasticity is increased, showing that a consistent description must take these correlations into account.

In this paper we revisit the one-dimensional system of $`N`$ point-like particles interacting via collisions that conserve momentum but dissipate kinetic energy. To fix notation, the particle velocities after a collision are given by

$$c_1^{}=qc_1+(1q)c_2,c_2^{}=(1q)c_1+qc_2,$$
(1)
where $`c_i`$ is the velocity of particle $`i`$ before a collision, and $`q=(1r)/2`$, being $`r`$ the restitution coefficient. For the elastic case ($`r=1`$) the particles simply exchange velocities. Since the particles are point-like, the system is then indistinguishable from a system in which the particles do not interact.

This one-dimensional system is interesting because dissipation is the first order correction to a free gas. Besides, results for the one-dimensional system have been found to have unexpected relevance for higher-dimensional problems. For example, in two-dimensions the particles involved in inelastic collapse lie roughly on a line. Also, the dissipation-induced temperature gradients calculated in  for the one-dimensional case inspired the authors to look for dissipation-induced Rayleigh-Bénard-like convection for a two-dimensional system without an externally imposed temperature gradient.

For a system with one thermal wall and open on the other side, under the influence of gravity, the quasielastic system may be kept fluidized: any cluster that starts to form is forced against the thermal wall, where it evaporates. In  the test-particle equation—which is the 1–D Boltzmann equation where the limit $`N\mathrm{}`$ is taken, but keeping $`qN`$ fixed—is successfully applied, with close matching of theory and simulations even at the level of the distribution function.

The one-dimensional system under study is left to evolve between two thermal walls at temperatures $`T_{\text{left}}`$ and $`T_{\text{right}}`$, with $`T_{\text{left}}T_{\text{right}}`$. We define the parameter

$$\mathrm{\Delta }\frac{T_{\text{right}}T_{\text{left}}}{T_{\text{right}}+T_{\text{left}}}$$
(2)
to quantify how far from the symmetrical case is the system. Under the same conditions, an elastic system has a perfectly bimodal velocity distribution with global homogeneous temperature equal to $`\sqrt{T_{\text{left}}T_{\text{right}}}`$. For the sake of comparison, we simulate systems with $`\sqrt{T_{\text{left}}T_{\text{right}}}=1`$. Here $`\mathrm{\Delta }=0`$ represents a symmetrical setting, while $`\mathrm{\Delta }=1`$ represents an infinitely strong temperature gradient.

As in , for $`\mathrm{\Delta }=0`$ a cluster unavoidably forms away from the thermal walls. After forming, the cluster performs an apparently random walk about the center, growing in size (and therefore mass) while the rest of the system grows more rarefied. With the decrease in density of the surrounding gas and the increased inertia of the cluster, an eventual collision with one of the walls is to be expected. When this happens part of the cluster evaporates, and what is left of it is expelled from the wall (see Fig. 1), thus restarting the growth process. Thus not only is the system highly clumped, but also in a non-steady state. Nevertheless, the gas that is far from the random-walk zone has a well-defined time average for the distribution function, as is seen in Fig. 2. A noteworthy feature of this distribution is that it exhibits apparently singular behavior for slow velocities.

Setting $`\mathrm{\Delta }0`$ the symmetry of the system is broken. For $`\mathrm{\Delta }1`$ the cluster performs a slightly asymmetric random walk, spending more time near the colder wall, and therefore colliding more often with it. Thus the cluster cannot grow as much as it did in the symmetric case before colliding with a wall. By increasing $`\mathrm{\Delta }`$ the cluster-wall collision frequency grows, even obtaining short “windows” in which the cluster completely evaporates. By further increasing $`\mathrm{\Delta }`$ these windows grow larger until a point is reached where no cluster forms. In this fashion a totally fluidized state is achieved, which may be tractable with the dissipative Boltzmann equation. However, the distribution function obtained from the molecular dynamics (MD) simulations exhibits a peculiar non-Gaussian feature for slow velocities. This feature is a smoothed version of the apparent singularity of the symmetric case. To discern whether this feature is due to correlations or is present before they settle in, we compared Newtonian molecular dynamics results with those obtained through direct-simulation Monte Carlo (DSMC), which neglects correlations. The results agree very well, except when the system approaches the clustering regime.

## II Simulation method

We simulate the system through event-driven molecular dynamics and through direct-simulation Monte Carlo. The direct-simulation Monte Carlo procedures use the null-collision technique where, overestimating the collision frequency (using the maximum relative velocity within a cell), the number of collisions to be attempted is calculated through a Poisson process. In the next step the collisions are attempted, choosing at random two particles within the cell, and making them collide with a probability proportional to their relative velocity. Most of the molecular dynamics and DSMC simulations were done with $`N=1000`$ particles.

In the MD simulations we detect clusters using a geometric criterion: we consider chains of particles that are nearer than a critical distance (in our case $`10^6`$ and $`10^5`$, to be certain that the conclusions are independent of the choice). The system length is one, and with a thousand particles the mean distance between neighbors for a homogeneous system is $`10^3`$. Thus we detect particles that are uncommonly near by three orders of magnitude. We discard chains of length three or less, since they may be random encounters. Measuring the total length of the cluster we have found that on average it is of the order of $`10^5`$; thus the choice of $`10^5`$ as link-link distance is much larger than the true distance between them.

As already stated, the boundary conditions are such that the (homogeneous) temperature of the corresponding elastic system (equal to $`\sqrt{T_{\text{left}}T_{\text{right}}}`$) is one.

## III Results

### A Clustering regime

Figure 1 shows the non-steady state of a granular system between two walls at the same temperature for $`qN=0.01`$. A cluster forms away from the walls, performing a random walk of varying amplitude. When the cluster reaches a wall, part of it evaporates, and the growth process begins anew.

As is usual for the quasielastic case, we relabel the particles when they collide. This enables us to visualize this system as a group of barely interacting particles passing through each other.

The picture for cluster evolution, as explained before, is the following: the cluster grows because the slowest particles, due to the asymmetry of the distribution function, drift towards the cluster. As it grows, the density of the gas surrounding it decreases, with the consequent saturation in growth. Thus we have a “Brownian particle” of increasing mass moving in an increasingly rarefied medium. This “particle” will be increasingly less affected by the surrounding medium, until it can no longer be kept away from the walls.

The cluster moves several orders of magnitude slower than the thermal speed (three orders of magnitude in Fig. 1). Upon reaching a wall, the front liners strike the wall and are expelled by it much faster than the other cluster members. These particles pass through the cluster, transferring momentum to it, as described in . Thus these fast particles push the cluster away from the wall, where it can absorb particles again. The fast particles, however, no longer belong to the cluster.

Since the slowest particles in the gas are the ones that will be absorbed by the cluster, it is the number of slow particles in the gas that will determine the cluster’s growth rate. After the cluster strikes a wall, the expelled particles will be fast particles, and they will not contribute to the growth of the cluster during the time it takes for the gas to cool down again: the only particles available for absorption are the ones that were available before the cluster-wall collision. This explains why the cluster keeps growing at approximately the same rate it did before the collision. After the gas has cooled down, the growth rate returns to its normal value. This is the end of the plateau seen in Fig. 1 after each cluster-wall collision.

To quantify the evaporation process we proceed as in : as soon as the first particle belonging to the cluster reaches a wall, it is expelled with a speed much higher than the cluster velocity; thus we may consider the ideal situation of a cluster of $`M`$ particles at rest being stricken by a fast particle with velocity $`v`$ (in this case $`v1`$). After colliding with the first particle in the cluster, the new fast particle’s speed will be $`(1q)`$. Thus, after traversing the cluster, the fast particle’s velocity will be $`(1q)^N`$. Since momentum is conserved in collisions, the center of mass of the cluster will have acquired a speed of

$$v_{\text{CM}}=\frac{1(1q)^N}{N}=\frac{1\left(1\frac{qN}{N}\right)^N}{N}.$$
(3)
Considering the case $`N1`$ with fixed $`qN`$, as in , we may simplify this expression to

$$v_{\text{CM}}\frac{1e^{qN}}{N}.$$
(4)
By further considering the case $`qN1`$ we get $`v_{\text{CM}}q`$. In this limit, if the cluster reaches the wall with velocity $`v_0`$ and is expelled from it with velocity $`v_1`$, the number $`m`$ of particles evaporated will satisfy $`|v_1v_0|mq`$. For the situation shown in Fig. 1, $`v_0`$ and $`v_1`$ are typically of the order of $`10^3`$ and $`q=10^5`$, hence the number of particles evaporated will be of the order of $`10^2`$.

Even if the system is in a non-steady state, the gas at the walls (far from the random-walk zone) has a well-defined time average for the distribution function. The distribution function at the left wall is shown in Fig. 2. There is an apparent singularity for slow velocities. The distribution is asymmetric as it should, since the particles leaving the wall ($`c>0`$) follow a Gaussian distribution. Figure 2 also shows the distribution multiplied by $`\sqrt{|c|}`$. Since the limit of $`\sqrt{|c|}f`$ for $`c0^{}`$ is finite and nonzero, we conclude that the distribution function exhibits a singularity that behaves like $`|c|^{1/2}`$ for slow velocities.

As shown in , when $`\mathrm{\Delta }=0`$ (the symmetric case) the distribution shown in Fig. 2 is not a solution of the steady-state test-particle equation:

$$c_xf(x,c,t)=qN_c\left[f(x,c,t)M(x,c,t)\right],$$
(5)
where

$$M(x,c,t)=_{\mathrm{}}^{\mathrm{}}f(x,c^{},t)(cc^{})|cc^{}|𝑑c^{}.$$
(6)
To establish this, let us study the behavior for small $`c`$ of a solution of this equation. Assume that, for small $`c`$, $`f(x,c0)f_0(x)c^\alpha `$, with $`\alpha <1`$ in order to have a finite density in the vicinity of $`x`$. Furthermore, assume that $`M(x,c0)M_0(x)c^\beta `$. Inserting this behavior in Eq. (5) we obtain

$$_xf_0qNM_0f_0(\beta \alpha )c^{\beta 2}.$$
(7)
Thus, in order to keep $`f_0`$ (the amplitude of the singularity) finite, we must have either $`\beta =\alpha `$ or $`\beta 2`$. Integrating the distribution of Fig. 2 we obtain $`\beta =1`$. Since $`\beta <2`$, we must have $`\alpha =1`$. But this corresponds to a nonintegrable distribution, and therefore the distribution cannot be steady.

### B Inhibition of cluster formation

Figure 3 shows the regions in $`(\mathrm{\Delta },qN)`$-space where clustering is inhibited for $`N=1000`$. As is to be expected, as the inelasticity increases, a stronger temperature gradient is necessary to inhibit cluster formation.

To discern whether the non-Gaussian features of the velocity distribution function are derived from correlations in the system we compared results from MD simulations (full Newtonian dynamics) with results from DSMC simulations (no velocity correlations assumed). Figures 4 and 5 show this comparison for a case far from the clustering threshold ($`\mathrm{\Delta }=0.6`$ and $`qN=0.1`$). The temperature of the left and right walls are chosen so that the global temperature for the elastic case (equal to $`\sqrt{T_{\text{left}}T_{\text{right}}}`$) is one. The curves match almost exactly.

At the level of the distribution function, the results also match closely. The peculiar non-Gaussian feature of the distribution function is clearly seen in Fig. 5. There is some slight mismatch near the peak.

For the fluidized case, since momentum is conserved and the system is stationary, the pressure is constant throughout the system. The number density and the granular temperature calculated here are related by $`p=nT`$ ($`T`$ in energy units). Thus when the normalized density $`n/N`$ varies little throughout the system ($`n/N1+ϵ(x)`$), the normalized temperature is

$$\frac{T}{T_0}=\frac{p/n}{p/n_0}=\frac{n_0/N}{n/N}=\frac{1}{1+ϵ(x)}1ϵ(x),$$
(8)
thus obtaining the nearly symmetric profiles seen in Figs. 4 and 6.

Figures 6, 7, and 8 compare the MD and DSMC results for cases near cluster formation. The non-Gaussian feature of the distribution function shows a systematic deviation for DSMC simulations: there is overpopulation for slow velocities. This is explained by considering that the DSMC method, like the Boltzmann equation, neglects correlations. When the system approaches the clustering regime, increased dissipation induces correlations which tend to make the particles collide less. In DSMC these correlations are neglected, with the corresponding systematic overestimation in the collision frequency. This overestimation results in a lower temperature of the system about the density peak.

## IV Conclusions

We have shown that a system not subject to gravity between thermal walls unavoidably reaches a non-steady state when the walls are at the same temperature. A cluster forms in the bulk, slowly roaming about the system while absorbing particles. As it grows, the amplitude of the random walk increases, until at last the surrounding gas cannot keep the cluster away from the walls. When the cluster reaches a wall, a part of it is ejected by the wall through the cluster (relabeling the particles on collisions), effectively pushing the cluster away from the wall, and leaving it to grow again.

Most of the time the cluster is far from the walls. Thus measuring the distribution function at a wall is measuring the distribution function of the gas that surrounds the cluster. This distribution function has a well-defined time average, and exhibits apparently singular behavior for slow particles, diverging like $`|c|^{1/2}`$.

Imposing an external temperature gradient forces the cluster against the colder wall, inhibiting its growth. Increasing the temperature difference leads to a system in which the cluster never forms: the system is completely fluidized. The distribution function of the gas exhibits peculiar non-Gaussian features: a smooth version of the aforementioned singularity. Therefore, any attempt at solving the Boltzmann equation through moment methods must consider this feature in the initial ansatz, as is done for the problem of an infinitely strong shock wave in  and . In fact, a solution for this problem was attempted using the four moment method of . As mentioned in , the fourth balance equation could not be freely chosen when the boundary conditions were symmetric: some choices gave undefined results. As is natural, by not including the non-Gaussian feature in the ansatz for this calculation we obtained absurd results, such as higher temperature in the middle of the system than near the walls.

We compared molecular dynamics with direct-simulation Monte Carlo. Agreement between these two methods shows that the non-Gaussian feature of the distribution function may be predicted by the dissipative Boltzmann equation. As the system approaches cluster formation, correlations settle in. These correlations reduce the collision frequency among particles. DSMC neglects these correlations, and thus overestimates the number of collisions. This exaggerates the effects of dissipation, producing steeper profiles.

###### Acknowledgements.

We thank Rodrigo Soto, Aldo Frezzotti, and Rosa Ramírez for helpful discussions. This work has been partially funded by Fundación Andes through a doctoral scholarship, Fondecyt through grants 2990108 and 1000884, and by FONDAP through grant 11980002.