text:
# Parallelization Strategies for ALI Radiative Transfer in Moving Media

## 1. Introduction

In general parallelization is a subject to be avoided (as nicely described by P. Höflich, this volume), however in order to take advantage of the enormous computing power and vast memory sizes of modern parallel supercomputers, allowing both much faster model calculations as well as more detailed models, we have implemented a parallel version of the general model atmosphere code PHOENIX (Hauschildt & Baron 1999 and references therein). Since the code uses a modular design, we have implemented different parallelization strategies for different modules in order to maximize the total parallel speed-up of the code. Our implementation allows us to change the load distribution onto different processor elements (PEs) both via input files and dynamically during a model run, which gives a high degree of flexibility to optimize the performance on a number of different machines and for a number of different model parameters.

Since we have both large CPU and memory requirements we have implemented the parallel version of the code on using the MPI message passing library (MPI Forum 1995).

We have chosen to work with the MPI message passing interface, since it is both portable (public domain implementations of MPI are readily available cf. Gropp et al. 1996), running on dedicated parallel machines and heterogeneous workstation clusters and it is available for both distributed and shared memory architectures. For our application, the distributed memory model is in fact easier to code than a shared memory model, since then we do not have to worry about locks and synchronization, etc. on small scales and we, in addition, retain full control over interprocess communication. This is especially clear once one realizes that it is fine to execute the same code on many PEs as long as it is not too CPU intensive, and avoids costly communication. We have added a few simple OpenMP directives to our code, but do not discuss SMP parallelization further here.

An alternative to an implementation with MPI is an implementation using High Performance Fortran (HPF) directives (in fact, both can co-exist to improve performance). However, the process of automatic parallelization guided by the HPF directives is presently not yet generating optimal results because the compiler technology is still very new. HPF is also more suited for problems that are purely data-parallel (SIMD problems) and would not benefit much from a MIMD approach. An optimal HPF implementation of PHOENIX would also require a significant number of code changes in order to explicitly instruct the compiler not to generate too many communication requests, which would slow down the code significantly. The MPI implementation requires only the addition of a few explicit communication requests, which can be done with a small number of library calls.

## 2. Basic numerical methods

The co-moving frame radiative transfer equation for spherically symmetric flows can be found in e.g., Mihalas & Mihalas (1984). The equation of radiative transfer is a integro-differential equation, since the differential operator on the left hand side acts on the specific intensity $`I_\nu `$, but the right hand side contains the emissivity $`\eta _\nu `$ which in turn contains $`J_\nu `$, the zeroth angular moment of $`I_\nu `$:

$$\eta _\nu =\kappa _\nu S_\nu +\sigma _\nu J_\nu ,$$

and

$$J_\nu =1/2_1^1𝑑\mu I_\nu ,$$

where $`S_\nu `$ is the source function, $`\kappa _\nu `$ is the absorption opacity, and $`\sigma _\nu `$ is the scattering opacity. With the assumption of time-independence $`\frac{I_\nu }{t}=0`$ and a monotonic velocity field the transfer equation becomes a boundary-value problem in the spatial coordinate and an initial value problem in the frequency or wavelength coordinate. The equation can be written in operator form as:

$`J_\nu =\mathrm{\Lambda }_\nu S_\nu ,`$ (1)
where $`\mathrm{\Lambda }`$ is the lambda-operator.

## 3. Parallel radiative transfer

### 3.1. Strategy and implementation

We use a modified version of the method discussed in Hauschildt (1992) for the numerical solution of the special relativistic radiative transfer equation (RTE) at every wavelength point (see also Hauschildt, Störzer, & Baron 1994 and Hauschildt & Baron 1999). This iterative scheme is based on the operator splitting approach. The RTE is written in its characteristic form and the formal solution along the characteristics is done using a piecewise parabolic integration (PPM, this is the “short characteristic method” Olson & Kunasz 1987). We use the exact band-matrix subset of the discretized $`\mathrm{\Lambda }`$-operator as the ‘approximate $`\mathrm{\Lambda }`$-operator’ (ALO) in the operator splitting iteration scheme (see Hauschildt et al. 1994). This gives very good convergence and high speed-up when compared to diagonal ALO’s.

The serial radiative transfer code has been optimized for superscalar and vector computers and is numerically very efficient. It is therefore crucial to optimize the ratio of communication to computation in the parallel implementation of the radiative transfer method. In terms of CPU time, the most costly parts of the radiative transfer calculation are the setup of the PPM interpolation coefficients and the formal solutions (which have to be performed in every iteration). The construction of a tri-diagonal ALO requires about the same CPU time as a single formal solution of the RTE and is thus not a very important contributor to the total CPU time required to solve the RTE at every given wavelength point.

In principle, the computation of the PPM coefficients does not require any communication and thus could be distributed arbitrarily between the PEs. However, the formal solution is recursive along each characteristic. Within the formal solution, communication is only required during the computation of the mean intensities, $`J`$, as they involve integrals over the angle $`\mu `$ at every radial point. Thus, a straightforward and efficient way to parallelize the radiative transfer calculation is to distribute sets of characteristics onto different PEs. As Auer (this volume) has emphasized, the solution of the radiative transfer equation is independent along a ray and in the sense that we treat each characteristic separately for the solution to the RTE, this parallelization step is one of domain decomposition. Figure 1 illustrates the characteristics for a typical supernova model calculation.

Within one iteration step, the current values of the mean intensities need to be broadcast to all radiative transfer PEs and the new contributions of every radiative transfer PE to the mean intensities at every radius must be sent to the master PE. The master radiative transfer PE then computes and broadcasts an updated $`J`$ vector using the operator splitting scheme, the next iteration begins, and the process continues until the solution is converged to the required accuracy. The setup, i.e., the computation of the PPM interpolation coefficients and the construction of the ALO, can be parallelized using the same method and PE distribution. The communication overhead for the setup is roughly equal to the communication required for a single iteration.

An important point to consider is the load balancing between the NRT radiative transfer PEs. The workload to compute the formal solution along each characteristic is proportional to the number of intersection points of the characteristic with the concentric spherical shells of the radial grid (the ‘number of points’ along each characteristic). Therefore, if the total number of points is $`N_\mathrm{P}`$, the optimum solution would be to let each radiative transfer PE work on $`N_\mathrm{P}/N_{\mathrm{RT}}`$ points. In general, this optimum cannot be reached exactly because it would require splitting characteristics between PEs (which involves both communication and synchronization). We therefore compromise by distributing the characteristics to the radiative transfer PEs so that the total number of points that are calculated by each PE is roughly the same and that every PE works on a different set of characteristics.

## 4. Parallelization of the line opacity calculation

The contribution of spectral lines by both atoms and molecules is calculated in PHOENIX by direct summation of all contributing lines. Each line profile is computed individually at every radial point for each line within a search window (typically, a few hundred to few thousand Doppler widths or about $`1000`$Å). This method is more accurate than methods that rely on pre-computed line opacity tables (the “opacity sampling method”) or methods based on distribution functions (the “opacity distribution function (ODF) method”).

In typical model calculations we find that about 1,000 to 10,000 spectral lines contribute to the line opacity (e.g., in M dwarf model calculations) at any wavelength point, requiring a large number of individual Voigt profiles at every wavelength point. The subroutines for these computations can easily be vectorized and we use a block algorithm with caches and direct access scratch files for the line data to minimize storage requirements. A block is the number of lines stored in active memory, and the cache is the total number of blocks stored in memory. When the memory size is exceeded, the blocks are written to direct access files on disk. Thus the number of lines that can be included is not limited by RAM, but rather by disk space and the cost of I/O. This approach is computationally efficient because it provides high data and code locality and thus minimizes cache/TLB misses and page faults.

There are several obvious ways to parallelize the line opacity calculations. The first method is to let each PE compute the opacity at $`N_\mathrm{r}/N_{\mathrm{PE}}`$ radial points, the second is to let each PE work on a different subset of spectral lines within each search window. A third way, related to the second method, is to use completely different sets of lines for each PE (i.e., use a global workload split between the PEs in contrast to a local split in the second method). We discuss the advantages and disadvantages of these three methods. All three methods require only a very small amount of communication, namely a gather of all results to the master PE with an MPI\_REDUCE library call.

The first and second methods, distributing sets of radial points or sets of lines within the local (wavelength dependent) search window over the PEs, respectively, are very simple to implement. They can easily be combined to optimize their performance: if for any wavelength point the number of depth points is larger than the number of lines within the local search window, then it is more effective to run this wavelength point parallel with respect to the radial points, otherwise it is better to parallelize with respect to the lines within the local window. It is trivial to add logic to decide the optimum method for every wavelength point individually. This optimizes overall performance with negligible overhead. The speed-up for this method of parallelization is very close to the optimum value if the total number of blocks of the blocking algorithm is relatively small ($`3\mathrm{}5`$).

However, if the total number of line blocks is larger (typically, about 10 to 20 blocks are used), the overhead due to the read operations for the block scratch files becomes noticeable and can reach 20% or more of the total wall-clock time. This is due to the fact that the ‘local parallelization’ required that each PE working on the line opacities needs to read every line block, thus increasing the I/O time and load to the I/O subsystem by the number of PEs themselves. We discuss I/O parallelization in § 8.

## 5. Parallelizing NLTE calculations

Our method for iteratively solving the NLTE radiative transfer and rate equations with an operator splitting method are discussed in detail in Hauschildt (1993), and in Hauschildt & Baron (1995, 1999), therefore, we summarize our approach here only a briefly. The method uses a “rate-operator” formalism that extends the approach of Rybicki & Hummer (1991) to the general case of multi-level NLTE calculations with overlapping lines and continua, i.e, “pre-conditioning”. We use an “approximate rate-operator” that is constructed using the exact elements of the discretized $`\mathrm{\Lambda }`$-matrix (these are constructed in the radiative transfer calculation for every wavelength point). This approximate rate-operator can be either diagonal or tri-diagonal. The method gives good convergence for a wide range of applications and is very stable. It has the additional advantage that it can handle very large model atoms, e.g., we use up to 11000 levels for NLTE model atoms in regular model calculations (Hauschildt et al. 1996, Baron et al. 1996, Short et al. 1999).

Parallelizing the NLTE calculations involves parallelizing three different sections: the NLTE opacities; the rates; and the solution of the rate equations. In the following discussion, we consider only a diagonal approximate rate-operator for simplicity. In this case, the computation of the rates and the solution of the rate equations as well as the NLTE opacity calculations can simply be distributed onto a set of PEs without any communication (besides the gathering of the data at the end of the iteration). This provides a very simple way of achieving parallelism and minimizes the total communication overhead. The generalization to a tridiagonal approximate rate-operator is in principle straightforward, and involves only communication at the boundaries between two adjacent PEs.

It would be possible to use the other two methods that we discussed in the section on the LTE line opacities, namely local and global set of lines distributed to different PEs. However, both would involve an enormous amount of communication between the NLTE PEs because each NLTE transition can be coupled to any other NLTE transition. These couplings require that a PE working on any NLTE transition have the required data to incorporate the coupling correctly. Although this only applies to the PEs that work on the NLTE rates, it would require both communication of each NLTE opacity task with each other (to prepare the necessary data) and communication from the NLTE opacity PEs to the NLTE rate PEs. This could mean that several MB data would have to be transferred between PEs at each wavelength point, which is prohibitive with current communication devices.

We have implemented the parallelization of the NLTE calculations by distributing the set of radial points on different PEs. In order to minimize communication, we also ‘pair’ NLTE PEs so that each PE works on NLTE opacities, rates, and rate equations for a given set of radial points. This means that the communication at every wavelength point involves only gathering the NLTE opacity data to the radiative transfer PEs and the broadcast of the result of the radiative transfer calculations (i.e., the ALO and the $`J`$’s) to the PEs computing the NLTE rates. The overhead for these operations is negligible but it involves synchronization (the rates can only be computed after the results of the radiative transfer are known, similarly, the radiative transfer PEs have to wait until the NLTE opacities have been computed). Therefore, a good balance between the radiative transfer tasks and the NLTE tasks is important to minimize waiting. The rate equation calculations parallelize trivially over the layers and involve no communication if the diagonal approximate rate-operator is used. After the new solution has been computed, the data must be gathered and broadcast to all PEs, the time for this operation is negligible because it is required only once per model iteration.

An important problem arises from the fact that the time spent in the NLTE routines is not dominated by the floating point operations (both the number and placement of floating point operation were optimized in the serial version of the code) but by memory access, in particular for the NLTE rate construction. Although the parallel version accesses a much smaller number of storage cells (which naturally reduces the wall-clock time), effects like cache and TLB misses and page faults contribute significantly to the total wall-clock time. All of these can be reduced by using Fortran-90 specific constructs in the following way: We have replaced the static allocation of the arrays that store the profiles, rates etc. used in the NLTE calculations (done with COMMON blocks) with a Fortran-90 module and explicit allocation (using ALLOCATE and SAVE) of these arrays at the start of the model run. This allows us to tailor the size of the arrays to fit exactly the number of radial points handled by each individual PE. This reduces both the storage requirements of the code on every PE and, more importantly, it minimizes cache/TLB misses and page faults. In addition, it allows much larger calculations to be performed because the RAM, virtual memory, and scratch disk space of every PE can be fully utilized, thus effectively increasing the size of the possible calculation by the number of PEs.

We find that the use of adapted array sizes significantly improves the overall performance, even for a small number of processes. The scaling with more PEs has also improved considerably. However, the overhead of loops cannot be reduced and thus the speed-up cannot increase significantly if more PEs are used. We have verified with a simple test program which only included one of the loops important for the radiative rate computation that this is indeed the case. Therefore, we conclude that further improvements cannot be obtained at the Fortran-90 source level but would require either re-coding of routines in assembly language (which is not practical) or improvements of the compiler/linker/library system. We note that these performance changes are very system dependent.

So far, (see also Hauschildt, Baron, & Allard 1997) we have described our method for parallelizing three separate modules: (1) The radiative transfer calculation itself, where we divide up the characteristic rays among processing elements and use an MPI\_REDUCE to send the $`J_\nu `$ to all the radiative transfer and NLTE rate computation tasks;111We define a task as a logical unit of code that treats an aspect of the physics of the simulation, such as radiative transfer. On the other hand, we use the term PE to indicate a single processing element of the parallel computer. Thus, a single PE can execute a number of tasks, either serial on a single CPU or in parallel, e.g., on an SMP PE of a distributed shared-memory supercomputer. (2) the line opacity which requires the calculation of about 10,000 Voigt profiles per wavelength point at each radial grid point, here we split the work amongst the processors both by radial grid point and by dividing up the individual lines to be calculated among the processors; and (3) the NLTE calculations. The NLTE calculations involve three separate parts: the calculation of the NLTE opacities, the calculation of the rates at each wavelength point, and the solution of the NLTE rate equations. So far our description has involved parallelization by distribution of the radial grid points among the different PEs or by distributing sets of spectral lines onto different PEs. In addition, to prevent communication overhead, each task computing the NLTE rates is paired on the same PE with and the corresponding task computing NLTE opacities and emissivities to reduce communication. The solution of the rate equations parallelizes trivially with the use of a diagonal rate operator.

In the latest version of our code, PHOENIX 13.0, we have incorporated the additional strategy of distributing each NLTE species (the total number of ionization stages of a particular element treated in NLTE) on separate PEs. Since different species have different numbers of levels treated in NLTE (e.g. Fe II \[singly ionized iron\] has 617 NLTE levels, whereas H I has 30 levels), care is needed to balance the number of levels and NLTE transitions treated among the PEs to avoid unnecessary synchronization problems.

## 6. Wavelength Parallelization

The division of labor outlined in the previous sections requires synchronization between the radiative transfer tasks and the NLTE tasks, since the radiation field and ALO operator must be passed between them. In addition, our standard model calculations use 50-100 radial grid points and as the number of PEs increases, so too does the communication and loop overhead.

Since the number of wavelength points in a calculation is very large and the CPU time scales linearly with the number of wavelength points, a further distribution of labor by wavelength points would potentially lead to large speedups and to the ability to use very large numbers of processors available on massively parallel supercomputers. Thus, we have developed the concept of wavelength “clusters” to distribute a set of wavelength points (for the solution of the frequency dependent radiative transfer) onto a different set of PEs. In order to achieve optimal load balance and, more importantly, in order to minimize memory requirements, each cluster works on a single wavelength point at a time, but it may consist of a number of “worker” PEs where the worker PEs use parallelization methods discussed above. In order to avoid communication overhead, the workers on each wavelength cluster are symmetric: each corresponding worker on each wavelength cluster performs identical tasks but on a different set of wavelengths for each cluster. This allows us to make use of communicator contexts, a concept which is built into MPI. The basic design is illustrated in Fig. 2. Each PE on a given wavelength cluster is assigned to a MPI\_GROUP and then the spatial parallelization occurs within that MPI\_GROUP so that, while the communicators in different wavelength clusters have identical names, they have different contexts and the code is much clearer and cleaner. For example, Worker 0 will be working on the same set of characteristic rays in the radiative transfer equation on all wavelength clusters. Therefore, Worker 0 on Wavelength cluster 0 needs to send only the intensities that it computed to its corresponding Worker 0 on Wavelength cluster 1 and so forth. The code has been designed so that the number of wavelength clusters, the number of workers per wavelength cluster, and the task distribution within a wavelength cluster is arbitrary and can be specified dynamically at run time.

In order to parallelize the spectrum calculations for a model atmosphere with a global velocity field, such as the expanding atmospheres of novae, supernovae or stellar winds, we need to take the mathematical character of the RTE into account. For monotonic velocity fields, the RTE is an initial value problem in wavelength (with the initial condition at the smallest wavelength for expanding atmospheres and at the largest wavelength for contracting atmospheres). This initial value problem must be discretized fully implicitly to ensure stability. In the simplest case of a first order discretization, the solution of the RTE for wavelength point $`i`$ depends only on the results of the point $`i1`$. In order to parallelize the spectrum calculations, the wavelength cluster computing the solution for wavelength point $`i`$ must get the specific intensities from the cluster computing the solution for point $`i1`$. This suggests a “pipeline” solution to the wavelength parallelization. Note that only the solution of the RTE is affected by this, the calculation of the opacities and rates remains independent between different wavelength clusters. In this case, the parallelization works as follows: Each cluster can independently compute the opacities, it then waits until it receives the specific intensities for the previous wavelength point, computes the solution of the RTE and immediately sends the results to the next wavelength cluster (to minimize waiting time, this is done with non-blocking send/receives), then proceeds to calculate the rates and the new opacities for its next wavelength point and so on.

In the simplest case that there are simply 2 PEs and the load balance has been set so that each PE acts as a wavelength cluster, we have one worker per wavelength cluster and two wavelength clusters. In this case, PE 1 and 2 begin by doing the various preprocessing required. PE 1 then begins working on wavelength point 1, PE 2 is working on wavelength point 2. Since PE 1 begins with the boundary condition, it calculates the opacities at this wavelength point and other necessary preprocessing. It then solves the RTE, using the boundary condition and then immediately sends the specific intensities off to PE 2. PE 1 then calculates the rates and various other post-processing at this wavelength point. At the same time PE 2 has finished calculating the opacities at its wavelength point, done the preprocessing for solving the RTE and then must wait until PE 1 has sent it the specific intensities for the previous wavelength point at which point it can then solve the RTE and immediately send its specific intensities back to PE 1, which will need them for the next wavelength point, PE 2 then calculates the rates at its wavelength point and the process repeats. Since PE 1 was busy doing post-processing, it may in fact have the specific intensities from PE 2, just in time and so minimal waiting may be required and the process proceeds in a round robin fashion. The basic pre- and post-processing required are illustrated in the pseudo-code of Figure 3.

This scheme has some predictable properties similar to the performance results for classical serial vector machines. First, for a very small number of clusters (e.g., two), the speedup will be small because the clusters would spend a significant amount of time waiting for the results from the previous cluster (the “pipeline” has no time to fill). Second, the speedup will level off for a very large number of clusters when the clusters have to wait because some of the clusters working on previous wavelength points have not yet finished their RTE solution, thus limiting the minimum theoretical time for the spectrum calculation to roughly the time required to solve the RTE for all the wavelength points together (the “pipeline” is completely filled). This means that there is a “sweet spot” for which the speedup to number-of-wavelength-clusters ratio is optimal. This ratio can be further optimized by using the optimal number of worker PEs per cluster, thus obtaining an optimal number of total PEs. The optimum will depend on the model atmosphere parameters, the speed of each PE itself and the communication speed, as well as the quality of the compilers and libraries.

The wavelength parallelization has the drawback that it does not reduce the memory requirement per PE compared to runs with a single wavelength cluster. Increasing the number of worker PEs per cluster will decrease the memory requirements per PE drastically, however, so that large runs can use both parallelization methods at the same time to execute large simulations on PEs with limited memory. On a shared-memory machine with distributed physical memory (such as the Origin), this scheme can also be used to minimize memory access latency.

## 7. Direct Opacity Sampling

In our model atmosphere code we have implemented and used direct opacity sampling (dOS) for more than a decade with very good results. During that time, the size of the combined atomic and molecular line databases that we used has increased from a few MB to $`>11`$GB. Whereas the floating point and memory performance of computers has increased dramatically in this time, I/O performance has not kept up with this speed increase. Presently, the wall-clock times used by the line selection and opacity modules are dominated by I/O time, not by floating point or overall CPU performance. Therefore, I/O performance is today more important that it was 10 years ago and has to be considered a major issue. The availability of large scale parallel supercomputers that have effectively replaced vector machines in the last 5 years, has opened up a number of opportunities for improvements of dOS algorithms. Parallel dOS algorithms with an emphasis on the handling of large molecular line databases are thus an important problem in computational stellar atmospheres. These algorithms have to be portable and should perform well for different types of parallel machines, from cheap PC clusters using Ethernet links to high performance parallel supercomputers. This goal is extremely hard to attain on all these different systems, and we compare two different parallel dOS algorithms and describe their performance on two very different parallel machines.

There are a number of methods in use to calculate line opacities. The classical methods are statistical and construct tables that are subsequently used in the calculation. The Opacity Distribution Function (ODF) and its derivative the k-coefficient method have been used successfully in a number of atmosphere and opacity table codes (e.g., Kurucz 1992). This method works well for opacity table and model construction but cannot be used to calculate detailed synthetic spectra. A second approach is the opacity sampling (OS) method (e.g., Peytremann 1974). This is a statistical approach in which the line opacity is sampled on a fine grid of wavelength points using detailed line profiles for each individual spectral line. In classical OS implementation, tables of sampling opacities are constructed for given wavelengths grids and for different elements.

In direct opacity sampling (dOS) these problems are avoided by replacing the tables with a direct calculation of the total line opacity at each wavelength point for all layers in a model atmosphere (Hauschildt & Baron 1999). In the dOS method the relevant lines (defined by a suitable criterion) are first selected from master spectral line databases which include all available lines. The line selection procedure will typically select more lines than can be stored in memory and thus temporary line database files are created during the line selection phase. The file size of the temporary database can vary, in theory, from zero to the size of the original database or larger, depending on the amount of data stored for the selected lines and their number. For large molecular line databases this can easily lead to temporary databases of several GB in size. This is in part due the storage for the temporary line database: its data are stored for quick retrieval rather than in the compressed space saving format of the master line databases. The number and identity of lines that are selected from the master databases depends on the physical conditions for which the line opacities are required (temperatures, pressures, abundances for a model atmosphere) and thus the line selection has to be repeated if the physical conditions change significantly.

The temporary line databases are used in the next phase to calculate the actual line opacity for each wavelength point in a prescribed (arbitrary) wavelength grid. This makes it possible to utilize detailed line profiles for each considered spectral line on arbitrary wavelength grids. For each wavelength grid point, all (selected) lines within prescribed search windows (large enough to include all possibly important lines but small enough to avoid unnecessary calculations) are included in the line opacity calculations for this wavelength point.

## 8. I/O Parallel Algorithms

There are currently a large number of significantly different types of parallel machines in use, ranging from clusters of workstations or PCs to teraflop parallel supercomputers. These systems have very different performance characteristics that need to be considered in parallel algorithm design. For the following discussion we assume this abstract picture of a general parallel machine: The parallel system consists of a number of processing elements (PEs), each of which is capable of executing the required parallel codes or sections of parallel code. Each PE has access to (local) memory and is able to communicate data with other PEs through a communication device. The PEs have access to both local and global filesystems for data storage. The local filesystem is private to each PE (and inaccessible to other PEs), and the global filesystem can be accessed by any PE. A PE can be a completely independent computer such as a PC or workstation (with single CPU, memory, and disk), or it can be a part of a shared memory multi-processor system. We assume that the parallel machine has both global and local logical filesystem storage available (possibly on the same physical device). The communication device could be realized, for example, by standard Ethernet, shared memory, or a special-purpose high speed communication network.

In the following description of the 2 algorithms that we consider here we will make use of the following features of the line databases:

* master line databases:

1. are globally accessible to all PEs
    2. are are sorted in wavelength
    3. can be accessed randomly in blocks of prescribed size (number of lines)
* selected line (temporary) databases:

1. the wavelength grid is known during line selection (not absolutely required but helpful)
    2. have to be sorted in wavelength
    3. can be accessed randomly in blocks of prescribed fixed size (number of lines)
    4. are stored either globally (one database for all PEs) or locally (one for each PE)
    5. are larger than the physical memory of the PEs

### 8.1. Global Temporary Files (GTF)

The first algorithm we describe relies on global temporary databases for the selected lines. This is the algorithm that was implemented in the versions of PHOENIX discussed above. In the general case of $`N`$ available PE’s, the parallel line selection algorithm uses one PE dedicated to I/O and $`(N1)`$ line selection PEs. The I/O PE receives data for the selected lines from the line selection PEs, assembles them into properly sorted blocks of selected lines, and writes them into the global temporary database for later retrieval. The $`(N1)`$ line selection PEs each read one block from a set of $`(N1)`$ adjacent blocks of line data from the master database, select the relevant lines, and send the necessary data to the I/O PE. Each line selection PE will select a different number of lines, so the I/O PE has to perform administrative work to construct sorted blocks of selected lines that it then writes into the global temporary database. The block sizes for the line selection PEs and for the global temporary database created by the I/O PE do not have to be equal, but can be chosen for convenience. The blocks of the master line database are distributed to the $`(N1)`$ line selection PEs in a round robin fashion. Statistically this results in a balanced load between the line selection PEs due to the physical properties of the line data.

After the line selection phase is completed, the global temporary line database is used in the line opacity calculations. If each of the $`N`$ PEs is calculating line opacities (potentially for different sets of wavelengths points or for different sets of physical conditions), they all access the temporary database simultaneously, reading blocks of line data as required. In most cases of practical interest, the same block of line data will be accessed by several (all) PEs at roughly the same time. This can be advantageous or problematic, depending on the structure of the file system on which the database resides. The PEs also cache files locally (both through the operating system and in the code itself through internal buffers) to reduce explicit disk I/O. During the line selection phase the temporary database is a write-only file, whereas during the opacity calculations the temporary database is strictly read-only.

The performance of the GTF algorithm depends strongly on the performance of the global file system used to store the temporary databases and on the characteristics of the individual PEs. This issue is discussed in more detail below.

### 8.2. Local Temporary Files (LTF)

The second algorithm we consider uses file systems that are local to each PE; such local file systems exist on many parallel machines, including most clusters of workstations. This algorithm tries to take advantage of fast communication channels available on parallel computers and utilizes local disk space space for temporary line list files. This local disk space is frequently large enough for the temporary line database and may have high local I/O performance. In addition, I/O on the local disks of a PE does not require any inter-PE communication, whereas globally accessible filesystems often use the same communication channel that explicit inter-PE communication uses. The latter can lead to network congestion if messages are exchanged simultaneously with global I/O operations.

For the line selection, we could use the algorithm described above with the difference that the I/O PE would create a single global (or local) database of selected lines. After the line selection is finished, the temporary database could then simply be distributed to all PEs, and stored on their local disk for subsequent use. This is likely to be slower in all cases than the GTF algorithm.

Instead, we use a “ring” algorithm that creates the local databases directly. In particular, each of the $`N`$ PEs selects lines for one block from the master database (distributed in round robin fashion between the PEs). After the selection for this one block is complete, each PE sends the necessary data to it next neighbor; PE $`i`$ sends its results to PE $`i+1`$ and simultaneously, receives data from the previous PE in the ring (the ring is periodic, i.e., PE $`N`$ sends data to PE $`0`$). This can be easily realized using the MPI\_SENDRECV call, which allows the simultaneous sending and receiving of data for each member of the ring. Each PE stores the data it receives into a buffer and the process is repeated until the all selected lines from the $`N`$ blocks are buffered in all $`N`$ PEs. The PEs then transfer the buffered line data into their local temporary databases. This cycle is repeated until the line selection phase is complete. The line opacity calculations will then proceed in the same way as outlined above, however, the temporary line databases are now local for each PE.

This approach has the advantage that accessing the temporary databases does not incur any (indirect) Network File System (NFS) communication between the PEs as each of them has its own copy of the database. However, during the line selection phase a much larger amount of data has to be communicated over the network between the PEs because now each of them has to “know” all selected lines, not only the I/O PE used in the first algorithm.

The key insight here is that low-cost parallel computers constructed out of commodity workstations typically have a very fast communication network (100 Mbs to 1 Gbs) but relatively slow NFS performance. This means that trading off the extra communication for fewer NFS disk accesses in the LTF algorithm is likely to give better performance.

### 8.3. Results

The performance of the GTF and LTF algorithms will depend strongly on the type of parallel machine used. A machine using NFS with fast local disks and communication is likely to perform better with the LTF algorithm. However, a system with a fast (parallel) filesystem and fast communication can perform better with the GTF algorithm. In the following we will consider test cases run on very different machines:

1. A cluster of Pentium Pro 200MHz PCs with 64MB RAM, SCSI disks, 100Mbs full-duplex Ethernet communication network, running Solaris 2.5.1.
2. An IBM SP system with 200MHz Power3 CPUs, 512MB RAM per CPU, 133 MB/s switched communication network, 16 PE IBM General Purpose File System (GPFS) parallel filesystem, running AIX 4.3.

We have run 2 test cases to analyze the behavior of the different algorithms on different machines. The small test case was designed to execute on the PPro system. It uses a small line database (about 550MB) with about 35 million lines of which about 7.5 million lines are selected. The second test uses a database which is about 16 times larger ($`9`$GB) and also selects 16 times more lines. This large test could not be run on the PPro system due to file size limitations and limited available disk space. The line opacity calculations were performed for about 21,000 wavelength points that are representative of typical calculations. The tests were designed for maximum I/O usage and are thus extreme cases. In practical applications the observed scaling is comparable to or better than that found for these tests and appears to follow the results shown here rather well.

The results for the line selection procedure on the PPro/Solaris system are shown in Fig. 4. It is apparent from the figure that the GTF approach delivers higher relative speedups that translate into smaller execution times for more than 2 PEs. For serial (1 PE) and 2 PE parallel runs the LTF line selection is substantially faster than the GTF algorithm. The reason for this behavior can be explained by noting that the access of the global files is done through NFS mounts that use the same network as the MPI messages. Therefore, $`n1`$ PEs request different data blocks from the NFS server (no process was run on the NFS server itself) and send their results to the I/O PE, which writes it out to the NFS server. In the LTF algorithm, each PE reads a different input block from the NFS server and then sends its results (around the ring) to all other PEs. Upon receiving data from its left neighbor, a PE writes it to local disk. This means that the amount of data streaming over the network can be as much as twice as high for the LTF compared to the GTF algorithm. This increases the execution time for the LTF approach if the network utilization is close to the maximum bandwidth. In this argument we have ignored the time required to write the data to local disks, which would make the situation worse for the LTF approach.

The situation is very different for the calculation of the line opacities (which uses only the local temporary database), as also shown in Fig. 4. Now the LTF approach scales well (up to the maximum of 8 available machines) whereas the GTF algorithm hardly scales to more than 2 PEs. The absolute execution times for the LTF approach are up to a factor of 4 smaller (more typical are factors around 2) than the corresponding times for the GTF algorithm (the GTF run with 8 PEs required roughly as much execution time as the LTF run with 1 PE!). The reason for this is clearly the speed advantage of local disk I/O compared to NFS based I/O in the GTF code. If more PEs are used in the GTF line opacity approach, the network becomes saturated quickly and the PEs have to wait for their data (the NFS server itself was not the bottleneck). The LTF approach will be limited by the fact that as the number of PEs gets larger, the efficiency of disk caching is reduced and more physical I/O operations are required. Eventually this will limit the scaling as the execution time is limited by physical I/O to local disks.

We ran the same (small) test on an IBM SP for comparison. The tests were run on a non-dedicated production system and thus timings are representative of standard operation conditions and not optimum values. The global files were stored on IBM’s General Purpose File System (GPFS), which is installed on a number of system-dedicated I/O PEs replacing the NFS fileserver used on the PPro/Solaris system. GPFS access is facilitated through the same “switch” architecture that also carries MPI messages on the IBM SP. For the small test the results are markedly different from the results for the PPro/Solaris system. The LTF algorithm performs significantly better for all tested configurations, however, scaling is very limited. The GTF code does not scale well at all for this small test on the IBM SP. This is due to the small size, i.e., processing is so fast (nearly 100 times faster than on the PPro/Solaris system) so that, e.g., latencies and actual line selection calculations overwhelm the timing. The IBM SP has a very fast switched communications network that can easily handle the higher message traffic created by the LTF code. This explains why the LTF line selection executes faster and scales better for this small test on the IBM compared to the PPro/Solaris system.

The line opacity part of the test performs distinctively different in the IBM SP compared to the PPro/Solaris system. In contrast to the latter, the IBM SP delivers better performance for the GTF algorithm compared to the LTF code. The scaling of the GTF code is also significantly better than that of the LTF approach. This surprising result is a consequence of the high I/O bandwidth of the GPFS running on many I/O PEs, the I/O bandwidth available to GPFS is significantly higher than the bandwidth of the local disks (including all filesystem overheads etc). The I/O PEs of the GPFS can also cache blocks in their memory which can eliminate physical I/O to a disk and replace access by data exchange over the IBM “switch”. Note that the test was designed and run with parameters set to maximize actual I/O operations in order to explicitly test this property of the algorithms.

The results of the large test case, for which the input file size is about 16 times bigger, are very different for the line selection, as shown in Fig. 5. Now the GTF algorithm executes much faster (factor of 3) than the LTF code. This is probably caused by the larger I/O performance of the GPFS that can easily deliver the data to all PEs and the smaller number or messages that need to be exchanged in the GTF algorithm. The drop of performance at 32 PE’s in the GTF line selection run could have been caused by a temporary overload of the I/O subsystem (these tests were run on a non-dedicated machine). In contrast to the previous tests, the LTF approach does not scale in this case. This is likely caused by the large number of relatively small messages that are exchanged by the PEs (the line list master database is the same as for the GTF approach, so it is read through GPFS as in the GTF case). This could be improved, e.g., by choosing larger block sizes for data sent via MPI messages, however, this will have the drawback of more memory usage and larger messages are more likely to block than small messages that can be stored within the communication hardware (or driver) itself.

The situation for the line opacities is also shown in Fig. 5. The scaling is somewhat worse due to increased physical I/O (the temporary files are about 16 times larger as in the small test case). This is more problematic for the LTF approach which scales very poorly for larger numbers of PEs as the maximum local I/O bandwidth is reached far earlier than for the GTF approach. This is rather surprising as conventional wisdom would favor local disk I/O over global filesystem I/O on parallel machines. Although this is certainly true for farms of workstations or PCs, this is evidently not true on high-performance parallel computers with parallel I/O subsystems.

## 9. Conclusions

While the parallelization of a large numerical code such as PHOENIX would ideally be performed at the compiler level, it is not possible with modern compilers to construct efficient code with good scaling properties. Thus, individual modules of the code have to be examined and parallelized separately. With MPI this task is manageable and allows the efficient use of both commodity parallel clusters and modern massively parallel supercomputers.

We found good speedup up to about 64 PEs for a typical supernova calculation. However, for more than 64 PEs the communication, synchronization, and loop overheads begin to become significant and it is not economical to use more than 128 PEs For static models and opacity tables we are able to use very large numbers of PEs with scaling at close to the theoretical maximum. Fully 3-D calculations will present yet another parallelization challenge.

#### Acknowledgments.

We thank our many collaborators who have contributed to the development of PHOENIX, in particular France Allard, Jason Aufdenberg, Andreas Schweitzer, Travis Barman, Eric Lentz, Ian Short, David Branch, Sumner Starrfield, and Steve Shore. This work was supported in part by NASA, the NSF, IBM, and the Pôle Scientifique de Modélisation Numérique at ENS-Lyon and the calculations were performed at NERSC and SDSC.

## References

Baron, E., Hauschildt, P. H., Nugent, P., & Branch, D. 1996, MNRAS, 279, 779

Baron, E., & Hauschildt, P. H. 1998, ApJ, 495, 370

Gropp, W., Lusk, E., Doss, N., & Skjellum, A. 1996, MPICH Model MPI Implementation Reference Manual. Technical report, Argonne National Laboratory, Argonne

Message Passing Interface Forum 1995, MPI: A Message-Passing Interface Standard, Version 1.1, (Knoxville, TN: Univ. of Tennessee)

Mihalas, D., & Mihalas, B. W. 1984, Foundations of Radiation Hydrodynamics, (Oxford: Oxford University)

Hauschildt, P. H. 1992, JQSRT, 47, 433

Hauschildt, P. H. 1993, JQSRT, 50, 301

Hauschildt, P. H., & Baron, E. 1995, JQSRT, 54, 987

Hauschildt, P. H., & Baron, E. 1999, Journal of Computational and Applied Mathematics, 102, 41

Hauschildt, P. H., Baron, E., & Allard, F. 1997, ApJ, 483, 390

Hauschildt, P. H., Baron, E., Starrfield, S., & Allard, F. 1996, ApJ, 462, 386

Hauschildt, P. H., Störzer, H., & Baron, E. 1994, JQSRT, 51, 875

Kurucz, R. L. 1992, Rev. Mex. Astron. Astrofis., 23, 181

Olson, G. L., & Kunasz, P. B. 1987, JQSRT, 38, 325

Peytremann, E. 1974, A&A, 33, 203

Rybicki, G. B., & Hummer, D. G. 1991, A&A, 245, 171

Short, C. I., Hauschildt, P. H., & Baron, E. 1999, ApJ, 525, 375


format: text

text:
# Comment on “Ruling out chaos in compact binary systems”

\[

## Abstract

In a recent Letter, Schnittman and Rasio  argue that they have ruled out chaos in compact binary systems since they find no positive Lyapunov exponents. In stark contrast, we find that the chaos discovered in the original paper under discussion, J.Levin, PRL, 84 3515 (2000), is confirmed by the presence of positive Lyapunov exponents.

04.30.Db,97.60.Lf,97.60.Jd,95.30.Sf,04.70.Bw,05.45 \]

Intuitively, black hole binaries are obvious candidates for chaotic dynamics. Firstly, they are highly non-linear systems and chaos is an expression of extreme non-linearity. Secondly, there are isolated unstable orbits around a Schwarzschild black hole, and unstable orbits are a red flag for the onset of chaos. Surely the orbits only become more unstable and more complex when there are two black holes.

There is no question that there is chaos in black hole binaries for some range of parameters . Yet the authors of Ref.  claim to rule out chaos by finding no positive Lyapunov exponents along the fractal of Ref. . Before carrying out a detailed calculation, there is reason to be suspicious of the absence of positive Lyapunov exponents. Even the Schwarzschild solution has positive Lyapunov exponents for unstable circular orbits . Though completely integrable and so non-chaotic, the orbits are unstable. As long as there are unstable circular orbits in the PN equations, there are positive Lyapunov exponents for these orbits. This is not to say that the orbits will definitely be chaotic, just that there will be instability and a positive Lyapunov exponent.

Chaos was detected in the Post-Newtonian (PN) expansion of the two-body problem for rapidly spinning bodies  using the method of fractal basin boundaries. One fully expects positive Lyapunov exponents at least very near the fractal basin boundaries, if not elsewhere in the phase space. We sampled some orbits near the boundary and find positive Lyapunov exponents as illustrated by the positive slope in Fig. 1. The reason for the discrepancy seems to be that the authors of Ref.  define the maximum exponent as “the Cartesian distance between the dimensionless 12-component coordinate vectors \[$`\stackrel{}{𝐫},\dot{\stackrel{}{𝐫}},\stackrel{}{𝐒}_\mathrm{𝟏},\stackrel{}{𝐒}_\mathrm{𝟐}`$\] and \[$`\stackrel{}{𝐫}^{},\dot{\stackrel{}{𝐫}^{}},\stackrel{}{𝐒}_\mathrm{𝟏}^{},\stackrel{}{𝐒}_\mathrm{𝟐}^{}`$\] of two nearby trajectories” . However, this is not a Lyapunov exponent. The Lyapunov exponent is obtained by linearizing the equations of motion about a given trajectory. An approximation to the Lyapunov exponent can be made using the Cartesian distance between two trajectories by continually rescaling the shadow trajectory so that it is always infinitesimally close to the original trajectory as was done by . However, the result does depend on the rescaling and can give false answers. A more thorough treatment of the Lyapunov exponents and their interpretation is given elsewhere .

Chaos has not been ruled out rather it has been confirmed by positive Lyapunov exponents. In fact, an even stronger claim can be made. The two-body problem has not been solved in relativity and it would not be too daring to conjecture that it will never be solved. The system shows every indication of being fully non-integrable.

Regardless of how we’d like nature to behave to make our lives as observers easier, she may not be so accommodating. Understandably there is a sense of urgency to the removal of any obstacles to the detection of gravitational waves. However, chaos need not be a terrible obstacle. There are ways that chaos can enhance signals or amplify salient features in the sky  to aid our quest to observe gravitational waves.


text:
# Quasiequilibrium sequences of synchronized and irrotational binary neutron stars in general relativity. III. Identical and different mass stars with 𝛾=2

## I Introduction

Coalescing binary neutron stars are expected to be one of the most promising sources of gravitational waves that could be detected by the ground based, kilometer size laser interferometers such as the Laser Interferometric Gravitational wave Observatory (LIGO), VIRGO, GEO600, and TAMA300 111For TAMA300, data taking has started in summer 1999, and the first results of the data analysis have been published Tagoshi01 ., and also are considered as one of candidates of gamma-ray burst source NaraPP92 .

Due to the emission of gravitational radiation, binary neutron stars decrease their orbital separations and finally merge. When discussing the evolution of the system, it is convenient to separate it into three phases. The first one is the inspiraling phase in which the orbital separation is much larger than the neutron star radius, and the post-Newtonian (PN) expansion constitutes an excellent approximation. Recently, two groups succeeded in deriving the 3PN equation of motion of point-mass binary systems DamourJS ; BlanchetG , and the equivalence of the results between these groups is shown in DamourJS01 ; AndradeBF01 . The second stage is the intermediate phase in which the orbital separation becomes only a few times larger than the radius of a neutron star, so that hydrodynamics as well as general relativity play an important role. In this phase, since the shrinking time of the orbital radius due to the emission of gravitational waves is still larger than the orbital period, it is possible to approximate the state as quasiequilibrium. The final stage is the merging phase in which the two stars coalesce dynamically. As in the intermediate phase, since hydrodynamics as well as general relativity play an important role, fully relativistic hydrodynamical treatments are required in this phase which therefore pertains to the field of numerical relativity. The first successful computations of the evolution of binary neutron stars from their innermost stable circular orbits (ISCO) to black hole or massive neutron star formation have been performed by Shibata and Uryu Shiba99 ; ShibaU00 . Other efforts are presented in OoharN99 ; FontGIMRSSST02 .

The present article belongs to the series of works BonazGM99a , GourGTMB01  (hereafter Paper I), TanigGB01  (hereafter Paper II), TanigG02 , devoted to the intermediate phase. This stage is interesting because we may get informations about equation of state of neutron stars through the ISCO FaberGRT02 . Furthermore, it is important from a numerical point of view since it provides initial data for the dynamical simulation in the merging phase Shiba99 ; ShibaU00 . Then numerous theoretical efforts are devoted in this phase including (semi-)analytic Newtonian LaiRS93 ; TanigN00  and post-Newtonian LombaRS97 ; TanigS97 ; ShibaT97 ; Tanig99  approaches and numerical Newtonian TanigGB01 ; TanigG02 ; HachiE84a ; HachiE84b ; UryuE98 , post-Newtonian Shiba96 , and general relativistic BonazGM99a ; GourGTMB01 ; BaumgCSST97 ; MarroMW98 ; MarroMW99 ; UryuE00 ; UryuSE00 ; UsuiUE99 ; UsuiE02  ones. Among these studies, we concentrate on numerical computations in the general relativistic framework. Kochaneck Kocha92  and Bildsten and Cutler BildsC92  have shown that the gravitational-radiation driven evolution is too rapid for the viscous forces to synchronize the spin of each neutron star with the orbit as they do for ordinary stellar binaries. Rather, the viscosity is negligible and the fluid velocity circulation (with respect to some inertial frame) is conserved in these systems. Provided that the initial spins are not in the millisecond regime, this means that close binary configurations are well approximated by zero vorticity (i.e. irrotational) states. Formulations for irrotational binaries in general relativity have been developed by several authors BonazGM97 ; Asada98 ; Shiba98 ; Teuko98  (see Appendix A of Paper I for a comparison between them). In addition to the quasiequilibrium and irrotational assumptions, we use the approximation of a conformally flat spatial 3-metric, introduced by Isenberg IsenbN80  and Wilson and Mathews WilsoM89  (hereafter IWM approximation; see Sec. IV.C of FriedUS02  and Sec. III.A of Paper I for a discussion).

Until now, several groups have produced quasiequilibrium configurations of binary neutron stars, as listed above. Among them, there are results for synchronized BaumgCSST97 ; MarroMW98  and irrotational BonazGM99a ; GourGTMB01 ; MarroMW99 ; UryuE00 ; UryuSE00  rotation states in the general relativistic framework, within the IWM approximation, and those for synchronized one within some axisymmetric approximation UsuiUE99 ; UsuiE02 . However they all deal with identical star binaries and calculations for different mass binary system in general relativity have not been performed yet (in Newtonian theory they have been performed in the synchronized case HachiE84b ; TanigG02 , and recently in the irrotational one TanigG02 ).

We will present here the first numerical results of relativistic binary systems composed of different mass stars. In addition, we give results for systems of identical stars which extend those already presented in Paper I. As a first step in the study of different mass systems, we restrict our computations to a polytropic equation of state, with the adiabatic index $`\gamma =2`$.

The plan of the article is as follows. A brief overview of our method is given in Sec. II. In Sec. III, new tests of the numerical code are shown, which were not given in Papers I and II. We present the numerical results in Sec. IV, and discuss them in Sec. V. Finally Sec. VI is devoted to the summary.

Throughout this article, we adopt geometrical units: $`G=c=1`$, where $`G`$ and $`c`$ denote the gravitational constant and speed of light, respectively.

## II Method

The reader is referred to Secs. II, III, and IV of Paper I for the complete set of equations governing perfect fluid binary stars under the assumption of quasiequilibrium and conformally flat spatial metric, as well as the numerical method which we use. We simply give here some outline of the method.

Let us summarize the assumptions expressed in Sec. I. The first one is quasiequilibrium, resulting from the fact that the time scale of orbital shrinking is larger than that of the orbital revolution until the ISCO. This assumption is taken into account by demanding that the spacetime is endowed with a helical Killing vector

$$𝒍=\frac{}{t}+\mathrm{\Omega }\frac{}{\phi },$$
(1)
where $`\mathrm{\Omega }`$ is the orbital angular velocity and the vectors $`/t`$ and $`/\phi `$ are associated with the time coordinate $`t`$ and the azimuthal coordinate $`\phi `$ of an asymptotical inertial observer. The second assumption regards the matter stress-energy tensor, which we assume to have the perfect fluid form:

$$T_{\mu \nu }=(e+p)u_\mu u_\nu +pg_{\mu \nu },$$
(2)
where $`e`$ denotes the fluid proper energy density, $`p`$ the fluid pressure, $`u_\mu `$ the fluid 4-velocity, and $`g_{\mu \nu }`$ the spacetime metric. The third assumption concerns the rotation state of the binary system. Although the realistic rotation state will be an irrotational one (cf. Sec. I), we consider both cases of synchronized and irrotational states in order to exhibit their differences. The fourth assumption is on the equation of state. For simplicity, we use a polytropic one,

$$p=\kappa n^\gamma ,$$
(3)
where $`n`$ is the fluid baryon number density and $`\kappa `$ and $`\gamma `$ are some constants. Here we assume that the two neutron stars have the same equation of state, i.e., the constants $`\kappa `$ and $`\gamma `$ are identical for both stars. The fifth and last assumption is the IWM one, namely a conformally flat spatial metric. Then the full spacetime metric takes the form:

$$ds^2=(N^2B_iB^i)dt^22B_idtdx^i+A^2f_{ij}dx^idx^j,$$
(4)
$`N`$ being the lapse function, $`B^i`$ the shift vector, $`A`$ the conformal factor, and $`f_{ij}`$ the flat spatial metric. Under these assumptions, the equations for fluid motion and gravitational field are presented in Paper I.

## III Code tests

In this series of studies BonazGM99a ; GourGTMB01 ; TanigGB01 ; TanigG02 , we use a multidomain spectral method with surface-fitted coordinates BonazGM98 ; BonazGM99b ; GrandBGM01 . This method is implemented by a numerical code constructed upon the C++ library Lorene lorene . Numerous tests of the code have been already presented in Papers I and II. In particular, the comparison with results from other groups BaumgCSST97 ; UryuE00  has been given in Sec. V.D of Paper I for the identical mass case. We will give here new tests based on the relative error in a general relativistic generalization of the virial theorem.

In Newtonian gravity, the virial theorem has proved to be useful to check the global error in numerical solutions for stationary fluid systems. A general relativistic version of the virial theorem has been obtained by Gourgoulhon and Bonazzola GourgB94  for stationary spacetimes. It has been recently extended to binary star spacetimes within the IWM approximation by Friedman, Uryu and Shibata FriedUS02  222The spacetime generated by a binary system is not stationary, due to gravitational radiation. However, within the IWM approximation, the gravitational radiation is neglected in the global spacetime dynamics, so that one is able to recover the virial theorem.. The virial relation is equivalent to

$$M_{\mathrm{ADM}}M_{\mathrm{Komar}}=0,$$
(5)
where $`M_{\mathrm{ADM}}`$ is the ADM mass:

$$M_{\mathrm{ADM}}=\frac{1}{2\pi }_{\mathrm{}}\overline{}^iA^{1/2}𝑑S_i,$$
(6)
and $`M_{\mathrm{Komar}}`$ is a Komar-type mass, defined by

$$M_{\mathrm{Komar}}=\frac{1}{4\pi }_{\mathrm{}}\overline{}^iN𝑑S_i.$$
(7)
The virial theorem obtained by Friedman, Uryu and Shibata FriedUS02  (see also Eq. (5.7) of Ref. ShibaU01 ) writes

$`VE(FUS)`$ $`=`$ $`{\displaystyle \left[2NA^3S+\frac{3}{8\pi }NA^3K_i^jK_j^i+\frac{1}{4\pi }NA(\overline{}_i\beta \overline{}^i\beta \overline{}_i\nu \overline{}^i\nu )\right]}`$ (8)
$`=`$ $`0.`$
By a straightforward manipulation, this integral can be recast in the form of the virial theorem as obtained by Gourgoulhon and Bonazzola GourgB94 :

$`VE(GB)`$ $`=`$ $`{\displaystyle \left[2A^3S+\frac{3}{8\pi }A^3K_i^jK_j^i+\frac{1}{4\pi }A(\overline{}_i\beta \overline{}^i\beta \overline{}_i\nu \overline{}^i\nu 2\overline{}_i\beta \overline{}^i\nu )\right]}`$ (9)
$`=`$ $`0.`$
Let us stress that the above identity has been derived by Gourgoulhon and Bonazzola GourgB94  only for stationary spacetimes and that its validity for IWM spacetimes with helical Killing vector has been obtained by Friedman, Uryu and Shibata FriedUS02 .

As error indicators of our numerical solutions, we have evaluated the quantities

$$\left|\frac{M_{\mathrm{ADM}}M_{\mathrm{Komar}}}{M_{\mathrm{ADM}}}\right|,$$
(10)
$$\left|\frac{VE(FUS)}{M_{\mathrm{ADM}}}\right|,$$
(11)
$$\left|\frac{VE(GB)}{M_{\mathrm{ADM}}}\right|.$$
(12)
The results are presented in Fig. 1 as a function of the orbital separation. One can see from this figure that the typical virial error is $`1\times 10^5`$ through the sequence.

## IV Numerical results

In this section, we present the numerical results for quasiequilibrium sequences of constant baryon number binaries (evolutionary sequences), for both synchronized and irrotational cases and for the adiabatic index $`\gamma =2`$. Although the interesting case is the irrotational one, as mentioned in Sec. I, we have computed configurations for both rotation states in order to compare their properties. The number of computational domains is 4 or 5 for each star and the space around it. The domains have an onion-like structure, being centered on each star: the innermost domain covers the star, two or three domains (shells) are placed outside it, and the external domain extends to infinity thanks to some compactification (cf. Sec. IV.A of Paper I for all details). The total number of domains is therefore $`2\times 4=8`$ or $`2\times 5=10`$. We use 8 domains for close orbital separation and 10 domains for large orbital separation. The numbers of collocation points in each domains are $`N_r\times N_\theta \times N_\phi =33\times 25\times 24`$ or $`25\times 17\times 16`$, where $`N_r`$, $`N_\theta `$, and $`N_\phi `$ denote the number of collocation points of the radial, the polar, and the azimuthal directions, respectively.

We parametrize the different mass systems by the values of $`(M/R)_{\mathrm{},\mathrm{star}1}`$ vs. $`(M/R)_{\mathrm{},\mathrm{star}2}`$, where $`(M/R)_{\mathrm{}}`$ denotes the compactness of the isolated spherical star with the same baryon number.

First of all, the relative changes in central energy density for the cases of $`M/R=0.14`$ vs. 0.14 and 0.14 vs. 0.16 are depicted in Fig. 2 as a function of the orbital separation. (Hereafter we abbreviate $`(M/R)_{\mathrm{},\mathrm{star}1}`$ vs. $`(M/R)_{\mathrm{},\mathrm{star}2}`$ as $`M/R`$ for simplicity.) One can see from these figures that the central energy density decreases monotonically in both synchronized and irrotational cases 333In Paper I (Table III), as well as in BonazGM99a  (Fig. 1), there was a slight increase in the central energy density for irrotational binaries at the large and medium separations. However, from Paper II we have improved the method to determine the surface of the star, which leads to a better accuracy in the present work.. The decrease of the central energy density for the synchronized case is about 10 % or more at the end of the sequence, while for the irrotational case, it is only about 1 %.

Next, the ADM mass of the binary system \[Eq. (6)\] is shown in Figs. 35 as a function of the orbital angular velocity, for different mass ratios and different states of rotation. One can clearly see the turning point in the synchronized cases for identical mass binaries. For more compact stars the turning points appear relatively earlier than those of less compact ones. On the other hand, we do not find any turning point in irrotational binaries. Moreover, it becomes more difficult for the turning point to appear even in the synchronized case when the mass difference in the binary system increases.

In Fig. 6, we show the total angular momentum of the binary system as a function of the orbital angular velocity for the compactness $`M/R=0.18`$ vs. 0.18. For synchronized binaries, the turning point of the total angular momentum coincides with that of the ADM mass within the numerical accuracy.

The equatorial to polar ratios of the radial derivative of the enthalpy are shown as a function of the orbital separation normalized by the radius of a star to the companion star in Figs. 7 and 8. As explained in Papers I and II, this quantity

$$\chi :=\frac{(H/r)_{\mathrm{eq},\mathrm{comp}}}{(H/r)_{\mathrm{pole}}}$$
(13)
is a good dimensionless indicator to detect the appearance of cusp at the stellar surface, which corresponds to the mass shedding limit. It is worth to plot the variation of $`\chi `$ with respect to the orbital separation normalized by the sum of the radius of each star toward its companion: $`d/(a_1+a_1^{})`$, where $`d`$ denotes the orbital separation between the two coordinate center (maximum of density), $`a_1`$ the radius of the less massive star in the direction of its companion, and $`a_1^{}`$ that of the more massive star. Note here that the separation $`d`$ and the radii $`a_1`$ and $`a_1^{}`$ are coordinate lengths. If the quantity $`\chi `$ becomes zero, we can conclude that the sequence ends by mass shedding. On the other hand, if the quantity $`d/(a_1+a_1^{})`$ becomes unity, the two stars will contact each other. It is found from Fig. 7 that for synchronized identical mass binaries the sequences seem to terminate by the contact between the two stars (the quantity $`\chi `$ will become zero at the same time). The important feature of the sequences is that they take almost the same track regardless of the compactness. On the contrary, one can see from Fig. 7 that for irrotational binaries the sequences seem to terminate by a cusp point (mass shedding). Such a behavior is enhanced for more compact cases. In Fig. 8, we compare the quantity $`\chi `$ of identical mass binary of the compactness $`M/R=0.14`$ with that of different mass binaries of the compactness $`M/R=0.14`$ vs. 0.15 and 0.14 vs. 0.16. It is found that the sequences seem to terminate by a cusp point of the less massive star for both synchronized and irrotational binaries when we extrapolate the lines, while the more massive star keeps its shape near spherical one (see the next paragraph). When we increase the difference in mass between two stars, this tendency is enhanced.

It is worth to mention that as soon as the two masses differ, even slightly, the sequence of a synchronized binary system terminates by a cusp point, because the line for the less massive star is always located below that for the identical mass stars which will reach $`\chi =0`$ at $`d/(a_1+a_1^{})=1`$. This implies that the line for the less massive star seems to always reach $`\chi =0`$ at $`d/(a_1+a_1^{})>1`$. Note here that the sequences for identical mass binaries and less massive stars in different mass binaries should reach the point $`\chi =0`$. However, in our numerical method, we cannot treat a cuspy figure, because we adapt the innermost numerical domain to the surface of the star and the multidomain spectral method that we use assume that the domain boundaries are smooth surfaces. For very close configurations, it becomes therefore difficult to make the iterative procedure converge to a sufficient level, i.e. to have the relative difference $`\delta H`$ in the enthalpy field between two steps lower than $`10^5`$. This explains why we stop the sequences at $`\chi 0.4`$. We use $`\delta H=10^{10}`$ for very large separations, $`\delta H=10^810^7`$ for medium ones, and $`\delta H=10^610^5`$ for very close ones.

In Figs. 9 and 10, isocontours of the baryon density of various binary systems are shown. There are three panels in each figure. The left panel shows the result for identical mass binary with the compactness $`M/R=0.14`$ vs. 0.14. The center and right ones show the results for different mass binaries with the compactness $`M/R=0.14`$ vs. 0.15, and 0.14 vs. 0.16, respectively. In these panels, the left hand side star is the less massive star. It is clearly seen that the less massive star is tidally deformed and elongated while the more massive star relatively does not deviate from the spherical shape so much (see the right panel). The orbital separations of these isocontours are listed at the last line of each case in Tables below which corresponds to the closest figure we can calculate, except for the case of $`M/R=0.14`$ vs. 0.15 which we do not show in the present paper as a table.

Finally, our numerical results for constant baryon number sequences are presented in Tables 14. We represent the orbital separation by two quantities:

$`\overline{d}_G`$ $`=`$ $`{\displaystyle \frac{d_G}{R_{\mathrm{poly}}}},`$ (14)
$`\overline{d}`$ $`=`$ $`{\displaystyle \frac{d}{R_{\mathrm{poly}}}}.`$ (15)
$`d_G`$ is the coordinate separation between the centers of mass of each star (same definition as Eq. (107) of UryuE00  or Eq. (128) of Paper I), $`d`$ that between the two stellar centers defined as the maxima of the density field, and $`R_{\mathrm{poly}}`$ the length constructed from the polytropic constants $`\kappa `$ and $`\gamma `$ \[Eq. (3)\]

$$R_{\mathrm{poly}}:=\kappa ^{1/2(\gamma 1)}.$$
(16)
The dimensionless orbital angular velocity is defined by

$$\overline{\mathrm{\Omega }}:=\mathrm{\Omega }\kappa ^{1/2(\gamma 1)},$$
(17)
the dimensionless ADM mass by

$$\overline{M}:=M_{\mathrm{ADM}}\kappa ^{1/2(\gamma 1)},$$
(18)
the dimensionless baryon mass by

$$\overline{M}_\mathrm{B}:=M_\mathrm{B}\kappa ^{1/2(\gamma 1)},$$
(19)
and the dimensionless total angular momentum by

$$\overline{J}:=J\kappa ^{1/(\gamma 1)}.$$
(20)
The virial errors $`\overline{VE(M)}`$, $`\overline{VE(GB)}`$ and $`\overline{VE(FUS)}`$ are defined by Eqs. (10), (12) and (11) respectively. $`a_1`$ denotes the coordinate radius toward the companion star, $`a_2`$ the radius perpendicular to it in the orbital plane and $`a_3`$ the radius perpendicular to the orbital plane. $`\delta e_\mathrm{c}`$ is the relative change in central energy density: $`\delta e_\mathrm{c}:=(e_\mathrm{c}e_{\mathrm{c},\mathrm{}})/e_{\mathrm{c},\mathrm{}}`$, and the cusp indicator $`\chi `$ is defined by Eq. (13). In the Tables, the superscript “ ’ ” denotes the values for the more massive star. and the symbol $``$ denotes the turning point in the ADM mass along the sequence.

## V Discussion

In this section, we discuss about the behaviors of the axial ratios, the relative change in central energy density, and the turning point of the synchronized sequences.

In Figs. 11 and 12, the axial ratios are shown as a function of the cusp indicator $`\chi `$. In Fig. 11, we show the axial ratios $`a_2/a_1`$ and $`a_3/a_1`$ for identical mass binaries with the compactness $`M/R=0.12`$ vs. 0.12 and 0.16 vs. 0.16. One can see from this figure that each axial ratio has almost the same track regardless of the difference of the compactness, in particular for synchronized binary systems. Moreover, such a coincidence of the sequence appears in the cases of different mass binary systems (see Fig. 12). In those cases, the sequences almost coincide with each other regardless of the compactness of the companion star. Note that the sequences of the axial ratios for more massive stars end before $`\chi =0`$, while those for less massive stars should reach the point $`\chi =0`$. However we cannot treat a cuspy figure with our numerical method, because we adapt the innermost numerical domain to the surface of the star and the domain boundaries are assumed to be smooth (see discussion in Sec. IV). Therefore we had to stop the sequences at around $`\chi =0.4`$ for the less massive star.

The relative change in central energy density is presented as a function of $`\chi `$ in Fig. 13. It is clearly seen that the sequences for the synchronized case with the same compactness coincide with each other, while those with the different compactness split into curves with respect to each compactness. For example, the relative change for the identical binary system with the compactness $`M/R=0.14`$ coincides with that for the more massive star of the different mass binary with $`M/R=0.12`$ vs. 0.14 and that for the less massive star with $`M/R=0.14`$ vs. 0.16. However, the sequences for stars with $`M/R=0.12`$, 0.14, and 0.16 are different from each other.

Finally, let us comment on the behavior of the turning point of the ADM mass (binding energy) along a sequence. We never found any turning points for irrotational binary systems with the polytropic index $`\gamma =2`$. On the other hand, it appeared clearly for the synchronized case. In Fig. 14, the cusp indicator $`\chi `$ at the turning point of synchronized identical mass binary systems is plotted as a function of the compactness $`M/R`$. Since the accuracy of the determination of the ADM mass is order of $`10^5`$ in our calculations, the relative errors on the orbital angular velocity and the quantity $`\chi `$ at the turning point become a few percents. Even when one is taking this fact into account, it is possible to conclude that the line in Fig. 14 increases proportionally to the compactness of the stars. This behavior agrees with the value of $`\chi _{\mathrm{turning}}`$ that we have found in Newtonian calculations ($`M/R=0`$), namely $`\chi _{\mathrm{turning}}=0.2862`$ (Table I of Paper II).

## VI Summary

In the present paper, we have computed constant baryon number sequences of binary neutron stars in quasiequilibrium, in both cases of synchronized and irrotational motion and in both cases of identical and different mass systems. We have performed a general relativistic treatment, within the IWM approximation (conformally flat spatial metric). We have used a polytropic equation of state with an adiabatic index $`\gamma =2`$.

The summary of our results is that (1) Among the quasiequilibrium sequences we have calculated in the present article, only the synchronized identical star binaries terminate by the contact between two stars, while all the other types of sequences (synchronized different star binaries, irrotational identical star and different star binaries) end at the mass shedding points. Note here that we cannot conclude about the end points of non-quasiequilibrium sequences (i.e. dynamical sequences), such that binary systems without viscosity but with intrinsic non-aligned spins, those with slight viscosity and deviating from both the irrotational state and the synchronized one, those with infalling radial velocity, and so on. (2) For identical mass binary systems there is no turning point of the ADM mass in the irrotational case, while there is clearly one before the end point of the sequence in the synchronized case. (3) It is more difficult to see the turning points of the ADM mass for different mass binary systems than for identical ones. (4) The deformation of the star is determined by the orbital separation and the mass ratio and is not affected much by its compactness. (5) The decrease of the central energy density depends on the compactness of the star and not on that of its companion.

###### Acknowledgements.

KT is grateful to all members of the numerical relativity group at Albert-Einstein-Institut leaded by E. Seidel for their hearty hospitality during his stay. KT also would like to thank K. Uryu for useful discussions. The code development have been performed on SGI workstations purchased thanks to a special grant from the C.N.R.S. The numerical computations were mainly carried out on a SGI Origin2000 at Albert-Einstein-Institut. KT acknowledges a Grant-in-Aid for Scientific Research (No. 14-06898) of the Japanese Ministry of Education, Culture, Sports, Science and Technology.


text:
# No observational proof of the black-hole event-horizon

## 1 Introduction

It is generally believed that the compact component in X-ray binary systems is either a star possessing a material surface (a neutron star or a quark star) or a black hole, i.e. an object whose surface is formed by an event horizon. The evidence of the presence of a material surface is obtained from two types of observations. First, stable periodic pulsations of the X-ray emission indicate the presence of a strong ($`10^910^{15}`$ G) magnetic field, which by the virtue of the “no-hair” theorem excludes the presence of a black hole in the system. Second, observations of X-ray bursts - thermonuclear explosions occurring in matter accumulated at the surface of the compact object are an obvious proof of the absence of an event horizon.

Although X-ray pulsations or X-ray bursts indicate the presence of a solid surface, their absence does not prove the presence of an event horizon. “Absence of evidence is not evidence of absence”. However, there is direct evidence that compact bodies in X-ray binaries form (at least) two types of objects: their masses show a bimodal distribution (see Miller, Shahbaz & Nolan msn98 (1998)). Neutron star masses are all concentrated around the “canonical” value of 1.4 M whereas the second class of bodies, usually called “black hole candidates” have higher masses in the range of $``$ 5 to 18 M(see e.g. Narayan, Garcia & McClintock ngm01 (2001) and Greiner, Cuby & McCaughrean gcmc01 (2001)). The reason for suspecting the more massive bodies of being black holes is that their masses are higher than the maximum mass of a neutron (or quark) star, which is never larger than $`3`$ M (see e.g. Salgado et al. setal94 (1994)).

In general, the maximum mass of a compact body can be expressed as $`8.4\left(\epsilon _0/10^{14}\mathrm{g}\mathrm{cm}^3\right)^{1/2}`$M, where $`\epsilon _0`$ is the fiducial density above which the equation of state is taken to be described by a causality-limit equation of state (Rhoades & Ruffini rr74 (1974); Friedmann & Ipser fi87 (1987)). Bahcall, Lynn & Selipsky (bls90 (1990)) showed that stars with a material surface can have masses as high as 10 M, if one is willing to entertain configurations of sub-nuclear density. These are the so-called “Q-stars.” The mean-field description of nuclear interactions given by Bahcall et al. (bls90 (1990)) allows baryonic matter to have densities this low.

Although, as shown by Miller et al. (msn98 (1998)), it is unlikely that bodies with masses larger than 10 M are Q-stars—because this would require unrealistically low densities at which hadronic bulk matter would persist—“unlikely” is not a very satisfactory argument in favour of the black-hole existence. One would rather wish a “positive” proof of the event-horizon’s existence. This has been attempted by Ramesh Narayan and collaborators. The claim is that properties of Advection Dominated Accretion Flows (ADAFs; Abramowicz et al. aetal95 (1995); Narayan & Yi ny94 (1994), 1995a ,b ) can be used to prove the existence of event-horizons.

## 2 Proof by ADAFs

ADAFs describe accretion with very low radiative efficiency in which energy released by viscous torques removing angular momentum from the accreting matter is not radiated away but stored in the flow. If an ADAF forms around a black hole, the stored energy will be lost forever under the event horizon, whereas if the accreting body is a “star” this energy must be radiated away once matter lands on its surface. Therefore, the argument runs, black holes should be dimmer than neutron stars, quark stars, etc., if in both cases an ADAF is present.

The best systems in which this hypothesis could be tested are the so-called Soft X-ray Transients (SXTs) which are close binary systems undergoing rare and powerful outbursts but spending most of their life in a low luminosity quiescent state (see Tanaka & Shibazaki ts96 (1996) for a review). In SXTs, like in Low-Mass X-Ray Binaries (LMXBs) in general, a compact body accretes matter lost by a Roche-lobe filling low-mass stellar companion. The accreting matter forms a disc whose instabilities trigger outbursts (see Lasota l01 (2001) for a review of the instability model). Narayan, McClintock & Yi (nmy96 (1996); see also Lasota, Narayan & Yi lny96 (1996) and Narayan, Barret & McClintock nbm97 (1997)) proposed that quiescent SXT discs are truncated and that the inner accretion flow forms an ADAF. This hypothesis has been recently vindicated from the theoretical point of view by Dubus, Hameury & Lasota (dhl01 (2001)) and is supported by observations (see Done cd02 (2002) for a review).

Narayan, Garcia & McClintock (ngm97 (1997)) compared quiescent luminosities of SXTs supposed to contain black holes with those of neutron-star SXTs and realized that, in accordance with the prediction of the ADAF model, systems containing black-hole “candidates” are dimmer. They came to the conclusion that they found evidence for the presence of event horizons.

This conclusion has been challenged by Chen et al. (chetal98 (1998)) who asserted that the relative dimness of black-hole candidate systems was due solely to Narayan et al. (ngm97 (1997)) comparison method. Things were clarified by Lasota & Hameury (lh98 (1998)) who suggested comparing systems with similar orbital period on the assumption such systems would have similar accretion rates – the ADAF model asserting only that accreting black holes should be dimmer than neutron stars for the same accretion rate. The new method showed, however, the same effect (Lasota & Hameury lh98 (1998); Menou et al. metal99 (1999)), recently confirmed by Garcia et al. (getal01 (2001)): black holes (candidates) are dimmer than systems known to contain neutron stars, or at least stars with surface.

This is a very strong argument in favour of the presence of event horizons, in fact this is the most conservative conclusion. However, it is not a proof.

## 3 Arguments against evidence based on relative dimness of black hole candidates

The arguments against the claim that the relative dimness of black-hole candidates is the proof of existence of event horizons are of two, not unrelated, types. First, it has been argued that the accretion flow in quiescent SXTs are not represented by ADAFs.

Narayan & Yi (1995a ) and Blandford & Begelman (bb99 (1999)) argued (see however Paczyński bp98 (1998) and Abramowicz, Lasota & Igumenshchev ali00 (2000) for criticism of the argument) that ADAFs are subject to mass loss and therefore the dimness of quiescent SXTs could result from the low accretion rate onto the compact object - most of the matter being lost with the wind. However, as shown by Menou et al. (metal99 (1999)), such wind models do not offer an explanation of the luminosity difference between neutron-star systems and those presumed to contain black holes. In fact, these authors also pointed out that the quiescent luminosity of neutron-star binaries is not consistent with the assumption of a $`10\%`$ radiative efficiency. Since the attempt to apply to these systems the windy-ADAF model of Quataert & Narayan (qn99 (1999)) failed, they proposed that the action of a magnetic propeller could be answer. However, a compelling signature of this effect has yet to be found.

Despite of this, Abramowicz & Igumenshchev (2001) suggested that the observed differences between quiescent luminosities of accreting black holes and neutron stars is well explained by the occurrence in such systems of a CDAF (Convection Dominated Accretion Flow; see Narayan, Igumenshchev & Abramowicz 2000) instead of an ADAF. They found that for low viscosities accretion flows around compact bodies form ADAFs only in their innermost regions but are convectively dominated at radii $`R\mathrm{}>10^2R_\mathrm{S}`$ (where $`R_\mathrm{S}=2GM/c^2`$ is the Schwarzschild radius). In such flows emission comes mostly from the convective region; the radiative efficiency is independent of accretion rate and equals $`\epsilon _{\mathrm{BH}}=10^3`$. Assuming that the efficiency of accretion onto a neutron star is $`\epsilon _{\mathrm{NS}}0.1`$ one obtains the observed ratio between black-hole and neutron-star luminosities. Unfortunately this cannot be the correct explanation of the luminosity difference (Lasota l02 (2002)) because, as mentioned above, neutron stars in quiescent transient systems do not seem to accrete with a 0.1 efficiency.

Another class of argument asserts that X-rays in quiescent SXTs are not emitted by the accretion flow.

Brown, Bildsten, & Rutledge (bbr98 (1998)) suggested that, in neutron-star systems, most (or all) of the quiescent X-ray luminosity is not due to accretion but results from cooling of the neutron-star crust heated by nuclear reactions. This crust-cooling model does not seem to be in perfect agreement with observations showing two spectral components and a variable flux (see Rutledge et al. retal02 (2002) and references therein). If the crustal-cooling model were right it would imply different X-ray emission mechanisms for the two classes of quiescent SXTs. However, luminosity variations observed also in quiescent black-hole systems (see e.g. Garcia et al. getal01 (2001)) would rather suggest a common origin. Attempts to ascribe quiescent X-ray luminosity in black-hole systems to active stellar companions (Bildsten & Rutledge br00 (2000)) are not based on a sound theoretical foundation (Lasota l01 (2001)) and have been refuted by observations (Garcia et al. getal01 (2001)).

Menou (m02 (2001)) presented an argument based on the settling-flow model of Medvedev & Narayan (medna01 (2001)) in which the accretion flow arrives with very low angular momentum at the surface of a rapidly rotating compact object. The X-ray luminosity is then due to rotation-energy loss by the accreting body. This requires viscous contact between this body and the accreting matter. Menou (m02 (2001)) pointed out that if black-hole candidates had, contrary to neutron stars, radii smaller than the inner-most stable orbit the accretion flow would be supersonic and viscous contact impossible. Black-hole candidates would be dimmer because unable to lose their rotational energy.

Finally, we note that very compact objects with a surface would be dimmer than less compact objects, simply because of redshift and light bending. If the surface is below the photon orbit, the fraction of “outward moving” photons which escape to infinity is in the Schwarzschild metric

$$\frac{\mathrm{\Delta }\mathrm{\Omega }}{2\pi }=1\left[1\frac{27}{4}\frac{(1R_S/R)}{(R/R_S)^2}\right]^{1/2}.$$
(1)
For the lowest possible value for a causality-limit equation of state $`R/R_S=9/8`$, this factor and the redshift squared yield a luminosity at infinity which is equal to only 0.040 of the luminosity at the source.

## 4 Absence of X-ray bursts

Three of the SXTs show millisecond pulsations, and two of them are X-ray bursters. They all have very short orbital periods, 2 hr in the case of SAX J1808.4-3658 (Wijnands and van der Klis wvdk (1998); Chakrabarty and Morgan cm98 (1998)), 43.6 min for XTE J1751-305 (Markwardt et al. marketal (2002)), and 42 min for XTE J0929-314 (Galloway et al. galetal (2002)). It is perfectly well understood that occurrence of coherent pulsations or of type I X-ray bursts is incompatible with the presence of an event horizon, so none of these sources can be found on the list of black hole candidates, even though their masses are unknown.

However, it is true, as pointed out by Narayan & Heyl (nh02 (2002)), that none of the longer (binary) period SXTs, with a measured mass function greater than $`3M_{}`$ is a type I burster. Narayan & Heyl (nh02 (2002)) compute instability of accretion onto a hypothetical $`10M_{}`$ star with a surface of radius between $`(9/8)R_S`$ and $`3R_S`$, and report that for a range of accretion rates compatible with observations of X-ray novae, the star is expected to give rise to an X-ray burst if the accreted column density is $`10^9\mathrm{g}/\mathrm{cm}^2\mathrm{\Sigma }10^{11}\mathrm{g}/\mathrm{cm}^2`$. From this, the authors conclude that black hole candidates cannot have a surface, as they do not exhibit X-ray bursts.

One concern is that the authors do not present the results separately for the lowest column density considered, $`10^9\mathrm{g}/\mathrm{cm}^2`$, and the higher values $`10^{10}\mathrm{g}/\mathrm{cm}^2`$ and $`10^{11}\mathrm{g}/\mathrm{cm}^2`$—for a $`10M_{}`$ star with a $`3R_S`$ radius, the mass transferred in the transient outburst $`6\times 10^{24}\mathrm{g}/\mathrm{cm}^2`$ corresponds to $`6\times 10^9\mathrm{g}/\mathrm{cm}^2`$, so the X-ray burst expected at one of the higher column densities may, in fact, not occur during a SXT outburst. However, there is a more fundamental doubt as to the relevance of the result.

Since the minimum radius of Q-star is 1.4 $`R_S`$ (Miller et al. msn98 (1998)), Narayan & Heyl (nh02 (2002)) consider not only objects composed of matter whose properties have been described by Bahcall et al. (bls90 (1990)) but also more compact configurations whose microscopic properties are not known at all. Therefore there is no reason to assume that the surface of such objects is composed of ordinary matter and is in the temperature range required for X-ray bursts to occur. The stellar surface could be too cold to support a thermonuclear runaway. As a matter of fact, the accreted matter could be converted right away to a more exotic form, as it would be on contact with quark matter in the color-locked phase (Alford, Rajagopal & Wilczek arw (1998), Rapp et al. 1998), or with the skin of a gravastar (Mazur & Mottola MM01 (2001), see below). This could happen even at zero density, contrary to the hypothesis advanced by Narayan & Heyl (nh02 (2002)). No nuclei, no bursts.

## 5 Gravastars

Mazur and Mottola (MM01 (2001)) have recently found a new static, spherically symmetric, solution of Einstein’s field equations. A gravastar, as it is called, has the standard vacuum Schwarzschild exterior, and an interior filled with matter that has the equation of state $`\rho =p`$. The interior is described by the de Sitter solution, and is matched to the exterior vacuum solution in a very thin shell of thickness on the order of the Planck length, $`\lambda _P=1.6\times 10^{33}`$ cm. The gravastar has no horizon or singularity. Its rigid surface is located at a radius just slightly greater than the gravitational radius, $`R_{}=R_S+f\lambda _P`$, $`f2`$.

There are several purely theoretical objections that one could raise against gravastars, none of them conclusive. For example, stellar-mass gravastars have entropy smaller than ordinary stars with the same masses and this would require extremely efficient cooling before gravastars could form during stellar collapse.

There is no observational way to distinguish what may seem to be a Schwarzschild black-hole from a gravastar. To see this, let us denote the surface redshift by

$$\epsilon =\left(1\frac{R_S}{R_{}}\right)^{1/2}=\left(\frac{f\lambda _P}{R_{}}\right)^{1/2}.$$
(2)
For astrophysically interesting gravastars, with mass greater than $`M_{}`$, i.e., $`R_S>3\times 10^5`$ cm, this quantity is very small,

$$\epsilon <10^{19}1.$$
(3)
The power of any radiation emitted by the surface of a gravastar is greatly reduced because only the radiation within the solid angle $`27\epsilon ^2/4`$ around the normal to the surface escapes to infinity. Further, because of gravitational redshift, the power of radiation received by a distant observer is only $`\epsilon ^2`$ of what was emitted at the gravastar’s surface. Therefore, the power emitted from the surface is reduced by

$$\epsilon ^4<10^{75}$$
(4)
by the time it reaches a distant observer. One should conclude that a gravastar with mass greater than $`M_{}`$ is to a distant observer as black as a black hole.

## 6 Conclusions

We have shown that it is fundamentally impossible to give an observational proof for the existence of a black-hole horizon. One could argue that it is not necessary to give such a proof – a black hole is a specific space-time metric, whose properties can, in principle, be determined through observations. If so, no ‘direct’ proof would be necessary, assuming one could determine that the spacetime around a compact object corresponds to the Kerr solution of the Einstein equations.

One way to distinguish a black hole from a rotating star is through the study of orbital and other frequencies (e.g., epicyclic) of accreting matter moving in strong-field gravity (Kato sk01 (2001); Wagoner wso (2001); Abramowicz & Kluźniak ak01 (2001) ak02 (2002)). Another method of determining the space-time geometry is by observation of the energy spectrum reflected from an accretion disc deep in the gravitational well of a compact object (Fabian et al. frsw (1989), firy (2000)). Finally, the capture of stellar-mass compact objects by supermassive black holes in galactic nuclei produces gravitational radiation whose properties reflect the structure of black-hole space-time (Ryan r95 (1995); Hughes sc01 (2001)). Such gravitational radiation could be observed by the gravitational-wave antenna LISA.

The last method is the most powerful since it uncovers the compact body’s multipolar structure. Like the other methods, however, it suffers from a fundamental weakness because it assumes that only a black hole can be the ‘source’ of the ($`a0`$) Kerr metric. Although it was shown that it is very unlikely that other sources exist (Abramowicz, Lasota & Muchotrzeb alm (1976)), such a possibility cannot be excluded.

Nevertheless, the case for the existence of black holes in the Universe is very strong and the evidence very convincing. We think, however, that a shadow of doubt will always cast its pall on our certainty in this matter. But it is a fertile doubt: it has already inspired new ideas and will surely continue to do so.