text:
# Peculiar Extended X-ray Emission around the “Radio-Loud” Black Hole Candidate 1E1740.7-2942

## 1 Introduction

Discovered by the Einstein satellite (Hertz & Grindlay 1984), 1E1740.7-2942 appears as a persistent X-ray source in the sky, only about 50′ away from the center of the Galaxy. It is in fact the brightest persistent source at high energies ($``$ 30 keV) within a few degrees of the Galactic center (e.g., Skinner et al. 1987; Sunyaev et al. 1991). At low energies ($``$ 3 keV), however, the X-ray emission from the source is strongly attenuated by the intervening matter along the line of sight. The inferred hydrogen column density is around $`10^{23}\text{ }cm^2`$ (see Sheth et al. 1996, Churazov et al. 1996, and Sakano et al. 1999 for most recent measurements), which is unusually high even for sources physically located near the Galactic center. Additional absorption can be attributed to the material associated with a molecular cloud in which the source is likely embedded (Bally & Leventhal 1991; Mirabel et al. 1991; Phillips et al. 1995; Yan & Dalgarno 1997).

1E1740.7-2942 has a very hard X-ray spectrum (measured up to $``$300 keV; e.g., Skinner et al. 1991; Cook et al. 1991; Kuznetsov et al. 1997), which is characteristic of black hole candidates (BHCs) such as Cyg X-1. The source also shows strong aperiodic and quasi-periodic X-ray variability (Smith et al. 1997; Lin et al. 2000), again similar to BHCs. The observed X-ray properties, therefore, make 1E1740.7-2942 a BHC, although dynamical evidence for such a candidacy is still lacking. The high foreground extinction precludes the identification of the optical counterpart. Extensive searches for the infrared counterpart have also been unsuccessful (e.g., Leahy et al. 1992; Marti et al. 2000). However, 1E1740.7-2942 shows up prominently at radio wavelengths (Mirabel et al. 1992). In fact, the most remarkable property known of the source is the presence of persistent, two-sided radio jets, reminiscent of extragalactic radio-loud quasars. Such comparisons support the notion that 1E1740.7-2942 is an X-ray binary that contains an accreting black hole.

Shortly after its discovery, 1E1740.7-2942 generated a tremendous amount of excitement for being a possible source of electron-positron annihilation radiation (Bouchet et al. 1991; Sunyaev et al. 1991), which had been detected from the general direction of the Galactic center by earlier balloon-borne and space-borne experiments (see Harris 1997 for a historical account). Suggestions were subsequently made that the jets are perhaps made of electron-positron pairs ejected from the central black hole at relativistic speeds (Mirabel et al. 1992) and that the particles will eventually be slowed down and annihilated in the surrounding cold, dense molecular cloud. However, the reality of 1E1740.7-2942 being the “great annihilator” in the Galaxy has been brought into fierce debate, mostly because of the failure of detecting any annihilation features by the instruments aboard the Compton Gamma-ray Observatory, even in some of the simultaneous or contemporaneous observations with those of SIGMA/GRANAT on which the claims of detections were based (see review by Harris 1997). This remains a highly controversial issue.

To shed more light on the unusual environment surrounding 1E1740.7-2942, we observed the source with the Chandra X-ray Observatory (Chandra; Weisskopf & O’Dell 1997), taking advantage of the unprecedented spatial resolution that Chandra offers (van Spybroeck et al. 1997). In this paper, we report the discovery of peculiar extended X-ray emission around 1E1740.7-2942, and present the first high-resolution X-ray spectra of the source.

## 2 Observation

1E1740.7-2942 was observed for about 10 ks on September 25, 1999 with the High-Energy Transmission Grating Spectrometer (HETGS; Canizares et al. 2000) on board Chandra, shortly after the Orbital Activation and Checkout phase of the mission. The grating spectra are read out by the spectroscopic array of the Advanced CCD Imaging Spectrometer (ACIS; Garmire et al. 2000) at the focal plane of the X-ray telescope. For the observation, we adopted an alternating exposure data mode with three primary (or long) frames (with a readout time of 3.3 s) followed by one secondary (or short) frame (with a readout time of 0.3 s). The reason for using the short frames is to obtain some zeroth-order data that is essentially free of the effects of “pile-up” (i.e., multiple photons hit the same event detection cell of the ACIS within one readout frame). This is important for observations of bright sources, such as 1E1740.7-2942, where serious pile-up is expected to distort the zeroth-order image. Although an unpiled-up spectrum can be obtained from the higher-order grating data in this case, such data is of little use for imaging analysis. However, the use of alternating exposure modes reduces the observing efficiency. In our case, the effective exposure time is only about 6.7 ks and 0.2 ks for the long-frame data and short-frame data, respectively.

## 3 Data Analysis and Results

We received data products of the standard processing pipeline from the Chandra X-ray Center (CXC), in which the long-frame data and short-frame data were conveniently separated. We checked the (Level 2) products against all known caveats (listed at the CXC web site) for possible problems and found none (except for a known offset in the aspect solutions; see below). We then analyzed the Level 2 data with the Chandra Interactive Analysis of Observations software (version 1.1) which we obtained from the CXC. For some aspects of the spectral analysis we also used custom software as well as XSPEC version 10.0.

As expected, the global image (not shown) shows a compact X-ray emission region in the zeroth-order data and a characteristic “X” pattern of higher-order (or dispersed) spectra from two grating assemblies, the High Energy Grating (HEG) and Medium Energy Grating (MEG). The intrinsic energy resolution of the ACIS is relied upon for sorting out the overlapping orders of the dispersed spectra. In addition, a faint trace is visible in the long-frame image that runs through the zeroth-order image along the direction of frame transfer. This is an artifact resulting from photons hitting the CCD chip during the frame transfer operation. Therefore, such a “readout trace” is present only in the observations of bright X-ray sources. No other point sources were significantly detected over the entire field of view (roughly 8′ $`\times `$ 50′).

### 3.1 Extended X-ray Emission

Zooming in on the zeroth-order image (shown in Fig. 1), we see that the innermost emission region appears to be elongated in the direction of northeast and southwest (but mostly east and west). In Fig. 1, we displayed an image made from the short-frame data and overlaid it with contours from the long-frame data. The image looks reasonably smooth and clearly shows a maximum at, presumably, the location of the point source. In contrast, the contours appear quite “lumpy” and actually show a local minimum at the source position. This is almost certainly due to the distortion of point spread function (PSF) caused by significant pile-up. The pile-up effects are obviously most severe at the position of the source and can cause an apparent “valley” there in the X-ray image. The effects can be crudely quantified by deriving the expected count rate of source in the zeroth-order image from that in the higher-order image, since the latter does not suffer from any pile-up. We determined that the expected zeroth-order count rate of the source should be about 1.9 count/s (or 6.3 count/frame), which explains the severe distortion of the PSF observed in the long-frame image. Actually, the situation is made worse by the fact that 1E1740.7-2942 has an unusually hard spectrum, which makes it more likely to cause additional “lost and undetected” events (Allen et al. 1998). Indeed, the measured count rate (from the zeroth-order image) is only about 0.03 count/s (compared to about 0.2 count/s for the extended emission surrounding the source). For the short-frame data, we estimated the count rate by adding up all the charge in the 3x3 event detection island around the central pixel in the 0th order image and determined a count rate of 1.28 count/s (or 0.38 count/frame), which corresponds to a $``$30% pile-up (compared to the expected 1.9 count/s). For comparison, the measured count rate of the surrounding extended emission is about 0.3 count/s (which is significantly larger than that from the long-frame data, again indicating severe pile-up effects in that case). To quantify the effects of lost-and-undetected events for an unusually hard source like 1E1740.7-2942, we also ran PIMMs to estimate the pile-up rate. We obtained a significantly lower value, $``$20%, which can clearly be attributed to the fact that PIMMS does not take into account the lost-and-undetected events. Even at a pile-up rate of 30%, little distortion of the PSF is expected, based on the results from pre-flight ground calibration tests (Allen et al. 1998). However, the short-frame image suffers from the lack of statistics, due to the short exposure, which is why it does not show the faint “halo” emission (only about 0.04 count/s) surrounding the point source and the elongated feature (see outer contours in Fig. 1).

To quantify the elongation seen in the zeroth-order image, we ran MARX simulations (Wise et al. 1999) for the HETGS 0th order image of a point source which included the observed level of background. In Fig. 2., we fit the simulated MARX PSFs to the short-frame count profiles collapsed along the elongation and across it. The measured profile (solid line) across the elongation (top of Fig. 2) cannot be distiguished from a point source at a confidence level of 99.9$`\%`$. Here we determined the full-width-at-half-maximum (FWHM) of the profile to be 0.83″ (please note, the pointing of the observation was slightly off-axis to avoid the CCD node boundary). At roughly the same confidence level, we can exclude the possibility that the measured profile along the elongation is point-like. The width of the elongation is 1.35″ FWHM.

Looking at Fig. 1 more carefully, the long-frame contours seem to indicate the presence of another region of extended emission, which lies to the south of the main feature. Given the severe distortion of the image by pile-up effects, however, it is entirely possible that this and the main feature are the integral part of the same extended emission. The smoothness of the short-frame image appears to support that. In addition, the fainter and more extended emission at the outer part of the image is statistically significant. It is approximately spherically symmetric about the point source and can be measured up to about 4″ away from the source.

### 3.2 Position of 1E1740.7-2942

We derived the position of the point source in two ways. First, we constructed the zeroth-order count profiles along the right ascension (RA) and declination (Dec) from the the un-distorted short-frame image (in sky coordinates). We fitted the profiles with a simple Gaussian function to determine the position of the peaks. Although Gaussian functions are not a very good representation of the overall count profiles (especially around the wings), they are used only to fit a limited range around the peak of the profiles for determing the peak position. This allows us to measure the position to a sub-pixel accuracy, which is only limited by the signal-to-noise ratio (S/N) of the short-frame data.

Second, we utilized the fact that the point source should be at the intersection point of the readout trace and the dispersed MEG and HEG images. In this case, the problem of deriving the source position is over constrained, which allows a more realistic estimate of the uncertainties (including all systematic effects). We started by fitting a Gaussian function to the count profile across the readout trace (in the long-frame image) to determine its center. This gives the position of the source along one direction (nearly along the Dec) in sky coordinates. We then determined the position of the intersection point between the readout trace and the MEG or HEG images. The average of the two positions was used as the position of the source along the orthogonal direction (nearly along the RA). The difference between the them provides an estimate of the uncertainty in the source position derived using this technique.

For comparison, the results from both methods are summarized in Table 1. It is worth noting that the first method produced slightly smaller error bars, in spite of poorer statistics of the short-frame data. This can be understood in terms of systematic uncertainties in the dispersed HETGS traces and the faintness of the readout trace. We also note that there is a known offset ($``$2.09″) in the aspect solutions for our observation (CXC web page). The numbers shown in the table have already been corrected for the offset. It is reassuring that the two methods are in an excellent agreement.

### 3.3 X-ray Spectra

We constructed the first-order ($`\pm 1`$) spectra of 1E1740.7-2942 from the MEG and HEG data. The higher-order data is ignored because of poor statistics. The individual spectra are compared for identifying common features which may only be marginal statistically in each case. The general agreement between the two MEG or HEG spectra or between the MEG and HEG spectra is good, although the statistics of the HEG data is poorer. To improve S/N, we first added the plus and minus orders for the MEG and HEG data separately; we then co-added the MEG and HEG first orders, with special care in deriving the wavelength scale for each spectrum. Fig. 3 shows the result. The spectrum was binned to $``$0.02 Å (which is about the energy resolution of the MEG). Note that there are hardly any counts above $``$5.5 Å , because the source is known to be heavily absorbed. A crude fit to the continuum with a simple power law plus absorption yielded a hydrogen column density $`11.8\pm 0.6\times 10^{22}\text{ }cm^2`$ and a photon index $`0.90\pm 0.15`$, in reasonable agreement with the ASCA results (Sheth et al. 1996; Churazov et al. 1996; Sakano et al. 1999). The integrated observed 2-10 keV flux is $`2.2\times 10^{10}\text{ }ergs\text{ }s^1\text{ }cm^2`$. The flux error is of the order of 10$`\%`$, which is the current status of calibration in the HETGS at high energies. The best-fit continuum model is also shown in Fig. 3 (in solid histograms) but is arbitrarily shifted down by 15% to accentuate possible line features (see below).

Deviations from the continuum model are quite apparent, mostly in the form of weak, narrow emission lines (though detected only at a significance level of roughly $`2\sigma `$; see Table 2). The feature at $``$1.74 Å (or $``$7.1 keV) appears to be an absorption edge, perhaps due to cold Fe atoms associated with the interstellar medium. The feature centered at $``$1.78 Å (or $``$6.96 keV) is probably due to the recombination line of H-like Fe ions (Fe XXVI). Note that the feature appears quite broad. Curiously, there appears to be two lines at 1.70 Å and 1.96 Å, respectively, which cannot be easily identified with any candidate elements. We speculate that these lines might also originate in the Fe XXVI recombination line but be blue- and red-shifted by the approaching and receding jets, respectively. It is worth pointing out that in this scenario the asymmetry of the lines about the rest-frame energy (see Fig. 3) would be expected because of the transverse Doppler shift. If the interpretation is correct, the observed radial and transverse Doppler shifts would imply that the velocity of the jets is about 0.26 c and the inclination angle (with respect to the line of sight) about 60°, assuming the approaching and receding jets are intrinsically identical. Applying the scenario to the remaining line candidates, we found that about 90% of them could be identified with H-like or He-like ions of Ca, S, and Ar and that the inferred blue- and red-shifted counterparts would give roughly the same jet parameters. Table 2 shows a list of the emission lines that are marginally detected, along with possible identifications, as well as the proposed Doppler-shifted counterparts of the lines. In contrast, we could only manage to identify about 40% of the lines without invoking any Doppler effects. Note that for the line candidates we only consider those that are present in at least three out of four first-order ($`\pm 1`$) MEG and HEG spectra. Any real spectral features in the HEG spectra should in principle always be present in the MEG spectra because of the improved statistics. The only exception is that a feature could fall right in the gap between two CCD chips along the MEG arm of the “X” pattern.

The zeroth-order spectrum could, in principle, shed some light on the spectral properties of the extended emission. Unfortunately, severe pile-up causes significant distortion of the spectrum obtained from the long-frame data, while the short-frame data is not very useful due to its poor S/N. We could not tell with sufficient confidence whether there are any emission lines that are not present in the dispersed spectra (and thus would be associated with the extended emission). For comparison, when we binned the observed grating spectrum (Fig. 3) to the ACIS energy resolution (e.g., 0.05 Åat 1.85 Å), all features vanish, except for a merge between the ones at 1.71 Åand 1.78 Å. We chose not to pursue the subject any further in this work.

## 4 Discussion and Conclusions

For the first time, Chandra made it possible to image a possible black hole jet system, 1E1740.7-2942, with sub-arcsecond resolution in X-rays. Our primary motivation for conducting the observation was to search for X-ray emission from the jets themselves, both by direct imaging and by searching for Doppler-shifted emission lines as in the well-known case of SS 433 (Kotani et al. 1996).

We found no evidence for any jet emission in the X-ray image. Instead, we discovered a region (or regions) of extended emission. For comparison, a radio image of the system is shown in Fig. 4, with the inset showing an X-ray image of the region around the radio core that is overlaid by radio contours with comparable resolution. We can see that the main X-ray feature is elongated in the direction roughly orthogonal to the axis of the radio lobes. The extent of the elongation is about 3″, which corresponds to a linear size of $``$0.12 pc, assuming that the distance to the source is 8.5 kpc. It seems, therefore, much too large to be associated with an accretion disk. We note that similar features have been detected in a recent VLBA+VLA observation of SS 433 (Paragi et al. 2000), although the extended radio emission is about a factor of 30 smaller than the main X-ray feature seen here in 1E1740.7-2942. It was speculated that the extended radio features in SS 433 might originate in the non-thermal radiation from relativistic electrons accelerated by shocks that are formed due to the interaction between the ejected matter from the accretion disk and the ambient medium (Paragi et al. 2000). Such an “excretion disk” scenario (through the L2 Lagrangian point) was thought to be possible for SS 433 during episodes of enhanced mass accretion. Perhaps, this scenario is also relevant for the much larger extended feature in 1E1740.7-2942. In this case, the extended X-ray emission observed could be due to the reprocessed radiation from the outflowing material that is illuminated by the central X-ray source. Alternatively, there might exist an equatorial wind from the companion star which produces the extended X-ray emission in a similar manner. In either case, fluorescent lines would be expected in the spectrum of the extended emission. Unfortunately, our data is not of sufficient quality to verify that.

On the other hand, it is perhaps worth asking whether the extended X-ray emission could be associated with a previously unknown young supernova remnant (SNR). To see whether such a small SNR (of radius only 0.06 pc) is plausible, we adopt the Sedov-Taylor similarity solution (e.g., Osterbrock 1989):

$$R_s=2t_{100}^{2/5}E_{51}^{1/5}N_0^{1/5}\text{ }pc,$$
(1)
where $`R_s`$ is the radius of the shock front, $`t_{100}`$ is the time since the explosion in units of 100 years, $`E_{51}`$ is the initial kinetic energy in units of $`10^{51}`$ erg, and $`N_0`$ is the number density of the ambient medium. Clearly, for any reasonable $`t_{100}`$ and $`E_{51}`$, an enormous ambient density would be required to produce a SNR of the size of the observed extended X-ray feature. Besides, 1E1740.7+2942 and its vicinity has been observed with VLA, and no SNR has ever been detected (or reported). Therefore, we believe that the SNR scenario can be safely ruled out.

Although the identification of the radio counterpart of 1E1740.7-2942 is reasonably firm (Mirabel et al. 1993), there were always lingering doubts about the association of the X-ray source with the compact radio core. The sub-arcsecond spatial resolution of Chandra has now allowed us to measure the position of the X-ray source to an accuracy $``$0.06″. The results (see Table 1) show that the X-ray source is about 0.6″ east and 0.1″ north of the compact radio core position. The difference between the X-ray and radio measurements is well within the uncertainty in the absolute aspect solutions ($`<`$ 1″; CXC web page). Given that in this case only one X-ray source is significantly detected over an approximately 8′ $`\times `$ 50′ field of view, the probability of the source to fall, by chance, to within $``$0.6″ of the radio core is $`8\times 10^7`$. It is worth emphasizing that such a bright source would have been detectable anywhere inside the field of view. Therefore, we conclude that it is highly unlikely that 1E1740.7-2942 was mis-identified at the radio wavelength.

The high-resolution HETGS spectrum of 1E1740.7-2942 shows evidence for emission lines (but only at roughly $`2\sigma `$ levels). The lines appear to be intrinsically narrow and weak, so they could have escaped the detections by previous X-ray spectrometers of inferior energy resolution (e.g., the GIS and SIS aboard ASCA; Sheth et al. 1996 and Churazov et al. 1996). The lines could be identified with H-like or He-like ions of Fe, Ca, S, and Ar. This would imply that the line-emitting material is highly ionized, which is perhaps not so surprising for an X-ray binary. The lines might originate from the accretion disk or the base of the jets where material moves relatively slowly. There is also evidence for the emission from the jets themselves, which is manifested in the pairs of Doppler red- and blue-shifted lines, although it is a bit speculative at present given the quality of the data. It is interesting to note, however, that every pair of lines would indicate approximately the same amount of Doppler shifts and that the inferred velocity and inclination angle of the twin jets are not unreasonable for this system (see § 3.3; compared to SS 433). A definitive confirmation or rejection of the proposed scenario awaits data of much improved S/N from future Chandra observations.

The faint, spherically symmetric halo emission is likely due to the scattering of X-rays from 1E1740.7-2942 by interstellar dust. Such scattering halos are known to exist around many X-ray sources and have proven useful for studying the spatial and size distributions of the dust grains in our Galaxy (e.g., Mauche & Gorenstein 1986; Predehl & Schmitt 1995). Chandra now makes it possible, for the first time, to directly image the phenomenon. The observations should, in principle, allow us to derive spatially resolved energy distribution of photons in the halo and thus to test theoretical models at a more detailed level. Unfortunately, our observation does not provide sufficient statistics to warrant such an investigation.

We are grateful to Paul Plucinsky and his colleagues at CXC for clarifying several critical issues regarding the integrity of the pipeline products for our observation. We also thank Nancy Evans and her colleagues at the CXC help desk for providing a definitive answer to our inquiry regarding the uncertainty in the absolute aspect solutions. This work was supported in part by NASA through contract NAS8-38252. NSS and SCT wish to acknowledge support from NASA through the Smithonian Astrophysical Observatiory contract SAO SV1-61010, and IFM support from Conicet/Argentina.


format: text

text:
# Replica-exchange multicanonical algorithm and multicanonical replica-exchange method for simulating systems with rough energy landscape

## I Introduction

With the great advancement of computational science over the past decades, partly due to the rapid growth of computer technologies, scientists in many research fields are now attacking realistic problems with large degrees of freedom and complexity. When dealing with systems with rough energy landscape, simulations by conventional molecular dynamics (MD) or Monte Carlo (MC) methods are of little use, because at low temperatures they tend to get trapped in local-minimum energy states and sample only a very limited region of the configurational space. The development of new powerful simulation algorithms that can alleviate this difficulty is thus particularly important.

Multicanonical algorithm (MUCA)  may be one of the most popular methods that can sample a wide phase space (for a review, see, e.g., Ref. ). A simulation in multicanonical ensemble is based on a non-Boltzmann weight factor, which we refer to as the multicanonical weight factor, and a free random walk in potential energy space is realized. The random walk can overcome any energy barrier and allows the simulation to escape from local-minimum states. Monitoring the potential energy throughout a MUCA simulation, one can obtain not only the global-minimum energy state but also canonical-ensemble averages of any physical quantity as a function of temperature. The latter is accomplished by the single-histogram reweighting techniques . The multicanonical weight factor is usually determined by iterations of short trial simulations (the details of this process are described, for instance, in Refs. ). This iterative process can be very difficult and tedious for complex systems, and there have been attempts to accelerate the convergence of the iterative process . Despite these efforts, the determination of the multicanonical weight factor still remains to be non-trivial .

Another powerful algorithm is the replica-exchange method (REM) . (A similar method was independently developed earlier in Ref. . REM is also referred to as multiple Markov chain method  and parallel tempering .) In this method, a number of non-interacting copies of the original system (or replicas) at different temperatures are simulated independently and simultaneously by the conventional MD or MC methods. Every few steps, pairs of replicas are exchanged with a specified transition probability. This exchange process realizes a random walk in temperature space, which in turn induces a random walk in the energy space so that a wide configurational space can be sampled during the simulation. The multiple-histogram reweighting techniques  (an extension of which is also referred to as WHAM ) are used to calculate various canonical-ensemble averages as a function of temperature. In the replica-exchange method the weight factor is just the product of Boltzmann factors, and so it is essentially known. However, REM also has a weak point: It is difficult to study the first-order phase transitions, because the replica-exchange rates passing the first-order transition temperatures are greatly reduced.

In a previous work we have developed a molecular dynamics algorithm for REM  (see, also, Refs. ). We then developed a multidimensional REM  (see, also, Refs. ). In this Letter we try to combine the merits of multicanonical algorithm and replica-exchange method. In the first method, which we refer to as replica-exchange multicanonical algorithm (REMUCA), a short replica-exchange simulation is performed and the multicanonical weight factor is determined by the multiple-histogram reweighting techniques . (We remark that the multiple-histogram reweighting techniques are also used during the regular iterations of multicanonical weight factor determination in Refs. ). We then present a further extension of the new method, which we refer to as multicanonical replica-exchange method (MUCAREM), where a replica-exchange multicanonical simulation is performed with a small number of replicas. The effectiveness of these methods is tested with a short peptide, Met-enkephalin, in gas phase.

## II Methods

### A Replica-exchange multicanonical algorithm

In the multicanonical ensemble  each state is weighted by the multicanonical weight factor, $`W_{mu}(E)`$, so that a flat potential energy distribution, $`P_{mu}(E)`$, is obtained: $`P_{mu}(E)n(E)W_{mu}(E)\mathrm{const},`$ where $`E`$ is the potential energy and $`n(E)`$ is the density of states. We can then write

$$W_{mu}(E)e^{\beta _0E_{mu}(E;T_0)}=\frac{1}{n(E)},$$
(1)
where we have chosen an arbitrary reference temperature, $`T_0=1/k_B\beta _0`$ ($`k_B`$ is Boltzmann’s constant), and the “multicanonical potential energy” is defined by

$$E_{mu}(E;T_0)=k_BT_0\mathrm{ln}n(E).$$
(2)
Since the density of states of the system is usually unknown, the multicanonical weight factor has to be numerically determined by iterations of short preliminary runs .

Once the weight factor is given, multicanonical Monte Carlo simulation is performed, for instance, with the Metropolis criterion. The molecular dynamics algorithm in multicanonical ensemble also naturally follows from Eq. (1), in which the regular constant temperature molecular dynamics (with $`T=T_0`$) is performed by solving the following modified Newton equation:

$$\dot{𝒑}_i=\frac{E_{mu}(E;T_0)}{𝒒_i}=\frac{E_{mu}(E;T_0)}{E}𝒇_i,$$
(3)
where $`𝒒_i`$ and $`𝒑_i`$ are the coordinate vector and momentum vector of the $`i`$-th atom, and $`𝒇_i`$ is the regular force acting on the $`i`$-th atom. In this Letter the formulations of the new methods in terms of molecular dynamics algorithm are presented. The details of the Monte Carlo versions of the new algorithms have also been worked out and will be given elsewhere .

We now briefly review the original replica-exchange method (REM) - (see Ref.  for details). The system for REM consists of $`M`$ non-interacting copies (or, replicas) of the original system in the canonical ensemble at $`M`$ different temperatures $`T_m`$ ($`m=1,\mathrm{},M`$). Let $`X=\{\mathrm{},x_m^{[i]},\mathrm{}\}`$ stand for a state in this generalized ensemble. Here, the superscript $`i`$ and the subscript $`m`$ in $`x_m^{[i]}`$ label the replica and the temperature, respectively. The state $`X`$ is specified by the $`M`$ sets of coordinates $`q^{[i]}`$ (and momenta $`p^{[i]}`$).

A simulation of the replica-exchange method (REM) is then realized by alternately performing the following two steps -. Step 1: Each replica in canonical ensemble of the fixed temperature is simulated $`simultaneously`$ and $`independently`$ for a certain MC or MD steps. Step 2: A pair of replicas, say $`i`$ and $`j`$, which are at neighboring temperatures $`T_m`$ and $`T_n`$, respectively, are exchanged: $`X=\{\mathrm{},x_m^{[i]},\mathrm{},x_n^{[j]},\mathrm{}\}X^{}=\{\mathrm{},x_m^{[j]},\mathrm{},x_n^{[i]},\mathrm{}\}`$. The transition probability of this replica exchange is given by the Metropolis criterion:

$$w(XX^{})w\left(x_m^{[i]}|x_n^{[j]}\right)=\{\begin{array}{cc}1,\hfill & \mathrm{for}\mathrm{\Delta }0,\hfill \\ \mathrm{exp}\left(\mathrm{\Delta }\right),\hfill & \mathrm{for}\mathrm{\Delta }>0,\hfill \end{array}$$
(4)
where

$$\mathrm{\Delta }=\left(\beta _m\beta _n\right)\left(E\left(q^{[j]}\right)E\left(q^{[i]}\right)\right).$$
(5)
Here, $`E\left(q^{[i]}\right)`$ and $`E\left(q^{[j]}\right)`$ are the potential energy of the $`i`$-th replica and the $`j`$-th replica, respectively. In the present work we employ molecular dynamics algorithm for Step 1. We remark that in molecular dynamics simulations we also have to deal with the momenta. It was shown that by rescaling the momenta uniformly by the square root of the ratio of $`T_m`$ and $`T_n`$, the kinetic energy terms in the Boltzmann factors are cancelled out and that the same criterion, Eqs. (4) and (5), which was originally derived for Monte Carlo algorithm - is recovered .

The major advantage of REM over other generalized-ensemble methods such as multicanonical algorithm  lies in the fact that the weight factor is essentially a priori known, while in the latter algorithm the determination of the weight factor can be very tedious and time-consuming. A random walk in “temperature space” is realized for each replica, which in turn induces a random walk in potential energy space. This alleviates the problem of getting trapped in states of energy local minima.

In the replica-exchange multicanonical algorithm (REMUCA) we first perform a short REM simulation (with $`M`$ replicas) to determine the multicanonical weight factor and then perform with this weight factor a regular multicanonical simulation with high statistics. The first step is accomplished by the multiple-histogram reweighting techniques . Let $`N_m(E)`$ and $`n_m`$ be respectively the potential-energy histogram and the total number of samples obtained at temperature $`T_m=1/k_B\beta _m`$ of the REM run. The density of states, $`n(E)`$, is then given by

$$n(E)=\frac{{\displaystyle \underset{m=1}{\overset{M}{}}}g_m^1N_m(E)}{{\displaystyle \underset{m=1}{\overset{M}{}}}n_mg_m^1e^{f_m\beta _mE}},$$
(6)
and

$$e^{f_m}=\underset{E}{}n(E)e^{\beta _mE}.$$
(7)
Here, $`g_m=1+2\tau _m`$, and $`\tau _m`$ is the integrated autocorrelation time at temperature $`T_m`$. Note that the density of states $`n(E)`$ and the “dimensionless” Helmholtz free energy $`f_m`$ in Eqs. (6) and (7) are solved self-consistently by iteration .

Once the estimate of the density of states is obtained, the multicanonical weight factor can be directly determined by using Eq. (1) (see also Eq. (2)). Actually, the multicanonical potential energy, $`E_{mu}(E;T_0)`$, thus determined is only reliable in the following range:

$$E_1EE_M,$$
(8)
where $`T_1`$ and $`T_M`$ are respectively the lowest and the highest temperatures used in the REM run and

$$\{\begin{array}{cc}\hfill E_1& =<E>_{T_1},\hfill \\ \hfill E_M& =<E>_{T_M}.\hfill \end{array}$$
(9)
Here, $`<E>_T`$ stands for the canonical expectation value of $`E`$ at temperature $`T`$. Outside this range we extrapolate the multicanonical potential energy linearly as follows:

$$_{mu}^{\{0\}}(E)\{\begin{array}{cc}\frac{E_{mu}(E;T_0)}{E}|_{E=E_1}(EE_1)+E_{mu}(E_1;T_0),\hfill & \text{for }E<E_1\text{ ,}\hfill \\ E_{mu}(E;T_0),\hfill & \text{for }E_1EE_M\text{ ,}\hfill \\ \frac{E_{mu}(E;T_0)}{E}|_{E=E_M}(EE_M)+E_{mu}(E_M;T_0),\hfill & \text{for }E>E_M\text{ .}\hfill \end{array}$$
(10)
In the present work the multicanonical potential energy function, $`E_{mu}(E;T_0)`$, and its derivative, $`\frac{E_{mu}(E;T_0)}{E}`$, were obtained by fitting the logarithm of the density of states, $`\mathrm{ln}n(E)`$, (which was determined by the multiple-histogram reweighting techniques) by the cubic spline functions in the energy range of Eq. (8). A long multicanonical MD run is then performed by solving the Newton equations in Eq. (3) into which we substitute $`_{mu}^{\{0\}}(E)`$ of Eq. (10).

Let $`N_{mu}(E)`$ be the histogram of the potential energy distribution, $`P_{mu}(E)`$, that is obtained from this multicanonical production run. The expectation value of a physical quantity $`A`$ at any temperature $`T`$ is then calculated from

$$<A>_T=\frac{{\displaystyle \underset{E}{}}A(E)n(E)e^{\beta E}}{{\displaystyle \underset{E}{}}n(E)e^{\beta E}},$$
(11)
where the (unnormalized) density of states is obtained by the single-histogram reweighting techniques :

$$n(E)=\frac{N_{mu}(E)}{W_{mu}(E)}.$$
(12)
We remark that our multicanonical MD simulation here actually results in a canonical simulation at $`T=T_1`$ for $`E<E_1`$, a multicanonical simulation for $`E_1EE_M`$, and a canonical simulation at $`T=T_M`$ for $`E>E_M`$ (a detailed discussion on this point will be given elsewhere ). Note also that the above arguments are independent of the value of $`T_0`$, and we will get the same results, regardless of its value.

Finally, although we did not find any difficulty in the case of protein systems, a single REM run in general may not be able to give an accurate estimate of the density of states (like in the case of a first-order phase transition). In such a case we can still greatly simplify the process of the multicanonical weight factor determination by combining the present method with the previous iterative methods -.

### B Multicanonical replica-exchange method

In the previous subsection we presented a new generalized-ensemble algorithm, REMUCA, that uses the replica-exchange method for the determination of the multicanonical weight factor. Here, we present a further modification of REMUCA and refer the new method as multicanonical replica-exchange method (MUCAREM). In MUCAREM the final production run is not a regular multicanonical simulation (as in REMUCA) but a replica-exchange simulation with a few replicas, say $``$ replicas, in the multicanonical ensemble. We expect that this enhances the sampling further because a replica-exchange process can be considered to be a global update, while multicanonical simulations are usually based on local updates. We remark that replica-exchange simulations based on another generalized ensemble are introduced in Ref. .

We now give the details of MUCAREM. As in REMUCA, we first perform a short REM simulation with $`M`$ replicas with $`M`$ different temperatures (we order them as $`T_1<T_2<\mathrm{}<T_M`$) and obtain the best estimate of the density of states $`n(E)`$ in the whole energy range of interest (see Eq. (8)). We then choose a number $``$ ($`M`$) and assign $``$ pairs of temperatures ($`T_L^{\{m\}},T_H^{\{m\}}`$) ($`m=1,\mathrm{},`$). Here, we assume that $`T_L^{\{m\}}<T_H^{\{m\}}`$ and arrange the temperatures so that the neighboring regions covered by the pairs have sufficient overlaps. We also set $`T_L^{\{1\}}=T_1`$ and $`T_H^{\{\}}=T_M`$. We choose $``$ (arbitrary) temperatures $`T_m`$ ($`m=1,\mathrm{},`$). We then define the following quantities:

$$\{\begin{array}{cc}\hfill E_L^{\{m\}}& =<E>_{T_L^{\{m\}}},\hfill \\ \hfill E_H^{\{m\}}& =<E>_{T_H^{\{m\}}},(m=1,\mathrm{},)\hfill \end{array}$$
(13)
and assign the following multicanonical potential energies (see Eq. (10)):

$$_{mu}^{\{m\}}(E)=\{\begin{array}{cc}\frac{E_{mu}(E;T_m)}{E}|_{E=E_L^{\{m\}}}(EE_L^{\{m\}})+E_{mu}(E_L^{\{m\}};T_m),\hfill & \text{for }E<E_L^{\{m\}}\text{,}\hfill \\ E_{mu}(E;T_m),\hfill & \text{for }E_L^{\{m\}}EE_H^{\{m\}}\text{,}\hfill \\ \frac{E_{mu}(E;T_m)}{E}|_{E=E_H^{\{m\}}}(EE_H^{\{m\}})+E_{mu}(E_H^{\{m\}};T_m),\hfill & \text{for }E>E_H^{\{m\}}\text{,}\hfill \end{array}$$
(14)
where $`E_{mu}(E;T)`$ is the multicanonical potential energy that was determined for the whole energy range of Eq. (8). The Newton equation in Eq. (3) again implies that our choice of $`_{mu}^{\{m\}}(E)`$ in Eq. (14) results in a canonical simulation at $`T=T_L^{\{m\}}`$ for $`E<E_L^{\{m\}}`$, a multicanonical simulation for $`E_L^{\{m\}}EE_H^{\{m\}}`$, and a canonical simulation at $`T=T_H^{\{m\}}`$ for $`E>E_H^{\{m\}}`$ .

The production run of MUCAREM is a replica-exchange simulation with $``$ replicas with $``$ different temperatures $`T_m`$ and multicanonical potential energies $`_{mu}^{\{m\}}(E)`$. By following the same derivation that led to the original REM, one can show that the transition probability of replica exchange, $`w\left(x_m^{[i]}|x_n^{[j]}\right)`$, is given by Eq. (4), where we now have

$$\mathrm{\Delta }=\beta _m\left(_{mu}^{\{m\}}\left(E\left(q^{[j]}\right)\right)_{mu}^{\{m\}}\left(E\left(q^{[i]}\right)\right)\right)\beta _n\left(_{mu}^{\{n\}}\left(E\left(q^{[j]}\right)\right)_{mu}^{\{n\}}\left(E\left(q^{[i]}\right)\right)\right).$$
(15)
Note that we need to newly evaluate the multicanonical potential energy, $`_{mu}^{\{m\}}(E(q^{[j]}))`$ and $`_{mu}^{\{n\}}(E(q^{[i]}))`$, because $`_{mu}^{\{m\}}(E)`$ and $`_{mu}^{\{n\}}(E)`$ are, in general, different functions for $`mn`$. We remark that the same additional evaluation of the potential energy is necessary for the multidimensional replica-exchange method .

Let $`N_m(E)`$ and $`n_m`$ be respectively the potential-energy histogram and the total number of samples obtained at $`T_m`$ ($`m=1,\mathrm{},`$). The expectation value of a physical quantity $`A`$ at any temperature $`T=1/k_B\beta `$, $`<A>_T`$, is then given by Eq. (11), where the density of states is now obtained from the multiple-histogram reweighting formulae  of Eqs. (6) and (7) with the replacements: $`M`$ and $`\beta _mE\beta _m_{mu}^{\{m\}}(E)`$.

## III Results and discussion

We tested the effectiveness of the new algorithms for the system of a penta-peptide, Met-enkephalin, in gas phase. This peptide has the amino-acid sequence, Tyr-Gly-Gly-Phe-Met. The N and C termini of the peptide were blocked with acetyl and N-methyl groups, respectively. The force-field parameters were taken from the all-atom version of AMBER developed by Cornell et al. , and the dielectric constant was set equal to 1.0.

We have performed MD simulations in the Cartesian coordinates based on the replica-exchange multicanonical algorithm (REMUCA) and multicanonical replica-exchange method (MUCAREM). As a thermostat for both MD simulations in the replica-exchange method (REM) and multicanonical algorithm (MUCA), we have adopted the constraint method  in which the total kinetic energy is constant, following Ref.  where the constraint method is recommended for MUCA simulations. The computer code developed in Refs. , which is based on PRESTO , was used. The unit time step was set to 0.5 fs. The simulations (of all replicas) were started from an extended conformation.

In Table I we summarize the parameters of the simulations that were performed in the present work. As discussed in the previous section, REMUCA consists of two simulations: a short REM simulation (from which the density of states of the system, or the multicanonical weight factor, is determined) and a subsequent production run of MUCA simulation. The former simulation is referred to as REM1 and the latter as MUCA1 in Table I. A production run of MUCAREM simulation is referred to as MUCAREM1 in Table I, and it uses the same density of states that was obtained from REM1. Finally, a production run of the original REM simulation was also performed for comparison and it is referred to as REM2 in Table I. The total simulation time for the three production runs (REM2, MUCA1, and MUCAREM1) was all set equal (i.e., 5 ns). Before taking the data, we made regular canonical MD simulations of 100 ps (for each replica) for thermalization. For replica-exchange simulations an additional REM simulation of 100 ps was made for further thermalization. In REM1 and REM2 there exist 10 replicas with 10 different temperatures, ranging from 200 K to 1000 K as listed in Table I (i.e., $`T_1=200`$ K and $`T_M=T_{10}=1000`$ K). For the present purpose of the work, we believe that this range of temperature is sufficiently wide. The temperatures are distributed exponentially between $`T_1`$ and $`T_M`$, following the optimal distribution found in Ref . In REM1 and REM2 a replica exchange was tried every 20 time steps (or 10 fs) as in our previous work .

In REM1 each replica was simulated by a constant temperature MD of $`2\times 10^5`$ steps (or 100 ps). We first check whether the replica-exchange simulation of REM1 indeed performed properly. For an optimal performance of REM the acceptance ratios of replica exchange should be sufficiently uniform and large (say, $`>10`$ %). In Table II we list these quantities. It is clear that both points are met (the values vary only between 13 % and 16 %). Moreover, in Fig. 1 the canonical probability distributions obtained at the chosen 10 temperatures from REM1 are shown. We see that there are enough overlaps between all pairs of neighboring distributions, indicating that there will be sufficient numbers of replica exchanges between pairs of replicas (see Table II). We did observe random walks in potential energy space between low energies and high energies.

After REM1, we obtained the density of states, $`n(E)`$, by the multiple-histogram method  (namely, by solving Eqs. (6) and (7) self-consistently). For biomolecular systems the integrated autocorrelation times, $`\tau _m`$, in the reweighting formulae can safely be set to be a constant , and we do so throughout the analyses in this section. The density of states will give the average values of the potential energy from Eq. (11), and we found

$$\{\begin{array}{cc}\hfill E_1& =<E>_{T_1}=30\mathrm{kcal}/\mathrm{mol},\hfill \\ \hfill E_M& =<E>_{T_M}=195\mathrm{kcal}/\mathrm{mol}.\hfill \end{array}$$
(16)
Then our estimate of the density of states is reliable in the range $`E_1EE_M`$. The multicanonical potential energy $`_{mu}^{\{0\}}(E)`$ was thus determined for the three energy regions ($`E<E_1`$, $`E_1EE_M`$, and $`E>E_M`$) from Eq. (10). Namely, the multicanonical potential energy, $`E_{mu}(E;T_0)`$, in Eq. (2) and its derivative, $`\frac{E_{mu}(E;T_0)}{E}`$, were determined by fitting $`\mathrm{ln}n(E)`$ by cubic spline functions in the energy region of ($`30E195`$ kcal/mol). Here, we have set the arbitrary reference temperature to be $`T_0=1000`$ K. Outside this energy region, $`E_{mu}(E;T_0)`$ was linearly extrapolated as in Eq. (10).

After determining the multicanonical weight factor, we carried out a multicanonical MD simulation of $`1\times 10^7`$ steps (or 5 ns) for data collection (MUCA1 in Table I). In Fig. 2 the probability distribution obtained by MUCA1 is plotted. It can be seen that a good flat distribution is obtained in the energy region $`E_1EE_M`$. In Fig. 2 the canonical probability distributions that were obtained by the reweighting techniques at $`T=T_1=200`$ K and $`T=T_M=1000`$ K are also shown (these results are essentially identical to one another among MUCA1, MUCAREM1, and REM2, as discussed below). Comparing these curves with those of MUCA1 in the energy regions $`E<E_1`$ and $`E>E_M`$ in Fig. 2, we confirm our claim in the previous section that MUCA1 gives canonical distributions at $`T=T_1`$ for $`E<E_1`$ and at $`T=T_M`$ for $`E>E_M`$, whereas it gives a multicanonical distribution for $`E_1EE_M`$.

In the previous works of multicanonical simulations of Met-enkephalin in gas phase (see, for instance, Refs. ), at least several iterations of trial simulations were required for the multicanonical weight determination. We emphasize that in the present case of REMUCA (REM1), only one simulation was necessary to determine the optimal multicanonical weight factor that can cover the energy region corresponding to temperatures between 200 K and 1000 K.

From the density of states obtained by REMUCA (i.e., REM1), we prepared the multicanonical weight factors (or the multicanonical potential energies) for the MUCAREM simulation (see Eq. (14)). The parameters of MUCAREM1, such as energy bounds $`E_L^{\{m\}}`$ and $`E_H^{\{m\}}`$ ($`m=1,\mathrm{},`$) are listed in Table III. The choices of $`T_L^{\{m\}}`$ and $`T_H^{\{m\}}`$ are, in general, arbitrary, but significant overlaps between the probability distributions of adjacent replicas are necessary. The replica-exchange process in MUCAREM1 was tried every 200 time steps (or 100 fs). It is less frequent than in REM1 (or REM2). This is because we wanted to ensure a sufficient time for system relaxation.

In Fig. 3 the probability distributions of potential energy obtained by MUCAREM1 are shown. As expected, we observe that the probability distributions corresponding to the temperature $`T_m`$ are essentially flat for the energy region $`E_L^{\{m\}}EE_H^{\{m\}}`$, are of the canonical simulation at $`T=T_L^{\{m\}}`$ for $`E<E_L^{\{m\}}`$, and are of the canonical simulation at $`T=T_H^{\{m\}}`$ for $`E>E_H^{\{m\}}`$ ($`m=1,\mathrm{},`$). As a result, each distribution in MUCAREM is much broader than those in the conventional REM (see Fig. 1) and a much smaller number of replicas are required in MUCAREM than in REM ($`=4`$ in MUCAREM versus $`M=10`$ in REM).

The acceptance probabilities of replica exchange in MUCAREM1 are listed in Table IV. Relatively high probabilities of replica exchange in MUCAREM1 should be noted. In MUCAREM1, the exchange probability is given by Eqs. (4) and (15). Suppose that a pair of replicas $`i`$ (at temperature $`T_m`$) and $`j`$ (at temperature $`T_n`$) are respectively in the energy ranges between $`E_L^{\{m\}}EE_H^{\{m\}}`$ and between $`E_L^{\{n\}}EE_H^{\{n\}}`$. Then Eq. (15) together with Eqs. (2) and (14) gives

$`\mathrm{\Delta }`$ $`=`$ $`\beta _m\left(E_{mu}(E\left(q^{[j]}\right);T_m)E_{mu}(E\left(q^{[i]}\right);T_m)\right)\beta _n\left(E_{mu}(E\left(q^{[j]}\right);T_n)E_{mu}(E\left(q^{[i]}\right);T_n)\right)`$ (17)
$`=`$ $`\mathrm{ln}n\left(E\left(q^{[j]}\right)\right)\mathrm{ln}n\left(E\left(q^{[i]}\right)\right)\mathrm{ln}n\left(E\left(q^{[j]}\right)\right)+\mathrm{ln}n\left(E\left(q^{[i]}\right)\right)`$ (18)
$`=`$ $`0.`$ (19)
Thus, in this case the replica exchange is always accepted. This is the main reason for the high acceptance ratio in MUCAREM simulations.

In Fig. 4 the time series of potential energy for the first 500 ps of REM2 (a), MUCA1 (b), and MUCAREM1 (c) are plotted. Note that in MUCA1 (Fig. 4(b)) some extra MD steps are required to escape from the lowest-energy region (from about 300 ps to 450 ps) and the transition from the highest to the lowest energy occurs only once during the first 500 ps in MUCA1. This is presumably due to the nature of local updates. In REM2 and MUCAREM1, such transitions occur twice or three times. This is because global updates by the replica-exchange process make the random walk more efficient (more quantitative comparisons of REM, MUCA, and MUCAREM will be given elsewhere ).

To check the validity of the canonical-ensemble expectation values calculated by the new algorithms, we compare the average potential energy as a function of temperature in Fig. 5. In REM2 and MUCAREM1 we used the multiple-histogram techniques  (see Eqs. (6) and (7) for REM2 and the corresponding variants for MUCAREM1), whereas the single-histogram method  (see Eq. (12)) was used in MUCA1. We can see a perfect coincidence of these quantities among REM2, MUCA1, and MUCAREM1 in Fig. 5.

## IV Conclusions

In this Letter we have proposed two new algorithms for configurational sampling of systems with rough energy landscape. In the first method, which we refer to as replica-exchange multicanonical algorithm (REMUCA), the multicanonical weight factor is determined from the results of a short replica-exchange simulation, and then a regular multicanonical production run is made with this weight. In the second method, which we refer to as multicanonical replica-exchange method (MUCAREM), the multicanonical weight factor is determined as in REMUCA, and then a replica-exchange multicanonical production run is made. The number of required replicas in MUCAREM is much smaller than in the original replica-exchange method. The new methods were tested with the system of a small peptide in gas phase. This is one of the smallest models of protein systems. The multicanonical weight factor was indeed obtained by a single, short replica-exchange simulation. It may occur, however, that this first step may not be as easy in more complicated systems such as a protein with explicit solvent molecules. In such a case one can still simplify the multicanonical weight factor determination by iteratively combining the present algorithms (and also the previous methods).

Acknowledgements

The authors are grateful to Dr. A. Kitao of Kyoto University for useful discussions on fitting procedures in the determination of the multicanonical weight factor. Our simulations were performed on the Hitachi and other computers at the IMS Computer Center. Part of our simulations were also performed on the supercomputers at the National Institute of Genetics. This work is supported, in part, by a grant from the Research for the Future Program of the Japan Society for the Promotion of Science (JSPS-RFTF98P01101).

Figure Captions

* Fig. 1. Probability distribution of potential energy obtained from REM1 (see Table I for the parameters of the simulation).
* Fig. 2. Probability distribution of potential energy obtained from MUCA1 (see Table I). The dotted curves are the probability distributions of the reweighted canonical ensemble at $`T=200`$ K (left) and 1000 K (right).
* Fig. 3. Probability distributions of potential energy obtained from MUCAREM1 (see Tables I and III).
* Fig. 4. Time series of potential energy for one of the replicas in (a) REM2, (b) MUCA1, and (c) MUCAREM1 (see Tables I and III for the parameters of the simulation).
* Fig. 5. The average potential energy as a function of temperature. The solid, dotted, and dashed curves are obtained from REM2, MUCA1, and MUCAREM1, respectively (see Tables I and III for the parameters of the simulations).


text:
# The Galaxy Environment of Quasars in the 𝐳≃1.3 Clowes-Campusano Large Quasar Group

## 1. Introduction

Quasars have been used as efficient probes of high-redshift clustering because they are known to favour rich environments. Quasars may also trace large-scale structures at early epochs ($`0.4<z<2`$) in the form of Large Quasar Groups (LQGs), which have comparable sizes to the largest structures at the present epoch. The largest of these is the Clowes-Campusano LQG of 18 quasars at $`z1.3`$, with a maximal extent of 200$`h^1`$Mpc (Clowes & Campusano 1991).

To examine the galaxy environment of quasars in this LQG, we have conducted an ultra-deep optical (to $`V27,I26`$) study of a $`30^{}\times 30^{}`$ field containing 3 quasars from this LQG, with additional $`K`$ imaging obtained for subfields around 2 of them. Any associated passively-evolving galaxies can be identified from their extremely-red optical-IR colours ($`IK4,VK7`$).

## 2. Galaxy Clustering Around a z=1.226 LQG Quasar

We find a 3.5$`\sigma `$ excess of $`K<20`$ galaxies in the $`2.25^{}\times 2.25^{}`$ ($`1.3\times 1.3h^2\mathrm{Mpc}^2`$) field centred on the z=1.226 LQG quasar 104420.8+055739 (Haines et al. 2000). This excess is due entirely to a factor $``$11 overdensity of $`IK>3.75`$ galaxies, which must have $`z>0.8`$ to explain their colour. We observe clear red sequences at $`VK6.9,`$ and $`IK4.3`$ (Figure 1), comparable to red sequences observed in other $`z1.2`$ clusters, and indicating a population of 15–18 passively-evolving massive ellipticals at the quasar redshift. We have produced a density map (Figure 2) of galaxies with the same optical colours and magnitudes as the red sequence galaxies ($`VI>2.25,I<23`$) to estimate the full extent of any associated clustering beyond the $`K`$ image (shown as box). There is clearly significant substructure, with the quasar located between two groups of red sequence galaxies and further clustering can be seen extending over 2–3$`h^1`$Mpc. The whole structure is suggestive of being a cluster in the early stages of formation through the coalescence of subclusters. The quasar is located within a concentration of blue ($`VI<1`$) galaxies, with $`>50\%`$ clearly at $`z>1`$, having $`IK>3.5`$, indicating a region of enhanced star-formation, similar to those found around several $`z1.1`$ LQG quasars (Hutchings, Crampton & Johnson 1995). The concentration appears extended to form a band which bisects the two groups of red sequence galaxies. We propose that this band of recent star-formation, and the quasar itself, have been triggered by the merging of the two groups, as the galaxies interact with the colliding ICMs.

The two other LQG quasars are in poor environments, although a background $`z=1.426`$ quasar lies on a cluster periphery. Quasars at these redshifts ($`1<z<1.5`$) are found either on the peripheries of clusters (see also Sanchéz & González-Serrano 1999) or in relatively undistinguished environments.

## References

Clowes, R. G., & Campusano, L. E. 1991, MNRAS, 249, 218

Haines, C. P., Clowes, R. G., Campusano, L. E., & Adamson, A. J. 2000,
MNRAS submitted

Hutchings, J. B., Crampton, D., & Johnson, A. 1995, AJ, 109, 73

Sanchéz, S. F., & González-Serrano, J.I. 1999, A&A, 352, 383


text:
# Classification of the Nuclear Multifragmentation Phase Transition

## I Introduction

The phenomenon of heavy-ion collision induced nuclear multifragmentation (NMF) has been extensively studied over the past decade. Experimentally it is found that nuclei with excitation energies of $`E^{}/A=`$ 5-10 MeV, i.e. in the region of the binding energy, expand forming a number of intermediate mass fragments Nebauer:1999 . Using subtle and advanced experimental techniques from a number of collision experiments with different nuclei and at different collision energies an experimental caloric curve has been derivedPochodzalla:1995 ; Ma:1997 , which supports the interpretation that NMF is the small particle number onset of the nuclear liquid-gas phase transition. Theoretically NMF has been described by a number of different approaches ranging from simple percolation modelsCampi:1988 ; Bauer:1995 , dynamical models with different levels of sophistication, e.g. quantum molecular dynamics of different flavors Aichelin:1986 ; Hartnack:1989 ; Maruyama:1992 ; Bondorf:1995b ; Peilert:1992  and fermion molecular dynamics Feldmeier:2000 , and a large variety of statistical models Bondorf:1995 ; Gross:1997 ; Parvan:2000 ; Bhattacharyya:1999 ; Elliott:2000 ; Chase:1995 ; Lee:1992 ; DeAngelis:1989 . Among the statistical models two mainstream-theories, namely the statistical multifragmentation model (SMM) Bondorf:1995  using the heat bath or canonical ensemble and the micro-canonical multifragmentation model (MMM) Gross:1997  using the constant energy or micro-canonical ensemble should be mentioned. For different ensembles the shape of the caloric curve may differ significantly. However, because all ensembles are uniquely connected by simple integral transforms the qualitative features of the phase transition should be the same in all ensembles. We therefore omit at this point the (surely) interesting question which ensemble is more appropriate for the description of nuclear multifragmentation.

In this paper we use the canonical statistical multifragmentation model based on the works of Mekjian et al. DeAngelis:1989 ; Mekjian:1990 ; Lee:1992 ; Lee:1993 ; Mekjian:1991a ; Chase:1994 ; Chase:1995  to analyze the order of the multifragmentation phase transition using a recently proposed classification scheme for phase transitions in finite systems Borr99b ; Muelken:2000 . The classification scheme is based on an analysis of the distribution of zeros of the canonical partition function in the complex plane and is in the thermodynamic limit in complete accordance with the Ehrenfest classification scheme. This scheme is shortly described in Sec. II. In Sec. III we briefly review the statistical model used in this paper. The results of our analysis which clearly show that NMF has the signature of a first order phase transition are presented in Sec. IV.

## II classification scheme

The classification scheme used below and described in far more detail in Borr99b ; Muelken:2000  is based on the analytic continuation of the inverse temperature $`\beta `$ into the complex plane $`\beta =\beta +i\tau `$. For macroscopic systems an equivalent scheme was developed by Grossmann et al. Gross1967 . All thermodynamic information about a system can be extracted from the distribution of complex zeros of the canonical partition function $`Z()=dE\mathrm{\Omega }(E)\mathrm{exp}(E)`$.

To avoid some difficulties with the high temperature limit we write the partition function as a product $`Z()=Z_{\mathrm{lim}}()Z_{\mathrm{int}}()`$, explicitly separating the high temperature limit ($`T\mathrm{}`$) of the partition function $`Z_{\mathrm{lim}}()`$ imposing $`lim_T\mathrm{}Z_{\mathrm{int}}()=1`$. The partition function is an integral function with complex conjugate zeros $`_k=_k^{}=\beta _k+i\tau _k(k)`$ and can be rewritten as

$`Z()`$ $`=`$ $`Z_{\mathrm{lim}}()Z_{\mathrm{int}}(0)\mathrm{exp}(_{}\mathrm{ln}Z_{\mathrm{int}}(0))`$
$`\times `$ $`{\displaystyle \underset{k}{}}\left(1{\displaystyle \frac{}{_k}}\right)\left(1{\displaystyle \frac{}{_k^{}}}\right)\mathrm{exp}\left({\displaystyle \frac{}{_k}}+{\displaystyle \frac{}{_k^{}}}\right).`$
All thermodynamic properties like the internal energy or specific heat are derivatives of the logarithm of the canonical partition function. Thus, the zeros of $`Z()`$ are poles of these quantities. In the thermodynamic limit different regions of holomorphy separate different phases by dense lines of zeros. In finite systems the zeros do not become dense on lines which leads to a less sharp separation of different phases. We interpret the zeros as boundary posts between two phases.

The distribution of zeros close to the real axis can approximately be described by three parameters, where two of them reflect the order of the phase transition and the third merely the size of the system. We assume that the zeros lie on straight lines with a discrete density of zeros given by

$$\varphi (\tau _k)=\frac{1}{2}\left(\frac{1}{|_k_{k1}|}+\frac{1}{|_{k+1}_k|}\right),$$
(2)
with $`k=2,3,4,\mathrm{}`$, and approximate for small $`\tau `$ the density of zeros by a simple power law $`\varphi (\tau )\tau ^\alpha `$. Considering only the first three zeros the exponent $`\alpha `$ can be estimated as

$$\alpha =\frac{\mathrm{ln}\varphi (\tau _3)\mathrm{ln}\varphi (\tau _2)}{\mathrm{ln}\tau _3\mathrm{ln}\tau _2}.$$
(3)
The second parameter describing the distribution of zeros is given by $`\gamma =\mathrm{tan}\nu (\beta _2\beta _1)/(\tau _2\tau _1)`$ where $`\nu `$ is the crossing angle of the line of zeros with the real axis (see Fig. 1). The discreteness of the system is reflected in the imaginary part $`\tau _1`$ of the zero closest to the real axis. For macroscopic systems we always have $`\tau _10`$. In this case the parameters $`\alpha `$ and $`\gamma `$ match those defined by Grossmann et al. Gross1967 , who showed the connection between certain values of $`\alpha `$ and $`\gamma `$ and different types of phase transitions. They claimed that a first order transition occurs for $`\alpha =\gamma =0`$, for $`0<\alpha <1`$ with $`\gamma =0`$ or $`\gamma 0`$ a second order transition, and for $`\alpha >1`$ a higher order transition. In the thermodynamic limit $`\alpha `$ is restricted to values greater than or equal to zero because otherwise this would cause a diverging internal energy. For finite systems values of $`\alpha `$ less than zero are possible.

In our classification scheme phase transitions in finite systems are of first order if $`\alpha 0`$, while the definitions for higher order transitions coincide with those given by Grossmann. Additionally we define $`\tau _1`$ as an unambiguous parameter for the discreteness of the system. Since in small systems no thermodynamic properties diverge, the specification of the critical temperature is difficult and several definitions are possible which coincide only in the thermodynamic limit. We define the real part of the first zero $`\beta _1`$ as the critical temperature. Another possible definition would be the crossing point of the line of zeros with the real temperature axis.

## III statistical multifragmentation model

In the statistical models for NMF Bondorf:1995 ; Chase:1994 ; Chase:1995 ; Mekjian:1991a  it is assumed that the fragmentation process can be described in the following manner. After initial excitation the system expands and fragments into a number of clusters. During the expansion a thermalization of all clusters is assumed. The thermalization abruptly stops when the mean separation between the clusters exceeds the range of the nuclear force. The properties of the systems at this so-called freeze-out density, which is assumed to be 0.1-0.3 of the normal nuclear density $`\rho _0=0.15`$ fm-3, are than carried by all fragments until experimental detection.

The probability for a fragmentation with fragmentation vector $`\stackrel{}{n}=(n_1,n_2,\mathrm{})`$ ($`n_k`$ is the number of clusters of size $`k`$ and $`kn_k=A`$), reads

$$W(\stackrel{}{n})=\underset{k1}{}\frac{1}{n_k!}\left(\frac{xy^{k1}}{\alpha _k}\right)^{n_k},$$
(4)
where $`x=V_{\mathrm{ac}}\left(\frac{M}{2\pi \mathrm{}^2}\right)^{3/2}`$ is the partition function for the ideal monoatomic gas of mass $`M`$ in the accessible volume $`V_{\mathrm{ac}}`$. $`y=\mathrm{exp}\left(a_v()+\frac{T_0/ϵ_0}{1+T_0}\right)`$ is related to the cluster binding and internal excitations. The cluster size dependence of the probability is given by $`\alpha _k=k^{5/2}`$. We assume an accessible volume three times larger than the normal nucleus volume $`V_0`$. The canonical partition function is then given as

$$Z_A=\underset{m=0}{\overset{A}{}}Z_A^{(m)}x^my^{Am},$$
(5)
where $`m=_kn_k`$ is the multiplicity, i.e. the total number of clusters, and $`Z_A^{(m)}`$ is independent of the thermodynamic variables and given by the simple recursion

$$Z_A^{(m)}=\frac{1}{m}\underset{k=1}{\overset{A}{}}\frac{1}{\alpha _k}Z_{Ak}^{(m1)},$$
(6)
with $`Z_A^{(1)}=1/\alpha _A`$.

Following Ref. Bondorf:1995  the binding energy for a composite of $`k`$ nucleons is given by

$$a_v()=W_0k\sigma ()k^{2/3},$$
(7)
with $`\sigma ()=\sigma _0\left((T_c1/^2)/(T_c+1/^2)\right)^{5/4}`$ being the surface tension ($`W_0=16`$MeV, $`\sigma _0=18`$MeV, and $`T_c=18`$MeV).

Within this framework different moments of the multiplicity are calculated by

$$m^i=\frac{1}{Z_A}\underset{m=0}{\overset{A}{}}m^iZ_A^{(m)}x^my^{Am}.$$
(8)
The variance of the multiplicity is given by $`\mathrm{\Delta }m=m^2m^2`$. It is straightforward to express thermodynamic quantities like the internal energy or the specific heat in terms of different moments of the multiplicity Chase:1995 .

## IV results

Due to the unfavorable scaling behavior of the partition function it is very arduous to find its zeros directly with numerical methods. The zeros of $`Z`$ are the poles of almost all thermodynamic expectation values, e.g. the internal energy, the specific heat and so on. In this paper we utilize the absolute value of the variance of the multiplicity $`|\mathrm{\Delta }m|`$ to locate the zeros of $`Z`$ in the complex temperature plane. Fig. 2 shows a contour plot of $`|\mathrm{\Delta }m|/A^2`$ in the complex temperature plane for nucleon-numbers $`A=25,75,150`$, and $`225`$ (note the different scaling of the $`\tau `$-axis.). This figure directly shows that multifragmentation is a first order phase transition. The zeros lie on straight lines perpendicular to the real temperature axis with constant spacings between them. Even at particle numbers as small as 25 the nature of the phase transition can be clearly identified from the lines of zeros of the partition function.

In contrast, it is very difficult to extract the order of the phase transition from the canonical caloric curve (internal energy Uint per particle number A vs. temperature T) shown in Fig. 3.

The recursion formula (6) has a very similar structure to that for Bose-Einstein systems bf93 ; Borr99a ; Muelken:2000 . This is because both systems have an intimate relationship to the permutation group Mekjian:1991a . It is therefore remarkable, that the multifragmentation phase transition for nuclear matter is of first order while Bose-Einstein condensation is a second or third order phase transition depending on the single particle density of states.

A detailed analysis of the distribution of zeros reveals that for all nucleon numbers ($`A>25`$) the most important classification parameter $`\alpha `$ is almost equal to zero (see Fig. 4(a)). The slight deviations for A=50 and A=150 are most likely due to uncertainties within our numerical zero detection method. The second parameter $`\gamma `$ can directly be extracted from Fig. 2. Even for small nucleon numbers the distribution of zeros is perpendicular to the real axis, corresponding to $`\gamma =0`$. The finite size of the system is connected to the third classification parameter $`\tau _1`$. Fig. 4(b) shows the $`A`$-dependence of $`\tau _1`$. With increasing nucleon number $`\tau _1`$ decreases. We approximated the scaling behavior as $`\tau _1A^{0.96}`$, i.e. $`\tau _1`$ is simply volume dependent and the phase transition would approach a true first order phase transition in the Ehrenfest sense with $`1/A`$. The critical temperature $`T=1/\beta _1`$ (Fig. 4(c)) decreases and adopts at $`A200`$ the experimentally found value of about $`5`$ MeV. Of course, for larger system sizes the simple equation (7) for the binding energy used in this study should be refined and augmented by additional terms. Additional calculations for different particle densities reveal that there are no qualitative changes to the first order nature of the phase transition. Fig. 4(c) also shows calculations for accessible volumes which are twice, three, and five times as large as the normal nucleus volume $`V_0`$. There are only quantitative changes in the critical temperature.

It should be noted that one can expect that refinements of the model, e.g. the introduction of an additional Coulomb term and other values of the particle density, may influence the details of the distribution of zeros. Especially tuning the particle density to the expected supercritical phase is an interesting challenge, worth to be pursued in the future.

## V conclusion

In conclusion we have found that within the statistical model nuclear multifragmentation is undoubtedly a first order phase transition, which is in complete agreement with recent works based on micro-canonical and grand-canonical calculations Chomaz:2000 ; Bugaev:2000 . Of course, the first order type of the NMF phase transition is not yet sufficient to interprete NMF as the onset of the nuclear liquid-gas transition. However, this seems to be a necessary condition. The distribution of zeros of the partition function in the complex temperature plane proofs as a reliable tool for the classification of phase transitions, adding a clear, detailed, and unambiguous view on the thermodynamic properties of small systems.

###### Acknowledgements.

We wish to thank Jens Harting, Heinrich Stamerjohanns and E. R. Hilf for fruitful discussions.


text:
# DAMTP-2000-12 hep-ph/0009130 Problems with empirical bounds for strangelet production at RHIC

## I Introduction

### A Background

High energy collider experiments tend to explore new realms of physics, in which theory tends to be uncertain. Catastrophic phenomena clearly occur in nature, and with sufficient speculation one can imagine theoretical mechanisms by which catastrophe could ensue from a collider experiment. Such mechanisms, of course, tend to be discussed and investigated for plausibility whenever any new experiment is proposed.

The experiments considered here are relativistic heavy ion collisions scheduled to take place at Brookhaven (RHIC) and CERN (ALICE). Several speculative catastrophe mechanisms have been discussed. We examine here the one which seems at first sight the least implausible and which has attracted most attention, namely the possibility that, if negatively charged metastable strange matter exists, it could be produced in the experiments and go on to destroy the Earth in a runaway reaction.

This possibility was investigated by Busza et al. (BJSW), who were requested to review candidate mechanisms for catastrophe scenarios at RHIC following public concern. Their initial report concluded that the possibility of strangelet-induced catastrophe can be firmly excluded by compelling theoretical arguments, and that very large safety factors can also be derived from the survival of the Moon and other astrophysical evidence. This analysis partly relies on independent work by Dar et al. (DDH), who noted a theoretical hypothesis which would imply that the survival of the Moon would provide no useful risk bounds, but concluded that the observed rate of supernovae provides large safety factors even under this hypothesis.

Following criticisms, Busza et al. produced a revised version of their report. In this version, BJSW shift their ground considerably, explicitly disclaiming any attempt to decide what might be an acceptable bound on the risk of catastrophe, and acknowledging explicitly the possibility that all empirical bounds on that risk would be invalidated if the hypothetical strangelets had particular unfortunate properties. However, they underline their view that compelling theoretical arguments provide sufficient, albeit unquantifiable, reassurance that there will be no catastrophe.

### B Scope of this paper

Considering the possibility, however remote, of so large a catastrophe is particularly difficult, since normal scientific practice is in some important respects inadequate. In everyday research, a hypothesis which has perhaps at best a $`1\%`$ chance of being correct is often dismissed as not worth pursuing, and a scenario which relies on two or three such hypotheses is very unlikely to be taken seriously. Such judgements have to be made: science would progress much more slowly if disproportionate attention was paid to unlikely hypotheses and implausible explanations. However, when trying to exclude the possibility of global catastrophe, a $`10^2`$, $`10^4`$ or even $`10^6`$ probability of error is far from negligible, indeed alarmingly high. So every argument and every piece of conventional wisdom needs to be very carefully scrutinised in case it might have some neglected flaw.

A good way of trying to uncover such flaws, which I gather is often used in analysing risks of large disasters, is to give different roles to different teams of people. One team is required to make the case for safety. The other is required to try to find flaws in that case, to look out for the unexpected, to search for ways in which things could go wrong — in short, to play devil’s advocate. This paper is written in the spirit of a contribution from the second, contrarian, team.

The discussion here focusses on the empirically based arguments for bounds on the risk of catastrophe, which form an important part of the case made for safety at RHIC and have been widely publicized in an attempt to allay public fears. Its aim is to probe the robustness of those arguments, not to argue that the risk is in fact overly large.

I also consider a potential danger that could arise from the production of positively charged strangelets. While the risk of catastrophe arising this way appears pretty small — the strangelets would have to be somehow transported to the Sun, or introduced into an artificial fusion reaction — it deserves careful consideration.

This paper does not examine whether the reassurance given by theoretical arguments is sufficient in itself, as BJSW suggest. Nor does it consider what might be an acceptable bound on the risk of catastrophe. These questions will be addressed elsewhere.

### C Empirically derived risk bounds: general comments

Approximately $`2\times 10^{11}`$ high energy collisions of pairs of gold ions are scheduled to take place over the $`10`$ year lifetime of RHIC. Following BJSW and DDH, we consider the speculative “killer strangelet” scenario in which the experiment might lead to disaster.

The scenario requires the following hypotheses. Small atomic mass metastable negatively charged strangelets must exist. It must be possible to produce them in heavy ion collisions. They must be sufficiently stable to survive long enough after production to interact with matter in the surrounding environment. Then, strange matter must be stable in bulk, with a spectrum of stable strangelets such that a negatively charged metastable strangelet can begin a runaway reaction by fusing successively with positively charged nuclei, producing larger negatively charged strangelets, and eventually consuming the Earth.

Empirical bounds on the risk of such a catastrophe can be derived because the type of collisions planned at RHIC also occur naturally when high energy cosmic rays collide, either with each other or with heavy nuclei in celestial bodies. If the probability of dangerous strangelet production in such collisions is sufficiently large, natural catastrophes should be induced, and we should be able to observe the consequences. The fact that we see no evidence of strangelet-induced catastrophes means we can bound the probability of strangelet production per collision, and hence the overall probability of catastrophe at RHIC. As quantum processes are probabilistic, small bounds on the probability of catastrophe are the best we can hope for from this line of reasoning. Empirical bounds cannot show that the probability is zero.

Moreover, no bound can be established from empirical data alone. Any argument establishing a bound must also rely on some theoretical hypotheses. Hence, to establish a bound of the form $`p_{\mathrm{catastrophe}}<p_0`$, we need to have at least the same level of confidence in all the relevant hypotheses. It will immediately be seen that this is likely to pose a serious problem, since on the one hand no bound will be thought reassuring unless $`p_0`$ is very small indeed, and on the other the necessary hypotheses involve a physical regime where theory is fairly uncertain.

The failure to consider this problem carefully and systematically is one of the weakest and most troubling points of BJSW and DDH’s analyses of the empirical evidence. The discussion below highlights some gaps in their arguments.

## II Risk bounds derived from the survival of the Moon

BJSW first consider the value of the Moon’s survival as empirical evidence against the possibility of catastrophe at RHIC. They note that cosmic rays, including high energy heavy nuclei, are constantly bombarding the Moon. The cosmic ray flux is believed to have been at the present level or higher during the Moon’s lifetime, roughly $`5`$ billion years. Since the Moon, unlike the larger planets in the solar system, is unprotected by an atmosphere, these cosmic rays collide with its surface, which contains a significant component of heavy nuclei and, in particular, is iron-rich. To derive a risk bound from these observations requires a series of assumptions.

* We need to decide precisely which types of heavy nuclei collisions should be considered good models for those at RHIC (or other collider experiments). RHIC involves gold-gold collisions. BJSW suggest considering all nuclei with $`Z>70`$ as ‘gold’, that is, effectively equivalent to gold for the purposes of this analysis. They also suggest, as a less conservative alternative, taking iron (Z=26) as effectively equivalent to gold.

BJSW then produce estimates the rate of relevant collisions on the Moon by the following reasoning. The flux of iron nuclei in cosmic rays at the relevant energies is fairly reliably known, and the flux of ‘gold’ nuclei can be estimated by extrapolating from lower energy data. The iron abundance on the Moon can conservatively be estimated by taking the lowest value of the abundance measured at six Apollo spacecraft landing sites; BJSW round this value down from 4.2% to 4%.

Accepting BJSW’s estimates and extrapolations, and the equivalence of ‘gold’ to gold, then gives four possible bounds, derived by considering respectively iron-iron, iron-‘gold’, ‘gold’-iron and ‘gold’-‘gold’ collisions. BJSW consider the first and third of these. The fourth, which most closely parallels the RHIC experiments, turns out to give no useful bound at all. BJSW note this important caveat in their analysis, but neglect it when summarising their conclusions.
* Any strangelets produced with high speed in Moon centre-of-mass frame may break up while being slowed down by lunar matter; if so, they will not survive to begin a catastrophic runaway reaction. We thus need some assumption about the rapidity distribution of strangelets produced, and some assumption about how likely a strangelet of any given rapidity is to survive the slowing process and begin a runaway reaction. These assumptions should also apply to any strangelets produced at RHIC. They lead to estimates for “suppression factors” — the fraction of those strangelets produced which survive to begin a runaway reaction — in the two cases.

Here, BJSW and DDH take different approaches. DDH consider the worst case hypothesis: a rapidity distribution which implies that no strangelet produced on the Moon will survive to begin a runaway reaction, while all strangelets produced at RHIC will. BJSW propose the rapidity distribution

$$\frac{d\mathrm{\Pi }}{dy}=Npy^ae^{by},$$
(1)
where $`a,b`$ are free parameters, $`2p`$ is the strangelet production probability per collision, and $`N`$ is a normalisation constant. They model the survival probability as a function of rapidity, atomic number and charge as

$$P(y,A,Z)=\mathrm{exp}(4.85(1+1/3A^{1/3})^2(\mathrm{cosh}y1)A/Z^2).$$
(2)
This expression is derived on the assumption that the strangelet break-up cross-section is independent of energy.

BJSW proceed by calculating the various relative suppression factors for $`a=1,2,3,4`$, and conservatively using the largest of these factors (corresponding to $`a=4`$) in their calculations.

They conclude that, if iron-iron collisions at RHIC energies are a good model for RHIC collisions, the probability of producing a dangerous strangelet at RHIC energies can be bounded by $`p_{\mathrm{catastrophe}}<10^5`$. If only ‘gold’-iron (or ‘gold’-‘gold’) collisions at RHIC energies are taken as relevant, no bound can be derived.

BJSW also calculate bounds using collisions at the lower energies of the AGS experiment, which they believe are (slightly) more plausible for strangelet production. They obtain $`p_{\mathrm{catastrophe}}<2\times 10^{11}`$ if iron-iron is a good model, $`p_{\mathrm{catastrophe}}<2\times 10^6`$ if ‘gold’-iron is a good model, and no useful bound if only ‘gold’-‘gold’ collisions are relevant.

BJSW summarise the situation thus: “If, however, one insists on recreating exactly the circumstances at RHIC and insists on the worst case rapidity distribution, then lunar limits are not applicable.” They add, in the penultimate paragraph of Ref. , “The rapidity distribution necessary to wipe out lunar limits is bizarre.”

But no assumption about the rapidity distribution is necessary to wipe out the lunar limits. They would not apply if there were a significant difference between the dangerous strangelet production rates for collisions involving iron and those involving gold (or ‘gold’). If only ‘gold’-‘gold’ collisions contribute significantly to the rate of strangelets which are produced and survive slowing, then, even granting BJSW’s rapidity distribution, none of their bounds are valid.

Other questionable assumptions in BJSW’s analysis are whether the precise forms of their rapidity distribution or survival probability function are correct, and whether, if so, $`a`$ must necessarily lie in the range $`1a4`$. To give a convincing argument that $`p_{\mathrm{catastrophe}}`$ is less than $`<2\times 10^6`$ (let alone $`2\times 10^{11}`$), one would need to argue that we can be all but certain that each of these assumptions is true. BJSW offer no such argument.

## III Risk bounds derived from supernova observations

Dar, De Rujula and Heinz’s analysis of the killer strangelet scenario starts by noting that the bounds based on lunar survival may contain a potential loophole. They suggest that the rapidity distribution eq. (1) might be replaced by a distribution which approximates a delta function:

$$\frac{d\mathrm{\Pi }}{dy}=p\delta (yY/2),$$
(3)
$`Y`$ being the total rapidity interval. If so, no strangelets produced on the Moon by cosmic ray collisions would survive slowing and cause a catastrophic runaway reaction, but strangelets produced at RHIC could. BJSW, while finding this hypothesis impossible to justify theoretically, review DDH’s analysis and reconsider their proposed bounds.

As noted above, one need not adopt DDH’s rapidity distribution (3) in order to challenge the lunar bounds, which are open to question for a variety of reasons. So their independent analysis, which claims “to derive a fool-proof and stringent limit on the potential danger of the BNL experiments”, needs and deserves careful examination.

DDH’s idea is the following. High energy collisions between heavy nuclei cosmic rays occur in space. If these collisions are capable of producing killer strangelets, some of the strangelets produced in the past will have survived slowing interactions with ambient hydrogen and eventually will have been swept up into protostellar material. If killer strangelets are capable of destroying the Earth catastrophically, they will be equally capable of causing catastrophic supernova-like explosions, which we should see. Even if killer strangelets are capable of causing only a very slow destruction of the Earth, over tens or hundreds of millions of years, they would add significantly to the luminosity of stars. So, the observed rate of supernovae, and the observed luminosities, should allow us to derive bounds on the probability of strangelet-induced catastrophe at RHIC.

Though DDH claim their limit is fool-proof, it actually relies on some important assumptions. First, they assume that killer strangelets produced with relative speed $`v<v_{\mathrm{crit}}=0.1c`$ with respect to ambient hydrogen will survive slowing collisions until it reaches the local mean velocity. Second, they assume that killer strangelets will then survive (growing as they occasionally hit and absorb hydrogen nuclei) during the long time interval before they form part of a proto-star.

The first of these assumptions is justified by assuming that the strangelet will have binding energy per nucleon similar to that of a typical nucleus, and will behave similarly in collisions with other nuclei. At first sight, neither claim seems beyond reasonable doubt. Yet the assumption is crucial in deriving DDH’s bounds, since the rate of dangerous strangelet production depends strongly on the value assumed for $`v_{\mathrm{crit}}`$ (cf. equation (3) of Ref. ).

This seems to be the weakest point in DDH’s argument. However, their second assumption could also, in principle, be invalid, for at least three reasons.

First, killer strangelets could be metastable rather than stable. If their lifetime were somewhere between, say, $`10^6`$ seconds and $`10^6`$ years, strangelets produced at RHIC would still cause catastrophes, while almost no strangelets produced in space would survive long enough to produce observable effects.

Second, it could be that the spectrum of stable strangelets is such that the strangelets produced in collisions tend to interact with hydrogen nuclei so as to produce harmless end-products, although they tend to produce larger negatively charged strangelets when interacting with other nuclei. Again, if so, strangelets produced in space would not survive to cause catastrophes, while strangelets produced on Earth would.

Third, it could be that, in the process of absorbing ambient hydrogen nuclei and producing decay products, an initially slowly moving strangelet would tend to be re-accelerated with respect to the ambient rest frame by randomly directed impulses from the decays. If so, it might tend to reach speeds greater than the relevant $`v_{\mathrm{crit}}`$, and then break up harmlessly in subsequent collisions with further hydrogen nuclei, long before forming part of a proto-star. It is not obvious that, if this was the typical fate of a strangelet produced in space, a similar fate should necessarily be expected for strangelets produced at RHIC: the respective environments are very different.

## IV Possible threat from a positively charged strangelet?

Empirical evidence aside, one of the theoretical arguments against the “killer strangelet” catastrophe scenario is that the existence of negatively charged strangelets would be theoretically surprising. Such a strangelet would have to contain a proportion greater than $`1/3`$ of $`s`$ quarks, which seems unlikely, since the higher mass of the $`s`$ quark suggests their proportion should be lower.

For this reason, it is generally believed that, if any strangelets were produced in RHIC or other collider experiments, they would be positively charged. Neither BJSW nor DDH seriously consider the possibility of producing positively charged strangelets, since they state that positively charged strangelets would pose no threat whatsoever.

The reason for this is that a positively charged strangelet, embedded in solid matter — for instance, the wall of a collider — would simply behave like an exotic nucleus. Coulomb repulsion would prevent it from approaching other nuclei and fusing with them, and so there would be no runaway reaction.

This argument does not, however, apply in all environments. In particular, if a positively charged strangelet reached the interior of the Sun, it seems quite possible that it would initiate a catastrophic runaway reaction, in the same way that a negatively charged strangelet would. So long as the Coulomb repulsion is overwhelmed by gravitational pressure (as it certainly would be initially) the strangelet charge would be essentially irrelevant.

The probability of a strangelet produced on Earth reaching the Sun without active human intervention, before the Earth is consumed by the Sun (by when it presumably would not matter), is presumably negligible. However, the possibility of human intervention needs to be considered. We do send spacecraft towards and even into the Sun, and objects in the vicinity of particle physics colliders — for instance, reused detector components — are considerably more likely than most to be included in these probes. Laboratory material potentially contaminated by positively charged strangelets would thus pose a potential danger which, though small, would need careful handling.***Another possible concern might be that a malicious entity could, at some point in the future, obtain material contaminated by positive strangelets and then deliberately threaten to cause catastrophe. Presumably, though, any group with the technological resources to send a spacecraft to the Sun, or to cause a fusion reaction on Earth, could credibly threaten catastrophe by less exotic means, if so inclined.

It should be stressed that BJSW give two theoretical arguments which suggest that the production of positively charged strangelets at RHIC, while less implausible, is still unlikely. These are the fact that theory seems to weigh against the possibility of stable strange matter of any type, and the fact that theory argues against the possibility of producing any form of strangelet in heavy ion collisions. BJSW find the second of these arguments particularly compelling. If either conclusion is correct, then, of course, there is no cause at all for concern. Still, theorists’ opinions do seem rather unsettled on this question. As recently as July 1999, the possibility of producing a (presumably positively charged) strangelet at RHIC was described as “a speculative but quite respectable possibility”.

Interestingly, DDH’s empirical arguments against the possibility of producing negatively charged strangelets at RHIC apply equally well to positively charged strangelets — though, of course, the weaknesses noted above remain. BJSW’s empirical arguments do not apply, since a positively charged strangelet created on the Moon would cause no catastrophe. Weaker bounds using similar reasoning might, however, be derived given that asteroids are also potentially vulnerable to strangelet production via cosmic ray impact, and that the infall rate of asteroid matter into the Sun has been estimated.

It is worth noting here that stronger bounds could be derived by assuming that our solar system is not very atypical and then combining BJSW’s and DDH’s lines of argument. Suppose that a significant proportion of stars are orbited by asteroids, satellites or planetoids unprotected by atmosphere, that those bodies tend to have similar chemical compositions to those of our solar system, and that the stellar infall rate of these bodies is similar to that of our solar system. Quantifying these assumptions, we could estimate the rate at which strangelets, created by cosmic ray impact on asteroids, were introduced into stars, causing supernovae. The bounds derived would apply whether the strangelets created are positively or negatively charged, assuming that the infall rate of asteroids converted into strange matter would be at least as great as that of unaffected asteroids. Unfortunately, though it seems very likely that our solar system is indeed not atypical, we have no firm evidence. Moreover, the potential flaws in DDH’s arguments, noted above, also apply to this line of reasoning. Nonetheless, it would be interesting to pursue it further.

## V Conclusions

To reiterate: this paper is written from the point of view of a devil’s advocate trying to point out potential flaws in the empirical arguments of BJSW and DDH. Those arguments rely on hypotheses which seem reasonable, but which could be wrong. To place any confidence in BJSW’s and DDH’s quantitative bounds on the risk of catastrophe, we need to be confident that the probability of their hypotheses being wrong is very low. BJSW and DDH have not carefully addressed this question.

Among the potential flaws in BJSW’s lunar survival arguments are (i) the possibility that the behaviour of iron nuclei may not be a reasonable model for that of gold nuclei; (ii) the possibility that their best guess at the strangelet rapidity distribution may be significantly wrong (without necessarily taking the extreme form suggested by DDH); (iii) the possibility that, although their distribution is roughly correct, the parameter $`a`$ lies outside the range they consider; (iv) the possibility that their strangelet survival probability function may be significantly wrong.

Among the potential flaws in DDH’s arguments based on supernovae observations are (i) the possibility that the critical speed above which strangelets are broken up in collisions with hydrogen nuclei might be substantially lower than their guess; (ii) the possibility that “killer strangelets” might be metastable; (iii) the possibility that they might tend to produce harmless end-products after combining with hydrogen nuclei, but not after combining with larger nuclei; (iv) the possibility that the process of absorbing stray hydrogen nuclei and emitting decay products might tend to re-accelerate a slowed strangelet in interstellar space, leading to its eventual breakup in a collision with a hydrogen nucleus.

It would be reassuring if all of these possibilities could be convincingly rebutted.

## VI Acknowledgements

Friends and colleagues too numerous to list here have provided invaluable help. I am very grateful to them all. This work was supported by the Royal Society and PPARC.