text:
# Stochastic resonance in the driven Ising model on small-world networks

## I Introduction

It has been known that a periodically modulated bistable system may display stochastic resonance (SR), arising from the cooperative interaction between random noise and periodic modulation ref:SR:general . Namely, when a system with an energetic activation barrier is subject to periodic but weak driving, the inherent thermal stochastic noise can enhance the signal out of the system rather than weaken it. Such SR phenomena, which have various practical applications ref:SR:practical , have been extensively investigated in a variety of systems: In parallel to studies of systems with relatively small numbers of degrees of freedom, the research focus has been shifting to extended complex systems with many degrees of freedom, where interesting collective dynamic behavior can emerge ref:SR:complex ; ref:ORfirst . In those studies of SR in extended systems, the underlying connection topology of dynamic variables has usually been assumed to be regular. However, recent studies of computer networks, neuronal networks, biochemical networks, and even social networks, have revealed that many real systems in nature possess quite complex structures, which can be described neither by regular networks nor by completely random networks ref:network .

In this paper we consider an Ising model on Watts and Strogatz (WS) type small-world networks ref:WS  and study the SR behavior in the presence of temporally oscillating external magnetic fields. The WS network is characterized by a short characteristic path length and high clustering ref:WS , both of which are commonly observed properties of real networks in nature. Accordingly, we believe that the study of SR in an extended system can be made more realistic if one uses the WS network as the underlying topology.

There are five sections in this paper: Section II introduces the system, i.e., the Ising model, driven by oscillating magnetic fields, on small-world networks. Defined in Sec. III is the dynamic order parameter, which describes conveniently the dynamic phase transition. The relaxation behavior is examined and the relaxation time is computed. Section IV is devoted to the stochastic resonance phenomena, characterized by double resonance peaks in the occupation ratio. The effects of long-range interactions on the resonance are investigated, revealing opposite trends in the ferromagnetic and paramagnetic phases. Finally, the main result is summarized and discussed in Sec. V.

## II Ising Model on Small-World Networks

The WS network in the present paper is constructed following Ref. ref:WS : First, a one-dimensional regular network with only local connections (of range $`k`$) is constructed with the periodic boundary condition. Next, each local link is visited once, and with the rewiring probability $`P`$ removed and reconnected to a randomly chosen node. After the whole sweep of the entire network, the average number of shortcuts in the network of size $`N`$ is given by $`NPk`$. Throughout this paper, the interaction range $`k`$ is set equal to two for convenience; longer ranges ($`k>2`$) are not expected to lead to any qualitative difference. After the WS network is built as above, an Ising spin is put on every node, and an edge (or a link) connecting two nodes is regarded as the coupling between the two spins at the two nodes.

The Hamiltonian for the field-driven Ising model on the WS network thus reads

$$H=\frac{J}{2}\underset{i}{}\underset{j\mathrm{\Lambda }_i}{}\sigma _i\sigma _jh(t)\underset{i}{}\sigma _i,$$
(1)
where $`J`$ is the coupling strength, $`\sigma _i(=\pm 1)`$ is the Ising spin at node $`i`$, the neighborhood $`\mathrm{\Lambda }_i`$ of $`i`$ stands for the set of nodes connected to $`i`$ (via either local edges or shortcuts), and an oscillating magnetic field $`h(t)=h_0\mathrm{cos}\mathrm{\Omega }t`$ is applied with the driving amplitude $`h_0`$ and frequency $`\mathrm{\Omega }`$. We perform Monte Carlo (MC) dynamic simulations, employing the heat bath algorithm ref:Newmanbook  and measuring the time $`t`$ in units of the MC time step. For thermalization, we start from sufficiently high temperatures and lower the temperature $`T`$ slowly with the increment $`\mathrm{\Delta }T=0.02`$ (in units of $`J/k_B`$ with the Boltzmann constant $`k_B`$). The driving amplitude and frequency are mostly taken to be $`h_0=0.1`$ and $`\mathrm{\Omega }=0.1`$ although different frequencies are also considered. While simulations are performed at a given temperature, the data from the first $`4\times 10^4`$ MC steps are discarded, which turns out to be sufficient for stationarity, and measurements are made for next $`10^4`$ MC steps. Networks of various sizes, up to $`N=6400`$, are constructed as described above, and averages are performed over 100 different network realizations.

In the absence of the long-range interaction ($`P=0`$), the network structure reduces to that of the one-dimensional regular network with only local couplings. Accordingly, when $`P=0`$, the driven Ising model described by Eq. (1) as well as the undriven model ($`h_0=0`$) should not exhibit long-range order at finite temperatures. For $`P0`$, on the other hand, it has been found that the (undriven) Ising model displays ferromagnetic order at finite temperatures ref:smallIsing . This suggests that the driven Ising model in Eq. (1) should undergo a dynamic phase transition at a finite temperature unless $`h(t)`$ is too large. When all nodes are fully connected, Eq. (1) describes the infinite-range Ising model, where double SR peaks have been observed and argued to be a generic property of the system with a continuous dynamic phase transition ref:BJKim . In view of this, similar double SR peaks are naturally expected when $`P0`$; this contrasts with the case $`P=0`$, where the absence of a dynamic phase transition implies the emergence of just a single SR peak ref:Brey .

## III Dynamic Phase Transition

The dynamic phase transition in the system can be conveniently described by the dynamic order parameter. We first measure the magnetization

$$m(t)\frac{1}{N}\underset{i=1}{\overset{N}{}}\sigma _i,$$
(2)
and take the time average of $`m(t)`$ during the $`n`$-th period of $`h(t)`$, to obtain

$$Q_n\frac{\mathrm{\Omega }}{2\pi }_{t_n}^{t_{n+1}}𝑑tm(t)$$
(3)
with $`t_n2\pi n/\mathrm{\Omega }`$. The dynamic order parameter is then defined to be

$$Q\underset{n\mathrm{}}{lim}Q_n,$$
(4)
which takes different values below and above the critical temperature $`T_c`$: Figure 1 shows that $`Q0`$ in the ferromagnetic phase at $`T<T_c`$ and $`Q=0`$ in the paramagnetic phase at $`T>T_c`$, with $`T_c3`$ for $`P=0.5`$. In order to determine $`T_c`$ more precisely, we measure the dynamic fourth-order cumulant ref:Sides ; ref:Binder

$$U_N1\frac{[Q_{n}^{}{}_{}{}^{4}]}{3[Q_{n}^{}{}_{}{}^{2}]^2},$$
(5)
where $`\mathrm{}`$ and $`[\mathrm{}]`$ denote the time average and the average over different network realizations, respectively. At $`T=T_c`$, the cumulant $`U_N`$ should have a unique value regardless of the size of the system; this yields the estimation $`T_c3.13`$, as shown in the inset of Fig. 1.

To investigate the time relaxation behavior, we begin with the initial condition $`m(t=0)=1`$ (i.e., $`\sigma _i=1`$ for all $`i`$), and measure $`Q_n`$ in Eq. (3) as a function of $`t_n`$ during MC simulations. As time proceeds, $`Q_n`$ approaches the dynamic order parameter \[see Eq. (4)\] with the value $`Q=0`$ and $`Q0`$ in the high- and low-temperature phases, respectively. We find that the relaxation of $`Q_n`$ is very well described by the exponential form, $`Q_nQe^{t_n/\tau }`$ with $`t_n=2\pi n/\mathrm{\Omega }`$, which defines the relaxation time $`\tau `$ ref:relaxtime . Figure 2 shows the relaxation behavior for $`P=0.5`$ and $`0.8`$ at $`T=3.5`$; it is observed that the exponential decay form describes the simulation results very well. We repeat the same procedure at other temperatures and obtain the temperature dependence of the relaxation time, shown in the inset of Fig. 2. As expected from the existence of the finite-temperature dynamic phase transition, the relaxation time diverges near the dynamic transition point ($`T_c3.13`$ and $`3.20`$ for $`P=0.5`$ and $`0.8`$, respectively), and decreases as we move away from $`T_c`$ in both directions ref:BJKim .

## IV Stochastic Resonance

In this section we study how the long-range interactions influence SR by varying the rewiring probability $`P`$. The SR behavior is conveniently captured by the occupancy ratio $`R`$, first introduced in Ref. ref:ORfirst  and defined to be the average fraction of the spins in the direction of the external field ref:BJKim:PRB :

$$R\frac{\text{number of spins in the direction of}h(t)}{\text{total number of spins}}.$$
(6)
In other words, $`R`$ measures how many spins follow the oscillating magnetic field. It is easy to understood that $`R`$ has the value 1/2 in both low- and high-temperature limits (see, e.g., Ref. ref:BJKim ) and becomes larger near the SR temperature, reflecting that more spins follow the external driving.

In Fig. 3, $`R`$ in the system with $`P=0.5`$ is plotted as a function of the temperature for various sizes $`N`$. The double SR peaks are clearly exhibited, with the dip emerging in the vicinity of the dynamic phase transition temperature $`T_c`$ determined from the crossing of $`U_n`$. As the size $`N`$ is increased, the first SR peak at the lower resonance temperature $`T_{SR}^{(1)}`$ becomes sharper and apparently saturates for $`N800`$, displaying negligible finite-size effects. The SR temperature $`T_{SR}^{(1)}`$ and the height $`h_{SR}^{(1)}`$ of the first resonance peak are displayed in the insets of Fig. 3, manifesting that there are no substantial size effects for $`N800`$. For the second SR peak at the higher resonance temperature $`T_{SR}^{(2)}`$, finite-size effects are insignificant in both the height $`h_{SR}^{(2)}`$ and the SR temperature $`T_{SR}^{(2)}`$. We thus conclude that the double SR peaks are not merely finite-size effects, and obtain $`T_{SR}^{(1)}2.70`$ and $`T_{SR}^{(2)}4.46`$ from the positions of the two SR peaks in Fig. 3.

In Ref. ref:BJKim , the positions of SR peaks in the infinite-range Ising model have been analytically obtained from the time-scaling matching condition: The intrinsic time scale given by the relaxation time $`\tau `$ should match the extrinsic time scale of the external driving. Furthermore, since $`\tau `$ diverges only at $`T_c`$ and decreases as the temperature $`T`$ is raised or lowered from $`T_c`$, the presence of double SR peaks should be a general property of the system with a finite-temperature continuous phase transition. However, the specific form of the matching condition for the infinite-range Ising model, $`\tau =1+\sqrt{1+\mathrm{\Omega }^2}`$, found in Ref. ref:BJKim  may not hold for the Ising model on the WS network studied in this work. Here we reveal the relation between the positions of double SR peaks and the relaxation time in the following manner: We first draw a vertical line at $`T=T_{SR}^{(2)}(4.46)`$ in the inset of Fig. 2 and locate the crossing point with the plot for $`1/\tau `$, from which a horizontal line is drawn. This horizontal line corresponds to the time-scale matching condition and its crossing point with the lower temperature branch of the $`1/\tau `$ plot then yields the position of the first SR peak $`T_{SR}^{(1)}`$. This gives $`T_{SR}^{(1)}2.70`$, which is in excellent agreement with the value obtained from the first peak of $`R`$ in Fig. 3.

The behavior of $`R`$ with the temperature $`T`$ at various values of the rewiring probability $`P`$ is shown in Fig. 4. When long-range interactions are not present ($`P=0`$), the Ising model does not establish ferromagnetic order at finite temperatures, thus not displaying a finite-temperature phase transition. Then expected is a single SR peak instead of two peaks ref:BJKim , which is confirmed in Fig. 4. As $`P`$ is increased, the system behaves more similarly to the mean-field system, and exhibits clear double SR peaks, which, in the mean-field (infinite-range) limit, have equal heights ref:BJKim . Note that the maximally connected case ($`P=1`$) of the WS network has $`O(N)`$ long-range connections while the infinite-range system has $`O(N^2)`$ connections. This explains why the double peaks of $`R`$ at $`P=1`$ do not have equal heights in Fig. 4.

The SR peaks in Fig. 4 may be located from the condition $`dR/dT=0`$ except for the first resonance peak at $`P=0.1`$, where we have used the condition that $`dR/dT`$ has a minimum positive value. In the insets of Fig. 4, the obtained SR temperatures and peak heights are shown as functions of $`P`$. It is noteworthy that as $`P`$ is increased the height $`h_{SR}^{(1)}`$ of the first peak also grows, whereas $`h_{SR}^{(2)}`$ reduces. In other words, addition of long-range interactions tends to enhance SR in the ferromagnetic phase, while the same addition suppresses SR in the paramagnetic phase. To understand these conflicting effects, which appear somewhat counterintuitive, we notice that as $`P`$ is increased, the critical temperature $`T_c`$ shifts toward higher temperatures, reducing the relative (temperature) distance to $`T_{SR}^{(2)}`$ in comparison with that to $`T_{SR}^{(1)}`$. In the paramagnetic phase, the correlation length should increase with $`P`$ since $`T_{SR}^{(2)}`$ becomes closer to $`T_c`$; this leads more spins to be correlated and thus resistive to the change of the external driving field, resulting in the suppression of SR. In the ferromagnetic phase, on the other hand, a larger number of shortcuts (i.e., long-range connections) increases the distance to $`T_c`$, yielding a shorter correlation length. Accordingly, the spins are allowed to follow better the external driving, and the SR is enhanced. It is also to be noted that SR in the present system tends to saturate for $`P0.7`$ and that almost the same SR as the random network ($`P=1.0`$) can be achieved with relatively small numbers of shortcuts. Similar enhancing and saturating tendency has also been found in the synchronization behavior of the coupled oscillators on small-world networks ref:Hong  as well as in the standard (single-peak) SR behavior of the system of coupled bistable elements ref:SR:Hu .

We have also considered different driving frequencies and investigated how the occupancy ratio behaves. Figure 5 displays the behavior of $`R`$ depending on the temperature $`T`$ for various driving frequencies. The interval between the two resonance peaks is observed to grow as the driving frequency is increased, which is consistent with the obtained relation between the relaxation time and the driving frequency. On the other hand, the heights of the two peaks tend to reduce with the driving frequency; eventually, the lower peak in the ferromagnetic phase disappears, which reflects that at high driving frequencies there exists only one crossing point between the relaxation time and the driving frequency.

## V Summary

We have examined the stochastic resonance phenomena in the Ising model, driven by oscillating magnetic fields, on small-world networks. Double resonance peaks have been found to develop as the rewiring probability is increased, and explained in terms of the diverging time scale at the dynamic phase transition. It has been demonstrated that the resonance behavior essentially the same as that of the random network can be achieved with relatively small numbers of shortcuts. The most interesting finding in our system is that while stochastic resonance is enhanced with the number of long-range connections in the ferromagnetic phase, it is suppressed in the paramagnetic phase. As an interpretation based on analogy, we suggest to consider the opinion formation in a social system: an instructor and a group of students. More long-range connections in this analogy correspond to stronger interactions, i.e., more active discussions, between students far away, and thus help the group behave as a whole with a majority of students having the same opinion (this situation is closely related to the enhanced synchronization for larger $`P`$ in Ref. ref:Hong ). However, this opinion developed from active discussions among students may not necessarily be the opinion preferred by the instructor of the class. The enhanced SR peak in the ferromagnetic phase corresponds to the situation that more interactions among students help the class follow the instructor, with the students having already an identical opinion (whatever it is). On the other hand, the suppression of SR in the paramagnetic phase is interpreted as the situation that more interactions make the students, among whom an agreement is not reached, not so obedient to the instructor.

###### Acknowledgements.

H.H. thanks J. Lee for providing the privilege of using the computing facility Gene. B.J.K. acknowledges the support of Ajou Univ. for the year 2002. M.Y.C. thanks the Korea Institute for Advanced Study for hospitality during his visit, where this work was performed, and acknowledges the partial support from the Korea Research Foundation through Grant No. 2000-015-DP0138.


format: text

text:
# Biases in Cometary Catalogues and Planet X

## 1 Introduction

Since the discovery of Pluto, the idea that there may be an undiscovered tenth planet (Planet X) lying still further from the Sun has often been put forward. Attempts to predict the existence of undiscovered planets from studies of the aphelia of periodic comets have been made many times (e.g., Forbes 1880; Kritzinger 1957; Opik 1971; Guliev 1992). Most of these authors argued for a close planet at $`50`$ astronomical units (au). There is also work motivated by periodicities in the cratering and fossil records of $`30`$ Myr. Here, the idea is that the advance or the regression of the perihelion or aphelion of Planet X through the Kuiper Belt triggers periodic comet showers (e.g., Whitmire & Matese 1985; Matese & Whitmire 1986). Assuming Planet X is responsible for the periodicities, then it has a semimajor axis of typically $`100`$ au and a mass of a few Earth masses. A different hypothesis, motivated by the same periodicities, is a distant unseen companion to the Sun (Nemesis) on an eccentric orbit with semimajor axis $`90000`$ au and period $`27`$ Myr (Whitmire & Jackson 1984; Davis et al. 1984). The orbit has an assumed perihelion of $`20000`$ au, close enough to make incursions into the inner Oort Cloud and trigger comet showers (see Tremaine (1986), Vandervoort & Sather (1993), for reviews of the status of these theories).

The idea has most recently surfaced in this Journal in Murray (1999), who claimed first that the aphelia of long period comets between $`30000`$ and $`50000`$ show an excess and second that they are aligned along a great circle inclined at $`30^{}`$ to the ecliptic. Murray interpreted this as evidence of the imprint of a body orbiting the Sun at $`32000`$ au with a period of $`5.8\times 10^6`$ years. Planet X is presumably massive (to produce a detectable family of comets) and its orbit is possibly retrograde, suggesting a recent capture by the Sun. A similar analysis was carried out by Matese, Whitman & Whitmire (1999), who looked at the same dataset of aphelia of long period comets. Matese et al. (1999) identified an excess of aphelia on a different great circle, this time with galactic longitude $`135^{}`$, and argued that this must be produced by a bound perturber in the Oort Cloud. They suggested that there was a planet or brown dwarf with mass $`3M_{\mathrm{Jup}}`$ at $`25000`$ au on an orbit roughly normal to the Galactic plane. The distances suggested by Murray and Matese et al. for Planet X are in good agreement, but the comets that belong to the great circle streams are quite different in the two papers, as is the orientation of the orbit.

Even so, there is a some supporting theoretical evidence for Planet X. First, studies of the ejection and growth of planetesimals in the early Solar System do suggest the existence of a number of small planets beyond Pluto (e.g., Stern 1991, Kenyon & Luu 1999). Second, there are unexplained residuals in the motions of Uranus and Neptune, although it is difficult to evaluate the trustworthiness of the historical data (e.g., Anderson & Standish 1986, Hogg, Quinlan & Tremaine 1991).

In this paper, we analyse the evidence put forward by Murray (1999) and Matese et al. (1999). We begin with a thorough examination of the biases that afflict catalogues of long period comets in Section 2 and develop a mathematical model of the biases in Section 3 to illustrate our main points. Section 4 returns to the great circle alignments found by Murray (1999) and Matese et al. (1999) and assesses whether they are unexpected or not. Simulations are used to illustrate the possible signal of a massive perturber. In the concluding Section 5, we draw attention to the fact that the GAIA satellite will provide an all-sky catalogue of comets brighter than 20th magnitude and largely free of directional biases. We argue that this may finally bring to fruition the long tradition of looking for the effects of planets on the orbits of comets.

## 2 Observational Biases

Murray used the 1994 version, and Matese et al. the 1996 version, of the catalogue of long period comets by Marsden & Williams. Our data are extracted from the 1999 version. It is vital that the observational biases in this dataset be analysed before any great circle alignments can be accepted as real. In the past, a number of people have suggested possible effects (see e.g., Everhart 1967, Hughes 1983, Kresák 1983, Rüst 1984, Neslušan 1996), though a comprehensive treatment seems lacking.

### 2.1 North-South Biases

Marsden & Williams (1999) divide the long period comets into four classes, namely 1A, 1B, 2A and 2B, in order of diminishing certainty of the orbital elements. The upper panel of Figure 1 shows the aphelia of all 1A and 1B comets. This is an equal area Hammer-Aitoff projection using heliocentric ecliptic longitude and latitude coordinates. The bold line is the projection of the Earth’s equator onto the celestial sphere. It is evident from the figure that many more cometary aphelia lie to the south rather than the north of the projection of the Earth’s equator. Table 1 shows the raw statistics on perihelion passage as seen from Earth together with the hemisphere of discovery using comets from Marsden & Williams and data from Vsekhsvyatskii (1964). The table demonstrates that objects discovered in the north generally have perihelia in the north, and vice versa. Of course, there is a known excess of observers in the northern as against the southern hemisphere. So, this bias leads to a greater number of comets with aphelia in the southern rather than the northern hemisphere. Note that comets with perihelia within $`30^{}`$ of the projected terrestrial equator may be discovered from either hemisphere. Hence, this bias is most pronounced in the sample of comets with high latitude aphelia. The lower panel of Figure 1 shows all comets (1A, 1B, 2A and 2B) in Marsden & Williams (1999) catalogue discovered prior to 1940. The bias in favour of southern aphelia is very pronounced, especially at high latitudes. This is expected because the excess of northern hemisphere observers relative to the southern hemisphere was greatest in the past.

Matese, Whitman & Whitmire (1998) contend that the north-south bias is minimized for comets whose orbits are well enough known. They argue that these comets are primarily bright. However, it could also be argued that, in order for a comet be studied enough for its orbit to be well known, it must be observed over a reasonable period of time, and that the excess of northern hemisphere observers again means that this is more likely to happen for comets with northern rather than southern perihelia.

Figure 2 shows a histogram of the number of comets in Marsden & Williams (1999) discovered per decade from the northern and southern hemisphere between 1810 and 1956, using the information in Vsekhsvyatskii’s (1964) book. This clearly shows the substantial bias towards northern hemisphere observers, which is still the case as late as the 1950s (the records in Vsekhsvyatskii’s book end in 1956). It even seems that the effects of great upheavals like the Boer War, the First World War and the Great Depression can be discerned in the diagram. However, statistics are low and so the passing interests of individual observers like Barnard, Perrine and Giacobini may also be a contributory factor. At first sight, it is surprising that there is not a dip during the Second World War. Looking at all single apparition long period comets in Vsekhsvyatskii (1964) gives us a larger sample than in the Marsden & Williams catalogue. We see that the number of comets discovered in five-year intervals from 1930 to 1950 by observers based in America rises from 3 to 4 to 5 to 7. For Southern Hemisphere observers, it changes from 3 to 2 to 6 to 13, while for observers based in Europe, it changes from 7 to 2 to 4 to 5. This suggests that the rise in discoveries from the United States and the Southern Hemisphere masks the deleterious effects of the Second World War on European discoveries. Similarly, plotting the cumulative total of numbered asteroids against time shows a steady, almost exponential increase, but there is a prominent plateau coinciding with the duration of the Second World War and ending in the 1960s. All this suggests that the cometary (and asteroidal) catalogues are not just telling us about astronomy!

### 2.2 Diurnal and Seasonal Biases

There are a number of effects introduced into the catalogue of comets due to diurnal and seasonal bias. The time of day influences the area of the sky in which comets can be discovered. In general, a comet cannot be discovered until it is properly dark (after astronomical twilight), unless the comet is unusually bright or passes very close to the Sun. Similarly, if the comet is near opposition, it is further from the Sun than at quadrature, and is therefore harder to see. This suggests that the distribution of comet discoveries has a minimum towards the Sun and a secondary minimum towards midnight. All other things being equal, we expect there to be two equal maxima, one in the evening, and one in the morning, displaced towards midnight from the point of quadrature due to the effects of twilight. The seasons also impart a bias to the data, since the length of night varies, and hence there is more time to observe in winter, and therefore more objects will be discovered in winter than in summer. This can be seen to an extent in the raw numbers from the Marsden & Williams catalogue, which are plotted in Figure 3. However, the overall effect of this is probably reduced by the fact that the weather is generally worse in winter than summer at the temperate northern latitudes in which cometary observers predominate, which leads to a lowering in the number of observations made in winter. The most obvious feature in Figure 3 is a pronounced dip in mid-summer, along with peaks in April and August. Discoveries between August and April are roughly constant (bearing in mind that February is 10 % shorter than other months). The effects of the weather seem to counterbalance the effect of longer nights in these months.

Figure 4 shows plots of aphelia of comets which passed through perihelion in the months of May, August, November and February respectively. Ideally, these plots should show the comets discovered in each month. Unfortunately, the discovery months are scattered through the literature for post-1956 comets. So, we assume that the discovery took place close to perihelion. All the comets in Marsden & Williams (1999) catalogue are used. This seems reasonable because the direction of the perihelion in heliocentric coordinates should be fairly well-known, even if the detailed orbit is subject to uncertainties. The north-south bias is most pronounced in the northern winter months (see Figures 4 (c) and (d)). This is the coupled effect of the longer nights of the northern hemisphere winter, together with the excess of northern observers. By contrast, in Figures 4 (a) and (b), which correspond to the northern hemisphere summer months, there are very nearly equal numbers of comets with perihelia in each terrestrial hemisphere. The north-south bias of observers is almost cancelled by the longer nights of the southern hemisphere winter.

The panels in Figure 4 also show a zone of avoidance (marked with a cross) at longitudes $`180^{}`$ away from the position of the Sun as viewed from the Earth. Comets cannot be discovered too close to the Sun. This means that comets with perihelia close to the Sun are unlikely to be discovered near perihelion. Holetschek (1891; see also Hughes 1983) pointed out a similar effect – namely, that comets with perihelion the far side of the Sun are less likely to be discovered both because of the proximity of the Sun and because they are further away and hence fainter. In May, the zone of avoidance is centered at an ecliptic longitude of $`225^{}`$ and advances by $`90^{}`$ as we move forward by three months. The zone is obvious in the May, August and February plots, although not at all in the November plot. Figure 5 shows the aphelia of comets which passed through perihelion in the month of July. It is included because the zone of avoidance is very clearly seen opposite the Sun at ecliptic longitudes of $`285^{}`$. Even more strikingly, the region of sky above $`30^{}`$ north of the ecliptic is completely devoid of all cometary aphelia. This is strong evidence that even in summer, the smaller numbers of southern observers are affecting the results. Another effect that is visible in the July and August plots is a tendency for comets to cluster about the apparent position of the Sun. These are comets with perihelia outside the Earth’s orbit discovered near opposition and possibly near perihelion. A northern hemisphere observer at this time of year has a comparatively small area of sky to search, due to the short summer nights. This biases cometary discoveries to the area around opposition, which explains the apparent clustering.

There is another diurnal effect that has been suggested in the literature. As is well-known, many more meteoroids are observed in the morning as compared to the evening (e.g., McKinley 1961). We encounter many more meteoroids on the forward side of the Earth, and those encountered are also brighter as they have a higher relative velocity. By contrast, in the evening the meteoroids have to catch up the motion of the Earth and so are scarcer and fainter. By analogy, Kresák (1983) suggested that such a diurnal asymmetry may be present in cometary discoveries. This, however, is not clear, as comets are not discovered through encounter with the Earth and so relative velocities do not appear to be important. Kresák’s bias may have been important in the days of naked eye discoveries, as comets that have passed through perihelion are often brighter. This may have biased statistics towards the discoveries of retrograde comets in the morning.

### 2.3 Directional Biases

There is a bias against areas of the sky with high stellar densities or dense nebulosity. The crowded nature and high surface brightness of such areas makes it harder to discover a comet against the background. Hence, we might expect a dearth of cometary discoveries, for example, in the plane of the Milky Way. However, there may be other subtle effects operating, since some areas of the sky are of much greater interest to observers than others. If an area of sky is receiving twice the attention of a neighbouring area, we might assume that the better observed area will yield a greater number of comet discoveries. An interesting example of a recent comet discovered in a region with a dense stellar background is Comet Hale-Bopp (C/1995 O1). This was found independently by Hale and Bopp, who were observing the globular cluster M70 in the direction of Sagittarius (see IAU Circular 6187). The comet just happened to lie in the same field of view.

Another directional bias comes from groups of comets that fragmented at a perihelion passage prior to discovery. The aphelia of such objects will form a tight cluster. The most obvious example is the Kreutz family of sun-grazing comets (see e.g., Marsden 1989). There are 3 Kreutz family sun-grazers in the catalogue of osculating elements of Marsden & Williams (namely, C/1882 R1-B, C/1963 R1, C/1965 S1-A). These can be seen in the upper panel of Figure 1 as the three superposed symbols at an ecliptic longitude of $`130`$ degrees.

Finally, directional biases may also come from showers of comets produced by weak stellar perturbations of the Oort Cloud. For example, Biermann, Huebbner & Lüst (1983) pointed out that there is a surplus of cometary aphelia in the region bounded by $`180^{}<\mathrm{}<240^{}`$ and $`0^{}<b<30^{}`$, where ($`\mathrm{},b`$) are galactic longitude and latitude. They attributed this to this to an encounter with a roughly solar mass star 3 Myrs ago. Irrespective of whether the details of the explanation are correct, the shower does seem to be a real dynamical effect, as it is more prominent in the more tightly bound comets (e.g., Matese et al. 1999). Further weak showers may remain undetected in the dataset.

### 2.4 Orbital Biases

Of course, comets do not all have an equal chance of discovery, as some orbits offer more favourable opportunities. It is well known that most comet searchers concentrate their efforts close to the plane of the ecliptic, and so high inclination comets are under-represented in catalogues. The bias towards low inclinations in long period comet catalogues is still further enhanced by contamination of the sample with misclassified short period comets. Comets with very small perihelion distances spend less time in the inner solar system than those at large perihelion distances and cover a much smaller arc. They are hence more likely to be missed, unless they are bright enough to be seen when very close to the Sun. Similarly, objects with large perihelion distance will be observable on a large arc. Comets with perihelia lying noticeably beyond the Earth’s orbit will be most favourably placed for discovery at opposition rather than at perihelion passage. Comets with perihelia the far side of the Sun are less likely to be discovered both because of the closeness of the Sun and because they are further away and fainter (the Holetschek effect mentioned in Section 2.2).

Comets with both high inclination and small perihelia offer particular difficulties. For example, if the the perihelion is in the north, then the comet spends most of its observable time in the southern hemisphere and may well be too close to the Sun at perihelion to be observed from the north. Table 1 demonstrates that comets tend to be discovered in the same hemisphere as the perihelion passage, but this is not true for for comets with small perihelia ($`<0.2`$ au).

Finally, objects in outburst are unpredictable and may be discovered at large heliocentric distances. Unusually bright and large comets, such as Hale-Bopp, are likely to be discovered a long time before perihelion passage. Comets are known to be much more active when new and fade on subsequent apparitions, presumably because the most volatile components are substantially lost at the first perihelion passage leading to a dimming of the intrinsic brightness thereafter. So, first time entrants from the Oort Cloud over the past 200 years (the time over which records have been maintained) are more likely to have been discovered than comets which were returning for a second or third time during the same period.

### 2.5 Sociological Biases

The heart of the problem is the uneven coverage of the sky at different epochs by different generations of observers. We have already suggested in Section 2.1 that the effects of historic upheavals like the First World War and the Great Depression may be able to be discerned in cometary catalogues. There may be additional effects that are sociological in origin. For example, especially in observations from casual observers, there may be an additional time of day bias. Many casual observers prefer to observe in the evening before going to bed, rather than in the morning after rising. Due to the difference in weather in winter and summer, casual observers in the northern hemisphere would perhaps rather observe in the warm summer months, which may slew the seasonal observations towards the months after the summer solstice. This is in accord with Figure 3, which shows that August is the best month for cometary discoveries.

### 2.6 Biases with Survey Programmes

Since 1999, the overwhelming majority of comets have been discovered by search programmes, such as LINEAR (“http://www.ll.mit.edu/LINEAR/”) and NEAT (“http://neat.jpl.nasa.gov/”). For example, 14 of the last 20 comets in the catalogue of Marsden & Williams were discovered by LINEAR and a further 2 comets were found by NEAT. It is apparent, given the success of these programmes, that in the future such surveys will constitute the overwhelming bulk of discoveries. The days of comet discoveries by amateur observers are drawing to a close. The surveys evade some of the biases above, but they introduce their own selection effects. One bias immediately obvious is that both LINEAR and NEAT are situated in the northern hemisphere. Neither will discover comets in the far south of the sky. This is the same situation as was present for visual observers in the past – a north-south bias is introduced into the dataset due to the observers rather than any astronomical effect. In addition, it appears from the skyplots of observations made by LINEAR that the area around the North Celestial Pole is not observed and so this area will show a dearth of discoveries. Another bias is that the survey telescopes only operate during ’dark time’. So, the periods around full moon will not be covered by the telescopes. This means that comet discoveries will by bunched away from periods of moonlight to a greater extent than in the past.

### 2.7 Anisotropies

Finally, there have been repeated claims of anisotropies in the aphelia directions (e.g., Yabushita 1979). The preferred direction has typically been found to coincide closely with the solar antapex. This played an important role historically as it was one of the pieces of evidence that was thought to favour Lyttleton’s (1953) accretion theory of cometary origin. The strength and direction of any anisotropy can be quantified by harmonic analysis. This idea is familiar in cosmic ray physics (e.g., Linsley 1975; Evans, Ferrer & Sarkar 2002), although apparently has not been used before in cometary science. Let ($`\lambda ,\beta `$) be heliocentric ecliptic coordinates and

$$N=\underset{i}{}h(\lambda _i,\beta _i).$$
(1)
The Fourier coefficients are

$$a=\frac{2}{N}\underset{i}{}h(\lambda _i,\beta _i)\mathrm{cos}\lambda _i,b=\frac{2}{N}\underset{i}{}h(\lambda _i,\beta _i)\mathrm{sin}\lambda _i,$$
(2)
where the index $`i`$ runs along all the comets in the sample and $`h(\lambda ,\beta )`$ is the selection function. The amplitude $`R`$ and phase $`\psi `$ of the first harmonic is

$$R=(a^2+b^2)^{1/2},\psi =\mathrm{atan}\left(\frac{b}{a}\right).$$
(3)
The probability that an amplitude higher than that measured $`R`$ can arise from an underlying isotropic distribution is given by the formula (Linsley 1975)

$$P(>R)=\mathrm{exp}(k_0),k_0=\frac{1}{4}R^2N.$$
(4)
The selection function $`h(\lambda ,\beta )`$ is both important for analysing any alleged anisotropies and nearly impossible to model. It is reasonable to argue that the cometary data accumulated over long periods of time is reasonably independent of biases in ecliptic longitude. A crude way to cope with the north-south bias is to split the cometary sample into those with aphelia in the northern or southern ecliptic hemispheres, and to analyse each sample separately.

This is performed in Table 2 for all comets in Marsden & Williams catalogue, as the the direction of the aphelion in heliocentric coordinates should be fairly well-known, irrespective of the other orbital uncertainties. The final column of Table 2 shows that the probability of an isotropic distribution yielding such large values of the amplitude $`R`$ is small. From the results of Jaschek & Vabousquet (1994), the solar antapex lies in the direction $`\lambda =88.3^{}`$, $`\beta =49.9^{}`$ (using equinox 2000). The longitude lies within $`20^{}`$ of the phase found by the anisotropy analysis of the whole sample. However, the phase does not point in a consistent direction for the three samples (comets with northern aphelia, with southern aphelia and the entire dataset) and hence, it seems that any correlation of anisotropy with the antapex direction is untrustworthy.

## 3 A Model of Discovery Probability

In this Section, we illustrate some of the biases with a simple model. Let us denote the relative probability of discovery of a comet at any heliocentric position $`𝐫`$ as $`P(𝐫)`$. This is equal to the relative probability $`P_{\mathrm{loc}}(𝐫)`$ that there is a comet at the location $`𝐫`$ multiplied by the probability $`P_{\mathrm{obs}}(𝐫)`$ that a comet at $`𝐫`$ will actually be observed.

### 3.1 The Algorithm

We assume that all comets move on parabolae and that the cometary perihelia occur uniformly in space about the Sun. The orientation of the orbit is completely specified if the rotation angle of the parabola about the major axis is given. We assume that this angle is distributed uniformly between $`0^{}`$ and $`360^{}`$. We assume that comets are discovered at random time intervals within 10 months either side of the perihelion passage. For each cometary arc, we pick a random time $`t`$ and solve for the true anomaly $`v`$ using (e.g., Danby 1992, p. 130)

$$\frac{1}{3}\mathrm{tan}^3\left(\frac{v}{2}\right)+\mathrm{tan}\left(\frac{v}{2}\right)=\frac{kt}{\sqrt{2q^3}},$$
(5)
where $`k`$ is the Gaussian gravitational constant and $`q`$ is the perihelion distance. This is used to give the heliocentric distance $`r=|𝐫|`$ at time $`t`$ through

$$r=\frac{q}{1+\mathrm{cos}v}.$$
(6)
This heliocentric position corresponds to two points on the orbit, symmetrically placed before and after perihelion. The same probability is allocated to each of these points. By generating $`10^6`$ such trials in Monte Carlo fashion, we find that the relative probability of a comet being at a position $`r`$ from the Sun is

$$P_{\mathrm{loc}}(𝐫)=\mathrm{min}(0.15r^{3/2},1).$$
(7)
The numerical pre-factor ($`0.15`$) varies only slightly as the window of discovery about the perihelion point is changed. An intuitive explanation of the power-law dependence with radius is that it is a consequence of Kepler’s Third Law. The larger the heliocentric distance, the more slowly a comet moves, and so the relative probability $`P_{\mathrm{loc}}`$ eventually reaches a plateau.

Let $`\theta `$ be the angle between the comet and the Sun as seen from the Earth. If $`\theta <18^{}`$, then the comet has set before the end of astronomical twilight. We assume that objects so close to the Sun are undetectable. If $`\theta =102^{}`$, then the comet is overhead at the end of nautical twilight. When $`102^{}<\theta 180^{}`$, then the comet is directly overhead some part of the Earth’s globe during the night. Of course, the airmass is unity when the comet is overhead. If $`18^{}<\theta <102^{}`$, then the minimum extinction is calculated by working out the highest the object can be in the sky under dark conditions. This is the zenith angle $`z=102^{}\theta `$. The definition of airmass $`a`$ is

$$a=\frac{1}{\mathrm{cos}z}.$$
(8)
When $`z<30^{}`$, then the airmass is computed from the slightly more accurate formula (e.g., Kristensen 1998)

$$a=\sqrt{0.635\mathrm{cos}^2z+2.954}0.79\mathrm{cos}z.$$
(9)
The apparent magnitude of the comet is

$$m=H_0+5\mathrm{log}\mathrm{\Delta }+2.5n\mathrm{log}r+(a1)ϵ.$$
(10)
Here, $`r`$ is the heliocentric and $`\mathrm{\Delta }`$ the geocentric distance of the comet, while $`n`$ is the “power law exponent”, which determines how rapidly the object brightens when it comes close to the Sun, and $`H_0`$ is the absolute magnitude (the magnitude at 1 au from Sun and Earth and zero phase angle). The quantity $`ϵ`$ is the number of magnitudes of extinction per airmass, which is taken as $`0.25`$ in our calculations. From the apparent magnitude, we assign the probability of observation as follows.

$$P_{\mathrm{obs}}=\{\begin{array}{cc}1\hfill & \text{if }m<10\text{,}\hfill \\ (15m)/5\hfill & \text{if }10<m<15\text{,}\hfill \\ 0\hfill & \text{if }m>15\text{.}\hfill \end{array}$$
(11)
In other words, if the comet is fainter than 15th magnitude, we assume that it is not observed. If the comet is brighter than 10th magnitude, the probability is unity. Between 10 and 15, the probability falls off in a linear manner. Of course, these choices are somewhat arbitrary, but they are a reasonably realistic simplification of a complex issue. Obviously, the limiting magnitudes are chosen to represent typical conditions over the last 200 or so years of comet searching.

### 3.2 Results

The observation probability $`P_{\mathrm{obs}}`$ multiplied by the location probability $`P_{\mathrm{loc}}`$ gives the relative probability that a comet is discovered at a given location. Figure 6 shows the probability distributions for a typical comet ($`H_0=5`$). The figure is plotted in the plane of the ecliptic and can be rotated about the Sun-Earth line to give the full three-dimensional plot. The opening angle of the two lobes is $`36^{}`$, which is a consequence of the assumption that no comets are discovered if they are so close to the Sun that they set before the end of astronomical twilight. The curvature of the lobes in the Sun-ward direction is controlled by the amount of extinction – the larger the extinction, the less likely an object is to be discovered at high airmass and the squatter the distribution. The overall size of the contours is set by the absolute magnitude and power-law index of the comet. Figure 7 shows the probability integrated along the line of sight as a function of the observed angular separation of the comet from the Sun. The left panels show a polar and the right panels a graphical representation of the results for a very faint ($`H_0=15`$) and very bright ($`H_0=0`$) comet. For a comet with such characteristics, these diagrams show the most likely angle of discovery. Note that – as anticipated in Section 2.2 – the distribution has a minimum towards the Sun and a second minimum at midnight. This effect is most obvious in the very faint comets, as bright comets can be discovered at greater heliocentric distances.

These figures hold good for observers on or close to the Earth’s equator (at latitudes $`<30^{}`$). As we move off the equator, two things happen. First, the angle at which the Sun sets varies from the vertical at the equator to the horizontal at the poles. For an object to be discovered at the poles, it would have to lie at a declination $`18^{}`$ or more to the north of the Sun. The effect of this is to vary the opening angle of the lobes. Second, the length of night and hence the maximum distance of the Sun beneath the horizon also vary over the course of the year. Again, the extreme is at the poles, which experience $`183`$ days of daylight followed by $`183`$ days of night. This causes the overall size of the figure to shrink or expand depending on the season.

So far, our analysis is appropriate for comets with a given value of $`H_0`$ and $`n`$. To provide results for an ensemble, we must integrate over the luminosity functions (LFs). We examine two LFs, the first being a simple power-law magnitude distribution (e.g., Hughes 2001)

$$dN2^{H_0}dH_0,3<H_0<15.$$
(12)
Here, $`dN`$ is the number of comets with intrinsic brightnesses in the interval $`H_0`$ to $`H_0+dH_0`$. The second is the improvement suggested by Hughes (2001), namely

$$dN10^{2.607+0.359H_0}dH_0,3<H_0<15.$$
(13)
Of course, both these laws are over-simplifications and we have extended them to magnitudes fainter than their limits of known validity. Nonetheless, this is justifiable as we are interested in understanding the qualitative effects of the increasing number of objects at faint magnitudes. Comet Sarabat of 1729 is thought to have had an absolute magnitude of -3 (e.g., Vsekhsvyatskii 1964), whilst the faintest comets considered by Hughes have an absolute magnitude of $`15`$. So we use these as our magnitude limits. Figure 8 shows the result of integrating over the LFs. The results are very similar, irrespective of whether the LF given by eq (12) or eq (13) is used. Both LFs are dominated by faint comets and so the secondary minima are very pronounced. This result is robust against quite drastic alterations in the LF.

Finally, Figure 9 illustrates the effects of the seasons. The position of the Earth is shown relative to the Sun at the equinox and at the summer and winter solstices. The position of the lobes stays constant with respect to the Earth-Sun line, which causes them to appear to rock in declination against the background sky to an Earth-bound observer. From the diagram, we can see that for a northern hemisphere observer, the probability of discovering a comet at southerly declinations is greatest in winter rather than summer.

## 4 Great Circle Alignments

Both Murray (1999) and Matese et al. (1999) claimed the detection of (different) great circle alignments in the aphelia of the long period comets. In this section, we sift the evidence that they presented with an eye as to whether biases and selection effects could have caused the suggested alignments.

### 4.1 The Planet of Murray

First of all, Murray (1999) argued that there is an excess of aphelia of long period comets in the range 30000 to 50000 au and that this may be a piece of evidence for an undiscovered planet. An immediate worry is that for objects with aphelia $`>10000`$ au, the eccentricity is close to unity. Therefore, small errors in the eccentricity can cause substantial errors in the value of $`1/a`$ and hence in the aphelion distance. In the osculating elements catalogue of Marsden & Williams (1999), there are comets classified as 1A which are nominally hyperbolic prior to their entry into the Solar System. The most extreme example is C/1996 J1-a which has a values of $`1/a`$ of $`0.000511`$. Even allowing for the fact that this is an unusual object, it is clear that the typical errors in the value of $`1/a`$ are sufficient to move cometary aphelia by $`10000`$ au at these distances. Another worry is that a maximum in the distribution of aphelia seems to arise quite naturally. Just from considerations of spatial volume alone, the number of comets with aphelia in the range \[$`Q,Q+dQ`$\] is expected to be an increasing function of $`Q`$. In the inner areas of the Oort Cloud ($`10^310^4`$ au), it is reasonable to assume that very few comets will be ejected from the Cloud by the passage of nearby stars and perturbations from the Galaxy. The further from the Sun, the more powerful are the effects of these ejection mechanisms and so the lower the probability of any individual comet surviving for the current age of the Solar System. Hence, in the outer parts of the Cloud, the number of comets with aphelia in the range \[$`Q,Q+dQ`$\] is expected to be a decreasing function of $`Q`$. It is apparent from the two competing effects that a turnover or maximum in the number of long period comets as a function of aphelion distance is inevitable and that this maximum could well occur in the range 30000 to 50000 au.

Figure 10 shows all 1A and 1B comets fulfilling Murray’s criteria (namely, with aphelia between 30000 and 50000 au) in a Hammer-Aitoff projection and in the Cartesian format as presented by Murray. In fact, Murray used only 1A comets in his original paper, arguing that these alone have sufficiently accurate orbits to be trustworthy. Marsden & Williams (1999) have four classes (1A, 1B, 2A and 2B) in order of decreasing orbital certainty and so we feel that the 1B comets should not be entirely neglected. Figure 10 presents the aphelia of the 1A comets as squares and the 1Bs as triangles. Murray’s plot (pictured in the lower panel) is not an equal area projection, but uses Cartesian axes. This does have the advantage that the equator and any great circle normal to the equator remain straight lines. However, it does come at the price of serious area and direction distortion near the poles. Any projection that is not equal area can mislead the eye in attaching too much weight to concentrations near the equator and rarefactions near the poles. The comets satisfying Murray’s criteria are recorded in Table 3. Note that our table differs in two ways from Murray’s original. First, it incorporates 2 additional 1A comets present in the 1999, but not the 1994, version of Marsden & Williams catalogue. It also includes comets classified as 1Bs, on the grounds that the direction of aphelia should still be fairly reliable, even if there are uncertainties in the other orbital elements. Second, one of the original comets identified by Murray for great circle membership (C/1993 A1) is omitted because it does not appear in Marsden & Williams (1999). Both the addition of the two new 1A comets (C/1997 J2 and C/1993 F1) and the removal of C/1993 A1 weaken the case for a great circle alignment and hence a massive object in the orbit proposed by Murray. Of the 1B comets, which have less well-determined orbits than the 1As, 4 are positioned far from the proposed path, 2 appear in fair agreement while 1 is in excellent accord. Murray’s stream appears clearest when the 1A comets only are used.

There is evidence that the dataset pictured in Figure 10 has been affected by some of the biases discussed in Section 2. It is apparent that most of the comets have southern aphelia rather than northern, as a consequence of the north-south bias. In fact, there is only one comet with aphelion at a latitude more northerly than $`30^{}`$. It seems likely that the absence of such comets is caused by observational selection effects. We argued in Section 3 that comets are usually not discovered in a band $`36^{}`$ about the Sun. The position of the Sun on the sky at the vernal equinox (21st March) is $`0^{}`$ ecliptic longitude. As the figure shows aphelia rather than perihelia, this is inverted to $`180^{}`$. This point moves forward by roughly $`30^{}`$ per month. Cometary datasets are biased against this moving band. Comets are typically discovered 2-4 months before or after the passage of the band across an area of sky. The months of discovery of Murray’s comets are given in Table 3. None of Murray’s comets were discovered in this band. So, the dataset of comets in Murray’s great circle is not entirely free from seasonal effects.

Murray’s third piece of evidence is the longitudinal trend in cometary orbital period. This is an ingenious idea which offers the possibility of illustrating the motion of Planet X around its orbit by showing the time at which the comets are presumed to have encountered it. Figure 11 shows our re-working of his plot, together with Murray’s original fit as a dashed line. Note that the postulated straight lines are a simplification valid for a low inclination planetary orbit. The removal of C/1993 A1 from the data set detracts somewhat from the goodness of fit. Of the two new 1A comets, one lies roughly in the expected position, but the other lies far away. Once the 1B comets are added to the plot, the saw-tooth shape is even less obvious, though it must again be noted that the orbits of these objects are subject to more uncertainty than their 1A cousins. Nonetheless, the technique involved in these plots is ingenious and could be used in the future to analyse similar data for evidence of a massive body.

Lastly, Murray claims that the bunching of all 13 aphelia of the 1A comets into a band within $`30^{}`$ of a great circle – and so covering only $`50\%`$ of the sky – has a probability of happening by chance of $`6\times 10^4`$. This number refers to the data in Murray’s original paper. The newer data in Table 3 merely requires 14 1A comets to fall within $`45^{}`$ of any great circle. Using Monte Carlo simulations, this happens with a probability of $`0.044`$, which is unusual but not astonishingly so. Note too that the statistical argument takes no account of selection effects. Comets with aphelion direction close to the ecliptic are more likely to be discovered than those close to the poles and so the dataset is more commonplace than the statistical argument implies.

For these reasons (the alternative explanation of the bunching of aphelia between 30000 and 50000 au; the weakening of both the sine curve fit to the aphelia and the saw-tooth fit to the longitudinal trend in cometary periods using newer data; the indication of north-south and seasonal biases in the dataset), we feel that the case for Murray’s great circle alignment is no longer compelling.

### 4.2 The Planet of Matese, Whitman & Whitmire

Matese et al. (1999) argue that evidence for a massive ($`3M_{\mathrm{Jup}}`$) body moving in the Oort cloud can be seen in the aphelion positions of 1A and 1B comets, when plotted in galactic latitude and longitude. They claim that there is an excess of objects in a great circle band defined by $`120^{}\mathrm{}150^{}`$ or $`300^{}\mathrm{}330^{}`$. Figure 12 shows the great circle alignment in galactic and ecliptic coordinates. The stream is particularly obvious amongst dynamically new comets ($`10^6/a<100`$), i.e., comets that are passing through the Solar system for the first time. Matese et al. identify 30 out of the total of 82 dynamically new comets with orbits designated 1A or 1B in the 1996 version of Marsden & Williams, which lie within the great circle band. Of these, four have orbits that are apparently hyperbolic prior to entering the inner solar system. (There are seven which also appear in Murray’s original list of objects). The upper part of table 4 lists all comets originally identified by Matese et al. (1999). The lower part lists all comets in Marsden & Williams (1999) satisfying the cut in galactic longitude alone. Three of these (C/1997 BA6, C/1999 F2 and C/1996 E1) are also dynamically new. The comets in the stream show good coverage of both the northern and southern ecliptic hemispheres, which mitigates concerns about the north-south bias. However, Matese et al.’s comet alignment is close to normal to the solar apex, which is curious but perhaps coincidental.

Figure 13 shows histograms of the aphelia of the 1A and 1B comets as a function of galactic longitude. The upper two panels show all 1A and 1B comets with no cut placed on the value of energy. The lower two panels show the equivalent histograms for the dynamically new comets ($`10^6/a<100`$). Matese et al.’s (1999) Figure 2 is a combination of the lower two plots. Without a cut-off in energy, the excess is present but not prominent. However, in the sample of new comets, and especially those whose orbits are best known, the excess is very obvious. Considering only the dynamically new objects, the stream contains 24 1A comets out of 51 and 9 1B new comets out of 37 and covers a sixth of the sky. In order to determine the significance, Monte Carlo simulations are run to assess how easily such deviations may be attributable to chance. Of course, there is no a priori basis for assuming that the great circle arc passes through the galactic poles. Therefore, it is equally impressive if the overpopulated great circle is at any orientation. This suggests that we ask for the frequency with which 33 or more out of 88 points occur in any band covering a sixth of the sky. Using Monte Carlo simulations, this occurs with a probability of $`1.4\times 10^2`$. Suppose however the sample is contaminated by a weak comet shower (like the Biermann shower) and this is responsible for 6 of the objects in the great circle alignment, which are therefore not statistically independent. In this case, we only require $`28`$ comets or more out of 83 in the stream, which occurs with a probability of $`8.2\times 10^2`$. However, to reproduce the data on the 1A comets alone, we require 24 or more comets in the stream out of a sample of 51. This arises by chance only $`1.5\times 10^3`$ times. Hence, it seems that the pattern in the data discovered by Matese et al. (1999) is quite unusual. Perhaps the most remarkable thing is that the signature is strongest in the 1A comets. The more careful we are to ensure that the comets are first time entrants, then the greater the statistical significance despite the decreasing numbers.

We use simulations to study whether a massive object in the orbit proposed by Matese et al. could give such a signal. Using the Mercury program (Chambers 1999), Planet X with a mass of $`3M_{\mathrm{Jup}}`$ is placed on a circular orbit around the Sun at a radius of 25000 au. Massless particles are then distributed in this two-body system. Each particle starts with a semi-major axis of 25000 au, the same as the planet. The orbits of the objects are chosen to give one hundred clones in eccentricity, varying in steps of 0.01 from 0 to 0.99. This seems reasonable as objects in the Oort Cloud are able to attain high eccentricities both because they were originally ejected and because of the flexing of the Galactic tide. The inclinations of the clones are $`0^{},90^{},180^{}`$ and $`270^{}`$ with respect to the planet’s orbital plane. In other words, the clones are restricted to lie either in the same plane as, or perpendicular to, the planet’s orbit. Finally, the objects are scattered in longitude of pericentre so that there are 100 clones spread evenly around the full range of $`360^{}`$. This gives a total of 40000 massless particles in the simulation. The timestep chosen for the simulation is $`10\%`$ of the orbital period of Planet X, and the simulation is run for 30 Myr.

Figure 14 shows the cumulative number of clones reaching a perihelion distance less than a certain value during the duration of the simulation. Let us regard objects attaining a perihelion of less than 5 au as satisfying the requirements for discovery as long period comets. During the course of the simulation, $`19`$ clones – all initially in the same plane as the planet’s orbit – managed to reduce their perihelion to less than 5 au. However, no object starting out of plane ever came close enough to the Earth to be discovered as a comet. So, the rough probability that an in-plane clone becomes visible over the course of 30 Myr is $`0.1\%`$. Given the numbers in Matese et al.’s comet stream, a new member is discovered on average once every 10 years. In other words, $`3\times 10^6`$ comets in the stream must be generated over 30 Myr. This suggests that the parent population of comets on perturbable orbits is $`3\times 10^9`$ objects. This is quite reasonable, as it is only $`0.3\%`$ of the total number of comets ($`10^{12}`$) in the whole of the Oort Cloud. Although our simulation is crude, it does unambiguously show two things. First, the planet on its own can perturb comets in its orbital plane so that they become visible on timescales of tens of millions of years. Second, the parent population required to maintain the observed flux of comets in the stream is not wildly implausible. So, Matese et al.’s planet warrants further and serious consideration.

## 5 Conclusions

The catalogue of long period comets (Marsden & Williams 1999) is used to investigate the independent claims of Murray (1999) and Matese, Whitman & Whitmire (1999) of evidence for a massive undiscovered body (Planet X) in the Oort Cloud of comets. This leads to an examination of the biases present in the catalogue, both returning to work by earlier authors and also studying several other biases which have apparently not been discussed previously.

The biases considered are as follows. The first is the north-south bias. Especially in the past, the excess of observers in the northern hemisphere over the southern hemisphere has led to a noticeable excess in those comets which could be observed from that hemisphere. Hence, long period comets with northern perihelia (and hence southern aphelia) predominate in catalogues. A plot of the number of comets discovered per decade from both hemispheres in the years preceding 1956 shows dips at the times of a number of cataclysms (such as the First World War). This is a indication that the long period comet catalogue is telling us about more than science. Nowadays, most comets are discovered by automated searches. However, both LINEAR and NEAT are based in the northern hemisphere and so this selection effect still persists. Second, there are a number of diurnal and seasonal biases. These introduce selection effects into the catalogue according to time of day or time of year. The time of day has an effect on whereabouts in the sky an object can be discovered. Obviously, if the comet is too close to the sun, it will not be seen. In addition, if the comet is near opposition, it is likely to be further from the sun, and hence fainter and harder to observe. This leads to minima in the observed distribution of comet discoveries both toward the sun and, less so, towards the position of opposition on the sky. During the months around midsummer, the amount of dark time available to an observer is less than that around midwinter. If other things were equal, this would in turn mean that a greater number of comets are discovered during the local winter than summer. However, this effect is somewhat lessened by the fact that the weather tends to be far less clement in the winter, and hence the amount of time to observe is reduced. Both these effects are clearly present in the comet catalogue. They create changes in the distribution that act to accentuate or to lessen the effects of other biases, such as directional biases. Comets are, in principle, harder to discover in areas of the sky with a dense stellar or nebular background than in sparser areas. However, this may well be countered by the fact that such dense areas are observed more often, which in turn boosts the likelihood of the object’s discovery (e.g., comet Hale-Bopp). A number of orbital biases necessarily afflict the catalogue. Comets with perihelia at very small distances spend much less time in the inner solar system and hence may be more likely to be missed, unless they are bright enough to observe when close to the Sun. Also, comets with large perihelion distance, beyond that of the Earth, are most likely to be discovered near to opposition, when the Earth-comet distance is the smallest. A further type of bias is that imposed by the observers themselves – for example, many observers concentrate their attention close to the ecliptic, and hence objects moving on high inclination orbits may be under-represented in the comet catalogue. In order to understand better some of the selection effects, a mathematical model has been developed which allows some of these effects to be illustrated.

Murray (1999) argues that there is an accumulation of cometary aphelia between 30000 and 50000 au and this is a piece of evidence for Planet X. However, this can also be explained by the competing effects of greater spatial volume, yet increasing depletion from Galactic tides on moving outwards in the Oort Cloud, inevitably giving rise to a maximum. Murray also identifies a possible great circle alignment in the aphelia directions of comets and a saw-tooth trend of orbital period against ecliptic longitude. However, the most recent data does not maintain these patterns. There is also evidence that the subsample in the great circle stream is affected by seasonal and north-south biases.

The stream of Matese et al. (1999) appears to stand up reasonably well to critical examination. Considering only first time entrants to the Solar system, the stream contains 33 comets out of 88 and covers a sixth of the sky. This great circle excess of comets cannot be attributed to any known selection effects. Such an excess also seems unlikely on probabilistic grounds. The frequency with which 33 or more out of 88 points occur in any band covering a sixth of the sky is $`1.4\times 10^2`$. To reproduce the data on the 1A comets (i.e., those with the most accurate orbits) is a much stiffer task. There are 24 1A comets out of 51 in the great circle stream. This occurs with a probability of only $`1.5\times 10^3`$ by chance. Hence, it seems clear that the pattern in the data discovered by Matese et al. (1999) is unusual. There are obvious concerns regarding the small dataset involved. But, it is reassuring that the trend identified by Matese et al. has been maintained with the addition of newer data. It is also reassuring that the more careful we are to ensure that the comets are first time entrants, the better is the statistical significance of the stream. Matese et al. advocate the existence of an additional planet of mass $`3M_{\mathrm{Jup}}`$ possibly on a circular orbit with a radius of 25000 au. Numerical integrations demonstrate that such a planet by itself can perturb objects in its orbital plane and reduce their perihelia to sufficiently small values ($`Q<5`$ au) so that they could be discovered as long period comets from the Earth. To maintain the observed flux of comets in the stream requires a parent population of $`3\times 10^9`$ on perturbable orbits close to the planet’s orbital plane. This seems a plausible number, as it is only $`0.3\%`$ of the total number of comets in the Oort Cloud. We conclude that Matese et al.’s planet is a possible, perhaps even likely, explanation of the unusual pattern in the data. The only alternative is to ascribe the pattern to an unlucky artefact of the complex observational biases.

Clearly, there is a need for a sample of long period comets that is free from unknown or hard-to-model selection effects. The European Space Agency satellite GAIA (http://astro.estec.esa.nl/gaia/), targetted for launch in 2010, will repeatedly survey the whole sky complete to 20th magnitude during its 5 year mission lifetime. The scanning law is complicated but readily modelled with a computer. A slow-moving solar system object is typically seen a few hundred times by GAIA during the mission lifetime, provided it is brighter than 20th magnitude. We can estimate roughly the number of long period comets that will be discovered by GAIA by scaling the experience of LINEAR, which has a very similar limiting magnitude. LINEAR has found 59 comets since 1998, averaging $`20`$ per year. We estimate that over the course of 1 year, LINEAR covers $`<1/10`$ of the sky. GAIA covers the whole sky repeatedly over the year, so may be expected to discover $`20\times 10\times 5=1000`$ long period comets during the mission lifetime. Such a dataset would allow the existence of the stream of Matese et al. (1999) to be confirmed or denied. If such a stream does indeed exist in a comet sample with easy-to-model biases, then inferences may be drawn on the perturber, which may eventually lead to its discovery.

## Acknowledgments

JH acknowledges financial support from PPARC, while NWE thanks the Royal Society. Helpful conversations with Mark Bailey were much appreciated. Both John Matese and John Murray gave us useful comments on the draft manuscript.


text:
# Charged multiplicity density and number of participant nucleons in relativistic nuclear collisions

## I INTRODUCTION

The main focus of the Relativistic Heavy-Ion Collider (RHIC) at the Brookhaven National Laboratory (BNL) is to explore the phase transition related to the quark deconfinement and the chiral symmetry restoration. The first available experimental data were the energy dependence of charged particle pseudorapidity density in central $`Au+Au`$ collisions at $`\sqrt{s}_{nn}`$=56 and 130 GeV from the PHOBOS collaboration pho1 . Soon later, the PHENIX collaboration published their data of centrality dependence of the charged particle pseudorapidity density in $`Au+Au`$ collisions at $`\sqrt{s}_{nn}`$=130 GeV phe1 . Recently the charged particle pseudorapidity density in central $`Au+Au`$ collisions at $`\sqrt{s}_{nn}`$=200 GeV was also reported by PHOBOS collaboration pho3 . The BRAHMS collaboration published more recently their results of charged particle density from Au+Au collisions at $`\sqrt{s}_{nn}`$=130 GeV bra1 ; bra2 .

It has been predicted that the rare high charged multiplicity in the final state of relativistic nucleus-nucleus collisions might indicate the formation of the Quark-Gluon-Plasma (QGP) phase in the early stage of collisions hov ; gor ; kap1 . In Ref. wang1 , the centrality dependence of the charged multiplicity has been further proposed to provide information on the relative importance of soft versus hard processes in particle production and therefore to provide a means of distinguishing various theoretical models.

The pQCD calculation with assumption of gluon saturation esk1 ; esk2  (referred to as the EKRT model later) was first used to study the centrality dependence of the charged particle pseudorapidity density at RHIC. The conventional eikonal approach and the high density QCD (referred to as the KN model later) kha1  were also used to investigate the centrality dependence and these two methods surprisingly obtained almost identical centrality dependence. Recently, authors in cape  reported their study of the same issue from the Dual Parton Model. It was found that the experimental observation, the charged particle pseudorapidity density per participant pair slightly increasing with the number of participant nucleons, was reproduced by wang1 ; kha1 ; cape , but contradicted the results of esk2 .

In this paper a hadron and string cascade model, JPCIAE sa1 , was employed to study this issue further. Within the framework of this model the experimentally measured energy dependence of the charged particle mid-pseudorapidity density per participant pair both in relativistic $`p+\overline{p}`$ and $`Au+Au`$ collisions at RHIC energies was reproduced fairly well without retuning the model parameters. The predictions for the $`Pb+Pb`$ collisions at LHC energy were also given. Both the PHENIX phe1  and the PHOBOS pho2  observations that the charged particle mid-pseudorapidity density per participant pair slightly increases with the number of participant nucleons could be reproduced fairly well by JPCIAE. In studying centrality dependence the focus was put on the uncertainties in the definitions and calculations of the number of participant nucleons. However, this study indicated that it is not suitable to use the charged particle mid- pseudorapidity density per participant pair as a function of the number of participant nucleons to constrain theoretical models for particle production, because the number of participant nucleons is not a well defined physical variable both experimentally and theoretically. The effect of the different definitions in nuclear radius (diffused or sharp) on the number of participant nucleons, the charged particle mid-pseudorapidity density, and the charged particle mid-pseudorapidity density per participant pair as a function of the number of participant nucleons was discussed. A brief version of the part contents of this paper has been published elsewhere sa2 .

## II MODELS

The JPCIAE model was developed based on PYTHIA sjo1 , which is a well known event generator for hadron-hadron collisions. In the JPCIAE model the radial position of a nucleon in colliding nucleus $`A`$ (indicating the atomic number of this nucleus as well) is sampled randomly according to the Woods-Saxon distribution

$$\rho (r)\frac{\rho _0}{1+exp((rR_A)/\alpha )}$$
(1)
where $`\rho _0`$ refers to the normal nuclear density, $`\alpha `$ 0.54 fm stands for the nuclear diffusion edge, and $`R_A`$ is the nuclear radius of nucleus $`A`$. The solid angle of the nucleon is sampled uniformly in 4$`\pi `$. Each nucleon is given a beam momentum in z direction and zero initial momenta in x and y directions. The collision time of each colliding pair is calculated under the requirement that the least approach distance of the colliding pair along their straight line trajectory (mean field potential is not taken into account in JPCIAE) should be smaller than $`\sqrt{\sigma _{tot}/\pi }`$. Here $`\sigma _{tot}`$ refers to the total cross section. The nucleon-nucleon collision with the least collision time is then selected from the initial collision list to perform the first collision. Both the particle list and the collision list are then updated such that the new collision list may consist of not only nucleon- nucleon collisions but also collisions between nucleons and produced particles and between produced particles themselves. The next collision is selected from the new collision list and the processes above are repeated until the collision list is empty.

For each collision pair, if its CMS energy is larger than a given cut, we assume that strings are formed after the collision and PYTHIA is used to deal with particle production. Otherwise, the collision is treated as a two-body collision cugn ; bert ; tai1 . The cut (=4 GeV in the program) was chosen by observing that JPCIAE correctly reproduces charged multiplicity distributions in $`AA`$ collisions sa1 . An important feature of JPCIAE at relativistic energies is that QCD parton-parton scatterings are included through PYTHIA marie , which causes charged particle yields to increase with collision energy as well as centrality since mini-jet production rates increase with energy and number of collisions suffered by participant nucleons. It should be noted that the JPCIAE model is not a simple superposition of nucleon-nucleon collisions since the rescattering among participant nucleons, spectator nucleons, and produced particles is taken into account. We refer to sa1  for more about the JPCIAE model.

Since the number of participant nucleons, $`N_{part}`$, plays a crucial role in the presentation of PHOBOS or PHENIX data we first make a study on $`N_{part}`$. As the direct measurement of $`N_{part}`$ is not available, in the fixed target experiments the number of participant nucleons from the projectile nucleus with atomic number $`A`$, for instance, is estimated by

$$N_{part}^p=A(1\frac{E_{ZDC}}{E_{beam}^{kin}}),$$
(2)
where $`E_{ZDC}`$ refers to the energy deposited in the Zero Degree Calorimeter (ZDC), dominated by the energy deposition from the projectile spectator nucleons, and $`E_{beam}^{kin}`$ is the kinetic energy of beam ahle . However, in the collider experiments, in order to obtain $`N_{part}`$ one has to relate the measurements to the Monte Carlo simulations. In PHENIX, for instance, simulations for the response of the Beam-Beam Counter and the ZDC were used to calculate $`N_{part}`$ via the Glauber model phe1 . In PHOBOS, $`N_{part}`$ is derived by relating HIJING simulations to the signals in the paddle counter pho2 . Therefore, $`N_{part}`$ here is a model- dependent variable.

Theoretically, the number of participant nucleons in a collision of $`A+B`$ at impact parameter $`b`$ can be estimated in different ways:

1. In the geometry method sa0  $`N_{part}(b)`$ reads

$`N_{part}(b)=N_{part}^A(b)+N_{part}^B(b),`$ (3)
$`N_{part}^A(b)=\rho _A{\displaystyle 𝑑V\theta (R_A(x^2+(by)^2+z^2)^{1/2})\theta (R_B(x^2+y^2)^{1/2})},`$ (4)
$`N_{part}^B(b)=\rho _B{\displaystyle 𝑑V\theta (R_B(x^2+y^2+z^2)^{1/2})\theta (R_A(x^2+(by)^2)^{1/2})},`$ (5)
where $`\theta (x)`$=0 if $`x<`$0 and $`\theta (x)`$=1 otherwise, $`\rho _A`$ ($`\rho _B`$) is the nuclear density of the projectile (target) nucleus and the nuclear density is normalized to the atomic number.
2. In the Glauber model, $`N_{part}(b)`$ is calculated through

$$N_{part}(b)=d^2sT_A(\stackrel{}{b}\stackrel{}{s})[1\mathrm{exp}(\sigma _{in}T_B(\stackrel{}{s}))]+d^2sT_B(\stackrel{}{s})[1\mathrm{exp}(\sigma _{in}T_A(\stackrel{}{b}\stackrel{}{s}))]$$
(6)
where $`\sigma _{in}`$40 mb is the inelastic nucleon-nucleon cross section at RHIC energies and $`T_A`$ ($`T_B`$) refers to the nuclear thickness function of nucleus $`A`$ ($`B`$) and is normalized to $`A`$ ($`B`$) esk2 .
3. In the dynamical simulation of $`A+B`$ collisions, $`N_{part}(b)`$ can be estimated via counting the participant or spectator nucleons and averaging over events simulated at a given impact parameter $`b`$. However, there are multifarious in simulating models such as: FRITIOF pi1 , VENUS venu , HIJING wang2 , JPCIAE sa1 , UrQMD uqmd , and AMPT ampt , etc. . Not only is the theoretical uncertainty related to the calculation of $`N_{part}(b)`$ in each of the above models large but also the definitions of participant nucleons or spectator nucleons are different among each other. We give only a necessary description for the following dynamical simulations mentioned in this paper:

* In FRITIOF pi1  and HIJING wang2 , the wounded nucleons, i.e., nucleons which suffer at least one inelastic collision, are counted and identified as $`N_{part}(b)`$. It should be pointed out here that in FRITIOF and HIJING the produced particles from the string fragmentation do not have rescattering.
    * In JPCIAE, we have counted the nucleons involved in at least one inelastic nucleon-nucleon collision with string excitation and identified them as $`N_{part}(b)`$.
    * In AMPT ampt1 , the spectator nucleons, $`N_{spec}(b)`$, are counted in the final state without rescattering and $`N_{part}(b)`$ is then calculated through

    $$N_{part}(b)=(A+B)N_{spec}(b).$$
    (7)
    The spectator nucleons here refer to the nucleons with zero transverse momentum and beam energy in the final state of the AMPT simulation.

In the PHENIX or PHOBOS experiment, centrality bin was defined by the cut in the particle multiplicity distribution and expressed as the percentage of geometrical (total) cross section phe1 , $`g`$. However, in the theoretical calculation it is more convenient to define the centrality by impact parameter $`b`$. In order to compare the experimental data of centrality dependence with the theoretical results a relation between $`g`$ and $`b`$ is required, which is obtained from the definition of the geometrical cross section, that reads

$$b=\sqrt{g}\times b_{max},b_{max}=R_A+R_B,$$
(8)
where $`R_A=1.12A^{1/3}+0.54`$ fm, for instance, is the radius of nucleus $`A`$ with diffused edge. However, $`R_{Au}`$= 7.5 fm is taken in our calculations since this value (corresponding to $`\sigma _{geo.}`$=7.2 b) was used in phe1  to extract the number of participant nucleons in $`Au+Au`$ collisions. We give in Tab. 1 the mapping between the centrality bin ($`g`$ bin, selected from Tab. 1 in phe1 ) and the $`b`$ bin according to Eq. (8), and the averaged impact parameter, $`\overline{b}`$, over the $`b`$ bin according to the $`b^2`$ law.

As the experimental data were averaged over events in each $`g`$ bin the theoretical results, to be compared with the experimental data, are also averaged over the corresponding $`b`$ bin. We denote the number of participant nucleons after averaging over $`g`$ or $`b`$ bin as $`<N_{part}>`$ later.

## III RESULTS AND DISCUSSIONS

In Fig. 1(a) the experimental data of charged particle pseudorapidity density per participant pair at mid-pseudorapidity in relativistic $`p+\overline{p}`$ (open triangles and rhombuses with error bar) and in central $`A+A`$ collisions (open circles with error bar for $`Pb+Pb`$ at SPS and full circles with error bar for $`Au+Au`$ at RHIC) pho1 ; pho3  were compared with the results of the JPCIAE model (full stars for $`p+\overline{p}`$, open squares for $`Au+Au`$ collisions at RHIC, and full squares for $`Pb+Pb`$ at $`\sqrt{s}_{nn}`$=17.3 and 5500 GeV). In addition, the results from other models were plotted as follows: the dashed and dotted-dash curves are from the HIJING model wang1  with and without jet quenching, respectively, the solid curve are from the EKRT model esk2 . The EKRT results were obtained from wang1  directly, except that the EKRT result for $`Au+Au`$ collisions at $`\sqrt{s}_{nn}`$=5500 GeV was taken from esk3 . Fig. 1 (b) is the same as (a) but for $`A+A`$ collisions only and the vertical coordinate here is the charged particle pseudorapidity density itself. One knows from Fig. 1 that both the data of $`p+\overline{p}`$ and $`A+A`$ collisions at relativistic energies are also reproduced fairly well by the JPCIAE model.

In Fig. 2 the $`<N_{part}>`$ extracted by PHENIX phe1  and PHOBOS pho2  from $`Au+Au`$ collisions at $`\sqrt{s}_{nn}`$=130 GeV were compared with the model computations. The horizontal axis in Fig. 2 is the percentage of geometrical cross section, $`g`$, and each $`g`$ bin is represented by its middle point for convenience in plotting (the same for Fig. 3, 5, and 6). The solid and open circles with error bar in Fig. 2 were PHENIX and PHOBOS results, respectively. $`<N_{part}>`$ from the geometry method, open stars, were from the geometric $`N_{part}(b)`$ (Eq. 3 - 5) after averaging over $`b`$ sampled randomly in each $`b`$ bin due to the $`b^2`$ law. Similarly, the FRITIOF, JPCIAE, and AMPT ampt1  results, full stars, full triangle-ups, and full triangle-downs in Fig. 2, were obtained averaging over events simulated for $`b`$ sampled randomly in each $`b`$ bin due to the $`b^2`$ law, respectively. In the EKRT model esk2 , the curves of $`N_{part}(b)`$ vs. $`b/R_A`$ were given and $`N_{part}`$ is calculated by the Glauber method with $`R_A`$=1.12$`A^{1/3}`$-0.86$`A^{1/3}`$. The open triangle-ups in Fig. 2 are the EKRT model results taken from the $`\sqrt{s}_{nn}`$=130 GeV one of those curves according to the $`\overline{b}`$ in Tab. 1. Originally, HIJING $`N_{part}`$ wang3  were calculated for individual $`b`$, the HIJING points, open squares, in Fig. 2 were plotted after relating $`b`$ to $`g`$ according to Eq. (8). The full squares in Fig. 2 were the results of the KN model taken from Tab. 1 in kha1  under the centrality bins of 0 - 6, 10 - 20, 20 - 30, 30 - 40, and 40 - 50 %, respectively. In the KN model, $`<N_{part}>`$ was the result of $`N_{part}(b)`$ from the Glauber model after averaging over centrality bin. In this approach, the particle multiplicity and impact parameter were related by a Gaussian distribution with parameters fixed via fitting the PHOBOS charged multiplicity distribution. As proved in wojc , such kind of average is approximately equivalent to the average method based on Eq.(8). From Fig. 2 one knows that the discrepancies among the PHENIX or PHOBOS and the model results are visible.

The charged particle pseudorapidity density at mid-pseudorapidity in $`Au+Au`$ collisions at $`\sqrt{s}_{nn}`$=130 GeV as a function of the percentage of geometrical cross section was given in Fig. 3. In this figure, the full and open circles with error bar are the PHENIX phe1  and PHOBOS pho2  data, respectively. The full and open triangles, respectively, are the JPCIAE results with and without rescattering. One knows from Fig. 3 that the rescattering only leads to a few percent increase in the charged multiplicity although rescattering might enhance yields of strangeness, $`\mathrm{\Xi }^{}+\overline{\mathrm{\Xi }^{}}`$ for instance, by a couple of times.

In Fig. 4 (a) we compared the PHENIX data of charged particle mid-pseudorapidity density per participant pair (full circles with shaded area of systematic errors) phe1  with the results of the JPCIAE model (full triangles) and the results of other models (obtained from phe1  directly): HIJING (the dotted curve), the KN model (the solid curve), and EKRT (the dashed curve). One sees that except EKRT, three other models predict an increase of $`(dN_{ch}/d\eta |_{\eta =0})/(0.5<N_{part}>)`$ as a function of $`<N_{part}>`$ though the theoretical results seem to underestimate the PHENIX data. Such an increase can be understood in JPCIAE as a result of increasing hard parton scatterings per participant nucleon. Fig. 4 (b) compared the PHENIX data to the results of single $`dN_{ch}/d\eta |_{\eta =0}`$ from JPCIAE normalized by the $`<N_{part}>`$ from different models (taken from the corresponding curve in Fig. 2 at the middle point of $`g`$ bins for KN, HIJING, and EKRT models): full squares, $`<N_{part}>`$ from the KN model, open squares from HIJING, open triangles from EKRT, full stars from FRITIOF, and full triangles from JPCIAE. One sees from Fig. 4 (b) that starting from the charged particle mid-pseudorapidity density obtained in JPCIAE, but using $`<N_{part}>`$ from different models, the corresponding results of $`(dN_{ch}/d\eta |_{\eta =0})/(0.5<N_{part}>)`$ are different visibly among each other, in peripheral collisions especially. Therefore it might be inappropriate using $`(dN_{ch}/d\eta |_{\eta =0})/(0.5<N_{part}>)`$ as a function of $`<N_{part}>`$ to distinguish various theoretical models for particle production since $`<N_{part}>`$ is not a well defined physical variable.

If one uses the sharp nuclear radius, $`R_A=1.12A^{1/3}`$ fm ($`R_{Au}`$=6.50 fm), instead of diffused one used above, the corresponding mapping between the centrality bin ($`g`$ bin) and the $`b`$ bin from Eq. (8) is given in Tab. 2 . Fig. 5, 6, and 7 give, respectively, the comparisons between the JPCIAE results of diffused nuclear radius and the JPCIAE results of sharp nuclear radius in the $`<N_{part}>`$ as a function of centrality, the $`dN_{ch}/d\eta |_{\eta =0}`$ as a function of centrality, and the $`(dN_{ch}/d\eta |_{\eta =0})/(0.5<N_{part}>)`$ as a function of $`<N_{part}>`$. One knows from Fig. 5 and 6 that there are observable discrepancies between the results of diffused nuclear radius and the results of sharp nuclear radius , in peripheral collisions especially. However, in the $`(dN_{ch}/d\eta |_{\eta =0})/(0.5<N_{part}>)`$ vs. $`<N_{part}>`$ plot both results are close to each other as shown in Fig. 7 . That indicates again that it is hard using $`(dN_{ch}/d\eta |_{\eta =0})/(0.5<N_{part}>)`$ as a function of $`<N_{part}>`$ to distinguish various theoretical models for particle production since $`N_{part}`$ is not a well defined variable. As one sees from Fig. 5 and 6 that for most central collision both diffused and sharp nuclear radii work well in mapping centrality bin and $`b`$ bin therefore we do not redraw a figure as Fig. 1 for sharp nuclear radius. However, it should be pointed out that using sharp nuclear radius in mapping centrality bin and $`b`$ bin is inconsistent with that the nuclear diffusion edge is introduced in Woods-Saxon distribution initiating the nucleons in nucleus.

## IV CONCLUSIONS

In summary, we used the hadron and string cascade model, JPCIAE, to investigate the energy and centrality dependences of charged particle pseudorapidity density at mid-pseudorapidity in relativistic $`p+\overline{p}`$ and $`A+A`$ collisions. Both the relativistic $`p+\overline{p}`$ experimental data and the PHOBOS and PHENIX data of $`Au+Au`$ collisions at RHIC energies could be reproduced fairly well within the framework of the JPCIAE model without retuning any parameters. The JPCIAE model predictions for $`Pb+Pb`$ collisions at the LHC energy were also given. Both the PHENIX phe1  and the PHOBOS pho2  observations that the charged particle mid-pseudorapidity density per participant pair slightly increases with the number of participant nucleons could be reproduced fairly well by JPCIAE. This study shows that since $`<N_{part}>`$ is not a well defined physical variable both experimentally and theoretically it may be hard to use charged particle pseudorapidity density per participant pair at mid-pseudorapidity as a function of $`<N_{part}>`$ to distinguish various theoretical models for particle production. A discussion for the effects of the different definitions in nuclear radius (diffused or sharp) is given and it is indicatd that using diffused nuclear radius might be better.

ACKNOWLEDGMENTS

Finally, the financial supports from NSFC (19975075, 10135030, and 10075035) in China and DOE in USA are acknowledged.