text:
# Fast Monte Carlo simulations and singularities in the probability distributions of non-equilibrium systems

## Abstract

A numerical technique is introduced that reduces exponentially the time required for Monte Carlo simulations of non-equilibrium systems. Results for the quasi-stationary probability distribution in two model systems are compared with the asymptotically exact theory in the limit of extremely small noise intensity. Singularities of the non-equilibrium distributions are revealed by the simulations.

The understanding of fluctuations in systems away from thermal equilibrium is a problem of long standing in statistical physics Onsager:53 . Well known examples of physical phenomena in which non-equilibrium fluctuations play a particularly important role include e.g. switching of polarization in lasers Keay:95 , switching between different configurations in proteins Serpersu:83 , the transition to instability in Josephson junctions Kautz:96 , and chemical reactions Dykman:99a .

In non-equilibrium systems, where symmetries of detailed balance are broken, no general methods exist for the calculation of even basic quantities like the probability distribution. This is a case where numerical and asymptotic theoretical methods for investigating the probability distribution are of particular importance.

Theoretical approaches, such as WKB-like or path-integral methods, are available in the limit of small noise intensity, $`D0`$ Freidlin:70b ; McKane:89 ; Dykman:90a ; Dykman:99a . In particular the theory suggests that a solution to the problem of non-equilibrium fluctuations requires an understanding of the dynamics of deviations from the steady state Onsager:53  and an analysis of singularities in the non-equilibrium potential Graham:84 ; Smelyanskiy:97 . Some ideas about how to extend the existing ($`D0`$ limit) theory for still small but finite noise intensities have recently been suggested Dykman:99b ; Hanggi:00a ; Bandrivskyy:02 .

The main numerical technique used to verify theoretical predictions, and to analyse the behavior of the dynamical system under study, is Monte Carlo simulation. The theory gives an asymptotically exact solution in the $`D0`$ limit. In contrast, $`D`$ in the numerical simulations is necessarily finite. Typically, the time required for Monte Carlo simulations grows exponentially as $`D0`$. This meant that theoretical predictions of interesting singular structures, and of the non-equilibrium probability distribution Graham:84 ; Jauslin:87 , for long remained untested either experimentally or by numerical simulation. Moreover there was no clear understanding of how the picture changes for small but still finite noise intensities.

Approaches that have been tried to speed up the simulations have focused mainly on finding optimal fluctuational paths and rates of transition between stable states of a system (e.g. efficient transition path sampling Dellago:98  and dynamics importance sampling Woolf:98 , following the earlier suggestion of Pratt:86 ). In Crooks:01  the path sampling method was adapted for non-equilibrium systems. Based on the same idea, the umbrella sampling technique was suggested to estimate the probability of reaching any point in the phase space of an equilibrium system starting from a fixed initial state Dellago:98 . A technique for improving sampling in equilibrium systems by splitting up the probability packets was introduced in Huber:96 . So far, however, no general algorithm has been suggested, able to give both the whole probability distribution and dynamical information like the optimal fluctuational paths for small noise intensities for non-equilibrium systems.

In this Letter we introduce a numerical method that enables the time required for Monte Carlo simulations to be reduced by an exponentially large factor. It is applicable to generic two-dimensional non-equilibrium systems, does not require any a priori knowledge about the system apart from its dynamical equations of motion, and it allows the quasi-stationary probability distribution to be computed directly over the whole phase space. Using this method, we reveal for the first time singular behavior of the non-equilibrium distribution in numerical simulations, and we show that the results are in quantitative agreement with the asymptotic theory.

The central idea is to perform the simulations in sequential steps. We construct the quasi-stationary distribution, patching together intermediate results: we start from one of the steady states and gradually move away from it. We find that the time required for the simulations at each step is reduced by an exponentially large factor as compared to the standard technique: if the time necessary for a conventional Monte Carlo simulation technique is $`T`$, our modified method would require only time $`T_mNT\mathrm{exp}^{(N1)\frac{\mathrm{\Delta }\mathrm{\Phi }}{D}}`$, where $`N`$ is the number of steps involved and $`\frac{\mathrm{\Delta }\mathrm{\Phi }}{D}`$ is their separation in terms of the logarithm of the probability 111It is illuminating to compare simulations of the Duffing system by the fast and conventional techniques: the results of Fig. 3 took us $``$15 minutes to simulate with $`D=0.02`$, whereas for the same noise intensity it takes $``$4 days of standard Monte Carlo simulation to obtain comparable statistics close to the boundary of attraction..

We first explain the method on a very simple equilibrium stochastic system, and then we apply it to two much-studied non-equilibrium systems and compare the numerical results with theoretical predictions. To illustrate the technique, we consider an overdamped Brownian particle moving in a bistable Duffing potential $`U(x)=x^2/2+x^4/4`$

$`\dot{x}=U^{}(x)+\xi (t),`$ (1)
where $`\xi (t)`$ is zero-mean white Gaussian noise with intensity $`D`$ and moments

$`\xi (t)=0,\xi (t)\xi (0)=2D\delta (t).`$
The probability distribution is completely defined by the potential $`U(x)`$, and is of the Boltzmann form $`\rho (x)\mathrm{exp}(U(x)/D).`$ As in the case of a non-equilibrium system (where the probability distribution is not defined by a potential) a standard Monte Carlo technique can be used to deduce $`\rho (x)`$. Numerical integration rm:00  of the Langevin equation (1), assuming the system to be located initially at one of the potential minima $`x_m`$, gives the discrete probability distribution $`\rho (x)`$, peaked at $`x_m.`$ The potential can be deduced as $`\mathrm{\Phi }(x)D\mathrm{ln}\rho (x)`$. If the noise intensity is very small, the system fluctuates in a close vicinity of $`x_m`$ and large deviations from it are extremely rare. Accordingly, the conventional Monte Carlo technique cannot be used to study the dynamics of optimal escape paths, or the properties of the probability distribution far from the potential minima: for small noise intensities the statistics required cannot in practice be collected within a realistic time.

In order to overcome this problem, we start from the distribution already obtained near $`x_m.`$ We fix two probability levels $`\rho _i`$ and $`\rho _f`$, lying well within the region where the numerical $`\rho `$ is accurate, with $`\rho _f<\rho _i`$ corresponding to two levels in the potential $`\mathrm{\Phi }_i`$ and $`\mathrm{\Phi }_f`$, and two coordinates $`x_i`$ and $`x_f`$, as shown in Fig. 1. We require the levels $`\rho _i`$ and $`\rho _f`$ to be fairly different, such that the corresponding $`x_i`$ and $`x_f`$ are sufficiently separated: the distance between them must exceed $`\sqrt{Dh}`$, where $`h`$ is the integration time step used in the Monte Carlo simulation, and must also exceed the discretization step $`\mathrm{\Delta }x`$ in the discrete probability distribution.

The next step of the simulation is started from the level $`\mathrm{\Phi }_i`$ (putting the system at $`x=x_i`$ as its initial condition). If the system starts to evolve along a fluctuational trajectory (towards the boundary of attraction) we just follow its natural dynamics according to (1) and collect the statistics for building the probability distribution in the usual way. If the system starts with a relaxation trajectory (towards $`x_m`$), or when it crosses the boundary $`x_i`$ due to relaxation some time later, we stop the simulation and reinject the system back to the initial state $`x_i.`$ In this way we simulate the full dynamics of the system at higher levels of the potential $`\mathrm{\Phi }(x)>\mathrm{\Phi }_i`$ (in the region of coordinate space $`x>x_i`$ for this particular case). Thus, in the subsequent simulation step we follow only those fluctuations that have already attained a certain level in the potential $`\mathrm{\Phi }_i`$, without waiting for this exponentially slow event to happen. In this way, a new piece of the probability distribution is built with a time saving $`\mathrm{exp}\mathrm{\Phi }_i/D`$ compared to a simulation starting from the potential minimum $`x_m`$. The computed new piece of the potential $`\mathrm{\Phi }_2(x)`$ is shown as the upper curve in Fig. 1.

The two pieces of the inferred potential (the original $`\mathrm{\Phi }_1(x)`$ and the new $`\mathrm{\Phi }_2(x)`$) are then merged at $`x_f`$ by a simple shift. Continuing this procedure, the probability distribution and the corresponding potential can be built, step by step, for the whole region of interest. The inset in Fig. 1 shows the resultant potential, built from 13 such pieces between the minimum at $`x_m=1`$ and the maximum at $`x=0`$. It coincides closely with the Duffing potential $`U(x)`$ itself. The potential $`\mathrm{\Phi }(x)`$ is thus inferred within a region of coordinate space that is inaccessible in a conventional simulation (shown as bold curve for comparison). We stress that no a priori knowledge of the dynamics has been used in the simulations, and that the method is robust to choice of parameters.

In the case of a two dimensional system, the procedure remains essentially the same. The main difference is that, instead of identifying two points $`x_i`$ and $`x_f`$, we need to identify two closed lines of constant probability. One line is a boundary line for starting simulations from, and the other is a reference line for matching together different pieces of the probability distribution (see Fig.2 for clarification) 222In the case of an $`N`$-dimensional system, this would correspond to $`(N1)`$-dimensional surfaces of constant probability dividing the phase space into separate regions.. The crucial point of our technique is that, in starting the simulations from the boundary line, we must not perturb the natural dynamics of the system. This implies that we should consider the reinjection location probability (RLP) along the boundary line corresponding to $`\rho _i`$. Starting from the second step of the simulations, the system should be reinjected back according to the RLP after it relaxes across the boundary. We emphasize that the RLP is not the same as the probability distribution $`\rho (𝐱)`$, which is constant on the boundary line. The RLP is an additional important measure which describes local discrete dynamics in the neighborhood of the boundary line. It is a distribution along the boundary of how often the system crosses it.

In an equilibrium system, detailed balance provides a symmetry that can be used to reinject the system back at the boundary level, without any need to compute the RLP. For non-equilibrium systems, however, this procedure is inapplicable. The RLP should be considered separately (and calculated explicitly) for the particular system being investigated. It can be obtained from an analysis of the finite difference equation corresponding to the model. In the limit of small integration time step the probability to cross the boundary is proportional to the diffusion-related term in the finite difference equation. Then the RLP is simply proportional to the projection of the vector orthogonal to the boundary onto the coordinate affected by the noise $`\xi `$. It can also be computed numerically.

For non-equilibrium systems, the limit of small noise intensity is of particular importance. A sufficiently small $`D`$ gives rise to the possibility of revealing the non-equilibrium potential

$`\mathrm{\Phi }(𝐱)=\underset{D0}{lim}D\mathrm{ln}\rho (𝐱),`$
directly through a numerical experiment. Observations of the predicted singular shape of $`\mathrm{ln}\rho (𝐱)`$, and of its dependence on $`D`$, are thus of considerable interest.

As a first example of a nonequilibrium system, consider the periodically-driven overdamped Duffing oscillator

$`\dot{x}=U^{}(x)+A\mathrm{cos}\mathrm{\Omega }t+\xi (t).`$ (2)
We infer $`\mathrm{\Phi }(x,t)`$ as $`D\mathrm{ln}\rho (x,t)`$. This quantity corresponds to the theoretical “global minimum of the modified action” in the Hamiltonian theory of large fluctuations Bandrivskyy:02  and, in the limit $`D0`$, it becomes the non-equilibrium potential.

The complete $`\mathrm{\Phi }(x,t)`$, constructed from 12 such pieces, is shown in Fig. 2 and a time section of $`\mathrm{\Phi }(x,t)`$ calculated for different noise intensities together with the results of theoretical calculations (Hamiltonian theory including the prefactor) Bandrivskyy:02  is shown in Fig. 3. The RLP in the simulations can be taken as constant if a small enough integration time step is used in the scheme. A small difference between the theory and the simulations results appears for larger noise intensities then the asymptotic theory starts to break down.

As a second, more complicated, nonequilibrium example, consider the inverted Van-der-Pol oscillator

$$\ddot{x}+2\eta (1x^2)\dot{x}+\omega _0^2x=\xi (t)$$
(3)
Here, in order to be able to merge more easily the different pieces of $`\mathrm{\Phi }(x,y)`$, we apply a coordinate transformation from $`x`$ and $`y=\dot{x}`$ to amplitude $`A`$ and phase $`\varphi `$ ($`x=A\mathrm{cos}(\varphi ),y=A\omega _0\mathrm{sin}(\varphi )`$). The probability $`\rho (x,y)`$ can be then analyzed in the $`(A,\varphi )`$ coordinate space. This makes the problem very similar to the periodically driven Duffing oscillator: the only difference is the RLP which, in the case of the Van der Pol oscillator, turns out to be strongly modulated. It is essential for this modulation to be taken into account when reinjecting the system back to the boundary of constant probability. Two sections of $`\mathrm{\Phi }(x,y)`$, obtained from the simulations for different parameters $`\eta `$, are compared with the theory in Fig. 4. Again, the agreement between numerics and theory is excellent.

The non-equilibrium systems considered in this Letter share the same structure of singularities. Using the fast Monte Carlo simulations we reveal plateaus, the essentially flat regions in the probability distribution, which can be observed close to boundaries of attraction. They result from a purely dynamical effect that is not associated with the flatness of any potential. We have shown that its origin is related to switching between different types of optimal fluctuational path, and it is a general feature of non-equilibrium systems with metastable states Bandrivskyy:02 ; Bandrivskyy:02pre . The switching lines Smelyanskiy:97  are revealed as lines along which the “global minimum of the modified action” $`\mathrm{\Phi }(𝐱)`$ exhibits sharp bends – corresponding to the predicted line at which the non-equilibrium potential is non-differentiable. In the boundary region we found the oscillations of the probability distribution and their dependence on noise intensity (see the inset in Fig.3) discussed in the recent publications Dykman:99b ; Hanggi:00a ; Maier:01a . Using the simulations we demonstrated noise induced shift of the singularities and the optimal escape path, which has stimulated a new step in the development of the theory Bandrivskyy:02 .

We emphasize that the singularities can be confidently observed only in the limit of extremely small noise intensity, and therefore that the use of our new technique is crucial in that it reduces by an exponentially large factor the time required for Monte Carlo simulations. In addition to being fast, it preserves dynamical information, can be modified to analyse optimal fluctuational paths, is applicable to the energy-optimal control problem Khovanov:00 , and can be further extended to encompass higher dimensional systems and maps.

The work was supported by the Engineering and Physical Sciences Research Council (UK), the Joy Welch Trust (UK), the Russian Foundation for Fundamental Science, and INTAS.


format: text

text:
# Dense gas clouds and the Unidentified EGRET sources

## 1 Introduction

The nature of the UnIDentified (UID) EGRET sources (Hartman et al 1999) is one of the most interesting puzzles in contemporary high energy astrophysics. This puzzle persists primarily because of the dimensions – typically of order a degree – of the gamma-ray error boxes: large uncertainties in source position prohibit deep searches for counterparts at other wavelengths. We know, however, that this is only part of the problem: in some cases the EGRET source is bright enough that it can be located quite accurately, yet still no clear counterpart can be found (e.g. Mirabel et al 2000; Reimer et al 2001).

The majority of the discrete sources detected by the EGRET instrument are UID sources (Hartman et al 1999); it is not clear at present whether they are fundamentally similar to the known GeV emitters, or whether they represent an entirely new population. Because the UID sources are predominantly located at low to mid Galactic latitudes, they must be a Galactic population. This rules out one of the major classes of known discrete celestial GeV gamma-ray sources, namely the flat-spectrum radio quasars as these are extragalactic and hence isotropically distributed on the sky. Much attention has therefore focussed on the possibility that the UID sources are related to massive stars and their evolutionary end-points. In particular there has been intense interest in young pulsars as a contributor to the UID population (e.g. Bailes and Kniffen 1992; Kaaret and Cottam 1996; Roberts, Romani and Johnston 2001; Grenier and Perrot 2001); several of the identified GeV sources are young pulsars (Hartmann et al 1999). However, the observed latitude distribution of the UID sources is too broad to be readily compatible with young pulsars (Yadigaroglu and Romani 1997). Supernova remnants (Sturner, Dermer and Mattox 1996; Torres et al 2003), molecular clouds and massive star clusters (Montmerle 1979; Benaglia et al 2001) share this difficulty; but Gehrels et al (2001) have noted that the sky distribution of UID sources is consistent with a strong contribution from objects in Gould’s Belt, thus supporting a connection with massive stars.

Given the difficulty of detecting counterparts to the UID EGRET population, it is natural to consider the possibility that these sources may be connected with the dark matter problem. The currently popular cosmological model, the Cold Dark Matter model (e.g. Peebles 1993), stipulates that the dark matter is composed of weakly interacting particles; it is possible that these particles are massive – i.e. they are WIMPS – and could yield gamma-rays by annihilation or decay. This possibility has been explored by many authors, e.g. Calcáneo-Roldán and Moore (2000), Bergström, Edsjö and Gunnarsson (2001).

There are of course other types of dark matter. Of particular interest in the present context is baryonic dark matter in the form of cold, dense gas clouds: unlike diffuse gas, cold, dense clouds are difficult to detect directly (Pfenniger, Combes and Martinet 1994; Gerhard and Silk 1996; Combes and Pfenniger 1997; Walker and Wardle 1999). There are many astrophysical motivations for considering dark matter in this form, ranging from the properties of late-type galaxies (Pfenniger, Combes and Martinet 1994; Walker 1999), to observations of radio-wave scintillation (Henriksen and Widrow 1995; Walker and Wardle 1998), and the “Blank Field” SCUBA sources (Lawrence 2001). Indirect constraints – from the small amplitude of fluctuations in the Cosmic Microwave Background (CMB), and from Big Bang Nucleosynthesis – are generally supposed to exclude any significant amount of baryonic dark matter (e.g. Turner and Tyson 1999), but these arguments are not entirely free of loopholes (Hogan 1993; Walker and Wardle 1999) and direct constraints are desirable. There is not yet any direct observational evidence which excludes a substantial fraction of the Galactic dark matter being in the form of cold, dense gas clouds (Pfenniger, Combes and Martinet 1994; Gerhard and Silk 1996; Walker and Wardle 1999; Rafikov and Draine 2000; Ohishi, Mori and Walker 2003; Walker and Lewis 2003). But the key motivation here for considering this form of dark matter is that gas clouds are expected to emit gamma-rays if they reside in the cosmic-ray disk (e.g. Bloemen 1989), and any unexplained gamma-ray flux might therefore be a signature of dark matter in this form (De Paolis et al 1995; Kalberla, Shchekinov and Dettmar 1999; Sciama 2000). The Galactic gamma-ray halo reported by Dixon et al (1998) has, for example, been interpreted in these terms (De Paolis et al 1999). The low mass ($`10^4\mathrm{M}_{}`$) of the individual clouds renders them undetectable by EGRET when they exist in isolation (see §5). However, clustering is ubiquitous in gravitating systems, and we therefore consider the possibility that the UID EGRET sources are simply clusters of cold, dense gas clouds.

The present paper is organised as follows: we begin by describing the model we have employed for the distribution of cold gas clouds (§2), and the emission which is expected from the gas as a result of cosmic-ray interactions (§3); our model cosmic-ray distribution is given in §4; we present the results of our Monte Carlo simulation in §5, and then compare these results to the EGRET data in §6.

## 2 Dark matter model

To proceed with a calculation we need to specify a model for the clustering and distribution of the dark matter. The success of the modern structure formation paradigm, exemplified by the Cold Dark Matter (CDM) model (e.g. Blumenthal et al 1984; Davis et al 1985; Peebles 1993), argues that any acceptable model must possess clustering properties that are similar to those of CDM, in order to match the data on large-scale-structure. In CDM simulations it is found that the dark matter within individual galaxy halos is clustered into mini-halos, with roughly equal contributions to the total dark matter density being made by all mini-halo mass scales (Kauffmann, White and Guideroni 1993, Moore et al 1999; Klypin et al 1999) — i.e. $`\mathrm{d}n/\mathrm{d}MM^2`$, for mini-halos of mass $`M`$ and number density $`n`$. Following Walker, Ohishi and Mori (2003: WOM03 hereafter), we adopt this mass spectrum for our calculations, spanning the mini-halo mass range $`0.1M/\mathrm{M}_{}10^{10}`$. We further assume that all of the Galactic dark matter is in the form of cold gas clouds, clustered into mini-halos. Numerical simulations of structure formation typically show only $``$10% of the dark matter within a halo to be clustered into mini-halos (Ghigna et al 1998; Klypin et al 1999). However, these simulations do not have the resolution required to find low mass ($`M<10^8\mathrm{M}_{}`$) mini-halos, and are relevant to only the top 20% (mass fraction) of our adopted spectrum. Assuming that the true clustering spectrum does indeed extend to $`M10^8\mathrm{M}_{}`$, we thus expect that the total mass fraction in mini-halos is actually several times larger than is revealed by simulations. We have therefore adopted the extreme model in which all of the dark matter is clustered; this has the virtue of being a limiting case.

The adopted mini-halo mass spectrum is assumed to be the same throughout the Galaxy, but the normalisation (i.e. total number density of mini-halos) varies in direct proportion to the mean dark matter density of the Galaxy. Gas clouds are a collisional form of dark matter, so that an initially singular isothermal halo develops a core whose radius increases with time (Walker 1999). Simple calculations of this evolution lead to a preferred value for the column density of the individual gas clouds of $`\mathrm{\Sigma }140\mathrm{g}\mathrm{cm}^2`$, and hence a preferred model for the Galactic dark matter density distribution. We adopt this model, namely an isothermal sphere, with circular speed of $`220\mathrm{km}\mathrm{s}^1`$, and a core radius of 6.25 kpc (Walker 1999).

The density profile of the individual mini-halos themselves is also of some interest to us, in that it determines the apparent structure of the gamma-ray sources. At present there is little motivation for a detailed appraisal of these properties. The key point, however, is that the mini-halos are not expected to be point-like gamma-ray sources, and we are therefore interested in estimating their sizes. We follow WOM03 in assuming (see Ghigna et al 1998; Moore et al 1999) that the mini-halos are tidally truncated at perigalacticon, with the location of the latter estimated as one half of the Galactocentric radius of the mini-halo. This leads to gamma-ray sources with estimated angular sizes of order degrees.

## 3 Gamma-ray emissivity

For high column-density gas clouds, predicting the effects of cosmic-ray bombardment is more complicated than for diffuse gas because the cosmic-ray spectrum is substantially modified by transport through the cloud (Kalberla, Shchekinov and Dettmar 1999; Sciama 2000). Scattering and absorption of the emergent gamma-rays is similarly important. We have therefore used the Monte Carlo event simulator GEANT111http://wwwinfo.cern.ch/asd/geant to estimate the gamma-ray spectrum resulting from cosmic-ray bombardment of dense gas clouds (Ohishi, Mori and Walker 2003). By using GEANT we are able to study the transport of cosmic-rays into dense gas clouds, and their resulting gamma-ray emission, with all relevant interactions included. We employ the results of GEANT simulations for clouds of mean column-density $`100\mathrm{g}\mathrm{cm}^2`$, this being the closest, simulated case to our preferred value of $`\mathrm{\Sigma }`$ (§2). For the most part we shall need only the integrated gamma-ray emissivity of the gas clouds. In the solar neighbourhood we estimate the $`>100`$ MeV emissivity due to cosmic-ray hadrons (mainly protons) to be $`J_p=5.1\times 10^2\mathrm{ph}\mathrm{g}^1\mathrm{s}^1`$, while that due to cosmic-ray electrons is $`J_e=5.2\times 10^3\mathrm{ph}\mathrm{g}^1\mathrm{s}^1`$ (Ohishi, Mori and Walker 2003). Here we have adopted the median cosmic-ray proton spectrum of Mori (1997), and the cosmic-ray electron spectrum of Skibo and Ramaty (1993). The total emissivity is thus estimated to be $`J=J_p+J_e=5.6\times 10^2\mathrm{ph}\mathrm{g}^1\mathrm{s}^1`$.

## 4 Cosmic-rays

Most of our information on the Galactic cosmic-ray spectrum comes from direct observation of energetic particles in the solar neighbourhood, and we know very little about the spectrum at other locations in the Galaxy. Although it is possible to construct theoretical models of the cosmic-ray spectrum and distribution throughout the Galaxy, based on assumed sources, sinks and diffusive particle propagation (e.g. Porter and Protheroe 1997, Strong, Moskalenko and Reimer 2000), such models introduce an additional layer of complexity which is unwarranted in the present context. Here we adopt the simple assumption that the shape of the cosmic-ray spectrum is the same throughout the Galaxy.

It then remains to specify the cosmic-ray energy-density as a function of position in the Galaxy. Webber, Lee and Gupta (1992: WLG92) constructed numerical models of cosmic-ray propagation in the Galaxy; they did not give any analytic forms for their model cosmic-ray distributions, but an appropriate analytic approximation can be deduced from the results which they obtained. They found that the cosmic-ray radial distribution reflects, in large part, the radial dependence of cosmic-ray sources, with a modest smoothing effect introduced by diffusion. We have therefore adopted WLG92’s preferred model (their model #3) for the radial distribution of sources as our model for the radial distribution of cosmic-rays.

The various spectra of cosmic-ray isotope ratios considered by WLG92 favour models in which the diffusion boundaries are in the range $`24`$ kpc above and below the plane of the Galaxy. We adopt the midpoint of this range. WLG92 do not give a simple functional form for the vertical variation of cosmic-ray density within this zone, so we have simply assumed an exponential model: $`\mathrm{exp}(|z|/h)`$. We know that in WLG92’s models the cosmic-ray density is fixed at zero at the diffusion boundaries, and consequently the scale-height of the exponential should be approximately half of the distance to the diffusion boundary, i.e. $`h1.5`$ kpc.

These considerations lead us to the model cosmic-ray density distribution $`U(R,z)`$:

$$\frac{U}{U_{}}=\left(\frac{R}{R_0}\right)^{0.6}\mathrm{exp}[(R_oR)/\varrho |z|/h],$$
(1)
in terms of cylindrical coordinates $`(R,z)`$. Here $`R_o8.5`$ kpc is the radius of the Solar circle, while $`\varrho =7`$ kpc, and $`h=1.5`$ kpc; $`U_{}=U(R_o,0)`$ is the cosmic-ray density in the Solar neighbourhood. This distribution has the character of a disk with a central hole.

## 5 Results

Using the model described above, we simulated $`2\times 10^8`$ halos with masses $`10^2M/\mathrm{M}_{}10^{10}`$ located within 50 kpc of the Galactic Centre. The model mini-halo mass spectrum which we utilised extends to lower masses ($`0.1\mathrm{M}_{}`$); these very-mini-halos (pico-halos?) were not included in the simulation because a negligible fraction of them can be detected (see later). By the same token, the individual gas clouds themselves ($`M10^4\mathrm{M}_{}`$) are not expected to be detectable by EGRET. In particular, if all of the dark matter is assumed to be in a spherical halo of unclustered clouds, then the closest example to the Sun should lie at a distance of order 0.1 pc, and have a flux of order $`7\times 10^9\mathrm{ph}\mathrm{cm}^2\mathrm{s}^1`$ above 100 MeV; this is roughly an order of magnitude below the detection limit of EGRET.

Although the dark halo of the Galaxy extends to very large radii ($`50`$ kpc), mini-halos with Galactocentric radii greater than 50 kpc are typically undetectable because of the low cosmic-ray energy density at these radii. The total mass of dark matter simulated is approximately $`3.3\times 10^{11}\mathrm{M}_{}`$; this is only $`8/11`$ of the total amount of dark matter within the simulated volume, with the remaining $`3/11`$ being made up of clusters in the mass range $`10^1M/\mathrm{M}_{}<10^2`$.

We find that a fraction $`1.5\times 10^6`$ of the simulated mini-halos have a predicted photon flux above 100 MeV of more than $`F=7\times 10^8\mathrm{ph}\mathrm{cm}^2\mathrm{s}^1`$, thus placing them above the approximate EGRET threshold (cf. Gehrels et al 2001). In other words, the simulations yield $`280`$ sources bright enough that they should appear in the EGRET catalogue. The source counts, as a function of flux, are shown in figure 1; the relation is clearly flatter than usual Euclidean result ($`NS^{3/2}`$), and at low flux levels is well approximated by $`NS^1`$.

The locations of the synthetic sources on the sky are shown in figure 2, where it can be seen that although they populate the entire sky, they are more concentrated towards the disk of the Galaxy, and are preferentially found toward the central regions. These features are primarily a reflection of our adopted model of the distribution of cosmic-rays in the Galaxy (§4). The median Galactic latitude of the sources is $`12^{}`$. Dividing the population into two equally-sized groups of low and high latitude sources yields no significant difference in median flux between the two groups.

The mass distribution of the detectable mini-halos is shown in figure 3. At the high mass end the sources may be very luminous and can thus be detected to large distances, whereas only the closest of the low mass mini-halos are detectable. The competition between increasing mini-halo number density and declining luminosity, towards lower mini-halo masses, results in a characteristic (median) mass of $`5\times 10^6\mathrm{M}_{}`$ for the detectable sources. In the vicinity of the Sun, a mini-halo of this mass is expected to have a luminosity of approximately $`8\times 10^{34}\mathrm{erg}\mathrm{s}^1`$, in photons of energy above 100 MeV. This result is in accord with the analysis of the UID EGRET source population by Kanbach et al (1996), who concluded that their typical luminosity is $`10^{35}\mathrm{erg}\mathrm{s}^1`$. The median line-of-sight distance for the mini-halos with flux greater than the EGRET detection limit is 3.9 kpc.

Bearing in mind that both Active Galactic Nuclei (AGN) and pulsars are effectively point-like gamma-ray sources, the non-trivial angular size predicted for the gamma-ray emission from dark mini-halos is an important feature of the model. We characterise the size of the sources by $`\theta _t`$, the tidal radius of the mini-halo; this is the radius which encloses all of the flux from the model source. We compute $`\theta _t`$ using the same procedure as WOM03. For a singular isothermal sphere, the surface brightness is $`I1/\theta `$ at all radii, and we expect roughly half of the flux to be contained within $`\theta _t/2`$. The median value of the tidal radius is $`\theta _t=2.2^{}`$, for sources detectable by EGRET. The full distribution of $`\theta _t`$ is shown in figure 4.

Intrinsic source confusion, i.e. the regime where sources are so numerous that they overlap on the sky, becomes a minor problem in the region around the Galactic Centre. Within $`10^{}`$ latitude/longitude of the Galactic Centre (an area of roughly $`400\mathrm{deg}^2`$), we find 32 sources, covering about $`100\mathrm{deg}^2`$ (the sources are smaller than average). This figure is an overestimate in the sense that our model predicts too many sources — see §6. With more sensitive instrumentation, having a lower flux limit, intrinsic source confusion becomes a major issue, as the mini-halo population covers the entire sky (WOM03).

## 6 Comparison with observations

It is clear that the model we have presented predicts too many EGRET sources. The entire EGRET source list contains only about as many sources as our simulation predicts, and many ($`40`$%) of the catalogued sources are considered to be identified. Furthermore, we must expect that there will be many examples of known types of gamma-ray sources, such as pulsars, amongst the UID EGRET sources, implying that our model over-predicts the source count by a factor $`2`$. However, it must be acknowledged that the model we have presented involves substantial extrapolation into poorly charted territory and is therefore illustrative, not definitive. There are, for example, considerable uncertainties in the cosmic-ray distribution model, in the shape of the Galaxy’s dark halo, and in the mass spectrum of the clustering. In addition, we can immediately point to one aspect of the dark matter model which leads to an over-prediction of the source counts: we have assumed that all of the dark matter in the Galactic halo is bound into mini-halos (§2; WOM03), but this cannot be the case. Even if all the dark matter were initially in mini-halos, the process of tidal stripping will gradually unbind material, leading to a significant fraction in the form of dark matter “streams”. Considering these uncertainties, we recognise that the over-prediction of source counts is not a fundamental deficiency of the model.

Other basic features of the UID source population can be used to gauge the success of the model: the distribution on the sky; spectra; variability; angular structure; and observations made in other wave-bands. We now address each of these aspects in turn.

### 6.1 Sky distribution

The model correctly predicts a preponderance of sources at low- to mid-Galactic latitudes, with about the right median latitude ($`12^{}`$) for the population. However, it does not exhibit the two component (thin-disk plus thick-disk) source population which the data suggest (Gehrels et al 2001).

Figure 2 shows a relatively strong concentration towards the Galactic Centre region, compared to the data. However, in this region the predicted source density is sufficiently high that source confusion is expected to be a major problem with EGRET observations of a such a population. The predicted confusion is predominantly instrumental, because the sources are typically unresolved (see §6.4). We find approximately 100 sources within $`\pm 30^{}`$ of latitude and longitude of the Galactic Centre – an area of roughly $`3400\mathrm{deg}^2`$. By comparison the point-spread-function of EGRET has a Full-Width Half-Maximum (FWHM) of approximately $`4^{}`$ (see §6.4), implying that this entire region appears, to EGRET, to be “covered” with overlapping sources.

### 6.2 Spectra

The spectrum predicted by the present model is given in Ohishi, Mori and Walker (2003): it exhibits a peak power – i.e. a peak in $`E^2\mathrm{d}N/\mathrm{d}E`$ – at several-hundred MeV. At this point the spectrum rolls over from $`\mathrm{d}N/\mathrm{d}EE^1`$, at low energies, approaching $`E^{2.75}`$ at $`E1`$ GeV.

In comparing this prediction with the data (e.g. Merck et al 1996), there are two important points to bear in mind. First, the model gamma-ray emission spectrum is unique only by default: we don’t know the cosmic-ray spectrum elsewhere in the Galaxy. The simple fact that various gamma-ray spectra are observed should therefore not be used as an argument against the present model. We note that the diffuse Galactic plane emission at $`E1`$ GeV (Hunter et al 1997) is difficult to understand if the Galactic cosmic-ray spectrum is everywhere the same as in the solar neighbourhood, suggesting that the typical cosmic-ray spectrum may be harder than measured locally (Mori 1997).

Secondly, the sources are predicted to be extended, with low-intensity wings on the profile extending out to $`2^{}`$, typically. Although the spectrum should be uniform across the source, the fact that the point-spread function of the detector changes with energy, coupled with the low-level “wings” on the source profile, could lead to spurious estimates of spectral shapes. In particular, some fraction of the source flux will be absorbed into the estimate of the background intensity, and this fraction will vary with photon energy. These problems, coupled with the low signal-to-noise ratio of many of the UID sources, make it difficult to assess the success of the model spectral predictions.

### 6.3 Variability

The model we have presented involves emission which is intrinsically steady on observationally accessible time-scales; it is therefore not relevant to any sources which are known to vary significantly. Most of the UID EGRET sources are not bright enough to permit strong constraints on their variability. There is no consensus in the literature regarding the variability of the UID sources: McLaughlin et al (1996) find that only a small fraction, roughly one in six, of the UID EGRET sources are significantly variable (see also Wallace et al 2000); by contrast Torres et al (2001) suggest that the fraction may be as large as one third, in the case of low-latitude UID sources (see also Torres, Pessah and Romero, 2001, and Tompkins 1999). We note that large-amplitude variations do not sit easily with the multiple/extended source designation carried by $``$50% of the UID EGRET sources (Hartman et al 1999).

### 6.4 Angular Structure

Although the estimated angular sizes of the detectable mini-halos are large (of order degrees), the resolving power of EGRET is quite modest, with a point-spread-function of $`5^{}`$ Full-Width-Half-Maximum (FWHM) at 100 MeV (Thompson et al 1993). (The in-flight calibration data are consistent with the pre-launch calibration in respect of the point-spread-function — Esposito et al 1999.) The EGRET angular resolution improves with increasing energy, scaling roughly as $`E^{0.534}`$. However, if $`\mathrm{d}N/\mathrm{d}E`$ is close to $`E^2`$, half of the photons contributing to a source detection are within a factor $`2`$ of the low-energy threshold. We therefore adopt a FWHM of $`4^{}`$ as the relevant instrumental width; sources would thus need to be at least $`8^{}`$ across in order to be fully resolved. For an isothermal density distribution within each mini-halo (§5), half of the total flux is contained within a diameter approximately equal to $`\theta _t`$. Referring to figure 4 we then find that only a tiny fraction of the synthetic population could be fully resolved by EGRET.

If the instrumental FWHM is comparable to the source size then the source structure will not be resolved, but the data can nevertheless indicate that the source is extended, by virtue of the observed intensity profile being broader than the point spread function. Most of the synthetic sources fall into this category. It is notable that half of the UID EGRET sources are recorded as extended/multiple by Hartman et al (1999). (The two possibilities cannot be differentiated if the source structure is not fully resolved.)

### 6.5 Counterparts at other wavelengths

The basic criterion for an EGRET source to be classified as “UID” is that it should not have an obvious counterpart at other wavelengths. The population we have modelled clearly meets this requirement because the emission comes from “dark” matter. The detectability of cold, dense gas has been discussed by a number of authors: Pfenniger, Combes and Martinet (1993); Gerhard and Silk (1996); Combes and Pfenniger (1997); Walker and Wardle (1999); and Sciama (2000). Perhaps the simplest and most robust of expectations is that there will be thermal emission from the clouds, implying bright microwave sources coincident with the gamma-ray sources. An important point to note is that gamma-ray and microwave luminosities are both proportional to mini-halo mass and cosmic-ray density; consequently the microwave flux can be estimated simply by scaling the observed EGRET flux (Wardle and Walker 1999; Sciama 2000).

To determine the bolometric microwave flux we need only take the ratio of the thermal emissivity due to cosmic-rays ($`\mathrm{\Gamma }_010^4\mathrm{erg}\mathrm{s}^1\mathrm{g}^1`$; WOM03), to the gamma-ray emissivity ($`J=J_p+J_e5.6\times 10^2\mathrm{ph}\mathrm{s}^1\mathrm{g}^1`$ above 100 MeV; §3), yielding a microwave flux of $`S1.8\times 10^{10}F_7\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1`$. Here the gamma-ray flux above 100 MeV is $`10^7F_7\mathrm{ph}\mathrm{cm}^2\mathrm{s}^1`$. The atmospheric temperature is estimated to be in the range $`4.24.9`$ K, if the individual clouds have mass $`M10^410^5\mathrm{M}_{}`$ (WOM03). The principal remaining uncertainty then lies with the nature of the emitted spectrum: we have previously suggested a blackbody spectrum (WOM03), whereas Lawrence (2001) concludes that this is inconsistent with the observed spectrum of the Blank Field SCUBA sources (which objects he interprets in terms of cold, planetary-mass gas clouds). More generally we can take $`S_\nu \nu ^\alpha B_\nu `$, with $`\alpha =0`$ for blackbody emission, and $`\alpha =2`$ for emission from small particles whose absorption increases in proportion to $`\nu `$. Lawrence (2001) concludes that the latter spectrum is consistent with the data on the Blank Field SCUBA sources, providing $`4.7<T(\mathrm{K})<6.4`$. Here we restrict attention to two cases which adequately represent the plausible range of thermal emission spectra for the dense gas: a $`T=4.2`$ K blackbody, and a $`T=4.9`$ K “dusty” ($`\alpha =2`$) spectrum.

Our adopted spectra differ greatly in the expected flux at radio frequencies, where we are far below the peak thermal emission — see figure 5. At low frequencies the two spectral models can be approximated by $`S_\nu 4.2\nu ^2F_7`$ mJy (blackbody) and $`S_\nu 1.4\times 10^5\nu ^4F_7`$ mJy (dust), with $`\nu `$ in GHz. However, even in the blackbody case, where the radio emission is relatively strong, the predicted flux would be difficult to detect because it is so extended. Taking the bright UID EGRET source 3EG J1835+5918 as an example (Hartman et al 1999; Mirabel et al 2000; Reimer et al 2001), with $`F_77`$, we find a predicted flux at 1.4 GHz of roughly 60 mJy in the blackbody case. This estimate is well above the point source detection limit of the observations, reported by Mirabel et al (2000), of 2.5 mJy. However, a mini-halo is certainly not point-like. In fact it would fill the primary beam of the radio telescope, and would be resolved-out on all but the shortest interferometeric baselines. The very extended sources predicted by the present model are expected to be difficult to detect using radio interferometers configured for high resolution imaging.

## 7 Discussion

Although the predicted microwave counterparts to UID EGRET sources are not expected to have been detected in the counterpart searches to date, they could be revealed in the near future by some of the various experiments designed to study anisotropies in the Cosmic Microwave Background (CMB). In particular, the MAP satellite222http://map.gsfc.nasa.gov, has now completed a sensitive all-sky survey, at several frequencies, and the data are about to be released. Should MAP have detected the predicted microwave counterparts? The answer to this question depends on the spectrum of the emission, because all of the MAP frequencies are well below the thermal peak. Even considering the highest frequency channel (90 GHz), the two spectral models differ by a large factor in their predicted flux. If the emission is blackbody, then $`S_\nu 20F_7`$ Jy, but only $`S_\nu 0.6F_7`$ Jy for the “dust” model (see figure 5). The point-source (i.e. single pixel) detection limit for MAP at this frequency should be approximately 1.7 Jy (5$`\sigma `$, and we have taken the CMB contribution to be roughly equal to the thermal noise of the instrument). On this basis we do not expect the typical UID sources to be detected by MAP if they have a “dusty” spectrum.

Even if the mini-halos have a blackbody spectrum they might not be detected by MAP, because its $`0.3^{}`$ pixels are substantially smaller than the predicted mini-halo sizes. For an isothermal mini-halo density profile, the enclosed flux varies roughly in proportion to radius. With a predicted median source radius of $`2.2^{}`$, it is evident that the largest single-pixel flux expected from a typical mini-halo is not much above the MAP detection limit. Computing the peak single-pixel flux for each source in our simulation, we expect that only 40% of the predicted population ought to be detectable by MAP even if the spectra are blackbody. However, we note that sources which are individually undetected may still be useful for constraining the typical microwave/gamma-ray flux ratio of the UID EGRET population.

Other satellite and balloon-borne CMB experiments will provide robust constraints on the model we have presented, by virtue of observing close to the predicted thermal peak of the mini-halos. In particular, the High Frequency Instrument on the Planck satellite333http://sci.esa.int/home/planck and its prototype, the balloon-borne Archeops (Benoît et al 2002), will map the sky at 353 and 545 GHz, sandwiching the point – around 400 GHz (see figure 5) – where our model spectra cross. At these frequencies the predicted flux of a typical UID EGRET source ($`F_7=1`$) is of order 30 Jy, for either spectral model. By contrast, the limiting flux ($`5\sigma `$) for Planck, at 353 GHz, will be 100 mJy in a single (5 arcmin) pixel, and roughly 5 Jy for a source of $`2.2^{}`$ radius. (Oberving at these high frequencies has the additional advantage of less confusion from the degree-scale CMB anisotropies.) Planck should therefore detect microwave counterparts to any of the UID EGRET sources in which the gamma-rays arise from cosmic-ray interactions with dense gas. We note that a substantial fraction of the sky has already been mapped by the Archeops experiment (Benoît et al 2002), and UID EGRET source counterparts might be present in the existing data. Counterparts are best searched for, initially, well away from the Galactic plane as the latter region is likely to be confused at high frequencies. Scaling from the Planck sensitivity estimates, using the estimated duration and coverage – 2 year mission with full sky coverage, for Planck, versus 24 hour Archeops dataset covering 25% of the sky – we estimate that the flux limit of Archeops, for the extended microwave sources we predict, should be approximately 70 Jy. Hence UID EGRET sources brighter than $`3\times 10^7\mathrm{ph}\mathrm{cm}^2\mathrm{s}^1`$ ($`>\mathrm{\hspace{0.17em}100}`$ MeV) may have counterpart microwave sources found in the Archeops data.

We note that the balloon-borne MAXIMA and BOOMERANG experiments also covered frequencies close to the thermal peak of cold gas emission, but these experiments covered only $``$1% of the sky (Hanany et al 2000; Coble et al 2003), and these data are therefore of limited utility in the present context. Ground-based studies are similarly limited to small patches of the sky, when observing at 400 GHz, but would nevertheless be helpful if the UID EGRET source population is specifically targetted for observations.

If microwave counterparts are discovered, then we will be able to study the density profiles of dark matter mini-halos directly. Further to the properties noted in §2 for the individual source profiles, we can make some generic predictions for their structure: (i) the intensity should rise to a high central peak, but there should be a core (Walker 1999) in the surface-brightness profile; (ii) the limb of the source should exhibit a sharp cut-off due to tidal truncation; and (iii) there may be tidal streams extending along the mini-halo’s orbit.

In addition to tidal streams associated with identifiable mini-halos, it is expected that some mini-halos have been completely disrupted by the tidal fields they have experienced. In these cases the tidal streams are not associated with a bound cluster, and thus represent a distinct category of microwave/gamma-ray source predicted by the model. Shells with sharply-defined edges could also appear in the microwave/gamma-ray maps, as these are a common feature of tidal debris (Hernquist and Quinn 1988; Hernquist and Spergel 1992).

The Gamma-ray Large Area Space Telescope444http://glast.gsfc.nasa.gov (GLAST) will provide a significant advance in our understanding of the UID sources, because GLAST will be much more sensitive than EGRET and will have better angular resolution. These improved capabilities will permit powerful tests of the model we have presented. However, the data from MAP and the balloon-borne CMB experiments will be available on a much shorter time-scale than those from GLAST, and it should be possible to make considerable progress with the microwave data alone. In particular, our prediction of bright, extended, thermal microwave counterparts is unique amongst existing models of the UID EGRET population.

Finally we note that some of the UID EGRET sources are sufficiently bright that they may be detectable by ground-based TeV telescopes (Aharonian et al 1997), even if their spectra are as steep as $`E^{2.75}`$ between the GeV and TeV bands. Studying UID sources at TeV energies may permit their nature to be discerned. For example, the extended nature of the sources predicted by the present model is unusual for gamma-ray sources, and the slightly extended (6 arcminute) TeV source, in the vicinity of 3EG J2033+4118, reported by Aharonian et al (2002) is of interest in this context.

## 8 Conclusions

If our Galaxy contains a significant component of dark matter in the form of cold, dense gas clouds, clustered into large aggregates, some of those clusters should have been detected by EGRET. Using a CDM-like mass spectrum for the clustering, we have shown that the predicted gamma-ray source population has properties which are broadly similar to those of a large fraction of the UID EGRET sources. In particular: the Galactic latitude distribution and the source size distribution anticipated in the model find support in the EGRET data. Furthermore the intrinsically “dark” nature of the predicted sources naturally explains most of the difficulty in finding counterparts; however, the large angular size of the clusters also plays a role, because counterpart searches to date have had poor sensitivity to $``$degree-sized sources. The total number of UID sources predicted by the model is too large, particularly bearing in mind that many of the observed UID sources are likely to be examples of known types of gamma-ray emitters. Data returned by the various CMB experiments will provide a powerful test of the interpretation we have presented: counterpart thermal microwaves should be detected from the cold gas, and the source structure should be resolved.

MAW thanks Gordon Garmire and Eric Feigelson for emphasising the importance of the gamma-ray data, and for helpful discussions on the relevant physics. We thank Ben Moore for advice on modelling the mini-halo population, and Simon Johnston for helpful comments on the manuscript. The Referee’s comments also improved the paper.


text:
# 1 . Introduction

## 1 . Introduction

The Magellanic Clouds provide an ideal opportunity to study in detail the structure and evolution of stars. Rich populations of stars approximately at the same, relatively small distance with small interstellar reddening make the Large and Small Magellanic Clouds very important targets for observing surveys.

RR Lyr stars were first discovered in the SMC near the cluster NGC121 by Thackeray (1951). The stars were more than a magnitude fainter than expected, what was a crucial confirmation of the major revision of the extragalactic distance scale proposed by Baade (1952). The first systematic search for field RR Lyr variables in the SMC was conducted by Graham (1975). During that survey 76 RR Lyr stars were discovered in a $`1^{}\times 1\text{.}^{}3`$ outlying field centered on the cluster NGC121. Smith et al. (1992) used B-band photographic photometry of a field in the northeast arm of the SMC to identify additional 22 probable RR Lyr variables. The same outlying region of the SMC was observed by Sharpee et al. (2002). They presented V-band and B-band CCD photometry of a few RR Lyr stars.

Walker (1989) surveyed five SMC clusters for RR Lyr stars. Four already known RR Lyr variables in the NGC121 were rediscovered, but no such stars were found in the other clusters. Because the age of NGC121 was estimated to be $`12\pm 2`$ Gyr (Stryker et al. 1985), and Lindsay 1, the next oldest SMC cluster, $`10\pm 2`$ Gyr (Olszewski et al. 1987), Walker concluded that the minimum age of RR Lyr stars is about 11 Gyr.

In 1990s the number of known variable stars in the Magellanic Clouds dramatically increased, when the large microlensing searches began regular photometric monitoring of both galaxies (e.g., MACHO – Alcock et al. 1993, EROS – Aubourg et al. 1993). Natural by-product of the microlensing surveys are huge databases with precise photometry of millions of stars.

The SMC was also included to the list of targets of the second phase of the Optical Gravitational Lensing Experiment (OGLE-II; Udalski, Kubiak and Szymański 1997). About 2.4 square degrees of central part of the SMC were observed each night during the observing seasons 1997–2000. Photometry was obtained with the $`BVI`$ filters, closely resembling the standard system.

The OGLE-II survey has yielded a particularly rich harvest of variable stars from the SMC. In the previous papers we presented the catalog of eclipsing binary stars (Udalski et al. 1998b), the catalog of Cepheids from the SMC (Udalski et al. 1999), and general catalog of variable stars detected in the Magellanic Clouds (Żebruń et al. 2001). In addition BVI maps of the SMC were released providing precise photometry and astrometry of about 2.2 million stars (Udalski et al. 1998a).

In this paper we present a sample of 571 RR Lyr stars and several other pulsating objects, likely $`\delta `$ Sct stars, detected in the OGLE-II fields in the SMC. The stars were selected from the reprocessed OGLE-II photometry based on the Difference Image Analysis (DIA) technique – Woźniak’s (2000) implementation of Alard and Lupton (1998) and Alard (2000) optimal Point Spread Function matching algorithm.

Similarly to the previous catalogs, all data presented in this paper, including individual observations, are available to the astronomical community from the OGLE Internet archive.

## 2 . Observations and Data Reductions

Observations presented in this paper were collected during the second phase of the OGLE microlensing search with the 1.3-m Warsaw telescope at Las Campanas Observatory, Chile. The observatory is operated by the Carnegie Institution of Washington. The telescope was equipped with the “first generation” camera with a SITe $`2048\times 2048`$ CCD detector working in drift-scan mode. The pixel size was 24 $`\mu `$m giving the 0.417 arcsec/pixel scale. Observations of the SMC were performed in the “slow” reading mode of the CCD detector with the gain 3.8 e-/ADU and readout noise of about 5.4 e-. Details of the instrumentation setup can be found in Udalski, Kubiak and Szymański (1997).

Observations of the SMC were collected between June 26, 1997 and November 25, 2000. Eleven driftscan fields (SMC$`\mathrm{\_}`$SC1–SMC$`\mathrm{\_}`$SC11) covering about 2.4 square degrees of central parts of the SMC were observed. The majority of frames were taken in the $`I`$ photometric band (about 280–340 epochs depending on the field). Other images were collected through the V-band (typically about 30 epochs) and B-band (about 20 epochs) filters. The effective exposure time lasted 125, 174 and 237 seconds for the $`I`$, $`V`$ and B-band, respectively. The median seeing was about 1.′′3 for our dataset.

The I-band photometry was obtained using Difference Images Analysis (DIA) – image subtraction algorithm developed by Alard and Lupton (1998) and Alard (2000), and implemented by Woźniak (2000). We introduced several modifications compared to the DIA techniques employed in a catalog of variable stars in the Magellanic Clouds (Żebruń et al. 2001). For instance, we performed the DIA photometry for each star found in the reference image instead of the variable objects only. Further details of the DIA analysis and calibration of photometry may be found in the papers by Woźniak (2000) and by Żebruń, Soszyński and Woźniak (2001).

The frames in the $`V`$ and $`B`$ bands were analyzed using the DoPhot photometry program (Schechter, Saha and Mateo 1993). Transformation of the instrumental photometry to the standard system is described by Udalski et al. (1998a).

Equatorial coordinates of all stars were calculated in the identical manner as described in Udalski et al. (1998a). The internal accuracy of the determined equatorial coordinates, as measured in the overlapping regions of neighboring fields, is about 0.′′15–0.′′20 with possible systematic errors of the DSS coordinate system up to 0.′′7.

Fig. 1 presents the picture of the SMC from Digitized Sky Survey CD-ROMs with contours of the OGLE-II fields. Positions of RR Lyr stars are marked with black dots. One can easily notice that the surface density of RR Lyr variables is very uniform within the galaxy. This confirms the results of Graham (1975) and Smith et al. (1992) who noticed that RR Lyr stars are not strongly concentrated toward either the bar or the center of the SMC.

## 3 . Interstellar Reddening

Determination of the interstellar reddening to the SMC fields was performed by Udalski et al. (1999). In short, they used red clump stars for mapping the fluctuations of mean reddening in OGLE fields, treating their mean I-band magnitude as the reference brightness. Since reddening in the SMC is smaller and more homogeneous than in the LMC, it was determined only in 11 lines-of-sight – one per OGLE field. The zero point of the reddening map was derived on the basis of previous determinations toward two stars clusters: NGC416 and NGC330.

The final $`E(BV)`$ reddening in the SMC is listed in Table 1. The error of the map is equal to $`\pm 0.02`$ mag. Similar values of reddening were estimated by Sharpee et al. (2002) for a field in the northeast arm of the SMC located close to the OGLE fields SMC$`\mathrm{\_}`$SC10 and SMC$`\mathrm{\_}`$SC11. Interstellar extinction in the BVI bands can be calculated using the standard extinction curve coefficients (e.g., Schlegel et al. 1998):

$$A_B=4.32E(BV),A_V=3.24E(BV),A_I=1.96E(BV)$$

## 4 . Selection of RR Lyr Stars

Selection of RR Lyr stars was performed in two stages. First, preliminary search for variable stars was performed using the regular OGLE-II PSF (DoPhot) photometry. The mean I-band magnitude of objects was limited to $`I<20`$ mag. The minimal number of individual measurements was set to 50. Candidates for variable stars were selected based on comparison of the standard deviation of photometry with typical standard deviation for stars with similar brightness. Light curves of selected candidates were searched for periodicity using the AoV algorithm (Schwarzenberg-Czerny 1989). Light curves of all objects revealing statistically significant periodic signal were then visually inspected. Candidates for RR Lyr stars were extracted on the basis of light curve shapes and magnitudes of the stars.

The second, final search for RR Lyr variables was performed with the DIA photometry of all stars detected in the reference images. The quality of DIA photometry is improved by a factor of at least two, compared to the DoPhot photometry (Woźniak 2000), what enabled us to increase considerably the completeness of the catalog.

We examined with the AoV period search technique the DIA photometry of objects with the mean magnitude $`18.4<I<19.4`$ and standard deviation at least 0.02 mag larger than typical standard deviation for non-variable stars of similar brightness. Additionally we checked stars with the V-band magnitude between 18.9 and 20.0 and the standard deviation of V-band photometry 0.05 mag larger than typical. New candidates for RR Lyr stars were selected from a sample of stars with periods smaller than 1 day based on visual inspection of their light curves and location in the color-magnitude diagram (CMD).

The second stage of variability search increased the number of candidates for RR Lyr stars by about 30%. In total 571 candidates for RR Lyr stars were identified. Several objects with colors outside the range $`0.2<VI<0.8`$ and objects with no color information but with evident RR Lyr-type light curves were also included to this sample. They can be highly reddened or blended stars.

## 5 . Classification

We divided all objects into four groups: fundamental mode RR Lyr stars (RRab), first overtone (RRc), double mode RR Lyr stars (RRd) and other variable stars, which, in the majority of cases, are probably $`\delta `$ Sct stars.

Selection of RRab-type stars from our sample was not difficult, because these variables form well separated groups in various diagrams. We decided to use period–amplitude diagram to classify fundamental-mode pulsators, because periods and amplitudes are the observables that are measured with the highest precision.

There were more problems with separation of the first overtone RR Lyr stars and short-period variable stars. Both groups of objects in our sample have similar luminosities, amplitudes and shapes of light curve. Additionally, the range of periods of $`\delta `$ Sct stars overlaps with the range of periods of RRc variables.

The models of Bono et al. (1997) predict that the first-overtone RR Lyr stars with the shortest periods should have the lowest amplitudes. Indeed, in the period–amplitude diagrams (Fig. 4) one can notice a minimum of IV-band amplitudes near the period of about 0.26 days. Therefore, we decided to mark as the first overtone RR Lyr candidates with periods larger than 0.26 days. Objects with shorter periods were classified as $`\delta `$ Sct stars. However, we stress that it is possible that some of the RRc variables on our list are $`\delta `$ Sct or other types of variable stars, and vice versa. It is also possible that some second-harmonic RR Lyr stars (RRe) can be included in the short period group. Further observations, especially in the Strömgren uvby filters, could provide additional information for final classification.

We used two methods to search for multi-periodic variable stars. First, we fitted a fourth order Fourier series to each folded light curve from our sample and subtracted fitted function from the observational data. Then, the residuals were searched for periodic signal and, if detected, such a candidate was marked for visual inspection.

Second search for double-mode RR Lyr stars was performed using the Clean algorithm of period determination (Roberts, Lehár and Dreher 1987). All RR Lyr candidates from our list were subject to the Clean period analysis. In further analysis we selected objects for which the ratio of the highest peak in the power spectrum and one of the next four strongest peaks was close to 0.745.

Final list of the double-mode variable stars presented in this paper was obtained after careful visual inspection of power spectra and folded light curves of candidates. We identified 57 RRd stars, 2 double mode $`\delta `$ Sct stars and several dozen stars with two closely spaced frequencies, that are probably members of a new class of multi-periodic RR Lyr stars described by Olech et al. (1999). About 10% in both, RRab and RRc groups, exhibited secondary periodicity very close to the primary pulsation frequency, with period ratios in the range 0.98–1.02. This frequency pattern cannot be explained by a superposition of radial pulsations and is therefore believed to be related to non-radial modes. Non-radial oscillations have also been detected in some RR Lyr stars in the Galactic bulge by Moskalik (2000) and in the LMC by Kovács et al. (2000). Theoretical models of these stars have been proposed by Dziembowski and Cassisi (1999).

## 6 . Catalog of RR Lyr Stars

### 6.1 Single-Mode RR Lyr Stars

458 RRab and 56 RRc variable stars passed our selection criteria. They are listed in Tables 2 and 3. First two columns of both tables contain the star identification: field\_name star\_ID. The star\_IDs are simultaneously equatorial coordinates, RA and DEC (J2000), of the objects. In the next columns period in days and moment of the zero phase corresponding to maximum light are given. Finally, last columns show intensity mean IVB photometry.

More parameters (period errors, Fourier parameters of the light curve decomposition, IVB-band amplitudes) are available in the electronic form from the OGLE Internet archive:

http://www.astrouw.edu.pl/~ogle/

ftp://sirius.astrouw.edu.pl/ogle/ogle2/var\_stars/smc/rrlyr/

or its US mirror

http://bulge.princeton.edu/~ogle/
ftp://bulge.princeton.edu/ogle/ogle2/var\_stars/smc/rrlyr/

Also individual BVI observations of all objects and finding charts are included.

Tables 2 and 3 contain together 536 entries but only 514 objects, because 22 stars were detected twice – in the overlapping regions of adjacent fields. We decided not to remove twice-detected RR Lyr stars from the final list, because their measurements are independent in both fields and can be used for testing quality of data and completeness of the sample. In Table 6 we provide cross-reference list to identify stars in the overlapping regions.

In Fig. 2 we present I-band light curves of a few typical RRab, RRc stars and one RRd star. All the diagrams have the same magnitude range to compare the amplitudes and brightness of the stars. The light curves are arranged according to the periods.

### 6.2 Double-Mode RR Lyr Stars

RRd stars discovered in the SMC are listed in Table 4. 59 such objects were detected but 57 of them are unique. Two stars are located in the overlapping regions between fields and they were discovered independently. Table 4 lists the field name, star ID, period of the first overtone, period of the fundamental mode, period ratio, and finally $`IVB`$ intensity-mean magnitudes. Other parameters, individual BVI measurements and finding charts are available from the OGLE Internet archive.

### 6.3 Other Stars

Table 5 presents a list of 19 detected short-period variables and a few other variable stars. The table contains 20 entries because one object was identified twice in the overlapping region of fields SMC\_SC5 and SMC\_SC6. Consecutive columns of Table 5 represent the same data as in Tables 2 and 3. Most stars from the list are probably $`\delta `$ Sct stars. Two stars: OGLE005542.36–732105.3 and OGLE010118.99–723323.8, are double-mode pulsators, with the period ratios 0.76259 and 0.77098, respectively. In Table 5 we present longer periods of both stars.

## 7 . Basic Parameters of RR Lyr Stars

### 7.1 Period Distribution

Distribution of periods of RR Lyr stars is an indicator of the morphology of the horizontal branch of the population of old, metal-weak stars. The final period search was carried out with program Tatry by Schwarzenberg-Czerny (2002, private communication). The program works in two steps. First, it calculates a periodogram and identifies its peaks. Then the corresponding peak frequencies are examined in detail to find the best one and its corresponding ephemeris. Tatry provides also reliable estimation of period errors, what we could check comparing periods of stars located in the overlapping regions of the adjacent fields. Accuracy of periods for the vast majority of the RR Lyr stars is better than $`10^5`$ days. The period errors for each star are included in the OGLE Internet archive data.

Period distribution for the global sample is shown in Fig. 3. The mean period of ab-type RR Lyr stars, and simultaneously the most preferred period, is $`P_{ab}=0.589`$ days. This result is in the ideal agreement with Graham’s (1975) determination of the mean period of RRab variables near NGC121. $`P_{ab}`$ is believed to be an indicator of the mean metallicity – it increases with decreasing metallicity. However, one should remember that Clement et al. (2001), in a study of the properties of the RR Lyr stars in Galactic globular clusters, suggested that the correlation breaks down when the mean period is larger than 0.60 days.

The period distribution of the first overtone pulsators indicates two preferred periods. Most of the RRc stars prefer periods near $`P_c=0.357`$ days, but there is also a significant excess of stars around $`P=0.29`$ days. Such short-period RR Lyr stars also exist in other environments, e.g., in the LMC (Alcock et al. 1996) and in the globular cluster IC 4499 (Walker and Nemec 1996). It has been suggested (van Albada and Baker 1973) that these stars are second overtone (RRe) pulsators. The pulsation of RR Lyr stars in overtone modes was theoretically analyzed by Stellingwerf et al. (1987).

RRd variable stars (first-overtone periods) and RRc stars occupy different, but overlapping regions in the period distribution histogram. Periods of the first overtone in the double-mode RR Lyr stars range from about 0.354 to 0.435 days with the mean period $`P_d=0.382`$ days. Single-mode first-overtone pulsators have typically shorter or longer periods than double-mode RR Lyr stars (we found three candidates for RRc stars with periods between 0.428 and 0.446 days).

### 7.2 Intensity Mean Photometry

BVI intensity mean photometry of each object from our sample of RR Lyr candidates was derived by integrating the light curve converted to intensity units. The light curve was approximated by the Fourier series of fifth order and results were converted back to the magnitude scale. Accuracy of the mean I-band photometry is about 0.02 mag while that of the V-band and B-band about 0.03 mag and 0.08 mag, respectively, what is a consequence of smaller number of observations in the BV-band filters.

To obtain the mean brightness of the SMC RR Lyr stars, we prepared histograms of the BVI-band apparent magnitudes separately for RRab, RRc and RRd variable stars. In the next step we fitted a Gaussian to each histogram, and obtained the most preferred brightness of these stars. Then we repeated this procedure using magnitudes corrected for interstellar extinction. Results are presented in Tables 7 and 8.

### 7.3 Period-Amplitude Relation for RR Lyr Stars

We determined amplitudes of light curves by fitting Fourier series of fifth order and calculating difference between minimal and maximal values of the function. The period-amplitude diagram for the SMC RR Lyr stars is shown in Fig. 4. The symbols are as follows: black dots indicate the position of RRab and RRc stars, triangles and diamonds represent RRd stars (first-overtone and fundamental mode, respectively) and crosses mark other stars.

Well-known dependence of the amplitudes on periods for RRab stars is seen, although some stars, especially RRab with shorter periods, do not follow this relation. The first overtone RR Lyr variables show weak correlation of amplitudes and periods but, inversely to the fundamental mode pulsators, amplitude slightly increases with increasing periods.

Comparison of amplitudes of fundamental mode and first overtone of RRd stars confirms that the dominant mode of pulsation in the majority of double-mode RR Lyr stars is the first overtone. The mean ratio of the first-harmonic and fundamental mode amplitudes is about 1.6.

### 7.4 Fourier Parameters of Light Curve Decomposition

Fourier decomposition of light curves of RR Lyr stars has been widely used for analyzing their properties. Several empirical and theoretical relations between Fourier parameters and the physical parameters of RR Lyr variables have been proposed (e.g., Kovács and Walker 2001).

Fifth order Fourier series were fitted to the magnitude scale I-band light curve. Then, we derived $`R_{21}`$, $`R_{31}`$, $`\varphi _{21}`$ and $`\varphi _{31}`$ Fourier parameters, where $`R_{ij}=A_i/A_j`$, $`\varphi _{ij}=\varphi _ii\varphi _j`$. $`A_i`$ and $`\varphi _i`$ are the amplitudes and phases of ($`i1`$) harmonic of the Fourier decomposition of light curve.

Light curves of double-mode pulsators were decomposed to the sum of two Fourier series of fourth order corresponding to both periodicities. Then we calculated Fourier parameters for both pulsating modes.

In Fig. 5 we present $`R_{21}\mathrm{log}P`$ and $`\varphi _{21}\mathrm{log}P`$ diagrams constructed for our sample of RR Lyr stars. The first-overtone and fundamental mode pulsators are well separated in both diagrams. $`R_{21}`$ for RRab stars tends to become smaller as the periods become longer, what means that the light curves of long-period ab-type RR Lyr are more sinusoidal than the light curves of short-period RRab stars.

Fourier parameters of the first-overtone mode of RRd stars fall in the sequence of the single-mode first overtone RR Lyr stars, but the fundamental mode pulsation have $`R_{21}`$ smaller than corresponding single-mode RR Lyr variables. This means that fundamental-mode pulsations in RRd stars have not only smaller amplitude, but the light curves are more sinusoidal than the light curves of RRab stars.

### 7.5 Petersen Diagram

RR Lyr stars pulsating simultaneously in the fundamental mode and the first overtone are especially interesting as a test of structural and evolutionary models of the horizontal branch stars. A powerful tool for diagnosing stellar models and for determination of masses, absolute magnitudes and metallicities of stars is the so called Petersen diagram, in which the period ratio of two excited modes is plotted against the longer period. Position of stars in the Peterson diagram can be precisely determined, because periods are the observables that are measured with the highest precision.

In Fig. 6 we present Petersen diagram of all discovered RRd stars. For comparison, we also plotted in the same diagram double mode RR Lyr stars from the LMC (catalog of RR Lyr stars from the LMC will be published in the forthcoming paper).

The SMC and LMC double-mode RR Lyr stars evidently form the same relation $`P_1/P_0`$ vs. $`P_0`$, but RRd stars from the SMC prefer longer periods and larger period ratios than the LMC RRd stars. Pulsation models suggest that the position of a given RRd star in the Petersen diagram is a function of stellar mass and metallicity. We can expect that the shift in the Petersen diagram between RRd stars from the LMC and the SMC is caused by different mean metallicity of RR Lyr stars in both galaxies. More information about analysis of the Peterson diagrams can be found in Popielski, Dziembowski and Cassisi (2000).

## 8 . Completeness of the Catalog

We estimated completeness of the Catalog of RR Lyr stars in similar manner as for the catalog of Cepheids from the Small Magellanic Cloud (Udalski et al. 1999), i.e., by comparison of objects located in the overlapping regions of the adjacent fields. Ten such regions exist between our fields (Fig. 1) allowing to perform 20 tests of paring objects from a given and neighboring fields. In total, 57 stars from Tables 2–5 should be theoretically paired with counterparts in the overlapping fields. We found counterparts in 50 cases. Counterparts of four from seven unpaired objects had the number of observations smaller than 50 and thus they were not searched for variability. Smaller number of measurement points in the edge regions of the fields is due to non-perfect pointing of the telescope. The remaining three objects were missed because of severe blending with another stars or because of unusually small amplitude of variability.

Udalski et al. (1998a) estimated completeness of detection of stars in the SMC OGLE fields using artificial star tests. For stars as bright as RR Lyr the completeness was found to be between 90 and 95%, depending on the field density. In our sample completeness should be much higher, because the stars are being detected on the DIA reference images, obtained by co-adding 20 best frames for each field. We conclude that the completeness of our sample should be about 90%.

## 9 . Conclusions

The long term photometric surveys of dense stellar systems may provide large samples of all types of variable stars (Paczyński 2000). In this paper we present the catalog of RR Lyr variable stars – objects playing a major role in studies of the processes of stellar pulsation, post-main-sequence evolution of stars, calibration of extragalactic distances, and structure and age of galaxies. This is the largest sample of RR Lyr variables detected in the SMC so far.

We reported detection of 57 double-mode RR Lyr stars. We also found RR Lyr variables exhibiting two closely spaced frequencies, most probably related to non-radial pulsations.

About 10% of the detected RR Lyr variables are RRc stars. Next 10% of the

whole sample are RRd variables. Relatively small number of the first-overtone pulsators and the value of the average period of RRab stars (0.589 days) places the SMC RR Lyr stars near the Oosterhoff type I group.

Additionally, we conducted a search of RR Lyr stars in the SMC clusters. Discovery of RR Lyr variables in the clusters younger than NGC121 (hosting four RR Lyr stars), would have a significant consequences on our understanding of evolution of old, solar-mass stars. We used the OGLE catalog of the SMC clusters (Pietrzyński et al. 1998) containing 238 clusters, their coordinates and angular sizes. We scanned our sample of RR Lyr variables looking for stars in a distance smaller than 1.5 cluster radius from the cluster center. We obtained the list of 21 RR Lyr stars located close to the line-of-sight to, in total, 19 SMC clusters. Only in two cases (OGLE-CL-SMC0048 and OGLE-CL-SMC0052) we detected two RR Lyr stars in the cluster neighborhood. It is very probable that the majority of these cases are just optical coincidences, because more or less similar number of stars could be expected in the background of SMC clusters (comparing the area of the clusters to the area of all our fields). Nevertheless, because of potential importance of the discovery of RR Lyr variables in the SMC clusters, we provide in Table 9 the list of stars located close to the SMC clusters.

Acknowledgements. We would like to thank Prof. Bohdan Paczyński for many discussions and valuable suggestions. We are very grateful to Prof. Wojciech Dziembowski for helping us in the classification of stars. The paper was partly supported by the Polish KBN grant BST to Warsaw University Observatory. Partial support for the OGLE project was provided with the NSF grants AST-9820314 and AST-0204908 and NASA grant NAG5-12212 to B. Paczyński. We acknowledge usage of the Digitized Sky Survey which was produced at the Space Telescope Science Institute based on photographic data obtained using the UK Schmidt Telescope, operated by the Royal Observatory Edinburgh.

## REFERENCES

* Alard, C., and Lupton, R.H. 1998, Astrophys. J., 503, 325.
* Alard, C. 2000, Astron. Astrophys. Suppl. Ser., 144, 363.
* Alcock, C., et al. (MACHO team) 1993, Nature, 365, 621.
* Alcock, C., et al. (MACHO team) 1996, Astron. J., 111, 1146.
* Aubourg, E., et al. (EROS team) 1993, Nature, 365, 623.
* Baade, W. 1952, Trans. IAU, 8, 397.
* Bono, G., Caputo, F., Castellani, V., and Marconi, M. 1997, Astron. Astrophys. Suppl. Ser., 121, 327.
* Clement, C., et al.  2001, Astron. J., 122, 2587.
* Dziembowski, W., and Cassisi, S. 1999, Acta Astron., 49, 371.
* Graham, J.A. 1975, P.A.S.P., 87, 641.
* Kovács, G., et al.  2000, in: IAU Colloq. 176, “The Impact of Large-Scale Surveys on Pulsating Star Resaerch”, Ed. L. Szabados and D.W. Kurtz, ASP Conference Series, Vol. 203, 313.
* Kovács, G., and Walker, A.R. 2001, Astron. Astrophys., 371, 579.
* Moskalik, P. 2000, in: IAU Colloq. 176, “The Impact of Large-Scale Surveys on Pulsating Star Resaerch”, Ed. L. Szabados and D.W. Kurtz, ASP Conference Series, Vol. 203, 135.
* Olech, A., et al.  1999, Astron. J., 118, 442.
* Olszewski, E.W., Schommer, R.A., and Aaronson, M. 1987, Astron. J., 93, 565.
* Paczyński, B. 2000, P.A.S.P., 112, 1281.
* Pietrzyński, G., Udalski, A., Kubiak, M., Szymański, M., Woźniak, P., and Żebruń, K. 1998, Acta Astron., 48, 175.
* Popielski, B.L., Dziembowski, W.A., and Cassisi, S. 2000, Acta Astron., 50, 491.
* Roberts, D.H., Lehár, J., and Dreher, J.W. 1987, Astron. J., 93, 968.
* Schechter, P.L., Saha, K., and Mateo, M. 1993, P.A.S.P., 105, 1342.
* Schlegel, D.J., Finkbeiner, D.P., and Davis, M. 1998, Astrophys. J., 500, 525.
* Schwarzenberg-Czerny, A. 1989, MNRAS, 241, 153.
* Sharpee, B., Stark, M., Pritzl, B., Smith, H., Silbermann, N., Wilhelm, R., and Walker, A. 2002, Astron. J., 124, 2341.
* Smith, H.A., Silbermann, N.A., Baird, S.R., and Graham, J.A. 1992, Astron. J., 104, 1430.
* Stellingwerf, R.F., Gauchy, A., and Dickens, R.J. 1987, Astrophys. J., 313, L75.
* Stryker, L.L., Da Costa, G.S., and Mould, J.R. 1985, Astrophys. J., 298, 544.
* Thackeray, A.D. 1951, The Observatory, 71, 219.
* Udalski, A., Kubiak, M., and Szymański, M. 1997, Acta Astron., 47, 319.
* Udalski, A., Szymański, M., Kubiak, M., Pietrzyński, G., Woźniak, P., and Żebruń, K. 1998a, Acta Astron., 48, 147.
* Udalski, A., Soszyński, I., Szymański, M., Kubiak, M., Pietrzyński, G., Woźniak, P., and Żebruń, K. 1998b, Acta Astron., 48, 563.
* Udalski, A., Soszyński, I., Szymański, M., Kubiak, M., Pietrzyński, G., Woźniak, P., and Żebruń, K. 1999, Acta Astron., 49, 437.
* van Albada, T.S., and Baker, N.H. 1973, Astrophys. J., 185, 477.
* Walker, A. 1989, P.A.S.P., 101, 570.
* Walker, A.R., and Nemec J.M. 1986, Astron. J., 112, 2026.
* Woźniak, P. R. 2000, Acta Astron., 50, 421.
* Żebruń, K., Soszyński, I., and Woźniak, P.R. 2001, Acta Astron., 51, 303.
* Żebruń, K., Soszyński, I., Woźniak, P., Udalski, A., Kubiak, M., Szymański, M., Pietrzyński, G., Szewczyk, O., and Wyrzykowski, Ł. 2001, Acta Astron., 51, 317.