text:
# The radial velocity curve of HD153919 (4U1700-37) revisitedBased on observations obtained with the International Ultraviolet Explorer

## 1 Introduction

The bright Of star HD153919 is the optical counterpart of the 3.4 day eclipsing high-mass X-ray binary 4U1700-37, discovered by Jones et al. (Jo73 (1973)). About 2.5 million years ago, the binary system escaped from the OB association Sco OB2 after the supernova explosion that formed the compact companion (Ankay et al. AK01 (2001)). The X-ray source powered by the stellar wind of the Of star shows no regular pulsations and the highly variable X-ray flux is less than expected from standard wind accretion (e.g. Rubin et al. RF96 (1996), Kaper Ka98 (1998)). The decrease in orbital period determined from mid-eclipse observations and first mentioned by Haberl et al. (HW89 (1989)) was confirmed by Rubin et al. (RF96 (1996)) and is probably due to wind-driven angular momentum loss from the Of star. Because no direct measurements can be made of the orbital motion of the compact object, system parameters depend heavily on radial velocity studies of the optical companion. The most recent analysis of the physical parameters cannot resolve the nature of the compact object in this system (Clark et al. CG02 (2002)). Its likely mass is considerably lower than those found for black holes, but it is significantly in excess to fit the high density nuclear equation of state for neutron stars. A point to remember, however, is that also the Vela X-1 pulsar has a high mass, comparable to that of 4U1700-37 (van Kerkwijk et al. vK95 (1995)).

The most thorough study of radial velocities of lines in the visible spectral region has been made by Hammerschlag-Hensberge (HH78 (1978)). It is clear from the observations that the strongly variable stellar wind in this system complicates a straightforward determination of the system parameters. The results for the orbital analysis of all lines together in the visible blue spectral region ($`37004900`$ Å) yield a semi-amplitude $`K=19\pm 1`$ km s-1 and an orbital eccentricity $`e=0.16\pm 0.08`$.

Heap & Corcoran (HC92 (1992)) have examined the system parameters using the available ultraviolet spectra of HD153919. They find $`K=18\pm 3`$ km s-1 and a negligible eccentricity. For their investigation they cross-correlated a few selected narrow regions of the ultraviolet spectra with the same regions of one of the spectra. To bring their spectra to the same frame of reference they used the heliocentric correction. However, instrumental and processing shifts during the reduction of spectra obtained with the International Ultraviolet Explorer (IUE) require a more sophisticated approach (Stickland St92 (1992)).

Stickland & Lloyd (SL93 (1993)) re-examined the same data, now cross-correlating both stellar and interstellar spectra with reference to one of the UV spectra. They used the complete ultraviolet spectral range excluding P-Cygni lines whose profiles are affected by mass loss. Amazingly, they found a much smaller semi-amplitude $`K`$ of 10 km s-1 and an eccentricity of 0.17.

In a similar study of the high-mass X-ray binary HD77581/Vela X-1 Barziv et al. (BK00 (2001)) found a disagreement in the semi-amplitude of the radial velocity curve with the results of Stickland et al. (SL97 (1997)) for this system. This strengthened our decision to re-analyse the ultraviolet spectra of HD153919 again. In Sect. 2 we describe the observations and the data reduction procedures. In Sect. 3 we present the resulting radial velocity curve and in Sect. 4 we compare our results with the ones found by Heap & Corcoran and by Stickland & Lloyd. In Sect. 5 we summarize our conclusions and discuss the evolutionary consequences.

## 2 Observations and data reduction

A large number of ultraviolet spectra of HD153919 have been obtained with the International Ultraviolet Explorer. The spectra obtained with the Short Wavelength Prime Camera in high-resolution mode are listed in Table 1. A separate analysis of the same IUE spectra was performed by Stickland & Lloyd (SL93 (1993)). These authors cross-correlated the set of spectra with one particular spectrum (SWP25600). Because their results were in large discrepancy with results from optical data we decided to re-analyse the same IUE data. We cross-correlated each spectrum with all other spectra of the data set to reduce errors. We did not use a template spectrum of another star with identical spectral type as independent reference spectrum- this is not realistic for such hot mass-losing stars which each have unique spectra- nor did we use an average spectrum to avoid auto-correlation effects leading to a systematic decrease of the amplitude of the radial velocity curve.

First the spectra - initially mapped onto a uniform wavelength grid of $`0.05`$ Å - were transformed to a logarithmic wavelength scale such that the Doppler shift becomes a linear displacement along the spectrum. The relevant part of the cross-correlation function is the area close to the top; it is this upper part which is fitted with an analytic function (a Gaussian plus a linear function) whose maximum determines the velocity shift of the spectrum. Because the cross-correlation profiles of IUE spectra show small excess peaks (see Evans (Ev88 (1988)) and Fig. 4 in Van Kerkwijk et al. (VK95 (1995))) due to the so-called fixed-pattern noise in the detector, we corrected for those by including an extra Gaussian peak in the fit. For details on the cross-correlation method and instructive figures we refer to Van Kerkwijk (VK93 (1993)).

The spectral regions which were used for the cross-correlation are listed in Table 2. The regions were chosen in a way to exclude parts of the spectrum containing variable line profiles such as the well-known P-Cygni profiles (see for instance Hammerschlag-Hensberge et al. (HHK90 (1990))); also the spectral region including variable emission lines $`15401740`$ Å (Kaper et al. KH90 (1990)) was excluded. Table 2 lists the wavelength regions including interstellar lines which were used to place the individual spectra onto the same frame of reference, because IUE spectra have a variety of instrumental and processing shifts, as mentioned before.

## 3 Results

In Table 1 we list the radial velocities derived from the cross-correlations of all spectra with each other. In Fig. 1, we show the radial velocities we obtained from the IUE spectra as a function of orbital phase. Because the spectra were cross-correlated with each other, no $`\gamma `$-velocity has been determined for the system and the velocities in the figure have been shifted to zero system velocity. The $`\gamma `$-velocity has been determined elsewhere (e.g. Gies & Bolton GB86 (1986), Humphreys Hu78 (1978)): $`60`$ km s-1. It was not included in our figure because the value - which is rather uncertain - did not follow from our data. At first glance, it is already clear that the data show significant deviations from a smooth Keplerian curve. We believe these are intrinsic to the star, since earlier optical observations have already shown that the star has intrinsic variations (e.g., Kaper et al. KH94 (1994)). As we will find that the deviations are substantially larger than the measurement errors, we give each measurement equal weight in our fits below.

We first fit a circular orbit, shown by the long-dashed line in Fig. 1, and find $`K=18.7\pm 1.0`$km s-1. The root-mean-square (rms) residual is $`6.36`$ km s-1. This is a factor two larger than the mean error in the measured velocities (using the formal errors, one finds a reduced $`\chi ^2=3.6`$).

The residuals relative to the best-fit circular orbit are shown in the lower panel of Fig. 1. These seem to show a double-wave pattern, suggestive of an eccentric orbit. Fitting an eccentric orbit, we find eccentricity $`e=0.22\pm 0.04`$, periastron angle $`\omega =49\mathrm{°}\pm 11\mathrm{°}`$, and the radial-velocity amplitude is $`K=20.6\pm 1.0`$ km s-1. The rms residuals are slightly lower, at $`5.85`$ km s-1. Doing an F-test, this is not a significant improvement. The F-test, however, ignores the correlations in the residuals. In order to obtain a more robust estimate of the significance, we ran Monte-Carlo simulations, in which we made artificial data sets in which the velocity at each measurement phase was given by the sum of the velocity expected from the best-fit circular orbit and a random deviation. For the latter, we used two components, one to mimic the measurement error (i.e., a random number drawn from a normal distribution with standard deviation equal to the measurement error), and one to mimic the intrinsic deviations. For the intrinsic deviations, we again assumed a normal distribution, with a standard deviation of 6 km s-1, so that our artificial data sets have the same rms residuals as our data. Our simulations confirmed our error estimate on the eccentricity, and showed that the probability of obtaining an eccentricity of 0.22 by chance due to the deviations was negligible.

We should note that what is proven above is that for a circular orbit there are systematic variations with orbital phase over and above the significant random excursions. An eccentric orbit is a simple explanation for these systematic variations, but we cannot exclude other systematic effects. A cause for caution is that Barziv et al. (BK00 (2001)) found from their optical data an eccentricity which was larger than the value from the X-ray pulsar in Vela X-1 (see their Fig. 18). Furthermore, model calculations of Zuiderwijk et al. (vP77 (1977)) showed that the radial-velocity curve for a tidally distorted star in a circular orbit might show an apparent eccentricity. The shape of the expected curve, however, is different from that observed. Given this, we cannot dismiss the possibility of a circular orbit out of hand, although we believe that the simplest and most likely solution is that it is in fact eccentric.

It is obvious that for some spectra the individual radial-velocity measurements deviate significantly from the best-fit Keplerian orbit. This was already known from the optical spectra. The Of supergiant has a strong stellar wind which can cause irregular variations in the radial velocities of some lines. Also the presence of a photo-ionization wake (Kaper et al. KH94 (1994)) trailing the X-ray source may cause irregular small orbit-to-orbit variations in the radial velocities. Nevertheless we devoted extra attention to some of the spectra. Particularly striking are the three velocities near phase 0.25 that are much lower than the other spectra at the same phase interval. These three are for spectra SWP1982, 1983 and 1986, which were taken in sequence. This would suggest an instrumental problem, were it not that SWP1987, which was taken immediately after SWP1986 at the same night as the mentioned three spectra, has a much higher velocity comparable with that of SWP 4751, 4752 and 4753. We compared the line profiles of all those spectra but could not find significant differences, although SWP 1982 and 1983 are very noisy. Note that the modest spectral resolution ($`R10,000`$) and low S/N ($`20`$) does not allow for a detailed study of individual photospheric lines.

For HD 77581, the companion of Vela X-1, it was found that sequences of radial velocities showed deviations that were correlated over periods of about one day (Van Kerkwijk et al. VK95 (1995)). To see whether such deviations are present in HD 153919 as well, in Fig. 1 lines connect residual velocities for spectra that were taken as part of a sequence. Unfortunately, the result is inconclusive. Some groups seem to show correlated deviations, but for others the deviations seem random. It would be certainly worthwhile to monitor HD 153919 during some nights continuously to have a clue to the origin of those excursions and stellar wind variations that certainly influence and hamper an accurate determination of the orbital parameters.

Finally, we should mention that almost half of the set of spectra was taken in the small aperture of the SWP camera, the others in the large aperture (see Table 1). The spectra taken in the small aperture are the ones of the first years of IUE and have generally lower S/N. We determined orbital solutions for the small and large-aperture sets separately, but did not find significant differences.

## 4 Comparison with earlier results

From the present results we can conclude that the derived radial velocities are in good agreement with the results obtained much earlier from optical data (see Hammerschlag-Hensberge HH78 (1978)).

Our results do not differ much from those by Heap & Corcoran (HC92 (1992)) for the same data. The advantage of our analysis is that we use a larger part of the spectrum and correlate all spectra with each other, which, as a consequence, results in a smaller error on $`K`$. Remarkably, Heap & Corcoran (HC92 (1992)) find a negligible eccentricity ($`e0.003`$) with a very small error of the order of $`0.003`$. At first glance, the radial velocities in their Fig. 2 seem to fit equally well to an eccentric orbit. Also, from a statistical point of view, their errors on $`e`$ and $`\omega `$ are much too small. For small eccentricities the deviation from a circular orbit is a sinusoid with half the binary period and it is easy to see that $`eA(P/2)/A(P)=A(P/2)/K`$, where $`A`$ is the amplitude of the sinusoid and $`K`$ the radial velocity amplitude. The error on $`A(P/2)`$ will be comparable to the one on $`K`$ and as a consequence $`\sigma _e\sigma _K/K`$ . This implies $`\sigma _e0.17`$ for the error on the eccentricity, instead of $`\sigma _e=0.003`$ found by Heap & Corcoran.

On the other hand, Stickland & Lloyd (SL93 (1993)) find a stellar orbit with an eccentricity close to ours, but their results show a much smaller semi-amplitude $`K`$ of 10 km s-1. Contact by one of the authors with Dave Stickland for similar discrepant velocity values in the X-ray binary Vela X-1 has proven that the main cause for the discrepancy in the results is that their cross-correlation procedure sets to zero the regions around interstellar lines and wind lines. When this is combined with slight inadequacies in normalization, artificial lines are introduced at fixed wavelengths in both spectra to be cross-correlated. As a consequence the radial velocities are biased systematically towards zero and the radial-velocity amplitude will be underestimated. A discrepancy of about 4 km s-1 had been found between the radial-velocity amplitude for the ultraviolet spectra of Vela X-1 analyzed by both groups (see Barziv et al. BK00 (2001) and Stickland et al. SL97 (1997)). The reason for the much larger discrepancy of 10 km s-1 for HD153919 might be the weaker and broader photospheric lines in this star combined with the better signal-to-noise ratio which give more striking errors when interstellar and wind lines are set to zero.

## 5 Conclusions and evolutionary consequences

The two main conclusions of this paper are that (i) the radial-velocity amplitude $`K`$ of HD153919 determined from ultraviolet spectra is consistent with that derived from optical spectra; and (ii) the orbit most likely has a significant eccentricity. Note that the obtained value for $`K`$ depends on whether the orbit is eccentric ($`K=20.6`$ km s-1, $`e=0.22`$) or not ($`K=18.7`$ km s-1).

Clark et al. (CG02 (2002)) performed a thorough study of the system parameters of HD153919/4U1700-37 and determined the most likely masses of both components using Monte Carlo simulations. These suggest a mass for HD153919 and 4U1700-37 of $`\mathrm{M}_{}=58\pm 11`$ M and $`\mathrm{M}_X=2.4\pm 0.3`$ M, respectively. Clark et al. used the results described in this paper in their analysis, but prefered the orbital solution with zero eccentricity. As the mass function depends on the eccentricity:

$$f=1.038\times 10^7K^3P(1e^2)^{3/2}$$
(1)
and

$$(\mathrm{M}_{}+\mathrm{M}_X)^2=\frac{(\mathrm{M}_X\mathrm{sin}i)^3}{f},$$
(2)
they noticed that the masses of both components become significantly higher when taking $`e=0.22`$. They argue that such high values for $`M_{}`$ are inconsistent with $`\mathrm{log}g`$ derived from the analysis of the photospheric spectrum. A non-zero eccentricity also results in a higher number of solutions that are rejected due to the inclination constraints. However, Clark et al. took the value of the radial-velocity amplitude belonging to the eccentric orbit ($`K=20.6`$ km s-1, so that $`f=0.0031`$ M), rather than $`K=18.7`$ km s-1 for the circular orbit ($`f=0.0023`$ M); for the latter value the masses would have become even higher. If one takes $`K=20.6`$ km s-1 and $`e=0.22`$ then $`f=0.0029`$ M. Therefore, we conclude that the masses of both components are about 4 % higher than proposed by Clark et al. (CG02 (2002)) if one uses the correct value for $`f`$ in case of the eccentric orbital solution.

The short orbital period of the system, combined with its eccentricity, will induce strong tidal interactions between the two components. Quaintrell et al. (QN03 (2003)) present evidence for tidally induced non-radial oscillations in HD77581, the B-supergiant companion to Vela X-1, like 4U1700-37 an eccentric ($`e=0.09`$), wind-fed high-mass X-ray binary with an orbital period of 8.9 days. These non-radial oscillations may well be the cause of the radial-velocity excursions detected in Vela X-1. The radial-velocity excursions in HD153919 could have a similar origin.

The eccentricity of the orbit of both Vela X-1 and 4U1700-37 suggests that these systems have just entered the high-mass X-ray binary phase and are on their way to become Roche-lobe overflow systems. The latter systems have circular orbits, possibly the result of tidal interaction. Both HD77581 and HD153919 do not rotate synchronously with their orbit: $`P_{\mathrm{orb}}/P_{\mathrm{rot}}`$ is 0.67 (Zuiderwijk 1995) and 0.46 (taking $`v\mathrm{sin}i=150`$ km s-1, $`i=90^{}`$, $`R_{}=21.9`$ R), respectively, which may indicate that the supergiants are still expanding and that tidal forces have not (yet?) managed to synchronize the system and circularize the orbit.

Conti (Co78 (1978)) and Rappaport & Joss (RJ83 (1983)) suggested that the OB supergiants in HMXBs are too luminous for their mass. Kaper (Ka01 (2001)) showed that the discrepancy between the measured and ”spectroscopic” mass is most severe for the Roche-lobe overflow systems, but is only modest for wind-fed systems such as Vela X-1. The proposed high mass of HD153919 is consistent with its high luminosity, again suggesting that HD153919 has not yet reached the phase of Roche-lobe overflow.

###### Acknowledgements.

LK and MHvK acknowledge support from a fellowship of the Royal Netherlands Academy of Arts and Sciences.


format: text

text:
# Cosmological dark energy from the cosmic QCD phase transition and colour entanglement

## Abstract

Recent astrophysical observations indicate that the universe is composed of a large amount of dark energy (DE) responsible for an accelerated expansion of the universe, along with a sizeable amount of cold dark matter (CDM), responsible for structure formation. At present, the explanations for the origin or the nature of both CDM and DE seem to require ideas beyond the standard model of elementary particle interactions. Here, for the first time, we show that CDM and DE can arise entirely from the standard principles of strong interaction physics and quantum entanglement. Quantitative agreement with the present data obtains without the need for any adjustable parameters.

Shibaji Banerjee1, Abhijit Bhattacharyya2, Sanjay K. Ghosh3,

Ernst-Michael Ilgenfritz4, Sibaji Raha3 ***Corresponding Author Electronic Mail : sibaji@bosemain.boseinst.ac.in, Bikash Sinha5,

Eiichi Takasugi6 and Hiroshi Toki4

1Physics Department, St. Xavier’s College, 30, Park Street, Kolkata 700016, INDIA
2Physics Department, Scottish Church College, 1 & 3, Urquhart Square, Kolkata 700006, INDIA
3Physics Department, Bose Institute, 93/1, A. P. C. Road, Kolkata 700009, INDIA
4Research Center for Nuclear Physics, Osaka University, Ibaraki, Osaka 567-0047, JAPAN
5 Saha Institute of Nuclear Physics, 1/AF, Bidhannagar, Kolkata - 700 064, INDIA
6 Graduate School of Physics, Osaka University, Toyonaka, Osaka 560-0043, JAPAN PACS: 12.38.Mh, 12.90.+b, 14.80.Dq, 96.40.-z

The past few decades have been momentous in the history of science insofar as several accurate astrophysical measurements have been carried out and cosmology, often considered a fair playground of exotic or fanciful ideas, has had to confront the reality of experiments. Based on the knowledge gleaned so far, the present consensus  is that the standard model of cosmology, comprising the Big Bang and a flat universe is correct. The Big Bang Nucleosynthesis  (BBN), which forms one of the basic tenets of the standard model, shows that baryons can at most contribute $`\mathrm{\Omega }_B`$ ( $`\frac{\rho _B}{\rho _c}`$, $`\rho _c`$ being the present value of closure density $`10^{47}`$ Gev4 ) $``$0.04, whereas structure formation studies require that the total (cold) matter density should be $`\mathrm{\Omega }_{CDM}0.23`$. Matter contributing to CDM is characterized by a dust-like equation of state, pressure $`p`$ 0 and energy density $`\rho >0`$ and is responsible for clustering on galactic or supergalactic scales. Dark energy (DE), on the other hand, is smooth, with no clustering features at any scale. It is required to have an equation of state $`p=w\rho `$ where $`w<0`$ (ideally $`w=1`$), so that for a positive amount of dark energy, the resulting negative pressure would facilitate an accelerated expansion of the universe, evidence for which has recently become available from the redshift studies of type IA supernovae . For a flat universe $`\mathrm{\Omega }1`$, $`\mathrm{\Omega }_{DE}0.73`$  implies that $`\rho _{DE}`$ today is of the order of $`10^{48}`$ GeV4.

There is no consensus within the community about the origin or the nature of CDM or DE. It has been argued that given the BBN limit of $`\mathrm{\Omega }_B0.04`$, CDM cannot be baryonic; as a result, various exotic possibilities, all beyond the standard model $`SU(3)_c\times SU(2)\times U(1)`$ of particle interactions, have been suggested. The situation is even more complicated for DE. The most natural explanation for DE would be a vacuum energy density, which a priori would have the correct equation of state ($`w=1`$). This possibility however is beset with the insurmountable difficulty that for any known (or conjectured) type of particle interaction, the vacuum energy density scale turns out to be many orders of magnitude larger than the present critical density. A trivial, but aesthetically displeasing, way out could be to a priori postulate a small cosmological constant. On the other hand, there exists in the recent literature a large number, too many to cite in this letter, of speculative suggestions, all substantially beyond the standard model. Apart from being non-standard, all these pictures require large (perhaps unacceptable) amounts of fine-tuning, as these authors themselves admit. It thus appears that one needs a second set of magic, in addition to the one needed for CDM, to explain DE. In this letter, we show that there can indeed be Magic Without Magic; there is no need to go beyond the standard model $`SU(3)_c\times SU(2)\times U(1)`$ to understand the nature of CDM and DE. They can both arise from the same process of the cosmic quark-hadron phase transition occurring during the microsecond epoch after the Big Bang, provided we admit the existence of quantum entanglement for a strongly interacting system. Even though entanglement in quantum mechanics has been established both theoretically and experimentally, there is till date no evidence of such effects in strong interaction. We, however, postulate that being a quantum mechanical property, it would also occur in strong interactions too; more on this later on. Towards the end of this work, we shall a posteriori argue that it could have been anticipated that the Quantum Chromodynamic vacuum was the natural solution for the dark energy problem.

The role of phase transitions  in the early universe has been recognized to be of paramount importance. In the strong interaction (Quantum Chromodynamics) sector, there is expected to be a distinct phase transition separating the confined (hadronic) phase from the deconfined (quark-gluon plasma) phase. In the early universe, this phase transition is predicted to occur during the microsecond epoch after the Big Bang. The order of this phase transition is at present an open issue. While it may be of second order or even a continuous transition (cross-over) in the laboratory, the cosmic phase transition could most likely be of first order, as has been forcefully argued . In what follows, we shall tacitly assume it to be of first order. We shall return to this issue later on. Other crucial ansätze in our scenario include that the universe is overall colour neutral at all times, which we believe nobody would argue with, and that baryogenesis is complete substantially before the QCD transition epoch and that the baryon number is an integer. The significance of these requirements would hopefully become obvious from what follows. Together with the assumption of quantum entanglement, the above picture is capable of explaining both dark matter and dark energy, as we shall argue below.

A first order phase transition can be described through a bubble nucleation scenario, as depicted in figure 1. At temperatures higher than the critical temperature $`T_c`$, the coloured quarks and gluons are in a thermally equilibrated state in the perturbative vacuum (the quark-gluon plasma). The total colour of the universe is neutral (i.e., the total colour wave function of the universe is a singlet). It is strongly believed from QCD analysis that this thermally equilibrated phase of quarks and gluons exhibits collective plasma-like behaviour, referred to in the literature as the Quark-Gluon Plasma or QGP. Then, as $`T_c`$ is reached and the phase transition starts, bubbles of the hadronic phase, denoted by grey circles, begin to appear in the quark-gluon plasma, grow in size and form an infinite chain of connected bubbles (the percolation process). At this stage, the ambient universe turns over to the hadronic phase in which the constituents (protons and neutrons) are denoted by small open circles. (Mesons are omitted from the diagram, as they do not carry any conserved quantum numbers, as opposed to the baryons.) Within this hadronic phase, the remaining high temperature quark phase gets trapped in large bubbles (denoted by big coloured and shaded circles). As is well known, this process is associated with a fluctuation in the temperature around $`T_c`$; the bubbles of the hadronic phase can nucleate only when the temperature falls slightly below $`T_c`$. The released latent heat raises the temperature again and so on. It is thus fair to assume that the temperature of the universe remains around $`T_c`$ at least upto percolation.

Witten  argued some time ago that the net baryon number contained in these Trapped False Vacuum Domains (TFVD) could be many orders of magnitude larger than that in the normal hadronic phase and they could constitute the absolute ground state of strongly interacting matter. It has been shown  that if these TFVDs possess baryon number in excess of $`10^{4244}`$, they would be stable on cosmological time scales and would form the so-called strange quark nuggets (SQN). In such situations, they would have spatial radii $``$ 1 m, while their spatial separation would be $``$ 300 m. TFVDs with baryon number less than this would evaporate into normal baryons quite rapidly, much before BBN starts. (To distinguish the baryon number contained in SQNs from that participating in BBN, we denote the SQN matter as quasibaryonic.) The SQNs could evolve  to primordial structures of approximate solar mass, which should manifest themselves as the Massive Compact Halo Objects (MACHO), discovered  through gravitational microlensing in the Milky Way Halo in the direction of the Large Magellanic Cloud (LMC). If all the CDM is taken to be in the form of MACHOs, then their estimated abundance in the Milky Way Halo is seen to correspond to the observed optical depth .

In all these considerations, the explicit role of colour, the fundamental charge of strong interaction physics, has been glossed over. It has been tacitly assumed that in a many-body system of quarks and gluons, colour is averaged over, leaving only a statistical degeneracy factor for thermodynamic quantities. In this letter, we argue that such simplification may have led us to overlook a fundamentally important aspect of strong interaction physics in cosmology.

Considering the heterogeneous readership at whom this letter is aimed, it may be worthwhile to elaborate on the situation in some detail. In the QGP, all colour charges are neutralised within the corresponding Debye length, which turns out to be $`\frac{1}{g_s(T)T}`$, where $`g_s`$ is the strong coupling constant. If this Debye length is smaller than a typical hadronic radius, then hadrons cannot exist as bound states of coloured objects. The lifetime of the QGP may be roughly estimated by the temperature when the Debye screening length becomes larger than the typical hadronic radius, and formation of hadrons as bound states of coloured objects becomes possible.

Since the size of the universe is many orders of magnitude larger than the Debye length, the requirement of overall colour neutrality is trivially satisfied. Another condition required for the existence of the QGP would be the occurrence of sufficient number of colour charges within the volume characterised by the Debye length; otherwise the collective behaviour responsible for the screening would not be possible. For the cosmic QGP, these conditions are satisfied till temperatures $``$ 100 - 200 MeV, the order of magnitude for the critical temperature for QCD.

In the light of the above, one would argue that all the physics of the QGP is contained within a Debye length. A quick estimate would show that upto $`T_c`$, the Debye length is less than a fermi and the total number of colour charges (including quarks, antiquarks and gluons) within the corresponding volume is greater than 10. Note, however, our emphasis on the second ansatz above, that the baryogenesis has already taken place much before $`T_c`$ is reached and the value of $`\frac{n_b}{n_\gamma }10^{10}`$ has already been established. Firstly, it has to be realised that this baryon number is carried in quarks at this epoch, so that the net quark number ($`N_qN_{\overline{q}}`$) in the universe at all times is an exact multiple of 3. The net quark number within a Debye volume then turns out to be $`10^9`$. Thus, to ensure overall colour neutrality and an integer baryon number, one must admit of long-range correlations beyond the Debye length in the QGP, the quantum entanglement property , which was identified as an essential feature of quantum mechanics by Schrödinger in as early as 1927 but was experimentally established only recently.

Let us then consider the process of the cosmic quark-hadron phase transition from the quantum mechanical standpoint of colour confinement. As already mentioned, the colour wave function of the entire universe prior to the phase transition must be a singlet, which means it cannot be factorized into constituent product states; the wave functions of all coloured objects are completely entangled  in a quantum mechanical sense. (It may be mentioned at this juncture that the requirement of the entangled colour singlet wave function of the universe also ensures the integer baryon number for any non-zero baryon chemical potential. A quick group theoretic analysis would show that one must have net quark numbers in exact multiples of 3, irrespective of the number of quarks, antiquarks and gluons.) In such a situation, the universe is characterized by a vacuum energy corresponding to the perturbative vacuum of Quantum Chromodynamics (QCD). As the phase transition proceeds, locally colour neutral configurations (hadrons) arise, resulting in gradual decoherence of the entangled colour wave function of the entire universe. Note that the coloured objects within the hadrons are entangled among themselves but not with those in the rest of the universe. This amounts to a proportionate reduction in the perturbative vacuum energy density, which goes into providing the latent heat of the transition, or in other words, the mass and the kinetic energy of the particles in the non-perturbative (hadronic) phase. (It should be mentioned here that the vacuum energy of the non-perturbative phase of QCD is taken to be zero; more on this later.) In the quantum mechanical sense of entangled wave functions, the end of the quark-hadron transition would correspond to complete decoherence of the colour wave function of the universe; the entire vacuum energy would disappear as the perturbative vacuum would be replaced by the non-perturbative vacuum. It must be emphasized that this is an example of decoherence through a classical mechanism; the phase transition leads to a different quantum state with different entanglement property.

Combining these observations with the formation of TFVDs as discussed above, it is obvious that in order for the TFVDs to be stable physical objects, they must be colour neutral. This is synonymous with the requirement that they all have integer baryon numbers, i.e., at the moment of formation each TFVD has net quark numbers in exact multiples of 3. For a statistical process, this is, obviously, most unlikely and consequently, most of the TFVDs would have some residual colour at the percolation time, as shown in fig. 1. (This is not at all inconsistent with the screening of colour within a Debye length, even if the size of the TFVD is much larger than the Debye length. It is well known from the study of the elctromagnetic plasma that whenever a boundary occurs in the plasma, the charge neutrality condition is violated over a small length of the order of the Debye length at the boundary, the plasma sheath effect.) Then, on the way to becoming colour singlet (shaded large grey circles), they would each have to shed one or two coloured quarks, as depicted.

Thus, the end of the cosmic QCD phase transition corresponds to a situation where there would be a few coloured quarks, separated by spacelike distances. It has to be noted that such a large separation, apparently against the dictates of QCD, is by no means unphysical. The separation of coloured TFVDs occurs at the temperature $`T_c`$, when the effective string tension is zero, so that there does not exist any long range force. By the time the TFVDs have evolved into colour neutral configurations, releasing the few orphan coloured quarks in their immediate vicinity, the spatial separation between these quarks is already too large to allow strings to develop between them; see below. (It should be emphasized here that such a situation could not occur in the laboratory searches for quark-gluon plasma through energetic heavy ion collisions. There the spatial extent of the system is of the order of a few fermi and the reaction takes place on strong interaction time scales.) Therefore, the orphan quarks must remain in isolation. In terms of the quantum entanglement and decoherence of the colour wave function, this would then mean that their colour wave functions must still remain entangled and a corresponding amount of the perturbative vacuum energy would persist in the universe. In this sense, the orphan quarks are definitely not in asymptotic states and no violation of colour confinement is involved. If we naively assume that the entanglement among these orphan quarks imply collective behaviour like a plasma, then one can estimate the corresponding Debye screening length for this very dilute system of quarks, which would still be governed primarily by the temperature. At temperatures of $``$ 100 MeV, the scale would be a few fm, much smaller than the mutual separation between the orphan quarks; thus the formation of bound states of orphan quarks would be impossible.

There does not exist any way to calculate this persisting, or even the full, perturbative vacuum energy, from first principles in QCD. For the latter quantity, one may adopt the phenomenological Bag model  of confinement, where the Bag parameter B ($`(145MeV)^4`$) is the measure of the difference between the perturbative and the non-perturbative vacua. Thus we can assume that at the beginning of the phase transition, the universe starts out with a vacuum energy density B, which gradually decreases with increasing decoherence of the entangled colour wave function. A natural thermodynamic measure of the amount of entanglement during the phase transition could be the volume fraction ($`f_qV_{colour}/V_{total}`$) of the coloured degrees of freedom; at the beginning, $`f_q`$ is unity, indicating complete entanglement, while at the end, very small but finite entanglement corresponds to a tiny but non-zero $`f_q`$ due to the coloured quarks. Accordingly, the amount of perturbative vacuum energy density in the universe at any time is the energy density B times the instantaneous value of $`f_q`$; within the scenario discussed above, the remnant perturbative vacuum energy at the end of the QCD transition would just be B $`\times f_{q,O}`$, where $`f_{q,O}`$ is due solely to the orphan quarks. An order of magnitude estimate for $`f_{q,O}`$ can be carried out in the following straightforward manner. On the average, each TFVD is associated with 1 orphan quark so that the number $`N_{q,O}`$ of orphan quarks within the horizon volume at any time is about the same as the number $`N_{TFVD}`$ of TFVDs therein. It is well known from the study of percolating systems  that percolation is characterized by a critical volume fraction $`f_c`$ 0.3 of the high temperature phase. In the present case, this would require $`f_q`$ in the form of TFVD-s to be $``$ 0.3. Following the ansatz of Witten  that the most likely length scale for a TFVD is a few cm, one can estimate $`N_{TFVD}`$ (and hence $`N_{q,O}`$) within the horizon at the percolation time of about 100 $`\mu `$sec  to be about 1018-20. The inter-TFVD separation comes out to be $``$ 0.01 cm at that time. (It is obvious that the orphan quarks, separated by distances of 0.01 cm, cannot develop colour strings between them, even if there is some non-zero string tension generated at temperatures slightly lower than $`T_c`$.) Then, if we naively associate an effective radius of $`10^{14}cm`$ (estimated from $`\sigma _{qq}=\frac{1}{9}\sigma _{pp};\sigma _{pp}20mb`$ ) with each orphan quark, we obtain $`f_{q,O}N_{q,O}\times (v_{q,O}/V_{total})10^{42}10^{44}`$ (where $`v_{q,O}`$ is the effective volume of an orphan quark), so that the residual pQCD vacuum energy comes out to be in the range 10-46 to 10$`{}_{}{}^{48}GeV^4`$, just the amount of DE.

Admittedly, the preceding estimate is somewhat crude; nevertheless, the simple and natural assumptions that have led to this value ensure that it cannot be very wrong. An alternate treatment within a slightly more microscopic model also leads to the same estimate, as we show now.

The problem of finding an effective inter-quark potential at large $`r`$ from QCD is well known. In the spirit of the phenomenological string model for a quark-antiquark system (mesons), one expects a linear potential. It was shown by Richardson  some time ago that one can postulate an effective coupling constant of QCD which gives the usual leading order form at high $`Q^2`$ (small $`r`$) but can be carried over to small $`Q^2`$ (large $`r`$) limit without going through the singularity. The potential derived from such prescription yields a linear potential in the large $`r`$ limit. This ”standard” form of interaction is, unfortunately, not applicable in the context of the present analysis since here we are interested in the dynamics of interaction in a dilute many body quark system, where there may not be any antiquarks and certainly no mesons. It may however be recalled that a phenomenological prescription for incorporating the effect of confinement in such systems was proposed many years ago , through the use of a density dependent effective quark mass $`m_q1/n_q`$ in the low density sector. This prescription has been used extensively by a number of authors to study the properties of quark matter. In addition, Plümer and Raha  showed quite some years ago that such a behaviour of the effective mass was consistent with a running coupling constant which behaved like $`\alpha _s1/log(1+Q^4/\mathrm{\Lambda }^4)`$. Note that this expression agrees with the familiar perturbative form (leading order) for large $`Q^2`$, just like the Richardson case but has a very different large $`r`$ behaviour. Using this $`\alpha _s`$, one can find the effective inter-quark potential by taking the Fourier transform of $`V(q)=\alpha _s(q^2)/q^2`$ a la Richardson . This gives

$$V(r)\mathrm{\Lambda }\left((\mathrm{\Lambda }r)^3\frac{12}{\mathrm{\Lambda }r}\right)$$
(1)
For large $`r`$, only the cubic term needs to be retained. This (cubic) potential, rather than the linear one, is the most natural form of the potential to use in the case of a dilute many body system of quarks.

The inter-quark separation $`r`$ can be estimated as follows:

$`n_{q,O}N_{q,O}/V_H=(3/4\pi )N_{q,O}/R_H^3`$ (2)
$`r=\left({\displaystyle \frac{3}{4\pi }}n_{q,O}\right)^{1/3}=R_H/N_{q,O}^{1/3}`$ (3)
Then, the potential energy density for this inter-quark interaction is

$$\rho _V=\frac{1}{2}n_{q,O}V(r)\frac{3}{8\pi }\mathrm{\Lambda }^4$$
(4)
Note that this is independent of density and hence, time; once it comes into being, it will retain the same value through the subsequent evolution of the universe and can indeed be the ideal candidate for dark energy. All we need to do is to fix the value of $`\mathrm{\Lambda }`$, which would be appropriate for this dilute system of quarks. A natural measure of $`\mathrm{\Lambda }`$ would be the inverse of the length scale characterising this dilute system of quarks. For hadronic bags or hard QCD processes, the length scale is of the order of Fermi, which gives $`\mathrm{\Lambda }`$ 100 MeV. In our case, the appropriate length scale should be the size of the smallest TFVDs, as Kodama, Sasaki and Sato  argued many years ago. For the stable SQNs, with baryon density $`10^{38}cm^3`$, we had calculated the length scale for smallest TFVD  to be several cm. Now the baryon number density in the universe during the microsecond epoch is $`10^{30}cm^3`$. So, for smaller TFVDs, whose baryon density at the time of formation would be $`10^{30}cm^3`$ (corresponding to $`T_c`$ 100 MeV), the appropriate length scale must be smaller by about two orders of magnitude, $``$ 0.01 cm. (It is no surprise that the inter-TFVD separation, the other possible length scale in the problem, as estimated above comes out to be of the same order of magnitude.) Then, the corresponding $`\mathrm{\Lambda }`$ would be $`10^{12}`$ GeV, so that $`\rho _V10^{48}GeV^4`$, just the value required for dark energy.

Even though the DE component appears during the microsecond era in the history of the universe, it remains negligible in comparison to the matter density for most of the history. Since matter density decreases as $`R^3`$ (R being the scale size), while DE density remains constant, the latter can become dominant only at very late times ($`z0.17`$) and thus would not affect the galaxy formation scenarios to any extent.

The density of the orphan quarks in the present universe is exceedingly small; compared to $`10^{77}`$ baryons which took part in BBN, there would be $`10^{4445}`$ orphan quarks. Their flux in the cosmic rays would be negligibly small; non-observance of fractionally charged objects is thus not a detractor. A definitive detection of fractional electric charge would certainly clinch the issue, but as mentioned above, the flux of orphan quarks seems to make it highly urealistic. A conclusive detection of strange quark matter in the cosmic ray flux, however, would go a long way in providing support to the idea. Many such experiments are currently under way.

The picture presented here relies on a putative first order transition, although any transition resulting in finite size fluctuations may suffice. This issue requires further investigation.

A posteriori, it is tempting to mention that it should perhaps have been anticipated that if the cosmological constant does arise from a vacuum energy, then the QCD vacuum is the most natural candidate. It is a known (and accepted) lore of (renormalisable) quantum field theories that the divergent vacuum energy density is renormalised to the physical parameters of the theory. Within all the field theories now in vogue, it is only the QCD which has a distinct perturbative and a non-perturbative vacuum, which are separated by a finite energy density, irrespective of the renormalisation prescription. (That we are unable to estimate this as yet from first principles in QCD is a technical shortcoming, not a conceptual one.) It is thus most plausible that even if suitable renormalisation prescriptions remove all vacuum energy densities, the finite part of the perturbative vacuum energy density of QCD should survive and play the role of the cosmological constant.

We therefore conclude by reiterating that the emergence of both CDM and DE from the same mechanism entirely within the standard model of particle interactions, and that too in quantitatively correct amounts without any fine-tuning of parameters, is nothing short of remarkable. The picture presented here satisfies, in our opinion, the most stringent criterion of a scientific theory, that of naturalness. The strangest property of the universe, as it is said, is how simple it is.

The authors are indebted to Prof. Anthony J. Leggett for his critical comments on an earlier version of the manuscript. SR would also like to thank the Research Center for Nuclear Physics, Osaka University for their warm hospitality during his sojourn there.


text:
# Impact of densitized lapse slicings on evolutions of a wobbling black hole

## I Introduction

The accurate, long-term stable evolution of binary black hole spacetimes is one of the most important challenges in general relativity research. While the solution of the relativistic two-body problem is highly relevant from an academic point of view, its importance is further enhanced by the need of such simulations for predicting the resulting gravitational wave emission. This will not only improve the chances of detecting gravitational waves but also be essential when it comes to the astrophysical interpretation of the measured wave patterns.

Early black hole evolutions (see e.g. Anninos et al. (1995a); Cook et al. (1998); Brandt et al. (2000)) were based on the canonical ADM formulation Arnowitt et al. (1962); York (1979) of the Einstein field equations. The stability properties of these evolutions, however, were severely limited. It gradually became accepted among the community that the underlying problems were not merely numerical in nature, but originated at least in part from the structure of the equations themselves. The striking success achieved in characteristic evolutions of single black hole spacetimes Gómez et al. (1998); Gómez and et. al. (1998) supports this view. As a result, numerous alternative 3+1 formulations have been constructed in recent years. The majority of these formulations are modified versions of the original ADM scheme. These alternative formulations either recast the Einstein field equations in explicitly hyperbolic form (see Kidder et al. (2001) and references therein) or exhibit some “hyperbolic flavor” such as the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) system Shibata and Nakamura (1995); Baumgarte and Shapiro (1999).

In recent years, several of these new systems have demonstrated vastly improved stability properties, either in fully non-linear 3-dimensional black hole evolutions or suitably chosen test problems which isolate some of the particular difficulties faced in numerical relativity. Even though the question as to the optimal formulation of the Einstein equations remains as yet unanswered, the development of alternative 3+1 schemes has been one of the key driving forces underlying recent progress in numerical black hole evolutions.

In addition to the form used to recast the Einstein equations, considerable improvements have also been achieved in dealing with spacetime singularities inevitably present in such simulations. One approach, commonly known as singularity avoidance, foliates the spacetime in such a way that the evolution of the code is basically frozen near the singularity. As originally implemented, singularity avoidance methods enjoyed rather limited success, mainly due to the ‘grid-stretching’ effect — the increase of the proper separation between neighboring grid-points. The problem of grid-stretching has been alleviated to some extent by the introduction of a non-zero shift vector Alcubierre et al. (2003a). It is not clear, however, whether this approach will be able to sustain evolutions for the dynamical time-scales that are needed.

An alternative treatment, often attributed to a suggestion by Unruh (quoted in Thornburg (1987)), is based on the cosmic censorship conjecture. Cosmic censorship states that a spacetime singularity will be surrounded by a causal boundary, the black hole horizon. Physical information contained within this boundary cannot escape. Thus, in principle, one can get away with murder and excise a region inside the event horizon containing the black hole singularity, without affecting computations in the exterior. This technique, known as black hole excision or the apparent horizon boundary condition, has become very popular in recent years and has been implemented in a number of cases Seidel and Suen (1992); Anninos et al. (1995b); Marsa and Choptuik (1996); Scheel et al. (1997); Cook et al. (1998); Gómez et al. (1998); Brandt et al. (2000); Alcubierre et al. (2001); Yo et al. (2001); Alcubierre and Brügmann (2001); Shoemaker et al. (2003); Calabrese et al. (2003). It is important to note, however, that the causal boundary strictly applies to physical information only. Unless one can guarantee that all characteristic speeds of the evolved fields (physical and unphysical) lie within the light cone (for example by using an explicitly hyperbolic formulation), gauge modes may well leave the black hole interior. While these modes are not supposed to affect the physical results in the exterior evolution, they may have a significant impact on the numerical stability. The success of the excision technique as applied to non-manifestly-hyperbolic systems of equations is therefore by no means a foregone conclusion and needs to be demonstrated explicitly for the application at hand.

Earlier implementations of black hole excision were still based on the ADM scheme and frequently made use of causal-differencing and/or apparent horizon locking coordinates Cook et al. (1998); Seidel and Suen (1992); Anninos et al. (1995b). Aside from demonstrating the feasibility of the method, these simulations led to improved evolutions of single black hole spacetimes as compared with techniques based on singularity avoiding slicing conditions. These codes were, however, still subject to the fundamental limitations of the ADM formulation. In more recent times, excision has been combined with alternative formulations of the field equations and more general classes of gauge conditions. The resulting evolutions have been able to reach the first major goal on the way toward solving the binary black hole problem, namely the evolution of single, stationary black hole spacetimes in three dimensions for arbitrary times Alcubierre and Brügmann (2001); Yo et al. (2002); Gómez et al. (1998). In addressing dynamic black hole scenarios, one may attempt to absorb part or all of the black hole movement by choosing appropriate coordinates such as co-rotating coordinates in the case of an in-spiraling binary black hole Brady et al. (1998). However, it is unclear how such a treatment will be possible, in particular in the later stages just before the merger, without introducing a considerable degree of grid-stretching. It is unlikely, therefore, that such dynamic simulations can be achieved without an excision region that can adjust dynamically to mask the black hole singularity as it moves through the computational domain. In the words of Ref. Scheel et al. (1997), “a key milestone in three-dimensional black hole simulations is the ability to stably move a hole through the numerical grid”.

In characteristic formulations, long-term stable evolutions using a dynamic excision region have been achieved as early as 1998 for the case of a wobbling black hole Gómez and et. al. (1998). In 3+1 formulations, to our knowledge, the earliest attempts were performed by Cook et al. (1998), who managed to track a boosted black hole through the grid for 60 $`M`$ ( where $`M`$ is the mass of the black hole). Further applications of dynamic singularity excision include grazing collisions of black holes Brandt et al. (2000) and the evolution of a scalar field in a static or boosted Kerr-Schild black hole background Yo et al. (2001). In Shoemaker et al. (2003), a dynamic scenario was simulated by evolving a single black hole in ingoing-Eddington-Finkelstein (iEF) coordinates transformed to a coordinate system in which the coordinate location of the singularity traced out bouncing or circular paths on the coordinate grid. As the excision region is moved across the numerical grid in such evolutions, one will inevitably encounter previously excised grid points which reemerge into the evolved domain and thus need to be filled with valid data for the grid variables. In evolutions with dynamic excision, this is commonly done via extrapolation from neighboring grid points. The motivation for this procedure arises from the same causal considerations underlying the singularity excision mentioned above and is subject to the same reservations for non-manifestly-hyperbolic evolution systems. Several extrapolation techniques have been studied in Yo et al. (2001), where the stability properties of the ensuing evolutions are found to be rather insensitive to the precise details of the scheme.

In this paper, we expand directly on the work of Shoemaker et al. (2003) (referred to as Paper I from now on) where third-order extrapolation was used for the population of reemerging grid points. The resulting dynamic excision enabled these authors to evolve non-boosted black holes traveling on circular or bouncing trajectories for about $`130M`$. While such evolution times represent considerable progress on previous dynamic simulations, they do not meet the requirements for simulations of in-spiraling binary black holes. There remains, thus, the fundamental question of whether the dynamic excision techniques currently under consideration are in principle capable of handling long-term stable evolutions of moving black holes. Addressing this question is one of the main goals of this paper. We demonstrate how the lifetime of the dancing or wobbling black hole evolutions of Paper I can be extended to at least several thousands of $`M`$ by prescribing the slicing condition in terms of the densitized lapse as opposed to the lapse function itself. This also confirms the belief expressed in Paper I that the gauge conditions were a crucial factor in limiting the stability properties of the code. We further present a simple and very efficient mechanism by which our excision region was able to track the black hole movement by merely analyzing the updated numerical data, with no prior knowledge of the black hole trajectory. This mechanism analyzes the fall-off behavior of certain grid variables near the black hole singularity and provides us with an alternative at negligible computational cost to more commonly used apparent horizon tracking methods. This method for tracking singularities is suitable for situations in which the black hole is mildly distorted.

The numerical code used for the evolutions presented in this paper is the Maya code; details of the code can be found in Paper I. The code is based on the BSSN formulation Shibata and Nakamura (1995); Baumgarte and Shapiro (1999); we will not repeat the BSSN equations here — see e.g. Alcubierre et al. (2003a). Whenever we encounter conformally transformed quantities, we will distinguish them by a “hat” from physical quantities. Latin and Greek indices will run from 1 to 3 and 0 to 3 respectively. We define the order of an extrapolation scheme as the order of the polynomial used.

The paper is organized as follows. We describe and motivate in Sec. II the use of a densitized lapse for our numerical evolution. The numerical results together with convergence tests will be presented in Sec. III. In Sec. IV, we describe the simple scheme used in our evolutions for tracking the black hole singularity. We conclude in Sec. V with a discussion of our results.

## II Introducing a densitized lapse

Under the canonical 3+1 decomposition of the Einstein equations, the line element is given by

$$ds^2=\alpha ^2dt^2+\gamma _{ij}(dx^i+\beta ^idt)(dx^j+\beta ^jdt),$$
(1)
where the 3-metric $`\gamma _{ij}`$ describes the geometry of the space-like slices. The lapse function $`\alpha `$ relates the separation in proper time between two space-like slices to their difference $`dt`$ in coordinate time, and the shift vector $`\beta ^i`$ determines how one identifies points with identical spatial coordinates on neighboring spatial hypersurfaces. For recasting the Einstein equations in a 3+1 form, it is convenient to also introduce the extrinsic curvature

$$K_{ij}=\frac{1}{2\alpha }(_t_\beta )\gamma _{ij},$$
(2)
where $`_\beta `$ denotes the Lie derivative along the shift vector. The extrinsic curvature describes the embedding of the 3-dimensional spatial slices in 4-dimensional spacetime. The 3+1 Einstein equations fall into two categories, 12 evolution equations for $`\gamma _{ij}`$ and $`K_{ij}`$, and 4 constraint equations (the Hamiltonian and momentum constraints). The majority of formulations of the Einstein field equations, and in particular the BSSN formulation we use for our code, are based on this decomposition but introduce various modifications, such as conformal decompositions, promoting auxiliary variables to fundamental status and combining the evolution and constraint equations to obtain systems with different mathematical and numerical properties.

The gauge functions $`\alpha `$ and $`\beta ^i`$, on the other hand, represent the coordinate or gauge freedom of general relativity and can be chosen arbitrarily without affecting the physical spacetime. Their choice does, however, have a crucial impact on the performance of a numerical scheme used to evolve the Einstein equations. Since the evolution of the gauge variables is not determined by the Einstein equations, one needs to introduce a recipe to determine the the lapse and shift. Such gauge conditions can in principle be divided into three types (see for example Khokhlov and Novikov (2002)). Fixed gauge conditions are those for which $`\alpha `$ and $`\beta ^i`$ are given functions of the spacetime coordinates $`x^\alpha `$. Algebraic conditions, on the other hand, are such that the lapse and shift are allowed to also depend on other evolved variables, such as $`\gamma `$, $`K`$ or derivatives of $`\gamma _{ij}`$. Finally, one may construct $`\alpha `$ and $`\beta `$ from solutions of differential equations (elliptic, parabolic or hyperbolic). This type is sometimes referred to as differential gauge conditions. Examples are minimal distortion shift and maximal slicing. A subclass of algebraic gauge conditions will be of particular importance for this work.

For the present work, we will only consider a fixed gauge condition for the shift vector; that is, the shift vector will be given by the exact solution for a single, non-rotating black hole in iEF coordinates. As mentioned before, in Paper I we also used a fixed gauge condition for the lapse function. The objective now is to replace this with an algebraic condition. In particular the algebraic lapse we consider is

$$\alpha =\gamma ^{n/2}q=e^{6n\varphi }q,$$
(3)
where $`\varphi `$ is the conformal factor relating the physical metric $`\gamma _{ij}`$ to the conformal metric $`\widehat{\gamma }_{ij}`$ such that

$$\gamma _{ij}=\gamma ^{1/3}\widehat{\gamma }_{ij}=e^{4\varphi }\widehat{\gamma }_{ij}.$$
(4)
We emphasize that in Eq. (3), $`q`$ replaces the lapse $`\alpha `$ as the independent variable. Notice that because $`\alpha `$ is a true scalar, $`q`$ is a density of weight $`n`$, hence the name densitized lapse.

Such a densitized lapse or a generalized version thereof has played a crucial rule in the development of various hyperbolic formulations of Einstein’s field equations. For example Kidder et al. Kidder et al. (2001) find the use of a densitized version of the lapse necessary for obtaining strong hyperbolicity. In particular these authors obtain exclusively physical characteristic speeds for the case corresponding to $`n=1`$ in our Eq. (3). Similarly, Frittelli and Reula Frittelli and Reula (1996) require a densitized lapse with $`n>0`$ in order to obtain a symmetric hyperbolic system. In Ref. Calabrese et al. (2002) 1-dimensional non-linear evolutions are presented to highlight the markedly different convergence and stability properties of strongly hyperbolic, weakly hyperbolic and completely ill-posed problems. Note that these authors switch between the different schemes by describing the spacetime slicing in terms of a densitized lapse of weight $`n=+1,\mathrm{\hspace{0.17em}\hspace{0.17em}0},1`$ (translated into our notation). Our results are compatible with their findings in that we observe extremely poor behavior in evolutions using $`n=1`$, somewhat better — though not convincing — behavior with $`n=0`$, and (as is demonstrated below) long-term stability and convergence for dynamic scenarios in three dimensions on the same time scale as in their evolutions by using $`n=+1`$.

With respect to the BSSN system used in our calculations, it is of particular interest to note that Sarbach et al. Sarbach et al. (2002) managed to establish equivalence between the BSSN system and hyperbolic formulations under certain conditions. One of these conditions is again the use of a densitized lapse with $`n>0`$. They also find that the special case $`n=1`$ leads to exclusively physical characteristic speeds. In a more general analysis of the stability of gauge conditions carried out within the framework of short wavelength perturbations, Khokhlov et al. Khokhlov and Novikov (2002) demonstrate that the use of fixed gauge conditions yields ill-posedness in all cases except for vanishing shift and spatially constant lapse. Allowing the lapse to be a function of coordinates and the volume element $`\gamma `$ on the other hand they obtain a well-posed system provided that $`\alpha ^2/\gamma >0`$. For the restricted class of densitized lapse defined in Eq. (3) this is equivalent to the requirement that $`n>0`$. It is also interesting to note that hyperbolic systems have resulted in significantly longer evolution times for a stationary single black hole than comparable evolutions of the BSSN or ADM systems if the gauge was specified in the form of given functions of $`x^\alpha `$. Scheel at al. Scheel et al. (2002), for example, present lifetimes of 3-dimensional black hole evolutions about two orders of magnitude longer than those commonly achieved with non-hyperbolic systems where the lapse (as opposed to the densitized lapse) has been specified algebraically. To our knowledge, the only case where a densitized version of the lapse has been used in combination with a BSSN-like system is that of Laguna and Shoemaker Laguna and Shoemaker (2002), who also report substantially improved stability properties. The results presented below will indeed strongly indicate that the use of a densitized lapse is the key ingredient for obtaining the improved stability properties in these cases.

As mentioned before, in our numerical evolutions we will mainly focus on an algebraic slicing which requires us to prescribe $`q`$ as a given function of time and space and calculate the lapse $`\alpha `$ from Eq. (3). Because $`\varphi `$ in Eq. (3) is a numerically evolved quantity, in a qualitative manner of speaking, we enable the gauge conditions to respond in a more flexible way to the numerical evolution.

For our evolutions of static black holes, we will also consider the case of a differential slicing condition. In particular we choose the live 1+log slicing which has been used with great success in such scenarios Yo et al. (2002); Alcubierre and Brügmann (2001). This slicing condition is a special case of the so-called modified Bona-Massó slicing

$$_t\alpha =\alpha f(\alpha )\left(\alpha K_i\beta ^i\right)=\alpha f(\alpha )\mathrm{\hspace{0.17em}6}_t\varphi ,$$
(5)
where $`f`$ is an arbitrary positive function of $`\alpha `$. The second equality in Eq. (5) is due to the evolution equation for $`\varphi `$, namely

$$6_t\varphi =\alpha K+_i\beta ^i.$$
(6)
Alcubierre et al. Alcubierre et al. (2003b) showed how this family of slicing conditions can be considered a generalization of a fixed densitized lapse. For our purposes, it will be sufficient to reformulate Eq. (5) in terms of the densitized lapse. This leads to

$$_tq=q(fn)\mathrm{\hspace{0.17em}6}_t\varphi .$$
(7)
The 1+log slicing is obtained for the special choice of

$$f=\alpha ^1=q^1e^{6n\varphi }.$$
(8)
## III Numerical Results

For the 3-dimensional simulations presented here, we evolve the data according to the BSSN equations \[Eqs. (6)-(10) of Paper I\] and actively enforce $`\widehat{A}_i^i=0`$ after each time step. We always use a Courant factor of $`0.25`$. The other grid parameters as well as gauge and boundary conditions will be discussed separately for each case studied below.

### III.1 Evolutions of a single static black hole

We first consider the numerical evolution of spacetimes containing a single static black hole. As mentioned above, this problem can be considered essentially solved in non-linear 3+1 evolutions Alcubierre and Brügmann (2001); Yo et al. (2002). It is imperative, therefore, that any modification of the code preserves the capability of accurately evolving this scenario with long-term stability. For the evolutions presented here, the excision region has cubical shape within a radius of $`1.5M`$. The computational domain assumes octant symmetry and has size $`12M`$. Fig. 1 shows the $`\mathrm{}_2`$-norm of the time derivative of $`K`$ for three simulations. The results show that, as expected for static solutions, $`_tK0`$.

The solid and dashed lines are from runs with resolution $`2M/5`$ in which simple excision is used Alcubierre and Brügmann (2001). That is, the grid-functions on the excision boundary are updated by copying the time derivative from the exterior neighbor closest to the normal direction. For the dotted line, we follow Paper I and use instead a second-order extrapolation of the updated grid variables on the excision boundary instead of merely copying the time derivative. The resolution for this run is $`M/5`$

Besides the way updating at the excision boundary is handled, these three runs differ in the recipe for fixing the slicing. For the solid line, we use the classical 1+log condition in terms of $`\alpha `$; that is, the lapse is computed from Eq. (5). In the dashed line, we apply instead the densitized 1+log version given by Eq. (7). It is then clear in Fig. 1 that for simple excision the choice of classical versus densitized 1+log slicing condition does not affect the quality of the simulation.

The run represented by the dotted line was intended to address the question of whether such long-term stable evolutions can also be achieved if we use the purely algebraic gauge conditions given by Eq. (3), with $`q`$ obtained from the exact iEF solution. The answer is affirmative; that is, condition (3) in conjunction with second-order extrapolation on the excision boundary yields qualitatively similar results to those depicted by the solid and dashed lines. The motivation for switching for this run from simple excision to an excision with extrapolation of the grid-functions lies in anticipation of the wobbling black hole simulations. When the evolved black hole solution is time-dependent (as is the case of a wobbling black hole), simple excision is no longer suitable. The copying of time-derivatives at the excision boundary used in simple excision becomes effectively a boundary condition on the spatial derivatives of grid-functions. This boundary condition is incompatible with the outflow nature of the excision boundary.

Comparison of the dotted line with the solid and dashed lines in Fig. 1 shows that the algebraic densitized lapse run possesses long-term stability. The difference is mostly in the variation of the time when machine precision is reached. In summary, using an analytic densitized lapse $`q`$ or equivalently an algebraic lapse $`\alpha `$, we are able to evolve a single black hole with stability properties comparable to those obtained with differential conditions such as 1+log slicing and the $`\mathrm{\Gamma }`$-driver condition for the shift vector. While it is not clear how helpful such analytic slicing conditions will be for simulations of a merging binary black hole, the introduction of a densitized lapse makes them at least available for serious consideration in long-term stable simulations carried out with the BSSN scheme.

### III.2 Evolutions of single moving black holes

We now turn our attention to single black holes moving across the numerical grid. One way to obtain such a scenario is to evolve a single boosted black hole. The ensuing motion will, however, move the black hole off the computational domain. With the sizes of the computational domain currently restricted by available hardware resources, this will happen on time scales significantly shorter than those considered relevant for black hole orbits or for testing long-term stability of simulations, namely simulations lasting more than $`1000M`$.

In Paper I, we introduced therefore an alternative approach which facilitates the motion of single black holes with trajectories confined to the computational domain. We transformed the iEF black hole solution to coordinates $`\overline{x}^\alpha `$ such that

$`t`$ $`=`$ $`\overline{t},`$ (9)
$`x^i`$ $`=`$ $`\overline{x}^i+\xi ^i(\overline{t}).`$ (10)
In terms of the new coordinates, the line element becomes

$$ds^2=\alpha ^2d\overline{t}^2+\gamma _{ij}\left(d\overline{x}^i+\overline{\beta }^id\overline{t}\right)\left(d\overline{x}^j+\overline{\beta }^jd\overline{t}\right),$$
(11)
where

$$\overline{\beta }^i=\beta ^i+\dot{\xi }^i.$$
(12)
In Paper I, this method was used to move black holes on circular and bouncing trajectories. Using either a cubical or spherical excision region, these runs lasted for about $`130M`$, though the apparent horizon started to intersect the excision region at about $`90M`$ and the apparent horizon finder failed to give reasonable results afterwards. The fixed gauge conditions used in these evolutions were considered a crucial limiting factor in those evolutions.

Using the densitized lapse condition (3) with $`n=1`$, we performed similar simulations to those in Paper I. We fix the value of $`q`$ in Eq. (3) from the exact single iEF black hole solution transformed according to Eqs. (10), (12). In Fig. 2, we plot the $`\mathrm{}_2`$-norm of the Hamiltonian constraint (upper panel) and normalized Hamiltonian constraint (lower panel) for simulations of a black hole on a circular path with radius $`2M`$ and orbital frequency $`1/(2\pi )`$. Equatorial symmetry is assumed for all the simulations, and the outer boundary conditions consist of setting the values of the grid-functions to the exact analytic solution.

Fig. 2 depicts results from four simulations. In the upper panel from top to bottom, each data set corresponds to computational domains $`20\times 20\times 7M^3`$, $`20\times 20\times 10M^3`$, $`30\times 30\times 7M^3`$ and $`30\times 30\times 10.5M^3`$, respectively. In the lower panel, the correspondence is reversed. That is, from top to bottom at early times, each data set is for domains $`30\times 30\times 10.5M^3`$, $`30\times 30\times 7M^3`$, $`20\times 20\times 10M^3`$ and $`20\times 20\times 7M^3`$, respectively. The shortest scale is along the direction perpendicular to the orbital plane.

The results in Fig. 2 from the $`20\times 20\times 7M^3`$ run should be compared with Fig. 9 of Paper I where simulations lasted $`<130M`$. This comparison clearly demonstrates the tremendous improvement on the duration of the simulations when a densitized lapse is used. Some of the evolutions are stable at least for $`6000M`$, when the runs were terminated due to limitations of computational resources. Fig. 2 also demonstrates that, even though the durations of the simulations have been improved by at least an order of magnitude, the simulations continue to be affected by boundary effects. This is expected since it is well known that setting boundary conditions to the exact analytic solution is conducive to numerical instabilities.

In Fig. 3 we show snapshots of the evolution obtained on the domain $`20\times 20\times 7M^3`$ at $`t=0`$, $`2000M`$, $`4000M`$, $`6000M`$. Here we plot the variable $`K`$ on the $`xy`$-plane at $`z=0`$, namely the orbital plane. The excision region is clearly visible and has been checked lie within the apparent horizon. We find the apparent horizon area to remain within a few per cent of its analytic value $`16\pi M^2`$ and the deviation from spherical shape to be of the same order. The method we have used to enable the excision region to follow the motion of the black hole is based on an analysis of the asymmetry and fall-off behavior of grid variables near the excision region (see next section for details). This method is in general not coordinate invariant and its performance needs to be monitored by verifying that the excision region always remains inside the apparent horizon. We have verified this for all runs discussed in this paper Den .

By construction, our numerical code is second-order accurate. Fig. 4 shows the ability of our code to maintain second-order convergence. Because demonstrating long-term convergence requires considerable computational resources, we have only performed the convergence analysis for times $`<1000M`$. The jaggedness in these plots is due to the requirement that the center of the excision region always be located on a grid point. The excision mask will therefore move across the computational domain in a discontinuous fashion.

In order to further test the robustness of our excision infrastructure, we have performed the same type of evolution using a different black hole trajectory. Having in mind the eventual target of simulating in-spiraling binary black holes, it will be of particular interest to see whether the code is able to evolve a black hole on an in-spiraling trajectory. We consider such a scenario by fixing the location of the black hole singularity, $`\stackrel{}{\xi }=(\xi ^1,\xi ^2)`$ \[cf. Eq. (10)\], from the solution of the differential equation

$$_{tt}\stackrel{}{\xi }+D_t\stackrel{}{\xi }+\mathrm{\Phi }=0,$$
(13)
with $`\xi ^3=0`$ because the black hole orbit is confined to the plane $`z=0`$. The potential is given by $`\mathrm{\Phi }=C[(\xi ^1)^2+(\xi ^2)^2]`$ and $`D`$ and $`C`$ are constants.

We emphasize that none of these quantities has any connection to realistic astrophysical simulations. They merely serve to determine the components $`\xi ^i`$ in the coordinate transformation (10) and, thus, the artificial black hole trajectory. We obtain an in-spiraling trajectory by choosing $`D=C=0.01`$. The black hole is initially placed at a distance of $`2.4M`$ from the origin at $`z=0`$ and given an initial purely tangential velocity $`v=0.3`$. Because of the damping, the velocity and radius will decrease and the black hole would eventually approach a steady state at the origin. In order to keep the evolution dynamic for long times, we switch the sign of the damping term when the radius shrinks below $`0.25M`$ after which the hole starts spiraling outward. The damping constant is switched back to a positive value once a radius of $`2.5M`$ is reached and so on. The resulting trajectory in the $`xy`$-plane is displayed in Fig. 5. This path corresponds to an evolution lasting $`3500M`$ when we decided to terminate the evolution.

The apparent horizon finder had no difficulty calculating the outermost trapped surface every $`20M`$. The resulting horizon area is shown as a function of time in Fig. 6. The area remains constant with good accuracy and deviates by a few per cent from the analytic value of $`16\pi M^2`$. Similarly we found the error in the Hamiltonian constraint normalized with respect to its individual addends to be constant to within about $`6\%`$.

## IV The excision tracking method

For the dynamic evolutions of the previous section, we have used a simple mechanism by which the excision region tracks the movement of the black hole. We will describe this method in more detail now.

In order to make use of the causal disconnection of the black hole interior it is imperative that the excision region be located within the event horizon of the black hole. Because the location of the event horizon requires analysis of the global (in space and time) data and can thus not be carried out during the evolution itself, one commonly resorts to searching for apparent horizons. An apparent horizon (also called a marginally trapped outer surface) is defined locally and can thus be calculated at each time step (see for example Thornburg (1996)). The existence of an apparent horizon guarantees that there is a surrounding event horizon and thus a causal boundary. Ideally one would perform a search for an apparent horizon at each time step and adjust the location of the excision region correspondingly. Unfortunately the calculation of apparent horizons is computationally expensive although recent work by Thornburg Thornburg (2003) seems to overcome some of these difficulties. It appears, nonetheless, desirable to have available a reliable, computationally inexpensive method of dynamically adjusting the location of the excision regions to accommodate the movement of the black holes. One may then calculate the apparent horizon at regular intervals to monitor the causal consistency of the excision procedure.

For this purpose we have used a simple scheme which analyzes the fall-off and symmetry behavior of evolution variables near the black hole singularity to predict a reasonable value for the central position of the excision region. The excision technique used in our Maya code has been described in detail in Paper I, and we will focus here on those parts added for the tracking of the excision center. In Fig. 7 we have schematically plotted a 2-dimensional slice of a spherical excision region (the precise shape of that region is irrelevant for our method). The current excision center is marked by the symbol $`\times `$. Assuming that the black hole is slightly offset with respect to this position, we have marked by a $`+`$ the desirable location of the excision center. In a loose sense we might call this the current black hole position. We note, however, that a precise definition of a black hole center, as for example in terms of a point-like singularity, can only be given for a restricted subclass of black holes.

For concreteness, we now consider analytic data for a static black hole of mass $`M=1`$ in iEF coordinates located at $`x=0.105`$, $`y=0`$, $`z=0`$. We describe this data in terms of an excision region with radius $`0.7`$ centered at $`x=0`$, $`y=0`$, $`z=0`$. Quite obviously this is not the optimal position of the excision region and the data will show some asymmetry in their fall-off behavior in the vicinity of the excision boundary. The idea is then to construct some combination of the field variables which adequately exhibits this behavior. Such a combination will in general not be coordinate-invariant and thus depends on the scenario under investigation. For all runs presented in this paper, we have found the trace of the extrinsic curvature $`K`$ a perfectly adequate choice. A 1-dimensional plot in Fig. 8 of $`K`$ along the $`x`$-axis through the excision center reveals the asymmetry. The $`\times `$ symbols in this figure correspond to the values of $`K`$ at the 8 points marked by boxes on the $`x`$-axis in Fig. 7. At the end of each evolution step in our code we fit a Gaussian

$$ae^{\frac{(xb)^2}{c^2}}$$
(14)
through these points. The parameters $`a`$, $`b`$ and $`c`$ are obtained from $`\chi ^2`$-minimization and the value for $`b`$ gives us the $`x`$-coordinate of the updated excision center. We then proceed similarly for the $`y`$ and $`z`$ direction. If further excision regions are

present in the computational domain they are treated in the same way. The total number of points used for this method is a free parameter but we typically find 8 (as in this example) to be sufficient.

With regard to non-stationary scenarios it is, of course, possible that a certain asymmetry of the data might arise from reasons unrelated to the black hole location, as for example in the case of boosted black holes. We emphasize, however, that the purpose of this algorithm is not to provide as accurate as possible an estimate of a black hole center (which in many cases will not be well defined anyway), but to prescribe a recipe for centering the excision region. The only requirement for a healthy evolution is that the excision region be contained entirely within the apparent horizon. This may be monitored, for example, by regular calculations of the apparent horizon. In the numerical evolutions presented above we have used a few buffer zones (layers of grid points inside the apparent horizon which are not excised) and verified that the excision region indeed remains confined to the interior of the apparent horizon.

## V Conclusions

We have studied the impact of a densitization of the lapse function on the stability properties of BSSN evolutions using singularity excision. We have demonstrated that the use of a densitized lapse preserves the capability of stably evolving a single Schwarzschild black hole using the simple excision technique in combination with a differential slicing condition of the live 1+log type. These numerical simulations settle down to a configuration that remains constant in time up to machine precision.

We have next considered the evolutions of single static black holes in the case where the slicing condition is a prescribed analytic function of the spacetime coordinates. In such scenarios numerical codes using the BSSN scheme have hitherto exhibited significantly poorer stability properties than those based on hyperbolic formulations. As the use of a densitized lapse (or a generalized version thereof) is a necessary ingredient in the majority of hyperbolic formulations, it has naturally been used in these codes in place of the lapse function itself. We have bridged the gap between hyperbolic and BSSN-like formulations, by evolving a stationary single black hole spacetime using the BSSN formulation in combination with a densitized version of the lapse function. We have thus been able to obtain evolutions lasting many thousands of $`M`$ in which the time-variations decay to machine precision. To our knowledge this has so far not even been achieved with strictly hyperbolic formulations. Regarding the use of purely analytic gauge conditions for such scenarios, our results indicate that the use of a densitized lapse (i.e. algebraic as opposed to fixed slicing) has been the key advantage of hyperbolic systems over BSSN-like formulations of Einstein’s equations. It is not clear to what extent algebraic gauge conditions will be beneficial in the simulation of in-spiraling binary black holes, but the results presented here demonstrate how such conditions can also be used for long-term stable evolutions carried out in the framework of the BSSN formulation.

These ideas have also shown advantages in evolving black holes moving across the computational domain. By virtue of time-dependent coordinates, we are able to move single iEF black holes through the computational domain on essentially arbitrary paths. The duration of such evolutions has previously been limited to about $`100M`$. Because the motion of the singularity is prescribed in terms of the gauge conditions, these simulations inevitably used analytic gauge. By reformulating this in terms of the densitized lapse, we were able to extend the lifetime of the runs to at least $`1000M`$. This also demonstrates that the dynamic excision techniques used for our simulations are, in principle, capable of handling long-term stable evolutions of moving black holes. We believe that the ideas presented in this work will be helpful for the development of long-term stable evolutions of black hole collisions.

###### Acknowledgements.

Work partially supported by NSF grant PHY-9800973 to Penn State University and PHY-9800737, PHY-9900672 to Cornell University. We acknowledge the support of the Center for Gravitational Wave Physics funded by the NSF under PHY-0114375. P.L. and D.S. thank the support while visiting the KITP during which part of this work was completed. KITP is supported in part by the NSF under PHY-9907949. Parallel and I/O infrastructure provided by Cactus.