text:
# Determination of masses and radii of the massive eclipsing binary system HV 2543 in the Large Magellanic Cloud

## 1 Introduction

The analysis of light curves of eclipsing binaries, in addition to radial velocity data, provides fundamental knowledge about the masses and physical dimensions of the stars. Studying massive binaries in the Magellanic Clouds we can learn about the evolution of these systems at metallicities lower than that of our galaxy.

The Harvard variable HV 2543 ($`\alpha =5^\mathrm{h}27^\mathrm{m}27^\mathrm{s}`$, $`\delta =67\mathrm{°}11\mathrm{}54\mathrm{}`$, J2000) is a hot binary star in the LMC. It was catalogued as an OB star by Sanduleak (sk (1969)), who assigned it the identification –67$`\mathrm{°}`$117. The eclipsing nature of this binary was discovered by Gaposhkin (gap70 (1970)), who published a photographic light curve fitting a period of 4.829052 days (see also Payne-Gaposhkin pg (1971)). Photoelectric photometry was performed by Isserstedt (isser (1979)) who found $`V=12.92`$, $`(BV)=0.18`$ and $`(UB)=1.03`$. Radial velocity orbit was obtained by Niemela & Bassino (virpi (1994)) who derived physical parameters of the binary components and concluded that HV 2543 was a semidetached system with the less massive component filling its equipotential Roche surface. On the basis of their spectroscopic data, they classified this system as O8V:+O9III. Smith Neubig & Bruhweiler (neubig (1999)) have published an UV spectral classification of LMC OB stars based on IUE data, assigning to HV 2543 the type O9III.

In this paper we present a CCD $`V`$ light curve for HV 2543. By means of the combined analysis of these data and previously published radial velocities, we derive new values of the fundamental parameters of this system.

The paper is organized as follows: in Sect. 2 we describe the observations, reductions and calibration steps. In Sect. 3 we describe the photometric results and light curve fitting. In Sect. 4 we discuss the results and in Sect. 5 we present our conclusions.

## 2 Data Acquisition and Reductions

### 2.1 Observations

The CCD images here analysed were acquired with the 2.15-m telescope at CASLEO, during three runs in 1995, 1997 and 1998. About five frames of the HV 2543 field were obtained each night, except during eclipses, when the observation was more intensive. During each observing night, a series of $`1015`$ twilight flats and bias frames was also obtained. In addition, for the 1997 and 1998 runs, $`1015`$ dome flats were also obtained each night. Tables 1 and 2 list details of the observations and Fig. 1 shows a finding chart for HV 2543.

### 2.2 Reductions and Photometry

Reductions were performed using IRAF routines. For each night, the frames were reduced in the standard way (overscan and bias corrections, flat fielding) using properly combined bias and flat frames. The illumination differences between sky-flats and dome-flats were found to be less than 0.3 % at the edges of the frames, and negligible at the center. Whenever high quality sky-flats (that is, zero clouds) were available, these were preferred to dome-flats.

For each observing run, the following steps were performed:

1. A master image was made by aligning and combining the highest quality frames (best seeing and darkest sky). On this combined image, a profile fitting photometry was performed using a stand-alone version of DAOPHOT II (Stetson pbs87 (1987), Stetson pbs91 (1991)). The resulting output was used as master coordinate list in the following steps.
2. Aperture photometry was performed on each individual frame using a circular diaphragm of $`9.5\mathrm{}`$ diameter. The sky background was determined locally for each star using a sky annulus of radii 9.76 and $`20.33\mathrm{}`$.
3. A list of instrumental magnitudes for each individual frame and for each star in the master list was constructed.

Five faint neighbour stars were detected in the combined frame within a radius of $`9\mathrm{}`$ around HV 2543, being the brightest of these some 5 mag fainter than the binary. Light contamination from these stars slightly contributes to the measured magnitude of HV 2543 (differences of the order of 0.014 mag). Consequently, an adequate flux correction was applied to the aperture photometry.

### 2.3 Photometric calibration

No single objects were chosen as comparison and check stars. Instead, zero point corrections were computed for each frame using a group of $``$15 stars, covering widely the range of magnitudes of HV 2543. To put all the observations in an unique instrumental system, we proceeded as follows:

1. A synthetic reference system was defined using the mean instrumental magnitudes of the calibration stars.
2. For each frame, a zero point offset was computed and the magnitudes were transformed to the synthetic reference system.
3. After this transformation, the RMS residuals were computed for each star. Only the stars with the smallest residuals ($`<0.01`$) were conserved as calibration stars.
4. We repeated these steps with the new, cleaned calibration star list.

Since all nights were not really photometric, this procedure allowed us to have a direct estimate of the internal errors of each image.

The above enumerated steps were performed independently for each observing run. We proceeded in this way because the stars available as calibrators are not the same for each run, due to pointing differences. On the other hand, some stars remain perfectly stable during a run, but show luminosity variations from one year to another, being then inadequate as calibrators between different runs.

After this first stage of calibration was complete, a list of robust estimates of instrumental stellar magnitudes and their dispersions was constructed for each run. After cross identification between the three lists, the most stable stars were selected to serve as calibrators between the different runs.

During the 1997 run, 19 frames through the Johnson’s $`B`$ band were also acquired. These frames were used to determine the transformations to the standard $`V`$ system.

The final transformations to the standard system were performed by means of aperture photometry of 33 standard stars – ranging from $`(BV)=.004`$ to $`2.192`$ – from Landolt (lan (1992)) acquired during the photometric nights.

We used the following transformation equations:

$`b=B+b_1+b_2X+b_3(BV)`$

$`v=V+v_1+v_2X+v_3(BV)`$

We adopted the mean values of the extinction coefficients given by Minniti et al. (minni (1989)) and obtained the following values for the other coefficients:

Nov. 19, 1997:

$`b_1=4.2826\pm 0.0105`$

$`b_3=0.0208\pm 0.0115`$

$`v_1=3.4446\pm 0.0051`$

$`v_3=0.0456\pm 0.0055`$

Nov. 20, 1997:

$`b_1=4.2986\pm 0.0076`$

$`b_3=0.03116\pm 0.0086`$

$`v_1=3.4502\pm 0.0057`$

$`v_3=0.0535\pm 0.0067`$

The RMS residuals of the $`V`$ transformations were 0.017 mag. We first shifted our unified instrumental system to that corresponding to a $`V`$ frame acquired immediately before a series of B exposures, on Nov. 19, 1997. Then, we applied the standard transformations derived for that night. The Nov. 20 transformations were used only to check the quality of the first ones.

Our photometry of HV 2543 is presented in Table 3. In successive columns the heliocentric julian day, standard V magnitude, internal errors, seeing and airmass are given. The errors ($`\sigma _\mathrm{i}`$) account for the internal photometric errors estimated by DAOPHOT and those of the transformation to an unified instrumental system. They do not include the error of the transformation to the standard system.

## 3 Results

### 3.1 Ephemeris

We have observed one primary minimum during 1997, at HJD$`=\mathrm{2\hspace{0.17em}450\hspace{0.17em}768}\stackrel{d}{.}802\pm 0.001`$. Combining this value with that given by Payne-Gaposchkin (pg (1971)) we obtained a new ephemerides for HV 2543:

$`E_0=\mathrm{2\hspace{0.17em}450\hspace{0.17em}768}\stackrel{d}{.}802\pm 0.001`$

$`P=4\stackrel{d}{.}\mathrm{829\hspace{0.17em}046}\pm \mathrm{0.000\hspace{0.17em}004}`$

We have used this ephemeris in calculating the phases for the radial velocity and photometric data.

### 3.2 Light and radial velocity curve analysis

The $`V`$ light curve of HV 2543 is typical of near-contact binaries, with different depths of the minima and continuous variations of light in out-of-eclipse portions. Also noticeable is the O’Connell effect, i.e., the small difference in magnitude between both maxima (see Davidge & Milone damilo (1984)). We solved simultaneously the light and radial velocity curves using the Wilson & Devinney code (hereafter WD), which is very well suited for the study of close binary systems (Wilson & Devinney wd (1971), Wilson wil (1990)). Unit weight was given to photometric CCD observations since they all were collected with the same instrumental configuration and have comparable quality. A unit weight was also given to most of the radial velocity points except for those near eclipses or indicated by Niemela & Bassino (virpi (1994)) as less confident. We also corrected the data points corresponding to the phase $`\varphi =0.05`$, that are inverted (i.e., the O8 velocity corresponds to the O9 and vice-versa) in Niemela & Bassino (virpi (1994)). The relation between photometric and spectroscopic weights was given through the values of the mean standard deviations (sigma) of the data which were estimated to be 17 and 12 km s<sup>-1</sup> for the radial velocities of the primary and secondary components, respectively, and 0.02 mag for the light points. Only the ratios between these values are significant.

The radii that result from our preliminary light curve analysis indicate that the O8 component is more alike an O8III star. Chlebowski & Garmany (chg (1991)) give a temperature of 36000 K for an O8III star. Given that Schmidt-Kaler (1982) gives a lower temperature for such spectral type, we also performed an analysis assigning a 34700 K temperature for the primary. Standard values of bolometric albedos, $`A=1.0`$ (Rucinski ruc (1969)), and gravity darkening coefficients, $`g=1.0`$ (Lucy luc (1976)), for radiative envelopes were used. Linear limb-darkening coefficients were determined from tables by van Hamme (ham (1993)). These parameters were not adjusted. The adjustable parameters in our computations were: $`a`$ (semimajor axis), $`V_\gamma `$ (systemic radial velocity), $`i`$ (orbital inclination), $`T_2`$ (temperature of the secondary component), $`q`$ (mass-ratio), $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$ (modified potential of both components) and $`L_1`$ (luminosity of the primary). The first approaches to fit the observations were made by means of the light curve (LC) program until an acceptable fit to the $`V`$ light curve was obtained.

Then we proceeded with the differential corrections (DC) code. We first adjusted the parameters $`a`$, $`V_\gamma `$ and $`q`$ using only the out-of-eclipses radial velocity data. Thereupon, we left those parameters fixed and fitted $`i`$, $`T_2`$, $`\mathrm{\Omega }_1`$, $`\mathrm{\Omega }_2`$ and $`L_1`$, using only the light curve. Then we re-computed the first set of parameters leaving fixed the second set, and iterated this procedure until the corrections were smaller than their errors. A better solution was achieved allowing the fit of $`q`$ with the photometric data, that is, including $`q`$ within the second parameter set. Mode 2 (Leung & Wilson lw (1977)) which stands for detached systems was used at the beginning, but after a few runs of the DC program the computations clearly evolved towards a semidetached configuration with the secondary component filling its critical Roche lobe. Thus, mode 5 of the WD code was employed until the final solutions were found. In this mode, the parameter $`\mathrm{\Omega }_2`$ cannot be adjusted and it is set equal to the critical filling-lobe value. The solutions never converged to contact configurations.

Since some systematic differences were obtained between the observed and modelled light curves (due to the O’Connell effect), we tried to improve the binary modelling by including one or two hot spots on the surface of the stars. The spots were placed at the equator of the stars (colatitude $`=90\mathrm{°}`$) while the other spot parameters where adjusted with the WD code. The range of the adjusted temperature for the secondary star, $`T_2`$, depends on the adopted value for $`T_1`$. To estimate the errors of the other parameters we considered the differences between the values that arised from different solutions (i.e., using the “spectroscopic” value of $`q`$, 0.55, or leaving it to increase until 0.64, to allow a better fit to the photometric data). The adopted solution, with $`q=0.61`$, is a compromise between the optimal fit of the light curve and the tolerable deviation from the spectroscopic data. Table 4, lists the model parameters. Fig. 2 and 3 show the modelled light and radial velocity curves, respectively, derived from the one spot solution and plotted along with the observations, and their corresponding residuals (O-C).

The final value of the orbital inclination is incidentally high, very close to $`90\mathrm{°}`$. This fact should produce a total eclipse at the secondary minimum, but, as the dimensions of both components are quite similar, this is not a noticeable feature of the observed light curve.

Absolute values of the masses, dimensions and bolometric magnitudes of HV 2543 were computed from the one spot solution. They are listed in Table 5 where $`R_1`$ and $`R_2`$ correspond to mean values of the derived polar, back and side radii of the stars. The final configuration of HV 2543 shows two components of similar dimensions but different masses and temperatures, with the secondary less massive star filling its Roche lobe.

### 3.3 $`BV`$ Field photometry

We combined the $`B`$ frames acquired during 1997 to make a master $`B`$ image, on which we performed a profile fitting photometry. We used this photometry together with that performed on the $`V`$ master frame to generate a colour-magnitude diagram, which is shown in Fig. 4.

In this figure, the presence of an OB association covering the whole frame ($`9\mathrm{}`$ of diameter) is evident. Ostrov et al. (yoetal (1999)) found evidence of a stellar association surrounding the massive binary system Sk–67$°$105. The relatively small angular distance between the two stellar groups ($`9.5\mathrm{}`$) suggests that they are probably related.

We have obtained spectroscopic data for five stars in the neighbourhood of Sk-67$`\mathrm{°}`$105. For these stars, we have estimated the reddening by comparison between the intrinsic $`(BV)_0`$ colours given by Schmidt-Kaler (s-k (1982)) and FitzGerald (fg (1970)) and the observed colours. We derived $`E(BV)=0.17\pm 0.015`$, while for HV 2543 itself we obtained $`E(BV)=0.20`$. A detailed study of these OB associations will be presented in a forthcoming paper.

## 4 Discussion

We note that the radius derived for the O8 component is substantially larger than the one obtained by Niemela & Bassino (virpi (1994)). In fact, the spectroscopic data suggest that the luminosity of the O9 star is larger than that of the O8 star, which caused Niemela & Bassino to refer to the O9 component as “primary”, even though it is the one presenting the largest radial velocity amplitude, and consequently the less massive component of the binary system. Hence, we explored alternative solutions to the light curve that could yield a smaller radius for the O8 component. The values $`R_112.75`$, $`R_215.31`$, $`M_121.91`$ and $`M_219.34`$, in solar units, provide a good fit to the light curve, although they require a value of $`q`$ near 0.9, which is not compatible with the radial velocity data. Hence we have discarded this solution, which on the other hand, implies a rather low value of the distance modulus, namely 18.18, assuming $`A_V=0.62`$ (see below) and adopting the temperature scale of Chlebowski & Garmany (chg (1991)). This value would result even smaller if a lower temperature scale is adopted.

On the other hand, the temperature difference between the two components of HV 2543 resulting from the light curve analysis is larger than that suggested by the corresponding spectral types. This fact might be due to the difference between the spectral features of the unperturbed back sides and the heated inner sides of the stars. Also should we have in mind that the mass transfer history of the system might account for significant departures of the He abundances relative to those regarded as normal for LMC members, and this effect could influence some spectral lines of the secondary component.

To solve these puzzles, high dispersion spectroscopy would be desirable.

The semimajor axis and star dimensions of this system are alike those determined by Pritchard et al. (prit (1998)) for HV 2241, but the masses are somewhat smaller. We presume that these systems have experienced case A mass transfer, being now in the slow stage of mass exchange. In these cases, the mass gainer should be indistinguishable of a normal star, excepting that it would stand on an isochrone corresponding to a shorter age (see Vanbeveren et al. vbv (1998)). From the stellar models of Schaerer et al. (sch (1993)) (for single stars) we found that a $`25\mathrm{M}_{\mathrm{}}`$ star takes $`6.7`$ Myr until its radius grows to $`15\mathrm{R}_{\mathrm{}}`$, but such star would have an effective temperature of only 30000 K. The radius and effective temperature derived for the O8 star are consistent whit those of a single star of some 3.7 Myrs and $`40\mathrm{M}_{\mathrm{}}`$, completely out of the range of masses compatible with the radial velocity data. We note that a similar problem arises from the analysis of AB Crucis (Lorenz et al. lor (1994)).

It is clear that there are a series of phenomena that we do not fully understand, and consequently, our analysis can not be considered definitive. The errors given for the derived parameters must be considered with caution. True errors are not easy to estimate analytically, since they depend on the importance of phenomena that are not properly accounted for, such as wind shocks, uncertainty of the adopted temperature scale, radiation pressure effects, etc. In fact, the sizes and positions of the spots used to model empirically the O’Connell effect are rather arbitrary, and an equally satisfactory solution could be found with other parameters, but these details do not affect meaningly the derived star dimensions.

From our photometry we determine $`(BV)=0.11\pm 0.015`$ for HV 2543. This value is somewhat redder than $`0.18`$, obtained by Isserstedt (isser (1979)). However, given that Isserstedt does not detect the variability of HV 2543, we think that both measures are still in reasonable agreement. If we assume a $`(BV)_0`$ of -0.31 (Schmidt-Kaler s-k (1982)), $`R=3.1`$ (Koornneef koo (1982)) then it results an $`A_V=0.62`$ for HV 2543. This value does not depend on which temperature scale we adopt, since for the range of temperatures of the O-type stars the $`(BV)`$ colours are degenerated.

Estimating the bolometric corrections according to Massey & Hunter (mh (1998)), we derive a distance modulus of $`(mM)_0=18.31\pm 0.2`$ to $`18.40\pm 0.2`$, depending on the adopted temperature scale. The error accounts for the uncertainties in the estimates of $`R`$ and the bolometric corrections. This distance modulus must be considered with caution, since this system has experienced strong mass transfer and exhibits the above mentioned anomalies.

## 5 Conclusions

We have presented a CCD $`V`$ light curve of the eclipsing binary HV 2543 (Sk-67$`\mathrm{°}`$117) in the LMC. The light curve appearance is almost symmetric, with a slight ($`0.02V`$ magnitudes) O’Connell effect. From $`BV`$ photometry of the surrounding field, we found that HV 2543 probably belongs to an OB association not previously identified.

We have analysed this light curve and the published radial velocity data using the Wilson-Devinney code, finding that HV 2543 is a semidetached system with the less massive and less luminous member filling its Roche lobe. The primary minimum occurs when the O8 component is behind the O9. From our analysis, we derived fundamental parameters for this system, obtaining $`M_1=25.63\pm 0.7`$, $`M_2=15.63\pm 1.0`$, $`R_1=15.54\pm 0.4`$ and $`R_2=13.99\pm 0.4`$ in solar units. Some discrepancies between the present results and those derived from the previous spectroscopic analysis (Niemela & Bassino virpi (1994)), should be thoroughly addressed in future works.

We found that this system has experienced significant mass exchange, being the present most massive star the originally less massive.

###### Acknowledgements.

The authors acknowledge use at CASLEO of the CCD and data acquisition system supported under U.S. National Science Foundation grant AST-90-15827 to R. M. Rich. The focal reducer in use at CASLEO was kindly provided by Dr. M. Shara. This research has made use of the Astronomical Data Center catalogs. We are indebted to the staff of CASLEO for valuable help during the observing runs. We are grateful to V.S. Niemela for her vigorous comments made during the preparation of this paper. We would also like to thank to the referee for the suggestions that allowed to improve the presentation of the paper.


format: text

text:
# A MICROLENSING SEARCH FOR COLD MOLECULAR CLOUDS IN VIRGO

## 1 Rationale

The number of microlensing surveys underway is growing rapidly. Most of these surveys aim to detect compact objects in the dark halo of our own Galaxy or in the nearby M31 halo by monitoring stars in the LMC, SMC, or M31 itself. These projects have proved very fruitful detecting and placing limits on the dark matter in the Galaxy halo. The latest results indicate that less than $`20\%`$ of the Galaxy halo is in the form of MACHOs (see refs 1 and 2). There are, however, a number of shortcomings of the current surveys. Firstly, one cannot measure the distance to the MACHO and there is currently some controversy over whether the events seen towards the LMC are due to Galactic dark matter at all. Secondly, the objects have to be extremely compact to produce a microlensing effect at the distance of the LMC. Hence, current surveys provide no constraint on the hypothesis that the dark matter is composed of diffuse cold molecular hydrogen clouds (see e.g. refs 3 and 4). Thirdly, source blending greatly complicates the calculation of survey efficiencies.

We have begun a project which addresses these issues and provides completely new information on the nature of dark matter in galaxy clusters. We are monitoring 645 quasar candidates behind the Virgo galaxy cluster in order to detect dark objects in the cluster in the mass range $`5\times 10^7`$ $`M_{}`$$`5\times 10^3`$ $`M_{}`$. The upper mass limit is set by the length of the monitoring program and the lower limit by the frequency of the monitoring. Quasar source size can also have an effect on the lower mass limit. If the quasar’s angular size is larger than the angular diameter of the Einstein ring of the lens, the magnification drops significantly and the event would not be detected. If all the mass in the Virgo cluster were in dark objects within the above mass range then we can expect to see significant numbers of microlensing events during our monitoring program (see ref 5 for details).

Several features of our monitoring program deserve to be emphasised. Firstly, if we see a microlensing event we have a reasonably accurate measure of the distance to the MACHO, allowing us to calculate its mass. Secondly, we can detect objects down to surface mass densities of a few tens of $`\mathrm{g}\mathrm{cm}^2`$ as these would act as point mass lenses when placed in Virgo. Such objects, for example cold molecular gas clouds, are accessible to no other current microlensing survey. Thirdly, because we are monitoring $`600`$ objects over a 30 sq. deg. field, we do not have blending of source images and the associated problems that blending causes.

One potential worry for a survey such as ours is the fact that quasars are intrinsically variable. If quasar variability is large on the time scales over which we are monitoring we would have to set our microlensing magnification threshold so high that the expected number of events would become very low. To date, little systematic quasar monitoring has been carried out on time scales matching those of our project, thus we have been careful to examine the feasibility of our project. As reported in this article, our pilot project proves that intrinsic quasar variability is not a serious problem on these time scales. A by-product of this microlensing survey is the information provided on the variability of quasars over these timescales, which will lead to constraints on models of quasar fuelling mechanisms.

## 2 Observations and data analysis

During the period 1999 February - June we obtained 28 R-band Schmidt plates of the central 30 sq. deg. of the Virgo cluster. The data are complete to an R-band magnitude of $`20`$. Using previous plates of the same area in the U and B bands we identified $`645`$ candidate quasars according to their colour. Of these we expect $`\mathrm{}<10\%`$ to be contaminant stars. Figure 1 shows $`2`$ typical quasar light curves extracted from this data set.

It can be seen qualitatively that the quasar candidates vary little during the monitoring period. We have quantified this by comparing the level of variability in the quasar candidates with the variability of $`3000`$ stars in our field. The candidate quasars are statistically no more variable than the stars. This confirms the feasibility of using quasars for microlensing surveys for monitoring periods of a few months. Our simulations show that with this low level of variability we can set a low detection threshold, and identify microlensing events down to a peak magnification of 0.46 mag.

We have run a matched filter analysis on each of the $`645`$ candidate-quasar light curves. The filter used was a MACHO light curve of peak amplification $`0.6`$ mag. We used a range of values for the event time-scale and the time of peak magnification. For each quasar light curve the event with the highest $`\mathrm{\Delta }\chi ^2`$ for a match to a MACHO event was retained. To this best event for each quasar, we performed a four parameter MACHO light curve fit.

To define selection criteria to pick out likely lensing events, we ran the same set of procedures on the light curves of the 3000 control stars in the field – where we do not expect any real microlensing events. Figure 2 shows the distribution of points in the plane of $`\mathrm{\Delta }\chi ^2`$ for the four parameter MACHO fit against the fitted impact parameter $`b_{min}`$. Quasar candidates are shown by crosses and stars as stars. We define the following criteria for an event to be considered as a serious microlensing candidate:
$``$ $`\mathrm{\Delta }\chi ^2>100`$
$``$ $`\text{Impact parameter}b_{min}<0.8`$
These criteria leave two events in the candidate region to be checked. Neither of these proved to be real events: one is a known highly variable BL Lac object and the other a probable stellar RR Lyrae contaminant. From this null detection we can now put constraints on the nature of dark matter in Virgo.

## 3 Detection Efficiency and Predicted Number of Events

To calculate the detection efficiency of the survey we have generated artificial microlensing events and added them to the 3000 stellar light curves. For each event time-scale we have generated $`10,000`$ artificial light curves with impact parameters distributed uniformly between $`0`$ and $`1`$ Einstein ring radius and added them to the stellar light curves, placed at random dates. We then applied our matched filter analysis, the fitting routine, and the selection criteria to these artificial light curves and measured how many of the artificial events were recovered. Our detection efficiency is shown as a function of event time-scale on the left hand side of Figure 3.

Using the detection efficiency we can calculate the number of microlensing events we would expect to see in our experiment if the dark matter in Virgo were in the form of dark objects of a particular mass. We assume that the Virgo cluster can be modelled as an isothermal sphere with a one-dimensional velocity dispersion of $`673\mathrm{km}\mathrm{s}^1`$. Our expected number of events is shown on the right hand side of Figure 3 and peaks at $`5`$, if all the dark mass in Virgo were in the form of $`10^5`$$`M_{}`$objects.

## 4 Conclusions and future work

From our null detection, and the fact that we would have expected $`5`$ events for lenses of mass $`1\times 10^5`$ $`M_{}`$, we can use Poisson statistics to put constraints on the mass in Virgo in the form of dark objects of this mass. The result of this pilot project is that less than $`1/2`$ the mass in Virgo can be in the form of $`10^5`$ $`M_{}`$objects at $`90\%`$ confidence.

The most important thing we have learnt from the pilot project is that a quasar monitoring microlensing project is technically feasible and we hope, over the next 2 years, to acquire of order $`100`$ more plates of Virgo. Our simulations indicate that this will allow us to improve our constraints on the nature of dark matter in Virgo by a factor $`510`$. A second season of monitoring Virgo started in February 2000.

We would also like to extend the experiment by monitoring other clusters - especially the Perseus cluster in the North as this is the most massive nearby cluster. Alternatively, another strategy would be to monitor several more distant clusters (thereby avoiding the need for a very wide field of view). In the long term we would also like to monitor a control field of quasars with no intervening massive cluster in order to precisely quantify quasar variability.

## Acknowledgements

We are immensely grateful to the Schmidt observing team: Fred Watson, Malcolm Hartley, Russell Cannon, Paul Cass, Ken Russell and Delphine Roussiel for taking all the photographic plates of Virgo used in this analysis.

## References


text:
# Coulomb gap in a model with finite charge transfer energy.

## I INTRODUCTION

Doping of solids might lead to drastic qualitative changes in their properties. The metal-insulator transition (MIT) is a spectacular manifestation of this. The understanding of the driving forces of the MIT is a long-standing problem. In the early seventies, the prediction  was made that on the dielectric side of the MIT the long-range Coulomb interactions deplete the density of one-electron excitations (DOE) $`g(\epsilon )`$ near the Fermi energy $`\mu `$. Further, analytical calculations of $`g(\epsilon )`$ with Coulomb correlation taken into consideration have been performed on the metallic side of the MIT. Altshuler and Aronov  showed that for the metallic case $`g(\epsilon )`$ in three dimensions has a cusp-like dependence $`g(\epsilon )|\epsilon \mu |^{1/2}`$ near $`\mu `$. This was later confirmed in electron tunneling experiments for amorphous alloys  and granular metals .

On the insulating side of the MIT charge transport occurs via inelastic electron tunneling hopping between states localized on the impurity sites with one-electron energies close to $`\mu `$. Mott  demonstrated that at low temperatures electrons seek accessible energy states by hopping distances beyond the localization length, leading to a hopping conductivity $`\sigma (T)\mathrm{exp}(T_0/T)^\nu `$ with $`T_0`$ being a characteristic temperature depending on localization length and with the hopping exponent $`\nu =1/4`$ for the non-interacting case in three dimensions. Efros and Shklovskii  (ES) argued that the ground state of a system with long-range Coulomb interactions is stable with respect to one-particle excitations only if $`g(\epsilon )`$ in the vicinity of $`\mu `$ has the symmetric shape

$$g(\epsilon )|\epsilon \mu |^{D1}$$
(1)
with the universal exponent $`D1`$ depending only on the dimensionality $`D`$ of the system. In particular, ES predicted that in $`D=3`$ $`g(\epsilon )=\frac{3}{\pi }\left(\frac{\chi }{e^2}\right)^3(\epsilon \mu )^2`$, where $`\chi `$ is the dielectric constant and $`e`$ is the electron charge. Because $`g(\epsilon )`$ vanishes only at $`\epsilon =\mu `$, this is called a “soft” Coulomb correlation gap with a width $`\mathrm{\Delta }\epsilon e^3(N_0/\chi ^3)^{1/2}`$, where $`N_0`$ is the DOE far away from $`\mu `$. The power law (1) gives  a hopping exponent $`\nu =D/(D+3)`$ at low temperatures, so for three-dimensional system with long-range Coulomb interactions $`\nu =1/2`$.

The intriguing hypothesis about universality of (1) has stimulated further theoretical research, both analytical  and numerical . To establish the hypothesis (1) Efros used the ground-state stability conditions for localized electrons (LES) with respect to charge transfer

$$\epsilon _j\epsilon _i\frac{e^2}{\chi r_{ij}}>0,$$
(2)
where $`\epsilon _i`$ and $`\epsilon _j`$ are the one-particle energies of a neutral donor on a site $`i`$ and of a charged donor on a site $`j`$, respectively, and $`r_{ij}`$ is the distance between the sites $`i`$ and $`j`$. The conditions (2) were used to heuristically derive a non-linear integral equation for $`g(\epsilon )`$  and then assymptotic analysis of this equation leads  to the power law (1).

LES have been studied using the so-called classical donor-acceptor ($`d`$-$`a`$) model (see, e. g. Ref.). Within this model, the system considered is modeled by a continuous sample with randomly distributed $`k\times N`$ ($`k1`$) acceptor and $`N`$ donor sites. Each acceptor site is negatively charged whereas out of $`N`$ only $`k\times N`$ donors have a positive charge which leads to a large number of configurations of charged donors. Moreover, each of these configurations must obey not only conditions (2) but also more complicated conditions related to many-particles excitations (e. g., charge transfer involving four, six, etc. sites). Efros conjecture about the universality implies that $`g(\epsilon )`$ does not depend on peculiarities of the particular model and, as a consequence, further theoretical studies of LES  were confined to a lattice $`d`$-$`a`$ model proposed in Ref.. In this model, $`N`$ donors are localized on all the sites of a D-dimensional lattice and the negative charge from $`k\times N`$ acceptors is uniformly smeared over the lattice sites so that each site $`i`$ has a charge $`e(n_ik)`$, where $`n_i=1`$ if a donor on the site $`i`$ is ionized and $`n_i=0`$ if a donor is neutral. Disorder in this model is ensured by introducing randomly distributed one-site potentials. Monte Carlo simulations  on very large specimens of the lattice $`d`$-$`a`$ model, however, have given rise to doubts about the universality of the $`g(\epsilon )`$ behavior.

Another hint about possible non-universal behavior of $`g(\epsilon )`$ has come from the intriguing and still not completely unfolded problem whether the so called spin-glass phase does exist in the classical $`d`$-$`a`$ model (see, e. g. Ref.). Grannan and Yu  studied the classical three-dimensional $`d`$-$`a`$ model with $`k=0.5`$ but with the total acceptor charge uniformly distributed over donor sites as in the lattice $`d`$-$`a`$ model. In this case, the classical $`d`$-$`a`$ model is equivalent to a model of Ising spins, localized on randomly distributed sites, with pairwise Coulomb interactions, a model in which a transition into the spin-glass state was found  to occur at non-zero temperature. It was then concluded that such a transition should exist in all $`d`$-$`a`$ models (with and without smearing of negative charge, defined on a lattice or on a continuous sample) as well because of the Efros universality hypothesis. Voita and Schreber , however, have shown that the spin glass transition does not exist in the lattice $`da`$ model. Besides, in recent work by one of us  it was unequivocally demonstrated that the ground state of the classical $`d`$-$`a`$ model and that of the model studied in Ref. are qualitatively different. An analysis of histograms $`[Q_{\alpha \beta }]`$ of the so called overlaps $`Q_{\alpha \beta }=\frac{1}{N}_i\delta (n_i^\alpha ,n_i^\beta )`$ (here $`\alpha `$ and $`\beta `$ refer to different pseudo-ground states (PGS) obtained by direct descents) has revealed that, indeed, for the model studied in Ref. $`[Q_{\alpha \beta }]`$ has a symmetric Gaussian shape with the maximum at $`Q_{\alpha \beta }=0`$ and with the dispersion $`Q_{\alpha \beta }^2N^1`$. This means that a large number of microscopically different PGS’s does exist in the model and according to Parisi’s theory  this implies the existence of a spin-glass state at low temperatures. Further Monte Carlo simulations at finite temperatures  revealed the typical finite-size scaling of the spin-glass susceptibility. In the classical $`d`$-$`a`$ model, however, $`[Q_{\alpha \beta }]`$ has its maxima at $`Q_{\alpha \beta }=1`$ which means that all PGS’s generated are the same from microscopical point of view. The absence of microscopically different PGS’s in the classical $`d`$-$`a`$ model was explained  by the pinning of all PGS’s on the electric field created by the discretely distributed acceptor charges.

Therefore, it is highly desirable to study the properties of not only the classical $`d`$-$`a`$ model, but of its various modifications as well. In the present work we consider a modified classical $`d`$-$`a`$ model (MCDAM) in which acceptors can be neutral, so the energy $`\mathrm{\Delta }`$ of the charge transfer from a donor to an acceptor ($`d^0+a^0d^++a^{}`$, where $`d^0`$, $`d^+`$, $`a^0`$, $`a^{}`$ stand for a neutral donor, a charged donor, a neutral acceptor, a charged acceptor, respectively) has to be finite. The classical $`d`$-$`a`$ model might be then viewed as the limit of the MCDAM as $`\mathrm{\Delta }\mathrm{}`$. We have investigated the shape of the Coulomb gap (i. e. $`g(\epsilon )`$ both for donor and acceptors in the vicinity of the Fermi level) in two- and three-dimensional MCDAMs at T=0 and found that the behavior of $`g(\epsilon )`$ is in strong contradiction to the Efros conjecture about the universality of $`g(\epsilon )`$. The rest of the paper is organized as follows. In Section II we introduce the MCDAM and arrive at some rigorous results which follow from a symmetry of the MCDAM with respect to the exchange of donor and acceptor sites. Further, the algorithm of energy minimization for the MCDAM including a discussion about inherent finite size effects is presented in Section III. Section IV is devoted to a description of the main results obtained. In Section V we discuss possible causes of universality violation in the MCDAM, analyze experimental data available in the literature and predict possible experimental situations in which the non-universal behavior of $`g(\epsilon )`$ might be observed. And finally, a summary is presented in Section VI.

## II BACKGROUND

### A Model

We consider a $`D`$-dimensional system of volume $`L^D`$, in which an equal number of acceptor and donor sites $`N`$ are allocated according to the Poisson distribution with a density $`n=N\times L^D`$. It is convinient to choose energy unit $`E_0`$ as an energy of the Coulomb interaction between a pair of acceptors, say, localized on the average distance $`n^{1/D}`$, $`E_0=e^2n^{1/D}/\chi `$. In typical bulk semicondustors $`n10^{18}`$ cm-3 and $`\chi 10`$, so $`E_00.02`$ eV. Hereafter all expressions will be written in dimensionless units $`n^{1/D}`$ for length and $`E_0`$ for energies. A microscopic state of a particular spatial arrangement of the donor and acceptor sites (henceforth referred to as the sample R) is determined by a set of occupation numbers $`(n_a,n_d)\{n_a(i),n_d(k),i=1,2,\mathrm{},N`$, $`k=1,2,\mathrm{},N\}`$ determined in the following way. For the acceptors, $`n_a(i)=1`$ if an acceptor on an acceptor site $`i`$ has captured an electron and $`n_a(i)=0`$ if an acceptor is neutral. For the donors, $`n_d(k)=1`$ if a donor on a donor site $`k`$ is neutral and $`n_d(k)=0`$ if a donor has given an electron away. We investigate the LES from the dielectric side of the MIT, so $`\alpha _B<1`$ ($`\alpha _B`$ is the localization length of the electron on donor). The energy of the sample, assuming that all the interactions are of Coulomb origin, then is

$`E(n_a,n_d)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{ij}{}}{\displaystyle \frac{n_a(i)n_a(j)}{r_{ij}^{aa}}}+`$ (3)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{kl}{}}{\displaystyle \frac{(1n_d(k))(1n_d(l))}{r_{kl}^{dd}}}`$ (4)
$``$ $`{\displaystyle \underset{i,k}{}}{\displaystyle \frac{(1n_d(k))n_a(i)}{r_{ik}^{ad}}}\mathrm{\Delta }{\displaystyle \underset{i}{}}n_a(i),`$ (5)
where indices $`i,j`$ and $`k,l`$ number acceptor and donor sites, respectively, $`r_{ij}^{aa}`$, $`r_{kl}^{dd}`$ and $`r_{ik}^{ad}`$ are the distances between the acceptors on the sites $`i`$ and $`j`$, between the donors on the sites $`k`$ and $`l`$, and between the acceptor on the site $`i`$ and the donor on the site $`k`$, correspondingly, and $`\mathrm{\Delta }`$ is the energy of charge transfer between acceptor and donors. When charge transfer occurs in the system, the energy of the sample changes by

$`\delta E(n_a,n_d)`$ $`=`$ $`{\displaystyle \underset{i}{}}\epsilon _a(i)\delta n_a(i)+{\displaystyle \underset{k}{}}\epsilon _d(k)\delta n_d(k)+`$ (6)
$`+`$ $`{\displaystyle \underset{i,k}{}}{\displaystyle \frac{\delta n_a(i)\delta n_d(k)}{r_{ik}^{ad}}}+{\displaystyle \frac{1}{2}}{\displaystyle \underset{ij}{}}{\displaystyle \frac{\delta n_a(i)\delta n_a(j)}{r_{ij}^{aa}}}+`$ (7)
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{kl}{}}{\displaystyle \frac{\delta n_d(k)\delta n_d(l)}{r_{kl}^{dd}}},`$ (8)
where $`\epsilon _a(i)`$ is the one-electron excitation (OEE) energy for the acceptors

$$\epsilon _a(i)\frac{\delta E(n_a,n_d)}{\delta n_a(i)}=\underset{ji}{}\frac{n_a(j)}{r_{ij}^{aa}}\underset{k}{}\frac{1n_d(k)}{r_{ik}^{ad}}\mathrm{\Delta },$$
(9)
$`\epsilon _d(i)`$ is the corresponding OEE energy for the donors and $`\delta n_a(i)`$ ($`\delta n_d(k)`$) denotes the change of the occupation number on the acceptor (donor) site $`i`$ ($`k`$). If a microscopic state ($`n_a^0`$,$`n_d^0`$) of the sample is the ground-state of this sample then for any excitation the relation

$$\delta E(n_a^0,n_d^0)0$$
(10)
holds. The specific appearance of the conditions (10) depends on what excitations are allowed in the model system considered.

In the present paper we investigate the simplest case when only pairs of sites are involved in the charge transfer which, in turn, is allowed to occur in four different ways: (i) via electron hops between a pair of the acceptors $`\{n_a(i)=1,n_a(j)=0\}\{n_a(i)=0,n_a(j)=1\}`$; (ii) via electron hops between a pair of donors $`\{n_d(k)=1,n_d(l)=0\}\{n_d(k)=0,n_d(l)=1\}`$; (iii) via ionization process $`\{n_a(i)=0,n_d(k)=1\}\{n_a(i)=1,n_d(k)=0\}`$ and (iv) via recombination process $`\{n_a(i)=1,n_d(k)=0\}\{n_a(i)=0,n_d(k)=1\}`$. For each of those processes there is an unique set of $`\{\delta n_a(i),\delta n_a(j),\delta n_d(k),\delta n_d(l)\}`$. For instance, for the acceptor-acceptor hops

$$\delta n_a(i)=1,\delta n_a(j)=1,\delta n_d(k)=0,\delta n_d(l)=0.$$
(11)
Substituting (11) into (8) one obtains the ground-state stability relation with respect to the charge transfer between the pair of acceptors on the sites $`i`$ and $`j`$

$$\epsilon _a^0(j)\epsilon _a^1(i)\frac{1}{r_{ij}^{aa}}0,$$
(12)
where $`\epsilon _a^{1(0)}(i)`$ denotes $`\epsilon _a(i)`$ if $`n(i)=1(0)`$. The stability conditions with respect to the other three manners of the charge transfer are obtainable in the similar manner.

The relation (12) implies that $`\epsilon _a`$’s for the neutral acceptors are, in general, larger than $`\epsilon _a`$’s for the charged acceptors. Furthermore, the pair of neutral and charged acceptors might be located on any distance and therefore in the thermodynamic limit the chemical potential for the acceptors (i. e. an energy level which separates the energies of the neutral and charged acceptors) is determined as

$$\mu _a=\text{min}\{\epsilon _a^0(i)\}=\text{max}\{\epsilon _a^1(i)\}.$$
(13)
Alike, there exist the chemical potential $`\mu _d`$ for the donors as well. Moreover, the stability relations with respect to the ionization and recombination lead to

$$\mu _a=\mu _d=\mu .$$
(14)
Despite the finite size of samples we investigated, the relation (14) with the chemical potentials calculated from (13) is valid within the limits of accuracy of our calculations (see Sect. III).

A macroscopic state of the sample R is characterized by degree of acceptor ionization

$$C_a(𝐑)=\frac{1}{N}\underset{i}{}n_a(i),$$
(15)
by the DOE for acceptors

$$g_a(\epsilon _a,𝐑)=\frac{1}{N}\underset{i}{}\delta (\epsilon \epsilon _a(i))$$
(16)
and by the corresponding DOE $`g_d(\epsilon _d,𝐑)`$ for the donors. Note, that for the finite samples (especially for the relative small systems we were able to investigate) $`C_a(𝐑)`$, $`g_a(\epsilon _a,𝐑)`$ and $`g_d(\epsilon _d,𝐑)`$ depend essentially on the particular implementation R of the spatial distributions of the donor and acceptor sites (if a sample would be big enough all quantities would be self-averaging). Therefore, in order to obtain reliable results, one has to work with the quantities $`C_aC_a(𝐑)`$, $`g_a(\epsilon )g_a(\epsilon _a,𝐑)`$ and $`g_d(\epsilon )g_d(\epsilon _d,𝐑)`$, where $`\mathrm{}`$ denotes the average over a number of R’s. Note, that the values $`g_{a(d)}(\epsilon _{a(d)},𝐑)d\epsilon `$ obtained for independent R’s are scattered according to the Gaussian distribution with the mean $`g_{a(d)}(\epsilon )d\epsilon `$ and the standard deviation $`\sqrt{g_{a(d)}(\epsilon )d\epsilon }`$. In the region of the Coulomb gap $`g_{a(d)}(\epsilon )d\epsilon 10^4`$ and dispersion is several orders of magnitude larger than the mean. Therefore, in order to reduce the statistical noise in the final $`g_{a(d)}(\epsilon )`$ dependences an average is needed over a sufficient large amount of independent samples (we performed calculations with up to 104 samples).

### B Acceptor-donor symmetry

Let us rewrite the energy (5) in terms of the OEE energies (9)

$`E(n_a,n_d)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i}{}}\epsilon _a(i)n_a(i)`$ (17)
$``$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{k}{}}\epsilon _d(k)(1n_d(k)){\displaystyle \frac{\mathrm{\Delta }}{2}}{\displaystyle \underset{i}{}}n_a(i).`$ (18)
The system investigated is electrically neutral, i. e. for any sample

$$\underset{i}{}n_a(i)=\underset{k}{}(1n_d(k)).$$
(19)
Then, the energies of the microscopic states ($`n_a,n_d`$) and ($`n_a^{},n_d^{}`$) for a sample R and its “mirror” reflection R (when the donor and acceptor sites exchange places keeping the spatial arrangement of sites unchanged), are equal under the following conditions

$$\epsilon _a(i)+\epsilon _d^{}(i)=\epsilon _d(k)+\epsilon _a^{}(k)=\mathrm{\Delta }$$
(20)
and

$$n_a^{}(i)=(1n_d(i))n_d^{}(k)=(1n_a(k)).$$
(21)
The stability relations (12) for the ground-state ($`n_a^0`$,$`n_d^0`$) of the sample R transform into stability relations for the ground-state ($`n_a^0`$,$`n_d^0`$) of the sample R through the relations (20,21) as well.

Since averaging over samples includes all possible pairs R and R, it follows from the symmetry relations (20) and (21) along with the definition (16) that $`g_a(\epsilon )`$ can be mapped to $`g_d(\epsilon )`$ using the relation

$$g_d(\epsilon )=g_a(\epsilon \mathrm{\Delta })$$
(22)
The symmetry of the model imposes also a relation between the Fermi energy $`\mu `$ (13,14) and the parameter $`\mathrm{\Delta }`$ of the model. Expressing $`n_{a[d]}(i[k])`$ in terms of the Heaviside’s step functions $`n_{a[d]}=\theta (\mu \epsilon _{a[d]}(i[k]))`$, the quantity $`C_a`$ can be written in the form

$$C_a=_{\mathrm{}}^\mu g_a(\epsilon )𝑑\epsilon =_\mu ^{\mathrm{}}g_d(\epsilon )𝑑\epsilon $$
(23)
The symmetry relation (22) transforms (23) into an integral relation

$$_{\mu \mathrm{\Delta }}^{\mathrm{}}g_d(\epsilon )𝑑\epsilon =_\mu ^{\mathrm{}}g_d(\epsilon )𝑑\epsilon ,$$
(24)
which has a meaning only if

$$\mu =\frac{\mathrm{\Delta }}{2}.$$
(25)
Thus, the Fermi energy of our model system in the thermodynamic limit is a fundamental quantity depending only on the energy of charge transfer from an acceptor to a donor.

## III METHOD

### A Algorithm of energy minimization

We start from a random allocation of $`N`$ donor and $`N`$ acceptor sites in the continuous $`D`$-dimensional system (generate a sample R) with the density $`n=1`$, so that the system has a linear size $`L=N^{1/D}`$ and then charge randomly chosen $`C_a\times N`$ both donors and acceptors (usually we take $`C_a=0.7`$), i. e. generate an initial microscopic state (IMS) ($`n_a`$,$`n_d`$) of the sample R. Further, we search for such microscopic state ($`n_a^0`$,$`n_d^0`$) which obeys the stability conditions (12) with respect to the four mechanisms of the charge transfer allowed in our model. We used an algorithm which is an extension of the algorithm proposed in Ref. to the case $`\mathrm{\Delta }\mathrm{}`$. The algorithm consists of the three main steps.

In order to save computer time, first, we look for pairs $`a^0a^{}`$ ($`d^0d^+`$) for which the “crude” stability relation $`\mathrm{\Delta }\epsilon \epsilon _{a(d)}^0\epsilon _{a(d)}^1>0`$ is violated. Then, the energy of the system is decreased by transferring an electron between such pair of sites for which $`\mathrm{\Delta }\epsilon `$ has its minimal non-positive value. This process is repeated until a state is reached, in which $`\mathrm{\Delta }\epsilon >0`$ for all possible $`a^0a^{}`$ and $`d^0d^+`$ pairs (step I). In the similar manner, we further minimize the energy of the system with respect to the “true” stability relations (12) for the charge transfer between the $`a^0a^{}`$ and $`d^0d^+`$ pairs (step II). And, finally, in the step III we diminish the energy of the system with respect to the stability relations for ionization and recombination processes. Since ionization and recombination processes change the degree $`C_a`$ of the system ionization, each time after one of these processes takes place during calculations, we go back to the step II. Repeating the steps II and III, we finally arrive at a microscopic state ($`n_a^0`$,$`n_d^0`$) for which all four stability conditions are fulfilled. We name the procedure ($`n_a`$,$`n_d`$)$``$ ($`n_a^0`$,$`n_d^0`$) via above steps I,II and III as “a single descent”.

It should be noted, however, that the state ($`n_a^0`$,$`n_d^0`$) is not necessarily the ground state of the sample R since for the ground state, in general, not only the simplest relations (12) with only pairs of sites included, but the more complicated relations involving quadruplets, sextets, etc. of sites have to be fulfilled. Therefore, the state ($`n_a^0`$,$`n_d^0`$) (after Ref. ) hereafter will be referred to as the pseudo-ground state (PGS) of the sample R. Then, two questions naturally arise: How close the PGS and the ground state of the given sample are and how this may influence the output of our calculations? In order to answer the first question, we calculate and analyze the histograms $``$ for the so-called overlaps

$$Q_{\alpha \beta }=\frac{1}{N}\underset{i}{}\delta (n_a^\alpha ,n_a^\beta ),$$
(26)
where indices $`\alpha `$ and $`\beta `$ refer to PGS’s which are obtained by means of the single descent on the same sample but with different IMS ($`n_a`$,$`n_d`$). If two PGS’s are identical then $`Q_{\alpha \beta }=1`$. We calculated for the $`D=2`$ system with $`N=500`$ at $`\mathrm{\Delta }=0`$ the mean $`Q(𝐑)=Q_{\alpha \beta }_{\alpha \beta }`$ for the sequence of 100 PGS’s generated by single descents from the different IMS of the same sample R. We further acquire $`Q(𝐑)`$ for 100 different samples and obtain that the mean $`\overline{Q}Q(𝐑)_𝐑=0.96`$. It means that in PGS generated by the single descent only 20 acceptors out of 500 are, in average, in the “wrong” states compared to those in the true ground state of the sample.

In order to evaluate how the “erroneousness” of PGS influences the outcome of our calculations we perform an analysis of ground states obtained by means of the so called multirank descents. Descent of rank $`m`$ comprises of a consequence of the single descents on the same sample with different IMS when calculations are stopped after the lowest observed PGS energy repeats $`m`$ times. We calculate $`\overline{Q}`$ (all other parameters were the same as described in the previous paragraph, where actually the case $`m=0`$ was explored) for descents with different ranks $`m=5,10,15`$ and found that, for instance, for $`m=15`$ (which implies drastic increase in the computation time) $`\overline{Q}=0.990`$. $`g_a(\epsilon )`$ and $`g_d(\epsilon )`$ obtained from the PGS’s generated by means of the single descents and by means of descents with $`m=10`$, say, do not differ within the limits of statistical errors. So, we conclude, that reliable results can be obtained by means of single descents already, thereby saving a lot of computer time and resources.

### B Finite-size effects

Due to constraints in computer resources, the largest samples, we were able to deal with, comprise up to $`N=2000`$ donor and $`N=2000`$ acceptor sites ($`L45`$ for $`D=2`$ and $`L12`$ for $`D=3`$). Such relative small sizes of the samples investigated might influence the outcome of calculations. Detailed analysis of finite size effects on the results obtained will be presented in Section IV and here we want to make two remarks about inherent finite size effects in the model system considered.

First, as follows from (12), the energies $`\epsilon _a^0`$ for the neutral acceptors and $`\epsilon _a^1`$ for the charged ones in finite samples at $`T=0`$ cannot be further away than $`(L\times \sqrt{D})^1`$. This implies that $`g(\epsilon _a)=0`$ within the $`\epsilon _a`$ interval

$$|\epsilon _a\mu |<(2L\times \sqrt{D})^1$$
(27)
Of course, the same holds for donors as well. The relation (27) gives the estimation how close to $`\mu `$ data on the energy spectrum are, in principle, obtainable from the calculations on finite samples.

Secondly, as follows from (9) the energies $`\epsilon _a`$ and $`\epsilon _d`$ for the finite samples are sensitive to the location of the donor and acceptor sites. Therefore, the Fermi energy $`\mu `$ for finite samples does differ, in general, from sample to sample. A straightforward averaging of $`g(\epsilon )`$ over different samples might thus lead to a distortion of the $`g(\epsilon )`$ shape especially in the region where the Coulomb gap is observed. In order to avoid this undesired effect, we used a trick first proposed in Ref.. During accumulation of the results for $`g(\epsilon )`$ we added together $`g(\epsilon )`$ for the same values of $`\epsilon \mu (𝐑)`$ rather than for the same values of $`\epsilon `$. Here $`\mu (𝐑)`$ denotes the Fermi energy for a finite sample R calculated as

$$\mu (𝐑)=\frac{1}{2}\left(\text{min}\{\epsilon _a^0(i)\}+\text{max}\{\epsilon _a^1(i)\}\right),$$
(28)
Such way of doing $`g(\epsilon )`$ average entirely excludes the influence of the fluctuations of the Fermi energy in the finite samples on the shape of the Coulomb gap.

Finally, we remark that all the data presented below were obtained for the open boundary condition. In order to ensure that results obtained are not determined by the type of the boundary conditions used in calculations, we performed calculations of the two-dimensional MCDAM at $`\mathrm{\Delta }=0`$ with different $`N`$ and found that periodic boundary conditions only effectively reduce the linear size of a sample, leaving the qualitative shape of the parameters calculated unchanged.

## IV RESULTS

According to the symmetry relation (22) $`g_a(\epsilon )`$ and $`g_d(\epsilon )`$ can be easily mapped to each other for any values of $`\epsilon `$ and hence all the results presented below concern the acceptor sites only. One can expect that the width of the Coulomb gap $`\mathrm{\Delta }\epsilon `$ and the energy scale in our model $`E_0=e^2n^{1/D}/\chi `$ are of the same order of magnitude. Fig.1 shows $`g_a(\epsilon \mu )`$ in the vicinity of the Fermi energy $`\mu `$ obtained for the two-dimensional samples with $`N=1000`$ and various values of $`\mathrm{\Delta }`$. As it is seen, $`g_a(\epsilon \mu )`$ depends considerably on $`\mathrm{\Delta }`$ except for a narrow window $`|\epsilon \mu |0.05`$, where all data merge into some “universal” curve symmetric with respect to $`\mu `$, the curve which can be anticipated to obey the Efros universality hypothesis (1). However, a double-logarithmic plot of the “universal” $`g_a(\epsilon \mu )`$ (insert in the Fig.1), reveals that the behavior of $`g_a(\epsilon \mu )`$ in the “universality” region is not even a power law. The width of this “universality” region is comparable to the width of the region where $`g_a(\epsilon \mu )=0`$ due to the finite size effects (for the data presented in Fig.1 relation (27) gives $`|\epsilon \mu |<0.011`$), so it is plausible to suggest that the “universal” behavior of $`g_a(\epsilon \mu )`$ is governed by the finite-size effects. This is clearly demonstrated in Fig.2 where $`g_a(\epsilon \mu )`$ are shown for several sizes of the samples investigated.

The $`\epsilon `$ window where finite size effects are severe, shrinks considerably with increasing $`N`$ for all values of $`\mathrm{\Delta }`$ we investigated. For instance, $`g_a(\epsilon \mu )`$ for $`N=500`$ and $`N=1000`$ at $`\mathrm{\Delta }=0`$ (see Fig.2a,c) merge when $`|\epsilon \mu |0.2`$ while corresponding curves for $`N=1000`$ and $`N=1500`$ are indistinguishable already at $`|\epsilon \mu |0.1`$. The statistical noise observed for the curves in Fig.2 is quite small even close to $`\mu `$ and hence, the influence of insufficient large statistics on the results obtained is excluded. Note, that the “universal” behavior of $`g(\epsilon )`$ in the vicinity of $`\mu `$ obtained for the classical $`da`$ model (see Fig.3 in Ref.) is most likely due to the finite size effects as well.

In the region $`|\epsilon \mu |0.2`$, where the curves for all $`N`$ collapse into a single curve (and where we believe the thermodynamic limit is reached), the behavior of $`g_a(\epsilon \mu )`$ is described by a power law $`g_a(\epsilon \mu )|\epsilon \mu |^\gamma `$. The deviation from the power-law observed far away from $`\mu `$ ($`|\epsilon \mu |0.7`$) is due to the boundaries of the Coulomb gap which, as was mentioned above, are $`1`$ in units of $`E_0`$. One can see from a comparison of the data shown in Fig.2 for different $`\mathrm{\Delta }`$, that the exponent $`\gamma `$ depends considerably on $`\mathrm{\Delta }`$. Furthermore, values of $`\gamma `$ in the region $`\epsilon \mu >0`$ and those in the region $`\epsilon \mu <0`$ differ as well with this difference increasing with increasing $`\mathrm{\Delta }`$. The data for $`\gamma `$ obtained for the two-dimensional MCDAM are summarized in Fig.3 where a significant deviation of $`\gamma `$ from the value $`D1`$ predicted by the hypothesis (1) is observed at all values of $`\mathrm{\Delta }`$ investigated except for the case $`\mathrm{\Delta }=0`$ when $`\gamma 1`$ within the limits of statistical accuracy. Note, that the deviation of $`\gamma `$ from its predicted value grows monotonically with increasing $`\mathrm{\Delta }`$. At $`\mathrm{\Delta }=10`$ where the features of the MCDAM are expected to be nearly the same as those of the classical $`da`$ model with all the acceptors being ionized (indeed, the degree of the acceptor ionization $`C_a0.9`$ for the two-dimensional MCDAM at $`\mathrm{\Delta }=10`$, see Fig.6 below) the deviation from the Efros exponent is very large.

The main results for $`g_a(\epsilon \mu )`$ obtained for the three-dimensional MCDAM are summarized in Figs. 4 and 5. It is seen, that the behavior of $`g_a(\epsilon \mu )`$ in three dimensions does not differ qualitatively from the behavior of $`g_a(\epsilon \mu )`$ in two dimensions. Some quantitative differences observed arise from the fact that at given $`N`$ (the parameter which determines the amount of computer memory needed for the calculations) the linear size of a two-dimensional sample with a given density of sites is larger than that of a three-dimensional sample with the same density of sites and thereby, the finite size effects for three-dimensional samples with given $`N`$ are more pronounced compared to those for the two-dimensional samples with the same $`N`$. For example, the lower boundary of the region where $`g_a(\epsilon \mu )`$ can be described by the power law $`|\epsilon \mu |^\gamma `$ shifts towards larger $`|\epsilon \mu |0.4`$ values (see inserts in Fig.4). Remarkably, the exponent $`\gamma `$ does not reach the value $`D1`$ predicted by the universality hypothesis (1) even at $`\mathrm{\Delta }=0`$ (Fig.5).

Unlike $`g_a(\epsilon \mu )`$ in the vicinity of the Coulomb gap, the density of ionized acceptors $`C_a`$ (15) describes the state of the entire sample and therefore reaches the thermodynamic limit much faster than $`g_a(\epsilon \mu )`$. This allows us to obtain quite accurate results for $`C_a`$ from data on a relatively small amount of samples with $`N=500`$ only. Fig.6 shows the variations of $`C_a`$ with $`\mathrm{\Delta }`$ both for two and three dimensions. In three dimensions almost all acceptors become ionized ($`C_a1`$) rather soon while for two dimensions even for the largest $`\mathrm{\Delta }`$ investigated around 10 % of the acceptors remain neutral. So, one can say, that the three-dimensional MCDAM at $`\mathrm{\Delta }7`$ reduces already to the classical $`da`$ model. It is known that the classical $`da`$ model exhibits in three dimensions the, so called, Coulomb fluctuational catastrophe . For calculations on finite samples it implies that statistical fluctuations of $`\mu (𝐑)`$ grow dramatically with increasing $`\mathrm{\Delta }`$ which is the case in our calculations (see Table I). Therefore, in order to reduce the statistical noise in three dimensions, the average of $`g_a(\epsilon \mu )`$ over a much larger (compared to $`D=2`$) number of samples is needed. Note, that $`\mu (𝐑)`$ in both two and three dimensions are scattered according to the Gaussian distribution with the mean $`\overline{\mu }`$ obeying the relation (25).

## V DISCUSSION

The behavior of $`g_a(\epsilon \mu )`$ calculated within the region of the Coulomb gap for the model (5) is in strong contradiction to the universality hypothesis (1). Despite the fact that $`g_a(\epsilon \mu )`$ is indeed described by the power law $`|\epsilon \mu |^\gamma `$ in a wide range of $`\epsilon `$ inside the region of the Coulomb gap, the exponent $`\gamma `$ is considerably smaller than that predicted by the hypothesis (1) both for the two- and three-dimensional cases. Moreover, the exponent $`\gamma `$ depends significantly on $`\mathrm{\Delta }`$ and is different for the cases $`\epsilon >\mu `$ and $`\epsilon <\mu `$. It is believed that information about $`g(\epsilon )`$ might be directly obtained from tunneling and photoemission experiments  and recent experiments  on boron-doped silicon crystals have shown that the density of one-electron excitations at higher energies obeys a power-law with an exponent slightly less than $`0.5`$ which is in good agreement with our results for $`D=3`$ and $`\mathrm{\Delta }8`$. However, the non-metallic samples show around the Fermi energy a nearly quadratic Coulomb gap, so the question arises whether our results could be related to the intermediate asymptotic behavior observed? Here we want to make three remarks concerning this question:

First, the power law $`g_a(\epsilon \mu )|\epsilon \mu |^\gamma `$ is valid above a value $`\epsilon _0(N)`$ below which the finite size effects take over (Figs.2 and 4). It seems from our results, that $`\epsilon _0(N)\mu `$ when $`N\mathrm{}`$. In two dimensions we were able to obtain size-independent results down to $`\epsilon _00.1`$, i. e. for $`90\%`$ of the whole Coulomb gap, the halfwidth of which is $`1`$ in units of $`E_0`$.

Secondly, as follows from the ground-state stability relations (2), the distance $`r_{ij}`$ between a neutral donor, with an energy, say, $`\epsilon _i^1[\epsilon ,0]`$ ($`\epsilon `$ here is the halfwidth of a narrow band around $`\mu =0`$) and a charged donor with an energy $`\epsilon _j^0[0,\epsilon ]`$ should be not less that $`\frac{1}{2\epsilon }`$. I. e., sites with energies $`\epsilon _i^1[\epsilon ,0]`$ cannot be inside a D-dimensional sphere of radius $`R_{sp}=\frac{1}{2\epsilon }`$ and with the center in a site with the energy $`\epsilon _j^0[0,\epsilon ]`$. Assuming that all such spheres do not intersect, the total volume occupied by the spheres is

$$V_{sp}=N\times S(D)\left(\frac{1}{2\epsilon }\right)^D_0^\epsilon g(\epsilon ^{})𝑑\epsilon ^{}$$
(29)
where $`S(D)`$ is the volume of a D-dimensional sphere with the radius equal to unity. Since $`V_{sp}`$ cannot exceed the total volume $`V`$ of a sample ($`V=N`$ at $`n=1`$) we arrive at the inequality

$$_0^\epsilon g(\epsilon ^{})𝑑\epsilon ^{}\frac{(2\epsilon )^D}{S(D)},$$
(30)
which is valid for all $`\epsilon `$ if

$$g(\epsilon )\frac{D\times 2^D}{S(D)}|\epsilon |^{D1}$$
(31)
The universality hypothesis (1) then is a limit case of (31). The density of sites with energies $`\epsilon _i^1[\epsilon ,0]`$ indeed decreases when $`\epsilon 0`$, so the assumption (29) for the spheres with finite radii seems to be plausible. However, simultaneously $`R_{sp}\mathrm{}`$ and consequently the plausibility of the assumption (29) and thereby of the hypothesis (1) becomes questionable.

And finally, the universality hypothesis (1) can be also obtained as the asymptotic behavior of a non-linear integral equation for $`g(\epsilon )`$ as $`\epsilon 0`$, the equation which, in turn, is heuristically obtained from the stability condition (2). The derivation of this integral equation (given, for example, in Ref.) is based on the implicit assumption that the sites with charged donors are randomly distributed in space according to the Poisson statistics. However, it was unequivocally demonstrated in computer studies of the Coulomb gap  that charged donor sites with energies close to $`\mu `$ tend to form clusters (Ref., Fig. 6).

We conclude that $`g_a(\epsilon \mu )`$ in the region of the Coulomb gap in model (5) has a power law behavior for all energies down to $`\mu `$ and that the universality hypothesis of Efros (1) is questionable. Note, that our results are in contradiction not only to the universality hypothesis (1), but to the inequality (31) as well. Up to now, all exponents found are in good agreement with this inequality. E. g. in Ref. specimens of 40 000 and 125 000 sites for two- and three-dimensional samples were investigated in the Efros’ lattice model and the power law $`g_a(\epsilon \mu )|\epsilon \mu |^\gamma `$ was found with $`\gamma =1.2\pm 0.1`$ and $`\gamma =2.6\pm 0.2`$ for two and three dimensions, respectively. The main conclusion

of our results and those of Ref. is that Efros’ lattice model can not be used to as a reliable approximation to the classical $`da`$ model with Poisson impurity distribution. Energy levels of donor (acceptor) impurities are usually close to the bottom (top) of the conduction (valence) band. Since in the most common semiconductors the energy gap $`E_g10^4`$ K and $`E_020`$ K, $`\mathrm{\Delta }1`$ and one may ask what physical relevance does the model (5) with a finite $`\mathrm{\Delta }10`$ have, except for being a pure academic exercise? However, in the case of deep impurities the energy levels for some donor–acceptor pairs are extremely close to each other not excluding even the case $`\mathrm{\Delta }=0`$. Table II shows some donor–acceptor pairs with $`\mathrm{\Delta }10`$ in the most common semiconductors. The solubilities of these impurities are rather low, thereby reducing the temperature at which the Coulomb gap with features described by the model (5) can be observed. Fortunately, these temperatures are high enough ($`10÷20`$ K) for modern experimental techniques and hence experimental observation of the Coulomb gap in the semiconductors with deep impurities is possible to accomplish.

## VI SUMMARY

We have studied a model of impurities in semiconductors with infinite-range Coulomb interactions between donors, between acceptors and between donors and acceptors. A new parameter introduced in the model is the finite energy $`\mathrm{\Delta }`$ of charge transfer between donors and acceptors, a parameter which enables processes of ionization of neutral impurities and of recombination of charged impurities. In the particular case of equal amounts of donor and acceptor impurities, we derived rigorous relations for the symmetry of the model with respect to exchange of donor and acceptor sites. We also extended the previously known algorithm to find the ground state including the stability relations with respect to ionization and recombination processes and performed computer studies of the model proposed at zero temperature on a number of two- and three-dimensional samples with randomly distributed $`N`$ donors and $`N`$ acceptors. We explored the energy region around the Fermi energy $`\mu `$ where the Coulomb gap in the density of one-electron excitations $`g(\epsilon )`$ is observed. The analysis of the calculated histograms $`g(\epsilon )`$ revealed that the behavior of $`g(\epsilon )`$ obtained from the simulations on finite samples in the immediate neighborhood of $`\mu `$ is determined solely by the finite size effects. In the region where finite size effects become negligible $`g(\epsilon )`$ is described by a power law with an exponent considerably depending on the parameter $`\mathrm{\Delta }`$ and on the sign of $`\epsilon \mu `$. Our findings challenge the Efros universality hypothesis. Moreover, our results are in contradiction to the main inequality (31) of which Efros’ universality hypothesis is a particular case. We have reexamined the heuristic derivation of the Efros hypothesis and shown that some implicit assumptions which lead to universality are questionable. From the analysis of experimental data on admixtures in semiconductors we put forward possible experimental situations where one could observe the Coulomb gap with the features being the same as those of the model with a finite $`\mathrm{\Delta }`$.

## Acknowledgements

This research was supported by The Swedish Natural Science Council and by The Swedish Royal Academy of Sciences.


text:
# Black holes on the brane

## I Introduction

Recent developments in string theory have shown that if matter fields are localized on a 3-brane in $`1+3+d`$ dimensions, while gravity can propagate in the $`d`$ extra dimensions, then the extra dimensions can be large (see, e.g., ). The extra dimensions need not even be compact, as in the 5-dimensional warped space models of Randall and Sundrum . (See also  for earlier work.) In particular, they showed that it is possible to localize gravity on a 3-brane when there is one infinite extra dimension.

If matter on a 3-brane collapses under gravity, without rotating, to form a black hole, then the metric on the brane-world should be close to the Schwarzschild metric at astrophysical scales in order to preserve the observationally tested predictions of general relativity. Collapse to a black hole in the Randall-Sundrum brane-world scenario was studied by Chamblin et al.  (see also ). They gave a ‘black string’ solution which intersects the brane in a Schwarzschild solution.

Here we give an exact localized black hole solution, which remarkably has the mathematical form of the Reissner-Nördstrom solution, but without electric charge being present. Instead the Reissner-Nördstrom-type correction to the Schwarzschild potential can be thought of as a ‘tidal charge’, arising from the projection onto the brane of free gravitational field effects in the bulk. These effects are transmitted via the bulk Weyl tensor (see below). The Schwarzschild potential $`\mathrm{\Phi }=M/(M_\mathrm{p}^2r)`$, where $`M_\mathrm{p}`$ is the effective Planck mass on the brane, is modified to

$$\mathrm{\Phi }=\frac{M}{M_\mathrm{p}^2r}+\frac{Q}{2r^2},$$
(1)
where the constant $`Q`$ is a ‘tidal charge’ parameter, which may be positive or negative.

A geometric approach to the Randall-Sundrum scenario has been developed by Shiromizu et al.  (see also ), and proves to be a useful starting point for formulating the problem and seeing clear lines of approach. The field equations in the bulk are (modifying the notation of )

$$\stackrel{~}{G}_{AB}=\stackrel{~}{\kappa }^2\left[\stackrel{~}{\mathrm{\Lambda }}\stackrel{~}{g}_{AB}+\delta (\chi )\left(\lambda g_{AB}+T_{AB}\right)\right],$$
(2)
where the tildes denote bulk quantities. The fundamental 5-dimensional Planck mass $`\stackrel{~}{M}_\mathrm{p}`$ enters via $`\stackrel{~}{\kappa }^2=8\pi /\stackrel{~}{M}_\mathrm{p}^3`$. The brane tension is $`\lambda `$, and $`\stackrel{~}{\mathrm{\Lambda }}`$ is the bulk cosmological constant. The brane is located at $`\chi =0`$ (so that $`x^4=\chi `$ is a natural choice for the fifth dimension coordinate), and $`g_{AB}=\stackrel{~}{g}_{AB}n_An_B`$ is the induced metric on the brane, with $`n_A`$ the spacelike unit normal to the brane. The brane energy-momentum tensor is $`T_{AB}`$, and $`T_{AB}n^B=0`$. The brane is a fixed point of the $`Z_2`$ symmetry.

## II Field equations on the brane

The field equations induced on the brane arise from Eq. (2), the Gauss-Codazzi equations and the matching conditions with $`Z_2`$-symmetry, and they may be written as a modification of the standard Einstein equations, with the new terms carrying bulk effects onto the brane :

$$G_{\mu \nu }=\mathrm{\Lambda }g_{\mu \nu }+\kappa ^2T_{\mu \nu }+\stackrel{~}{\kappa }^4S_{\mu \nu }_{\mu \nu },$$
(3)
where $`\kappa ^2=8\pi /M_\mathrm{p}^2`$. The energy scales are related to each other and to the cosmological constants via

$$M_\mathrm{p}=\sqrt{\frac{3}{4\pi }}\left(\frac{\stackrel{~}{M}_\mathrm{p}^2}{\sqrt{\lambda }}\right)\stackrel{~}{M}_\mathrm{p},\mathrm{\Lambda }=\frac{4\pi }{\stackrel{~}{M}_\mathrm{p}^3}\left[\stackrel{~}{\mathrm{\Lambda }}+\left(\frac{4\pi }{3\stackrel{~}{M}_\mathrm{p}^3}\right)\lambda ^2\right].$$
(4)
Typically, the fundamental Planck scale is much lower than the effective scale in the brane-world: $`\stackrel{~}{M}_\mathrm{p}M_\mathrm{p}`$. Local bulk effects on the matter are transmitted via the ‘squared energy-momentum’ tensor $`S_{\mu \nu }`$, but since we will consider vacuum solutions, the precise form of $`S_{\mu \nu }`$ (see ) will not be needed. In vacuum, $`T_{\mu \nu }=0=S_{\mu \nu }`$, and we also choose the bulk cosmological constant to satisfy $`\stackrel{~}{\mathrm{\Lambda }}=4\pi \lambda ^2/3\stackrel{~}{M}_\mathrm{p}^3`$, so that $`\mathrm{\Lambda }=0`$ by Eq. (4). Then Eq. (3) reduces to

$$R_{\mu \nu }=_{\mu \nu },R_\mu {}_{}{}^{\mu }=0=_\mu {}_{}{}^{\mu },$$
(5)
where $`_{\mu \nu }`$ is the limit on the brane of the projected bulk Weyl tensor :

$$_{AB}=\stackrel{~}{C}_{ACBD}n^Cn^D.$$
(6)
The Weyl symmetries ensure that this is symmetric and tracefree ($`_{[AB]}=0=_A^A`$) and has no orthogonal components ($`_{AB}n^B=0`$, so that $`_{AB}_{\mu \nu }\delta _A{}_{}{}^{\mu }\delta _{B}^{}^\nu `$ as $`\chi 0`$). It carries the influence of nonlocal gravitational degrees of freedom in the bulk onto the brane, including the tidal (or Coulomb) and transverse traceless (gravitational wave) aspects of the free gravitational field. (See  for a fuller discussion of $`_{\mu \nu }`$ in the general case.) On the brane, in the vacuum case, this tensor satisfies the divergence constraint

$$^\mu _{\mu \nu }=0,$$
(7)
where $`_\mu `$ is the brane covariant derivative. In view of the Bianchi identities on the brane, this is an integrability condition for the field equation $`R_{\mu \nu }=_{\mu \nu }`$. For static solutions, Eqs. (5) and (7) form a closed system of equations on the brane .

A vacuum solution outside a mass localized on the brane must satisfy equations (5) and (7). This leads to a prescription for mapping 4-dimensional general relativity solutions to brane-world solutions in 5-dimensional gravity: a stationary general relativity solution with tracefree energy-momentum tensor gives rise to a vacuum brane-world solution in 5-dimensional gravity. The 4-dimensional general relativity energy-momentum tensor $`T_{\mu \nu }`$ (where $`T_\mu {}_{}{}^{\mu }=0`$) is formally identified with the bulk Weyl term on the brane via the correspondence

$`\kappa ^2T_{\mu \nu }_{\mu \nu }.`$

The general relativity conservation equations $`^\nu T_{\mu \nu }=0`$ correspond to the constraint equation (7) on the brane. In particular, Einstein-Maxwell solutions in general relativity will lead to vacuum brane-world solutions. This is the observation that led us to the Reissner-Nördstrom-type solution.

## III Solutions with tidal charge

Algebraic symmetry properties imply that in general we can decompose $`_{\mu \nu }`$ irreducibly with respect to a chosen 4-velocity field $`u^\mu `$ as

$$_{\mu \nu }=\left(\frac{\stackrel{~}{\kappa }}{\kappa }\right)^4\left[𝒰\left(u_\mu u_\nu +\frac{1}{3}h_{\mu \nu }\right)+𝒫_{\mu \nu }+2𝒬_{(\mu }u_{\nu )}\right],$$
(8)
where $`h_{\mu \nu }=g_{\mu \nu }+u_\mu u_\nu `$ projects orthogonal to $`u^\mu `$. Here

$`𝒰=\left({\displaystyle \frac{\kappa }{\stackrel{~}{\kappa }}}\right)^4_{\mu \nu }u^\mu u^\nu `$

is an effective energy density on the brane arising from the free gravitational field in the bulk—but note that this energy density need not be positive. Indeed, as we argue below, $`𝒰<0`$ is the natural case. The effective anisotropic stress from the free gravitational field in the bulk is the spatially tracefree and symmetric part, i.e.,

$`𝒫_{\mu \nu }=\left({\displaystyle \frac{\kappa }{\stackrel{~}{\kappa }}}\right)^4\left[h_{(\mu }{}_{}{}^{\alpha }h_{\nu )}^{}{}_{}{}^{\beta }\frac{1}{3}h_{\mu \nu }h^{\alpha \beta }\right]_{\alpha \beta },`$

where round brackets denote symmetrization. The effective energy flux from the free gravitational field in the bulk is

$`𝒬_\mu =\left({\displaystyle \frac{\kappa }{\stackrel{~}{\kappa }}}\right)^4h_\mu {}_{}{}^{\alpha }_{\alpha \beta }^{}u^\beta .`$

In a static vacuum, with $`u^\mu `$ along the static Killing direction, we have $`𝒬_\mu =0`$, and the effective conservation equation (7) reduces to the single spatial equation

$$\frac{1}{3}\mathrm{D}_\mu 𝒰+\frac{4}{3}𝒰A_\mu +\mathrm{D}^\nu 𝒫_{\mu \nu }+A^\nu 𝒫_{\mu \nu }=0,$$
(9)
where $`\mathrm{D}_\mu `$ is the projection (orthogonal to $`u^\mu `$) of the brane covariant derivative , and $`A_\mu =u^\nu _\nu u_\mu `$ is the 4-acceleration. Static spherical symmetry means that

$`A_\mu =A(r)r_\mu ,𝒫_{\mu \nu }=𝒫(r)\left[r_\mu r_\nu \frac{1}{3}h_{\mu \nu }\right],`$

for some functions $`A(r)`$ and $`𝒫(r)`$, where $`r`$ is the areal distance and $`r_\mu `$ is a unit radial vector. The Reissner-Nördstrom-type solution in Eq. (1) corresponds to the solution

$$𝒰=\left(\frac{\kappa }{\stackrel{~}{\kappa }}\right)^4\frac{Q}{r^4}=\frac{1}{2}𝒫$$
(10)
of Eq. (9).

We can verify that the solution in Eqs. (1) and (10) satisfies the field equations (5), using natural coordinates, for which the metric on the brane is

$$ds^2=A(r)dt^2+B(r)dr^2+r^2(d\theta ^2+\mathrm{sin}^2\theta d\phi ^2).$$
(11)
Then it may be verified that

$`A=B^1=1+{\displaystyle \frac{\alpha }{r}}+{\displaystyle \frac{\beta }{r^2}},`$ (12)
$`_t{}_{}{}^{t}=_r{}_{}{}^{r}=_\theta {}_{}{}^{\theta }=_\phi {}_{}{}^{\phi }={\displaystyle \frac{\beta }{r^4}},`$ (13)
where $`\alpha `$ and $`\beta `$ are constants. Equations (12) and (13) satisfy all the field equations (5), and hence also the divergence equations (7). By Eqs. (8) and (10), we see that $`\beta =Q`$, and the far-field Newtonian limit \[see Eq. (1)\] shows that $`\alpha =2M/M_\mathrm{p}^2`$.

In summary, we have shown that an exact black hole solution of the effective field equations on the brane is given by the induced metric

$$g_{tt}=(g_{rr})^1=1\left(\frac{2M}{M_\mathrm{p}^2}\right)\frac{1}{r}+\left(\frac{q}{\stackrel{~}{M}_\mathrm{p}^2}\right)\frac{1}{r^2},$$
(14)
where $`q=Q\stackrel{~}{M}_\mathrm{p}^2`$ is a dimensionless tidal charge parameter. The projected Weyl tensor, transmitting the tidal charge stresses from the bulk to the brane, is

$$_{\mu \nu }=\left(\frac{q}{\stackrel{~}{M}_\mathrm{p}^2}\right)\frac{1}{r^4}\left[u_\mu u_\nu 2r_\mu r_\nu +h_{\mu \nu }\right].$$
(15)
## IV Properties of the black hole

The 4-dimensional horizon structure of the brane-world black hole depends on the sign of $`q`$. For $`q0`$, there is a direct analogy to the Reissner-Nördstrom solution, with two horizons:

$$r_\pm =\frac{M}{M_\mathrm{p}^2}\left[1\pm \sqrt{1q\frac{M_\mathrm{p}^4}{M^2\stackrel{~}{M}_\mathrm{p}^2}}\right].$$
(16)
As in general relativity, both horizons lie inside the Schwarzschild horizon: $`0r_{}r_+r_\mathrm{s}=2M/M_\mathrm{p}^2`$, and there is an upper limit on $`q`$:

$$0qq_{\mathrm{max}}=\left(\frac{\stackrel{~}{M}_\mathrm{p}}{M_\mathrm{p}}\right)\left(\frac{M}{M_\mathrm{p}}\right)^2.$$
(17)
The intriguing new possibility that $`q<0`$, which is impossible in the general relativity Reissner-Nördstrom case, leads to only one horizon, lying outside the Schwarzschild horizon:

$$r_+=\frac{M}{M_\mathrm{p}^2}\left[1+\sqrt{1q\frac{M_\mathrm{p}^4}{M^2\stackrel{~}{M}_\mathrm{p}^2}}\right]>r_\mathrm{s}.$$
(18)
In the $`q<0`$ case, the (single) horizon has a greater area than its Schwarzschild counterpart, so that bulk effects act to increase the entropy and decrease the temperature of the black hole. In general relativity, the electric field in the Reissner-Nördstrom solutions acts to weaken the gravitational field, and the same is true for the brane-world black hole with $`q>0`$. This can be seen clearly from Eq. (1). By contrast, the $`q<0`$ case corresponds to the opposite effect, i.e., bulk effects tend to strengthen the gravitational field. By Eq. (10), we see that in this case, the effective energy density $`𝒰`$ on the brane contributed by the free gravitational field in the bulk is negative. This is in accord with the (Newtonian) notion that the gravitational field of an isolated mass has negative energy density. Furthermore, it agrees with the perturbative analysis by Sasaki et al.  and the nonperturbative analysis of Maartens . The gravitational field generated by a source on the brane tends to squeeze test matter in the 5th direction, thus acting as an attractive field with a negative energy contribution. The tidal acceleration measured by static observers along the 5th direction $`n^A`$ is

$`\stackrel{~}{R}_{ABCD}u^An^Bu^Cn^D=\left({\displaystyle \frac{\stackrel{~}{\kappa }}{\kappa }}\right)^4𝒰+\frac{1}{6}\stackrel{~}{\kappa }^2\stackrel{~}{\mathrm{\Lambda }},`$

where the right hand side follows from Eqs. (2) and (6). The negative bulk cosmological constant contributes to acceleration towards the brane, reflecting its confining role on the gravitational field. In order for $`𝒰`$ to reinforce confinement, it must be negative. In other words, negative tidal charge $`q<0`$ is the physically more natural case. (See also  for further discussion of negative energy density from bulk effects.) Furthermore, $`q<0`$ ensures that the singularity is spacelike, as in the Schwarzschild solution, whereas $`q>0`$ leads to a timelike singularity, which amounts to a qualitative change in the nature of the general relativistic Schwarzschild solution.

It is widely assumed that astrophysical black holes cannot exhibit macroscopic electric charge due to the presence of neutralizing plasma in their vicinity. Discussions of astrophysical black hole phenomena are therefore commonly restricted to the Kerr family of black holes. The solution presented here raises the possibility that an effective Reissner-Nördstrom metric emerges by geometric considerations and is hence not constrained by the above argument. The tidal charge $`q`$ affects the geodesics and the gravitational potential, so that indirect limits may be placed on it by observations. Current observational limits on $`|q|`$ are rather weak, since the correction term in Eq. (1) dies off rapidly with increasing $`r`$, and astrophysical measurements (lensing and perihelion precession) probe mostly (weak-field) solar scales. These measurements require the correction term in the gravitational potential to be much less than the Schwarzschild term, so that

$$|q|2\left(\frac{\stackrel{~}{M}_\mathrm{p}}{M_\mathrm{p}}\right)^2M_{}R_{}.$$
(19)
This still allows for large values of $`|q|`$, which would modify the spacetime geometry of a nonrotating black hole in the strong-gravity regime , with implications for example for the last stable circular orbit for compact binaries. Although the strong-gravity regime is not currently directly accessible to observations, indirect limits could emerge from the way in which tidal charge modifies general relativity strong-gravity effects. This deserves further investigation.

Further indirect limits could arise from the effect of tidal charge on primordial black holes, which also merits further analysis. An intriguing possibility is that strong tidal effects in the early universe could lead to black hole formation even in the absence of gravitational collapse of matter. If such matter-free tidal collapse is possible, then the endstate is a metric with $`M=0`$ and $`q<0`$:

$$g_{tt}=(g_{rr})^1=1+\left(\frac{q}{\stackrel{~}{M}_\mathrm{p}^2}\right)\frac{1}{r^2},$$
(20)
and the horizon on the brane is given by

$$r_\mathrm{h}=\frac{\sqrt{q}}{\stackrel{~}{M}_\mathrm{p}}.$$
(21)
The tidal charge parameter arises from purely scalar (Coulomb-like) effects of the free gravitational field in the bulk, given the static spherical symmetry. In the original Randall-Sundrum model, as well as in perturbative studies of it (see, e.g., ), the bulk is assumed to be exactly anti-de Sitter in the absence of any source on the brane. The geometric approach of , which we have adopted here, makes no assumptions about the bulk metric, other than that it satisfies the 5-dimensional Einstein equations with cosmological constant. Thus the bulk need not be conformally flat in the absence of sources on the brane. This means that in general, the tidal charge parameter $`q`$ will be determined by both the brane source, i.e., the mass $`M`$, and any Coulomb part of the bulk Weyl tensor that survives when $`M`$ is set to zero.

## V Conclusion

We have not investigated fully the effect of the brane-world black hole on the bulk geometry, and in particular the nature of the off-brane horizon structure. This has been done for solutions which reduce to the Schwarzschild black hole on the brane . In these solutions, the bulk metric is

$`d\stackrel{~}{s}^2=\left({\displaystyle \frac{6}{\stackrel{~}{\kappa }^2\stackrel{~}{\mathrm{\Lambda }}}}\right){\displaystyle \frac{1}{z^2}}\left[dz^2+g_{\mu \nu }(x^\alpha )dx^\mu dx^\nu \right],`$

where $`g_{\mu \nu }`$ is the Schwarzschild metric. We have adopted a different approach: instead of starting from an induced Schwarzschild metric, we have solved the effective field equations for the induced metric on the brane (which form a closed system, since the metric is static), and found a generalization of the Schwarzschild solution. It turns out that if $`ds^2`$ is given by our solution, as in Eqs. (11)–(14), then

$`d\stackrel{~}{s}^2=f(z)\left[dz^2+ds^2\right]`$

cannot satisfy the bulk field equations for any $`f(z)`$ if $`q0`$. Finding an exact form of the bulk metric that is consistent with our exact induced metric on the brane is more difficult than the case where the induced metric is Schwarzschild.

Perturbative studies, which start from an exactly anti-de Sitter background, show that the first weak-field correction of the Newtonian potential on the brane is proportional to $`1/r^3`$ (see ). Our solution has by contrast a $`1/r^2`$ correction, so that it is incompatible with the long-distance limit on the brane. (Such a correction can also arise in thick-brane models .) However, in the short-distance limit on the brane, the lowest order correction to the potential is proportional to $`1/r^2`$ . This term will dominate the $`1/r`$ term, which reflects the fact that gravity becomes effectively 5-dimensional at high energies. Thus our solution should describe well the strong-gravity regime on the brane. In the short-distance limit, the perturbative analysis shows that

$`q={\displaystyle \frac{M}{\stackrel{~}{M}_\mathrm{p}}},`$

so that the tidal charge is negative, as we argued above, and is determined by the black hole mass (which is to be expected if the background bulk Weyl tensor vanishes).

In order to pursue our nonperturbative analysis, we need to look at the off-brane equations for the curvature (see  for the general form of these equations). The induced field equations on the brane given in Eq. (5) are supplemented by off-brane equations which determine

$`_𝐧_{AB},_𝐧_{ABC},_𝐧R_{ABCD},`$

where $`_𝐧`$ is the Lie derivative along $`n^A`$, $`R_{ABCD}`$ is the 4-dimensional Riemann tensor, and

$`_{ABC}=g_A{}_{}{}^{D}g_{B}^{}{}_{}{}^{E}\stackrel{~}{C}_{DECF}^{}n^F,`$

with $`_{\mu \nu \sigma }=0`$ on the brane. These equations together with Eq. (5) form a closed system .

In conclusion, we have shown how a Reissner-Nördstrom-type metric satisfies the effective field equations on the brane in Randall-Sundrum-type gravity, which form a closed system because of staticity. There is no electric charge present, but instead a tidal charge, arising from the imprint of the free gravitational field in the bulk on the brane. The tidal charge correction to the Schwarzschild potential is negative, and the solution describes the strong-gravity regime on the brane. Negative $`q`$ leads to an horizon on the brane that is outside the Schwarzschild horizon, corresponding to lower temperature and greater entropy. Current observations place only weak limits on the tidal charge, and in principle significant tidal modifications could arise in the strong-field regime or in the early-universe case of primordial black holes. Further investigation is in progress to probe these modifications and any indirect limits that they may impose on the tidal charge, as well as to find the off-brane behaviour of the horizon.

Acknowledgements: VR was supported by a Royal Society grant while at Portsmouth, and thanks the Relativity and Cosmology Group for hospitality. RM thanks IUCAA, Pune for hospitality during a visit, which was partially supported by the Royal Society. We thank Bruce Bassett, Marco Bruni, Roberto Casadio, Roberto Emparan, David Langlois, Jose Senovilla, Carlo Ungarelli and David Wands for useful discussions and comments, and especially Tetsuya Shiromizu.