text:
# Complementary approaches to the ab initio calculation of melting properties

## Abstract

Several research groups have recently reported ab initio calculations of the melting properties of metals based on density functional theory, but there have been unexpectedly large disagreements between results obtained by different approaches. We analyze the relations between the two main approaches, based on calculation of the free energies of solid and liquid and on direct simulation of the two coexisting phases. Although both approaches rely on the use of classical reference systems consisting of parameterized empirical interaction models, we point out that in the free energy approach the final results are independent of the reference system, whereas in the current form of the coexistence approach they depend on it. We present a scheme for correcting the predictions of the coexistence approach for differences between the reference and ab initio systems. To illustrate the practical operation of the scheme, we present calculations of the high-pressure melting properties of iron using the corrected coexistence approach, which agree closely with earlier results from the free energy approach. A quantitative assessment is also given of finite-size errors, which we show can be reduced to a negligible size.

## 1 Introduction

The study of solid-liquid equilibrium by computer simulation has a long history, going back to the classic work of Alder and Wainwright  on the hard-sphere system in the 1950’s. Many techniques have been used to determine the pressure-temperature relation at equilibrium and other melting properties, such as the volume and entropy of fusion. In the last few years, there has been an upsurge of interest in the accurate ab initio treatment of the melting properties of real materials , which has focused attention again on the techniques used to locate the melting transition. Recent ab initio work has been based on two main approaches. The first locates the melting transition by requiring equality of the Gibbs free energies, which are calculated ab initio for liquid and solid ; we call this the ‘free energy’ approach. The second proceeds by fitting a potential model to ab initio calculations and using this to simulate a system containing liquid and solid in coexistence ; we call this the ‘coexistence’ approach. If appropriate measures are taken, the two approaches should clearly give the same results. The purpose of this paper is to analyze the relation between the two approaches and to propose what these ‘appropriate measures’ should be. We illustrate our analysis by presenting ab initio calculations on the high-pressure melting of iron performed using the coexistence approach, which we compare with earlier free-energy results from the same ab initio technique . We shall show that, once the necessary corrections have been applied, the two approaches yield the same results.

The free-energy approach has been well established for many years for calculations based on classical interaction models (see e.g. Refs. ). Typically the procedure for the solid has been to start from the harmonic approximation at low temperatures; the free energy of the high-temperature anharmonic system is then obtained by using the Gibbs-Helmholtz relation to integrate the internal energy given by molecular dynamics (MD) simulation . Alternatively, the anharmonic free energy has sometimes been obtained by thermodynamic integration starting from a reference model such as the Einstein solid . For the liquid, a common procedure has been to obtain the free energy at one thermodynamic state from the work done in reversible expansion  to low density or heating to high temperature. The Gibbs-Helmholtz relation is then used to obtain the free energy at other states. We also mention an important alternative approach, known as ‘Gibbs-Duhem’ integration, which allows the boundary between coexisting phases to be directly mapped out .

In early work, the interaction models were parameterized by fitting to experimental data, but recent years have seen a major shift towards the calculation of melting properties  from ab initio methods based on density-functional theory (DFT) . In DFT, the total energy function of a system is determined by the approximation used for exchange-correlation energy $`E_{\mathrm{xc}}`$. An important ambition then follows: the determination of melting properties with no statistical-mechanical or other approximations except those inherent in $`E_{\mathrm{xc}}`$ itself. This raises major new issues, because it is extremely costly to perform ab initio MD (AIMD) simulations on large systems long enough to reduce statistical errors to an acceptable level. The cost is particularly great for crystals, since extensive electronic $`k`$-point sampling may be needed. Long accepted methods like reversible expansion to low density are impracticable with AIMD. Instead, it is essential to use empirical interaction models that closely mimic the DFT system. The usual techniques are then employed to treat these models, and in some ab initio work thermodynamic integration is used to obtain the free energy difference between the model and DFT systems . In parallel with these developments, the advantages of avoiding intricate free-energy calculations by using the coexistence approach have made this route popular . Again, the simulations must be based on empirical models matched to DFT data.

The present work was stimulated by recent reports of ab initio or nearly ab initio calculations on the melting of a number of materials, including Si , Al , Cu  and Fe . Our own work on Al  and Fe  was based on the free-energy approach, which used an inverse-power reference model, with thermodynamic integration to calculate the difference between the ab initio and reference free energies. A major effort was made to ensure that finite-size and other systematic errors were reduced to an almost negligible level. The other ab initio work on Fe, by two independent groups, used the coexistence approach . The results for the melting curve differed substantially from each other and from our results, and we urgently need to understand the reasons for the disagreements. We aim to shed light on possible reasons here.

We shall explore a number of technical issues. The first concerns the correction of the coexistence method for errors due to the difference between the ab initio system and the empirical model that mimics it. A second issue is the fitting of models to ab initio data. We have already discussed this in depth for the free-energy approach, so here we focus mainly on coexistence. We want to study what physical quantities should be fitted, and how to tell if the fit is good enough. A third important issue concerns finite-size errors, which arise mainly from the limited system sizes that can be handled ab initio. It has been claimed  that the coexistence approach suffers less from size errors than the free-energy approach. We shall demonstrate that in fact this type of error affects both approaches in essentially the same way. One thing we shall not do is to pass judgment on which method is ‘better’, since we shall argue that the answer depends on what one wishes to achieve, and that both approaches are vital. However, we shall comment on the strengths and weaknesses of both.

In the next Section, we define the problem, summarize briefly the free-energy and coexistence approaches as they have been applied in practical ab initio calculations, analyze the corrections that need to be applied in the two approaches, and discuss size errors. Section 3 then presents our new calculations on the high-pressure melting properties of Fe using the coexistence approach. We shall compare with our earlier free-energy results and show the practical necessity of the corrections outlined in Sec. 2. The final Section gives further discussion and summarizes our conclusions.

## 2 Theory of the two approaches

### 2.1 Definition of the problem

We start by defining the problem to be addressed. At a given pressure $`p`$, two phases are in thermodynamic equilibrium when their Gibbs free energies $`G(p,T)`$ are equal. We regard $`G`$ as derived by the relation $`G=F+pV`$ from the Helmholtz free energy $`F(V,T)`$. The ab initio value of $`F`$ at any volume and temperature, denoted by $`F_{\mathrm{AI}}(V,T)`$, is given in classical statistical mechanics by:

$$F_{\mathrm{AI}}(V,T)=k_\mathrm{B}T\mathrm{ln}\left\{\frac{1}{N!\mathrm{\Lambda }^{3N}}𝑑𝐫_1\mathrm{}𝑑𝐫_Ne^{\beta U_{\mathrm{AI}}(𝐫_1,\mathrm{}𝐫_N)}\right\},$$
(1)
where $`U_{\mathrm{AI}}(𝐫_1,\mathrm{}𝐫_N)`$ is the ab initio total energy as a function of the positions $`𝐫_1,\mathrm{}𝐫_N`$ of the $`N`$ nuclei, which range over the volume $`V`$ of the cell containing the system. For the purpose of treating phase equilibria, we need the free energy per atom $`F_{\mathrm{AI}}/N`$ in the thermodynamic limit $`N\mathrm{}`$, $`V\mathrm{}`$ at constant number density $`N/V`$.

Even though the coexistence approach does not work directly with free energies, we shall show later that the errors that need to be overcome in this approach, as well as in the free-energy approach, can be formulated as free-energy errors. But the errors in calculating $`F_{\mathrm{AI}}`$ are of two kinds: electronic-structure errors, i.e. imperfections in the calculation of the ab initio total energy $`U_{\mathrm{AI}}`$ at each set of nuclear positions; and statistical-mechanical errors, i.e. errors in calculating $`F_{\mathrm{AI}}`$ from the given $`U_{\mathrm{AI}}`$, and in taking the thermodynamic limit.

For the purposes of this work, electronic-structure errors are irrelevant: our sole concern is the treatment of phase equilibrium using some given algorithm for computing $`U_{\mathrm{AI}}`$. The ‘free energy’ and ‘coexistence’ routes differ only in the way they address the statistical mechanics. Our problem is therefore to assess and compare the ways that the two approaches control the purely statistical-mechanical errors. When we come to the practical calculations in Sec. 3, we shall compare the results of the two approaches applied using exactly the same algorithm for $`U_{\mathrm{AI}}`$.

### 2.2 The free-energy approach

In the free-energy approach , we use a reference system with total energy function $`U_{\mathrm{ref}}(𝐫_1,\mathrm{}𝐫_N)`$, whose solid and liquid free energies $`F_{\mathrm{ref}}`$ are calculated for very large systems, so that finite-size errors are negligible. Then the only demanding problem is the calculation of the difference $`\mathrm{\Delta }FF_{\mathrm{AI}}F_{\mathrm{ref}}`$, which is accomplished by thermodynamic integration :

$$\mathrm{\Delta }F=_0^1\mathrm{\Delta }U_\lambda 𝑑\lambda ,$$
(2)
where $`\mathrm{\Delta }UU_{\mathrm{AI}}U_{\mathrm{ref}}`$, and the thermal average $`_\lambda `$ is taken in the ensemble generated by the total-energy function $`U_\lambda (1\lambda )U_{\mathrm{ref}}+\lambda U_{\mathrm{AI}}`$. In practice, $`\mathrm{\Delta }U_\lambda `$ is computed as a time average in an AIMD simulation whose dynamics is governed by $`U_\lambda `$. The main kinds of error are: statistical errors in the evaluation of $`\mathrm{\Delta }U_\lambda `$; integration errors due to inadequate numbers of $`\lambda `$ points; and finite-size errors. In the systems studied so far by this approach, all these errors can be brought under tight control.

It is a crucial feature of this approach that $`F_{\mathrm{AI}}`$, calculated as $`F_{\mathrm{ref}}+\mathrm{\Delta }F`$, does not depend on the choice of reference system, provided all technical tolerances are set so as to suppress the errors we have just mentioned. However, for reasons expounded in detail elsewhere , the choice of $`U_{\mathrm{ref}}`$ has a major influence on the overall computational effort, which is minimized by reducing as much as possible the strength of the fluctuations $`\delta \mathrm{\Delta }U^2_\lambda /N`$, where $`\delta \mathrm{\Delta }U\mathrm{\Delta }U\mathrm{\Delta }U_\lambda `$. In particular, minimization of this fluctuation strength is important in ensuring that $`F_{\mathrm{ref}}`$ accounts for almost all the free energy, so that size errors in the small residue $`\mathrm{\Delta }F`$ are negligible. This can be achieved by using a parameterized reference model whose parameters are adjusted so as to minimize $`\delta \mathrm{\Delta }U^2_\lambda /N`$, as in our work on the melting of Fe. In this approach, it is an advantage that different reference systems can be used for solid and liquid, since for many materials it may be difficult to create $`U_{\mathrm{ref}}`$ functions that mimic $`U_{\mathrm{AI}}`$ with high precision in both phases.

### 2.3 The coexistence approach

The coexistence approach is also based on a reference model that mimics the ab initio system. (This model is sometimes called by names such as ‘optimised potential model’ , but here we give it the same name as in the free-energy approach.) Various ways have been used to fit the reference model to the ab initio system. Since it will be relevant later, we note that recent coexistence work on Fe has used the ‘force-matching’ procedure of Ercolessi and Adams , in which the reference parameters are adjusted so that the ab initio atomic forces are reproduced as well as possible for representative sets of atomic positions.

There are also several ways of using the reference model to simulate coexisting phases, and hence to determine the phase boundary between them. In the work of Morris et al. , coexisting solid and liquid Al were simulated with the total number of atoms $`N`$, volume $`V`$ and internal energy $`E`$ fixed. They showed that, provided $`V`$ and $`E`$ are appropriately chosen, the two phases coexist stably over long periods of time, and the average pressure $`p`$ and temperature $`T`$ in the system give a point on the melting curve. The underlying concept is that the mean volume per atom $`\overline{v}V/N`$ is given by:

$$\overline{v}=(1x_l)v_s(p)+x_lv_l(p),$$
(3)
where $`x_l`$ is the fraction of the atoms in the liquid phase, with $`v_s(p)`$ and $`v_l(p)`$ the volumes per atom in the coexisting solid and liquid as a function of pressure. For fixed $`\overline{v}`$, the pressure $`p`$ traverses a certain range as $`x_l`$ goes from 0 to 1. As $`x_l`$ varies in this way, the mean internal energy per atom $`\overline{e}(1x_l)e_s(p)+x_le_l(p)`$ also traverses some range (here, $`e_s(p)`$ and $`e_l(p)`$ are the internal energies per atom in the two phases). Provided $`\overline{e}`$ lies in this range for the given $`\overline{v}`$, the simulation will yield stably coexisting solid and liquid. An alternative procedure would be to simulate at constant $`(N,V,T)`$. Then coexistence will be achieved for a given $`\overline{v}`$ provided $`T`$ is chosen so that the corresponding $`p`$ on the melting line curve falls in the range specified by Eq. (3). Yet another approach was used in the work of Laio et al.  on the high-pressure melting of Fe; this used constant-stress simulations, with enthalpy almost exactly conserved. The approach of Belonoshko et al.  is different again. Here, the $`(N,p,T)`$ ensemble is used. The system initially contains coexisting solid and liquid, but since $`p`$ and $`T`$ generally do not lie on the melting line, the system ultimately becomes entirely solid or liquid. This approach does not directly yield points on the melting curve, but instead provides upper or lower bounds, so that a series of simulations is needed to locate the transition point. Whichever scheme is used, some way is needed of monitoring which phases are present. In the $`(N,V,E)`$ method of Morris et al. , graphical inspection of particle positions appears to have been used, supplemented by calculating of radial distribution functions to confirm the crystal structure of the solid. In the $`(N,p,T)`$ method of Belonoshko et al. , the primary diagnostic is the discontinuity of volume as the system transforms from solid to liquid.

The coexistence calculations presented here on the high-pressure melting of Fe employ the $`(N,V,E)`$ method applied to reference systems consisting of the ‘embedded-atom model’  fitted to ab initio data. For geophysical reasons, we are interested in pressures near that at the boundary between the Earth’s inner and outer cores, namely 330 GPa . In this pressure region, the most stable crystal structure just below the melting curve is believed to be hexagonal close packed (h.c.p.) , and we assume here that melting occurs from this phase. We start with a simulation cell containing only the h.c.p. solid, with the basal plane parallel to one face of the cell, and the system is allowed to thermalize at a temperature where coexistence is expected. At some instant of time, the simulation is stopped, and the atoms in one half of the system are held fixed; the boundary between the two halves is taken parallel to the h.c.p. basal plane. The atoms in the other half are raised to a very high temperature, and dynamical evolution of these atoms is allowed to proceed so that melting occurs. With this half of the system molten, its temperature is now reduced to the original value, the atoms in the other half still being fixed. Finally, the atoms in the fixed half are given thermal velocities and released, and the whole system is allowed to evolve freely. The system is monitored by calculation of the average number density in slices of the cell taken parallel to the boundary between solid and liquid. As we shall show, the density in the solid part is a periodic function of slice number, while in the liquid it fluctuates rather weakly about its average value. The total energy, temperature and pressure are, of course, monitored throughout the simulation. A feature of the $`(N,V,E)`$ method is that at equilibrium the stress in the solid phase will generally not be hydrostatic, and manual adjustment of the cell parameters is needed to achieve hydrostatic stress. We shall show in Sec. 3 that this does not present a problem.

### 2.4 Correcting the coexistence approach

Both the free-energy and the coexistence approaches are subject to errors, which need to be assessed and corrected for. It may be too costly to perform all the ab initio calculations with the precision needed to obtain accurate melting properties, so that there are errors due to inadequate $`k`$-point sampling, incompleteness of the basis set, or other approximations. In the coexistence approach, inevitable differences between the reference and ab initio total energies will create errors in the melting properties, which need to be evaluated. The correction of errors in the free-energy approach has been extensively discussed elsewhere , so our main concern here is correction of the coexistence approach.

Assuming the ab initio energy is calculated with adequate precision, the main error comes from differences between $`U_{\mathrm{AI}}`$ and $`U_{\mathrm{ref}}`$. The key question posed here is therefore: how are predicted melting properties changed by small changes in the total energy function? (System-size errors will be discussed separately later.) Related questions have been discussed before (see e.g. Ref. ), so we give only a brief summary of the main points.

The difference $`U_{\mathrm{AI}}U_{\mathrm{ref}}`$ is denoted by $`\mathrm{\Delta }U`$, as in the free-energy approach. The change of any quantity resulting from the replacement of $`U_{\mathrm{ref}}`$ by $`U_{\mathrm{AI}}`$ will also be indicated by $`\mathrm{\Delta }`$; for example, the shift of melting temperature at a given pressure is called $`\mathrm{\Delta }T_\mathrm{m}`$. The melting temperature is shifted because the liquid and solid Gibbs free energies $`G^l(p,T)`$ and $`G^s(p,T)`$, and hence their difference $`G^{ls}(p,T)G^l(p,T)G^s(p,T)`$, are shifted. Working at the given pressure, we take the variable $`p`$ as read and express the ab initio value of $`G^{ls}`$ as:

$$G_{\mathrm{AI}}^{ls}(T)=G_{\mathrm{ref}}^{ls}+\zeta \mathrm{\Delta }G^{ls}(T),$$
(4)
where the parameter $`\zeta `$ is introduced so that the reference melting temperature $`T_\mathrm{m}^{\mathrm{ref}}`$ can be written as a power series:

$$T_\mathrm{m}^{\mathrm{AI}}=T_\mathrm{m}^{\mathrm{ref}}+\zeta T_\mathrm{m}^{}+\zeta ^2T_\mathrm{m}^{\prime \prime }+\mathrm{}.$$
(5)
Since the Gibbs free energies are equal in the two phases, this $`T_\mathrm{m}^{\mathrm{AI}}`$ is the solution of $`G_{\mathrm{AI}}^{ls}(T)=0`$, which is:

$$G_{\mathrm{ref}}^{ls}\left(T_\mathrm{m}^{\mathrm{ref}}+\zeta T_\mathrm{m}^{}+\zeta ^2T_\mathrm{m}^{\prime \prime }\mathrm{}\right)+\zeta \mathrm{\Delta }G^{ls}\left(T_\mathrm{m}^{\mathrm{ref}}+\zeta T_\mathrm{m}^{}+\zeta ^2T_\mathrm{m}^{\prime \prime }+\mathrm{}\right)=0.$$
(6)
Expanding in powers of $`\zeta `$ and equating powers, one obtains formulas for $`T_\mathrm{m}^{}`$, $`T_\mathrm{m}^{\prime \prime }`$, etc:

$`T_\mathrm{m}^{}`$ $`=`$ $`\mathrm{\Delta }G^{ls}\left(T_\mathrm{m}^{\mathrm{ref}}\right)/S_{\mathrm{ref}}^{ls}`$
$`T_\mathrm{m}^{\prime \prime }/T_\mathrm{m}^{}`$ $`=`$ $`{\displaystyle \frac{1}{S_{\mathrm{ref}}^{ls}}}\left[T_\mathrm{m}^{}C_{p,\mathrm{ref}}^{ls}/2T_\mathrm{m}^{\mathrm{ref}}+\mathrm{\Delta }S^{ls}\right],`$ (7)
where $`S_{\mathrm{ref}}^{ls}S_{\mathrm{ref}}^lS_{\mathrm{ref}}^s`$ is the reference entropy of fusion, $`C_{p,\mathrm{ref}}^{ls}C_{p,\mathrm{ref}}^lC_{p,\mathrm{ref}}^s`$ is the liquid-solid difference of the constant-pressure specific heats, and $`\mathrm{\Delta }S^{ls}`$ is the shift of the entropy of fusion. We note particularly the implication of the formula for $`T_\mathrm{m}^{}`$. Since entropies of fusion are on the order of $`k_\mathrm{B}`$ per atom, then a difference $`\mathrm{\Delta }G^{ls}`$ of 10 meV/atom implies a shift of melting temperature of ca. 100 K, so that substantial errors will need to be corrected for unless the reference total energy function matches the ab initio one very precisely. Although we have included a formula for $`T_\mathrm{m}^{\prime \prime }`$, it seems unlikely that this will be used in practice, except perhaps for a rough estimate of $`T_\mathrm{m}^{\prime \prime }/T_\mathrm{m}^{}`$, since $`\mathrm{\Delta }S^{ls}`$ would be difficult to compute without rather extensive free-energy calculations.

If simulations are done in the isothermal-isobaric ensemble, then for closely matching $`U_{\mathrm{AI}}`$ and $`U_{\mathrm{ref}}`$, with small fluctuations of $`\mathrm{\Delta }U`$, the Gibbs free energy shifts $`\mathrm{\Delta }G^l`$ and $`\mathrm{\Delta }G^s`$ can be evaluated using the well-known expansion:

$$\mathrm{\Delta }G=\mathrm{\Delta }U_{\mathrm{ref}}\frac{1}{2}\beta \delta \mathrm{\Delta }U^2_{\mathrm{ref}}+\mathrm{},$$
(8)
where $`\delta \mathrm{\Delta }U\mathrm{\Delta }U\mathrm{\Delta }U_{\mathrm{ref}}`$, and the averages are taken in the reference ensemble. The ab initio simulations presented later were performed in the isothermal-isochoric ensemble, so that the quantity that is readily evaluated is $`\mathrm{\Delta }F(V,T)`$, the change of Helmholtz free energy when $`U_{\mathrm{ref}}`$ is replaced by $`U_{\mathrm{AI}}`$ are constant volume $`V`$. This $`\mathrm{\Delta }F`$ is given by the same formula Eqn (8), but with the averages evaluated in the isothermal-isochoric ensemble. In this case, we need to consider the relation between $`\mathrm{\Delta }G`$ and $`\mathrm{\Delta }F`$, which is readily shown to be:

$$\mathrm{\Delta }G=\mathrm{\Delta }FV\kappa _T\mathrm{\Delta }p^2,$$
(9)
where $`\kappa _T`$ is the isothermal compressibility and $`\mathrm{\Delta }p`$ is the change of pressure when $`U_{\mathrm{ref}}`$ is replaced by $`U_{\mathrm{AI}}`$ at constant $`V`$ and $`T`$.

### 2.5 Size effects

Whichever approach is used, coexistence must be treated in the thermodynamic limit. In the free-energy approach, almost all the free energy is that of the reference system, for which size errors can be made negligible by doing simulations on very large systems. Appreciable size errors remain only in the difference $`F_{\mathrm{AI}}F_{\mathrm{ref}}`$ between the ab initio and reference systems. These can only be assessed and corrected for by explicitly calculating $`\mathrm{\Delta }F\mathrm{\Delta }FF_{\mathrm{AI}}F_{\mathrm{ref}}`$ for systems of increasing size, as we report elsewhere for the cases of Al  and Fe . Similarly, in the coexistence approach, explicit coexistence simulations are performed with the model system on very large systems, so that size errors are made negligible. However, the shift of melting temperature due to the difference $`U_{\mathrm{AI}}U_{\mathrm{ref}}`$ requires calculations of $`\mathrm{\Delta }F`$ of exactly the same kind as are needed in the free-energy approach. This means that the size errors are essentially the same in the two approaches. We give a quantitative assessment of size errors in the melting properties of Fe by the coexistence approach in the next Section.

## 3 The melting properties of iron

### 3.1 Technical details

Our coexistence calculations on high-pressure Fe use precisely the same DFT techniques used in our free-energy work , so we give only a brief summary here. The exchange-correlation functional $`E_{\mathrm{xc}}`$ is the generalized gradient approximation known as Perdew-Wang 1991 . We use the projector-augmented-wave (PAW) implementation of DFT , a technique that shares the properties both of all-electron methods such as full-potential linearized augmented plane waves (FLAPW)  and the ultrasoft pseudopotential method . The calculations were done using the VASP code . Details of the core radii, augmentation-charge cut-offs, etc. are exactly as in our PAW work on liquid Fe . Our division into valence and core states is also the same: the 3$`p`$ electrons are treated as core states, but their response to the high compression is represented by an effective pair potential, with the latter constructed using PAW calculations in which the 3$`p`$ states are explicitly included as valence states.

The reference model for our coexistence simulations is the embedded-atom model (EAM) recently used by Belonoshko et al.  to calculate the high-pressure melting curve of Fe. This EAM has the standard form , in which the total energy $`E_{\mathrm{tot}}=_iE_i`$ is a sum of energies $`E_i`$ of atoms $`i`$, with each $`E_i=E_i^{\mathrm{rep}}+F(\rho _i)`$ consisting of two parts: first, a purely repulsive energy $`E_i^{\mathrm{rep}}`$ represented as a sum of inverse-power pair potentials $`E_i^{\mathrm{rep}}=_j^{}ϵ(a/r_{ij})^n`$, where $`r_{ij}`$ is the interatomic separation and the sum excludes $`i=j`$; second, an ‘embedding’ part $`F(\rho _i)`$ which accounts for the metallic bonding mainly due to partial filling of the $`d`$-bands. The embedding function $`F(\rho )`$ is represented as $`ϵC\rho ^{1/2}`$, and the density $`\rho _i`$ for atom $`i`$ is given by the sum over neighbours $`\rho _i=_j^{}(a/r_{ij})^m`$. The parameters in this EAM were determined in Ref.  by fitting to full-potential linearized muffin-tin orbital (FPLMTO) energies for typical configurations of liquid iron, and they are: $`n=8.137,m=4.788,ϵ=0.0173`$ eV, $`a=3.4714`$ Å, and $`C=24.939`$.

The aim of our calculations is not to generate the entire melting curve, but to obtain a point on this curve at a pressure close to the value of 330 GPa at the boundary between the Earth’s solid inner core and liquid outer core . Our main simulations have been done on cells containing 8000 atoms, constructed as a $`10\times 10\times 40`$ h.c.p. supercell. The long axis is perpendicular to the basal plane of the crystal, and we refer to it as the $`z`$-axis. We have tested the adequacy of this system size by performing coexistence calculations on systems of up to 20,480 atoms. Within the statistical errors, we were unable to detect any difference between results for coexistence pressure and temperature with this system size and those with the 8,000-atom system. The mean volume per atom for the calculations reported here was $`v=7.12`$ Å<sup>3</sup>, which is near the volumes of the solid and liquid in the pressure region of interest. Since our calculations are performed at constant volume and constant cell shape (see Sec. 2.3), we checked that non-hydrostatic stresses do not affect the melting properties of the system. To do this, we performed the simulations with the two $`c/a`$ values 1.64 and 1.66, which give the $`P_{zz}`$ diagonal component of the stress tensor slightly larger or smaller than the $`P_{xx}`$ and $`P_{yy}`$ components in the two cases. We find that the effect on the melting temperature is undetectably small.

In preparing the coexisting solid and liquid as described in Sec. 2.3, the system was initially equilibrated at 6100 K, and the high temperature used to ensure complete melting of half the system was $`5\times 10^4`$ K. Once melting had been achieved, the liquid part was re-equilibrated at 6100 K before free evolution of the whole system was started.

### 3.2 Results

In Fig. 1, we show the temperature of the system as a function of simulation time for $`c/a=1.64`$. We also show the three diagonal components of the stress tensor; the off-diagonal components fluctuate around their average value of zero, so there is no shear stress on the cell. After an initial equilibration period, one sees that the temperature and the pressure settle around the values $`T=6550\pm 100`$ K and $`p=305\pm 1`$ GPa. In Fig. 2, we display the density profile, calculated by dividing the simulation cell into 400 slices parallel to the liquid-solid interface and counting the number of atoms present in each slice. The profile shown corresponds to the last configuration in the simulation with $`c/a=1.64`$; a similar profile is observed in the simulation with $`c/a=1.66`$. It is easy to identify the solid and the liquid regions in the system: in the solid region the density is a periodic function of slice number, but in the liquid it fluctuates randomly around its average value. This confirms that we do indeed have coexisting solid and liquid in the cell. At 305 GPa, the melting temperature reported by Belonoshko with precisely the same EAM is 6680 K, which agrees with our value within the combined statistical errors.

We now correct the EAM reference melting temperature to obtain the fully ab initio melting temperature using the techniques presented in Sec. 2.4. In order to assess possible system-size errors, the free-energy corrections of Eqs. (8,9) were calculated with systems of 64, 150 and 288 atoms. To do this, we used the EAM reference model to generate long simulations for the solid and liquid separately. From these simulations we extract typically 50 and 100 statistically independent configurations for solid and liquid respectively, for which we calculate the DFT total energies $`U_{\mathrm{AI}}`$. The differences $`\mathrm{\Delta }UU_{\mathrm{AI}}U_{\mathrm{ref}}`$ are then used to compute the free-energy corrections and hence the shift of melting temperature (see Eq. (7)). The temperature was set equal to the value of 6550 K that emerges from the coexistence simulations, and volumes per atom of $`v_s=7.05`$ and $`v_l=7.218`$ Å<sup>3</sup>/atom were also deduced from the coexistence simulations; the pressure in both cases was 303 GPa. (The very small difference from the pressure of 305 GPa mentioned above is of no consequence.) From these two volumes we can extract the volume change on melting, which is $`\mathrm{\Delta }V2.5\%`$, the same value reported as by Belonoshko et al. . (The value from our ab initio free-energy work was 1.8 % at 305 GPa .) We also calculated the entropy of melting from the relation $`TS_{\mathrm{ref}}^{ls}=E_{\mathrm{ref}}^{ls}+pV_{\mathrm{ref}}^{ls}`$; the energy change on melting $`E_{\mathrm{ref}}^{ls}`$ was calculated using the two separate simulations of solid and liquid at $`p=303`$ GPa and $`T=6550`$K, and we obtained $`S_{\mathrm{ref}}^{ls}=0.88`$ k<sub>B</sub>/atom; a similar value of $`S^{ls}`$ can be deduced from the work of Belonoshko et al.  by using their value for $`V^{ls}`$ and the Clausius-Clapeyron relation on the reported melting curve. (The value from our ab initio free-energy work was 1.07 $`k_\mathrm{B}`$/atom at 305 GPa .) The DFT energies were carefully checked for electronic $`k`$-point errors, as in our free-energy work . As expected from that work, $`\mathrm{\Gamma }`$-point sampling for the 64-atom system underestimates the fully converged energies of liquid and solid by ca. 10 and 50 meV/atom respectively. With 150 atoms, $`\mathrm{\Gamma }`$-point sampling gives negligible errors for the liquid, but an error of ca. $`8`$ meV/atom for the solid. With 288 atoms, we have used only $`\mathrm{\Gamma }`$-point sampling, but the indications are that $`k`$-point errors should be negligible for both phases.

We report in Table 1 the results for $`\mathrm{\Delta }U_{\mathrm{ref}}`$ and $`(\delta \mathrm{\Delta }U)^2_{\mathrm{ref}}`$ for the three system sizes, the results being already corrected for electronic $`k`$-point errors. An important feature of these results is that there is no discernible system size effect on these corrections within the statistical errors of 5 meV or less, so that systems of 64 atoms are adequate for calculating the corrections. From Eq. (7), this imples that the shift of $`T_\mathrm{m}`$ due to size errors will not more than $`50`$ K. In order to obtain the corrections $`\mathrm{\Delta }G`$ to the Gibbs free energy, we need to include the term in $`\mathrm{\Delta }p^2`$ in Eq. (9). We find that the pressure differences between the EAM and ab-initio systems are only 22 and 12 GPa for solid and liquid respectively, which give this term values of ca. 5 meV. From Table 1, we see that the free-energy differences between the ab initio and EAM systems have the effect of stabilizing the liquid with respect to the solid by about 35 meV/atom compared with the EAM. In order to obtain the correction to the melting temperature from Eq. (7), we use the value for the entropy of melting $`S_{\mathrm{ref}}^{ls}`$ of $`0.88k_\mathrm{B}`$/atom quoted above. This gives the first-order correction $`T_\mathrm{m}^{}=450`$ K, so that our corrected $`T_m`$ at 305 GPa is $`6100\pm 100`$ K. This result is in very close agreement with the free-energy approach, which at $`p=305`$ GPa gives $`T_m=6100`$ K. (This value is somewhat lower than the preliminary result of our free-energy calculations ; the downward revision came from a careful reanalysis of the anharmonic free energy of the h.c.p. crystal, as reported in Refs. .)

Although the EAM of Belonoshko et al.  mimics our ab initio systems reasonably well, we have found that the model can be still further improved by refitting it to our ab initio energies of the solid and liquid. In doing this, we allowed only the repulsive potential and the strength of the embedding energy to change, so that only the parameters $`n`$, $`ϵ`$ and $`C`$ are allowed to vary. These were adjusted to minimise the $`\delta \mathrm{\Delta }U`$ fluctuations for the solid, while also maintaining the correct pressure. The new parameters are $`n=5.93,ϵ=0.1662`$ eV, $`C=16.55`$, the parameters $`a=3.4714`$ Å and $`m=4.788`$ retaining their original values. We then repeated the coexistence simulation using this new EAM and obtained a coexistence temperature of $`T_m=6200\pm 100`$ K at a pressure of $`p=323.5\pm 1`$ GPa. The very small free-energy corrections (Table 2) slightly stabilize the solid, resulting in an increase of melting temperature of ca. 50 K to give $`6250\pm 100`$ K. The differences with the ab-initio pressures are negligible in this case, and the term in $`\mathrm{\Delta }p^2`$ in Eq. (9) does not contribute to the free energy difference. At this pressure, the melting temperature from our free energy approach was $`T_m=6290`$ K , so that once again we find very good agreement between the two approaches.

## 4 Discussion and conclusions

We have advocated the aim of calculating the melting properties that follow from a chosen approximation for the exchange-correlation functional $`E_{\mathrm{xc}}`$, with the errors due to all other approximations being made negligible. We have stressed that in the free-energy approach the total energy of the reference model does not have to agree exactly with the ab initio total energy, since this approach includes the calculation of the free energy difference between the two systems. Up to now, no allowance has been made for such differences in the coexistence approach, but we have shown how corrections can be made so that the final results are again independent of the reference model. Our practical results for high-pressure Fe demonstrate that the two approaches then give the same melting temperatures, as would be expected.

Our analysis has implications for the way in which reference systems are constructed. In the free-energy approach, the crucial requirement of a reference model is that the fluctuations of the difference $`\mathrm{\Delta }U`$ between the ab initio and reference energies be as small as possible; we have the freedom to use different reference models for the solid and liquid phases, and a constant offset between the ab initio and reference energies is of no consequence. In the coexistence approach, it seems clear that the model should be chosen so that the $`\mathrm{\Delta }G`$ corrections, and the resulting shift of melting temperature away from that of the reference system, be as small as possible. This demands more of the reference system than in the free-energy approach. It is necessary but not sufficient that the $`\mathrm{\Delta }U`$ fluctuations be small, in order that $`\mathrm{\Delta }G`$ be small (see Eq. (9)). In addition, a constant energy offset will shift $`T_\mathrm{m}`$, if the offset is not identical in the two phases, since the same model must simultaneously reproduce the energetics of both. This implies that force matching  may not always be a reliable way of constructing reference models: even if the reference and ab initio forces match precisely in both phases, there may still be a non-cancelling $`\mathrm{\Delta }U`$ offset that will shift $`T_\mathrm{m}`$. In this sense, energy matching is safer than force matching.

Concerning finite-size errors, we have shown that these affect appreciably only the free-energy differences between the ab initio and reference systems, since in both approaches the calculations on the reference model will always be done on systems large enough to make size errors negligible. If anything, these errors may be more troublesome in the coexistence approach, because of the need to make the reference model fit both phases at once, so that the free-energy corrections may sometimes be larger. However, in our coexistence calculations on high-pressure Fe, we have shown that these errors are unlikely to shift $`T_\mathrm{m}`$ by more than $``$ 50 K.

In commenting on the strengths and weaknesses of the two approaches, we note that they differ mainly in the way of treating the thermodynamic properties of the reference model. In both approaches, corrections are then needed to obtain the properties of the ab initio system, and these involve ab initio simulations to perform thermodynamic integration or to compute the free-energy differences of Eq. (8). The melting properties of the reference model are probably more straightforward to calculate by the coexistence approach, since the free-energy approach requires rather intricate thermodynamic integrations. However, the heaviest computational effort comes in the calculation of the free-energy corrections, and here we believe the free-energy approach may have the advantage, since the effort can be reduced by using different reference systems for the two phases. In practice, we have found it very helpful to use both approaches, as in the work on Fe reported here.

## Acknowledgments

The work of DA is supported by a Royal Society University Research Fellowship, and the work of MJG by Daresbury Laboratory and GEC. The computational facilities of the UCL HiPerSPACE Centre (JREI grant JR98UCGI) and the Manchester CSAR Service (grant GST/02/1002 to the Minerals Physics Consortium) were used. The authors thank A. Laio and S. Scandolo for detailed discussions about their coexistence calculations.


format: text

text:
# The X-ray – radio alignment in the 𝑧=2.2 radio galaxy PKS 1138–262

## 1 Introduction

Radio galaxies are the most massive galaxies known at high redshift, and have been shown to reside in regions of high ambient gas and galaxy density (Röttgering, Best, & Lehnert 1999; McCarthy 1993). As such, these sources provide unique laboratories for the study of the formation of galaxies and clusters of galaxies, acting as beacons to biased galaxy formation within large scale structure at high $`z`$ (Carilli et al. 2001).

One of the more extensively studied high redshift radio galaxies in a dense environment is the source PKS 1138–262 at $`z=2.156`$ (Pentericci et al. 1997). The 1138–262 radio source has the most disturbed morphology of any $`z>2`$ radio galaxy yet identified, appearing as a string of bright knots with a total linear extent of 15′′ (Carilli et al. 1997). The radio source is enveloped in a Ly$`\alpha `$ emitting halo with a luminosity of $`2.5\times 10^{44}`$ erg s-1, spatially extended on a scale of 20′′ and with a major axis oriented along the radio source axis (Pentericci et al. 1997). There are clear spectroscopic signatures of strong interaction between the expanding radio source and the ambient gas (Kurk et al. 2001). The polarized radio emission from 1138–262 shows the most extreme values of Faraday rotation of any source at $`z>2`$, with (rest frame) rotation measures (RMs) up to 6250 rad m-1 (Carilli et al. 1997). These RMs are comparable to values observed toward low $`z`$ powerful radio galaxies in dense, X-ray emitting cluster atmospheres. For the low $`z`$ sources the origin of the large RMs is thought to be magnetized cluster gas, with magnetic fields of a few $`\mu `$G or more (Carilli & Taylor 2001).

High resolution optical images show that the parent galaxy of 1138–262 is comprised of many compact ‘knots’ of emission distributed on a scale of 10′′. Wide field imaging and spectroscopy with the VLT has revealed many compact Ly$`\alpha `$ emitting galaxies distributed on arcminute scales around 1138–262 (Pentericci et al. 2000), with an implied over-density of objects similar to $`z3`$ proto-clusters recently discovered in searches for UV dropout galaxies by Steidel et al. (1998). Pentericci et al. (2000) interpret these results on 1138–262 as indicating a forming giant elliptical galaxy within a proto-cluster.

Carilli et al. (1998) searched for X-ray emission from 1138–262 using the HRI detector on ROSAT. While the source was clearly detected, the sensitivity and spatial resolution of the ROSAT observations were insufficient to determine the relative contribution to the X-ray luminosity from the AGN and possible extended emission. In this paper we present Chandra observations which clarify the X-ray emission processes in 1138–262. While the dominant X-ray source is the AGN, we also detect an extended soft X-ray emission component with an extent and orientation similar to the radio source. We consider mechanisms for this extended emission, and discuss the origin of the X-ray – radio alignment in 1138–262. We use Ho = 65 km s-1 Mpc-1, $`\mathrm{\Omega }_M=0.3`$, and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, leading to DL = 18300 Mpc and Dθ = 1844 Mpc (1′′ = 9 kpc). We defined spectral index, $`\alpha `$, as a function of frequency, $`\nu `$, and flux density, Sν, as: $`\mathrm{S}_\nu \nu ^\alpha `$.

## 2 Observations and Results

PKS 1138–262 was observed for 40 ksec on June 6, 2000, with the back-illumiated ACIS-S CCD detector array on the Chandra Observatory (Weisskopf, O’Dell, & Van Speybroeck, 1996; Garmire 1997). The source was positioned one arcmin from the standard aim point so as to avoid the node boundary on the ACIS-S3 chip. Standard ACIS settings were used: TE mode with 3.2 s readout and ’faint’ telemetry format. The data were reprocessed in 2001Jan. For analysis of extended regions the event file was filtered to reject times of high background, resulting in a a live-time of 28 ksec. For unresolved sources we use the entire observation (live-time of 39.5 ksec). Data analysis was performed using the CIAO 2.1 software package for spectroscopic imaging. Image analysis was also performed using the Astronomical Image Processing System.

### 2.1 Images

Figure 1 shows the Chandra total X-ray image (energy range = 0.3 to 10 keV), convolved with a Gaussian of Full Width at Half Maximum (FWHM) of 2′′, along with the VLA radio image at 1.4 GHz, 2′′ resolution in greyscale (from Carilli et al. 1997). A total of 880 counts are detected in a $`17^{\prime \prime }\times 22^{\prime \prime }`$ box centered on the radio galaxy. We estimate a background contribution in this box of 102 counts.

Two components are evident in Figure 1: (i) compact emission at the location of the AGN, and (ii) extended emission on a scale of about 20′′. Comparison of the compact AGN component to the Chandra point spread function, and Gaussian fitting, shows that the emission is unresolved with a size $`<1^{\prime \prime }`$. The total emission from the nucleus is 663 counts. There is also a point source about 5′′ northwest of the nucleus which contains 15 counts at a position of 11h 40m 47.93s, –26o 29 06.5′′ (J2000). This source is also detected in the Ly$`\alpha `$ and the HST images, but not in the 5 GHz image to a 3$`\sigma `$ limit of 90$`\mu `$Jy (Figure 2a). Subtracting the background and the two point sources leads to 98$`\pm `$30 counts for the extended emission from 1138–262, or 13$`\%`$ of the total.

Figure 2a shows the image of total X-ray emission from 1138–262 convolved with a Gaussian of FWHM = 1′′, along with the radio image at 5GHz, 0.5′′ resolution (2a), the HST F606W image (2b), and the VLT Ly$`\alpha `$ image (2c). The position of the nucleus of 1138–262 as determined in near IR images, H$`\alpha `$ line images, and the HST image, is 11h 40m 48.37s, –26o 29 08.9′′ $`\pm `$ 0.1′′ (J2000) using USNO astrometry (Monet et al. 1996). The radio nucleus is located at 11h 40m 48.35s, –26o 29 08.8′′ $`\pm `$0.1′′, implying that the radio and optically determined positions for the nucleus are within 0.2′′ of each other. The nominal astrometry from Chandra places the X-ray nucleus at 11h 40m 48.31s, –26o 29 08.6′′, or $`0.5\pm 0.1^{\prime \prime }`$ and $`0.2\pm 0.1^{\prime \prime }`$ away from the radio position in RA and Dec, respectively. This is consistent with the offset seen for the $`z=2.183`$ QSO in the 1138–262 field (see below), and with the mean offsets of $`0.4^{\prime \prime }`$ and $`0.0^{\prime \prime }`$ in RA and Dec, respectively, between radio and X-ray positions of sources in the Abell 370 field found by Barger et al. (2001). In the analysis below, the images at different frequencies are aligned using the nucleus.

Figure 3 shows the X-ray emission in three energy bands: 0.3 to 1.2 keV (3a = soft), 1.2 to 2.4 keV (3b = mid), and $`>2.4`$ keV (3c = hard). The 5 GHz radio continuum emission is shown in greyscale in Figure 3a. The extended X-ray emission is softer than the emission from the nucleus, while the emission from the point source 5′′ northwest of the nucleus is harder. For the soft image the extended emission contributes $`21\pm 3.5\%`$ to the total source counts. This percentage drops to $`7\pm 1\%`$ for the mid energy range, and $`<11\%`$ in the hard X-ray image.

Figure 4 shows the surface brightness profiles averaged azimuthally in rings centered on the nucleus of 1138–262 for quadrants oriented parallel and perpendicular to the radio axis. The profile along the radio axis shows extended emission to about 12′′ radius. For illustrative purposes we have fit a standard cluster $`\beta `$ model to the on-axis profile, as shown with the solid line. The formal results for the model parameters are $`\beta =2.5\pm 1`$ and $`\mathrm{r}_\mathrm{c}=12\pm 4^{\prime \prime }`$. This value of $`\beta `$ is much steeper than what is expected for a typical cluster atmosphere, even for high redshift clusters, for which $`\beta \frac{2}{3}`$ is observed typically (Stanford et al. 2001). Such a steep $`\beta `$ may simply reflect the fact that the extended emission region is clearly not circularly symmetric, but is highly elongated. An elongated morphology will lead to a decrease in the fractional emitting area intercepted with increasing annular radius, and hence an artificially steep radial profile. Perhaps a more interesting result in Figure 4 is that the profile perpendicular to the radio source axis shows very little emission beyond about 4′′ distance from the nucleus. Averaging the signal out to 12′′ radius, excluding the inner 5′′, shows that the luminosity of the extended emission perpendicular to the radio source major axis is less than 25$`\%`$ of that seen along this axis.

We have searched for X-ray emission from the Ly$`\alpha `$ emitting galaxies in the (proto) cluster around 1138–262 (Pentericci et al. 2000). There are no detected X-ray sources within 1′′ of any of the cluster galaxies, outside of the Ly$`\alpha `$ halo of 1138–262 itself. The 3$`\sigma `$ upper limit to the 2 to 10 keV X-ray luminosities of these sources is $`1.2\times 10^{43}`$ erg s-1, assuming a standard powerlaw spectrum of spectral index $`\alpha =0.8`$. X-rays are detected from a $`z=2.183`$ QSO located 2 west of 1138–262 (Pentericci et al. 2000), at 11h 40m 39.73s –26o 28 44.9′′ (J2000). A total of 56 counts are detected from this source, implying an X-ray luminosity of $`1.8\times 10^{44}`$ erg s-1 for $`\alpha =0.8`$.

### 2.2 The X-ray Spectra

The X-ray spectrum of the nuclear emission from 1138–262 is shown in Figure 5. The spectrum is reasonably fit by a power-law with $`\alpha =0.8\pm 0.2`$, plus absorption at the redshift of the source of N(HI) = $`2.6\pm 0.6\times 10^{22}`$ cm-2. The reduced $`\chi ^2`$ for the fit is 0.52, and the error bars for the fit parameters represent changes in reduced $`\chi ^2`$ by $`\pm 1`$. The Galactic value of N(HI) in this direction is $`4.5\times 10^{20}`$ cm-2 (Stark et al. 1992). The unobscured 2 to 10 keV luminosity for the nucleus is $`\mathrm{L}_\mathrm{X}^{\mathrm{nuc}}=4.0\times 10^{45}`$ erg s-1. The implied extinction toward the nucleus is $`\mathrm{A}_\mathrm{V}=16`$, assuming a Galactic dust-to-gas ratio. This large extinction may explain the lack of strong Ly$`\alpha `$ emission from the nucleus of 1138–262. For completeness, we note that the nuclear spectrum can be equally well fit by a thermal spectral model (Raymond-Smith) with $`\mathrm{kT}=16_4^{+7}`$ keV.

Considering the radio properties of 1138–262, the luminosity between rest frame frequencies of 1 to 10 GHz for the nucleus of 1138–262 is: $`\mathrm{L}_{110\mathrm{G}\mathrm{H}\mathrm{z}}^{\mathrm{nuc}}=1.5\times 10^{43}`$ erg s-1. The radio luminosity for the entire source for rest frame frequencues of 0.1 to 1.0 GHz is: $`\mathrm{L}_{0.11\mathrm{G}\mathrm{H}\mathrm{z}}^{\mathrm{tot}}=1.8\times 10^{45}`$ erg s-1 (Carilli et al. 1997).

The archetype powerful radio galaxy Cygnus A has a highly absorbed X-ray emitting nucleus with N(HI) = $`3.8\times 10^{23}`$ cm-2 and an unobscured 2 to 10 keV luminosity of $`\mathrm{L}_\mathrm{X}^{\mathrm{nuc}}=1\times 10^{45}`$ erg s-1 (Ueno et al. 1994). The radio nucleus of Cygnus A has $`\mathrm{L}_{110\mathrm{G}\mathrm{H}\mathrm{z}}^{\mathrm{nuc}}=1.3\times 10^{42}`$ erg s-1, and the radio luminosity for the entire source from 0.1 to 1.0 GHz is $`\mathrm{L}_{0.11\mathrm{G}\mathrm{H}\mathrm{z}}^{\mathrm{tot}}=3.8\times 10^{44}`$ erg s-1 (Carilli et al. 1991). Hence, to within a factor two for both 1138–262 and Cygnus A: $`\frac{\mathrm{L}_\mathrm{X}^{\mathrm{nuc}}}{\mathrm{L}_{110\mathrm{G}\mathrm{H}\mathrm{z}}^{\mathrm{nuc}}}500`$, and $`\frac{\mathrm{L}_\mathrm{X}^{\mathrm{nuc}}}{\mathrm{L}_{0.11\mathrm{G}\mathrm{H}\mathrm{z}}^{\mathrm{tot}}}2.5`$. Ward et al. (1996) show that the obscuration corrected properties of the X-ray nucleus of Cygnus A are typical of broad line radio galaxies.

Due to the paucity of counts, the spectrum of the extended X-ray emission from 1138–262 is very poorly constrained, with an allowed range 0.8 to 15 keV for a Raymond-Smith thermal plasma emission model. In the analysis below we adopt a value of 5 keV. The implied 2 to 10 keV (rest frame) luminosity for the extended emission is then $`3\times 10^{44}`$ erg s-1. Further observations of 1138–262 are required to constrain this important parameter.

## 3 Discussion

The most interesting aspect of the extended X-ray emission from 1138–262 is that it is oriented along the radio axis and has a similar spatial extent. There are a number of possible mechanisms to explain this emission, and its alignment with the radio axis.

### 3.1 Inverse Compton, synchrotron, and scattering

One mechanism for the extended X-ray emission from 1138–262 is synchrotron radiation from the relativistic electrons in the radio lobes. The total extended X-ray emission corresponds to $`1\times 10^8`$ Jy at 2 keV. The radio spectral index for the integrated emission from 1138–262 at cm wavelengths is –1.1, and the integrated flux density from the source at 5 GHz is 160 mJy (Carilli et al. 1997). Extrapolating the radio spectrum to 2 keV ($`4.8\times 10^{17}`$ Hz) implies an expected flux density of $`3\times 10^{10}`$ Jy. Hence the X-ray emission cannot be a simple extrapolation of the radio synchrotron spectrum. The spectrum would have to flatten at frequencies above 15 GHz. Such a flattening can be ruled-out by the fact that the predicted I-band ($`3.7\times 10^{14}`$ Hz) emission would be $`6\times 10^6`$ Jy, while the upper limit is $`1\times 10^8`$ Jy.

A second possibility is Thompson scattering of X-ray emission from the nucleus. This is easily ruled out since the spectrum of the extended X-ray emission is different than that of the nucleus.

A third possible mechanism is inverse Compton (IC) emission, ie. up-scattering of the cosmic microwave background (CMB) by the relativistic electrons in the radio source. In this case comparing the radio synchrotron and X-ray IC emission constrains the magnetic fields within the radio source, since both depend on the relativistic electron population. Using equation 11 in Harris and Grindlay (1979), and using the integrated flux densities of the extended emission from 1138–262 in the radio and X-ray, the resulting magnetic field is 7$`\mu `$G assuming $`\alpha =1`$. For comparison, the minimum energy magnetic field is about 200 $`\mu `$G in the radio source (Carilli et al. 1998). Whether the fields and particles conspire to reach a minimum configuration remains unknown, however, a field value a factor 30 lower than the minimum energy value raises the internal pressure by a factor 150. If the fields are at the minimum energy level, then IC emission constitutes 0.1$`\%`$ of the total extended X-ray emission.

Brunetti et al. (2001) have shown that in some sources the energy density in the photon field in the radio lobes due to optical/IR photons from the AGN is larger than that in the CMB. In principle, this would allow for larger magnetic fields in the radio lobes. However, assuming a luminous (obscured) AGN of optical through IR luminosity $`10^{46}`$ erg s-1 in 1138–262 (as adopted by Brunetti et al. for 3C 295), and correcting the energy density in the CMB for redshift, it is clear from their Figure 2 that the energy density in the photon field will be dominated by the CMB beyond about 10 kpc from the nucleus. The X-ray emission in 1138–262 extends well beyond this radius.

Perhaps the most telling argument against IC being the dominant emission mechanism for the X-rays seen from 1138–262 is that there are off-sets between the positions of high surface brightness radio and X-ray emitting features (see Figures 1 and 2). In all radio galaxies for which IC X-ray emission has been seen, it is coincident with high surface brightness radio features (Harris et al. 2000, Harris, Carilli, & Perley 1994). If the dominant photon field is the microwave background, the observed X-rays are predominantly from electrons with Lorentz factors $`\gamma _{\mathrm{IC}}=1000`$, corresponding to radio emission at 270 MHz for $`\mathrm{B}=200\mu \mathrm{G}`$. These data imply that if the X-ray emission is IC, then the radio morphology of the source would have to change significantly from low to high frequencies. High resolution images at frequencies $`500`$ MHz are required to test this hypothesis, although 1.4 GHz images show no dramatic change in the source morphology relative to that seen at 5 GHz.

### 3.2 Hot Gas

A final mechanism we consider is thermal emission from hot gas. The extended emission from 1138–262 has a 2 to 10 keV luminosity of $`3\times 10^{44}`$ erg s-1. For comparison, the cluster atmosphere enveloping the powerful radio galaxy Cygnus A has a 2 to 10 keV luminosity of $`4\times 10^{44}`$ erg s-1 (Arnuad et al. 1987). The Cygnus A cluster atmosphere is roughly azimuthally symmetric, with a FWHM $`400`$ kpc (Reynolds and Fabian 1996). The spatial extent of the X-ray emission from 1138–262 is confined to a region of $`180\times 90`$ kpc, and is aligned along the radio axis. So while the extended X-ray luminosities are within a factor two for Cygnus A and 1138–262, the relative volumes differ by a factor of about 30. The implication is that the density of hot gas in 1138–262 must be a factor of five or so larger than that seen in Cygnus A, implying gas densities $`0.05`$ cm-3. The mass in hot gas is $`2.5\times 10^{12}`$ M.

The pressure in the hot gas is $`8\times 10^{10}`$ dynes cm-2, assuming $`\mathrm{T}=6\times 10^7`$ K. The minimum pressure in the radio source is $`6\times 10^{10}`$ dynes cm-2, while the pressure in the Ly$`\alpha `$ emitting gas is about $`1\times 10^9`$ dynes cm-2 (Pentericci et al. 1997). The pressure in the hot gas appears to be roughly adequate to confine both the radio source and the low filling factor, line emitting clouds. Note that this confinement mechanism applies only to the Ly$`\alpha `$ emitting regions within the maximum radius of the X-ray source.

The cooling time for such gas is $`1\times 10^9`$ years, which is more than an order of magnitude larger than the dynamical timescale for the system, but a factor three shorter than the age of the universe at that redshift. A possible heating mechanism of this gas is the expanding radio source, ie. the X-ray emission is from the ‘cocoon’ of shocked material enveloping the radio source, as predicted by dynamical models of supersonically expanding jets (Begelman, Blandford, & Rees 1986). Assuming that the radio luminosity is $`10\%`$ of the jet kinetic luminosity implies a lower limit to the jet kinetic luminosity of $`5\times 10^{46}`$ ergs s-1. The total energy stored in the X-ray gas is about $`2\times 10^{61}`$ ergs. So it appears that the radio jet could have heated the medium over a timescale of $`2\times 10^7`$ years, which is reasonable for a radio galaxy of the size of 1138–262 (Blundell & Rawlings 1999).

Our search for X-ray emission from 1138–262 was originally motivated by the extreme values of Faraday rotation, with the idea that the source is embedded in a dense, magnetized cluster atmosphere. While we have detected extended X-ray emission, the morphology of the emitting regions is distinctly different than that expected for a normal cluster atmosphere. So whence the large rotation measures? If the X-ray emission is from the coccoon of shocked gas around the radio lobe, then large rotation measures will occur if the shocked material contains magnetic fields of order 20 $`\mu `$G, ordered on scales of a few kpc (Bicknell, Cameron, & Gingold 1990).

The lack of extended emission perpendicular to the radio source axis allows us to set an upper limit to the 2 to 10 keV luminosity of an extended (relaxed) cluster atmosphere around 1138–262 of $`1.5\times 10^{44}`$ erg s-1, or $`<40\%`$ of the X-ray luminosity of the Cygnus A cluster. So while the overdensity of galaxies in the region suggests that the system will evolve into a massive cluster, the lack of an extended X-ray atmosphere is consistent with the idea that it is not yet a dynamically relaxed system. This is also suggested by the irregular velocity distribution of the Ly$`\alpha `$ emitting galaxies (Pentericci et al. 2000).

The lack of a hot cluster atmosphere also implies that the ambient medium into which the radio source is expanding in 1138–262 is very different than the diffuse cluster atmospheres at $`10^8`$ K seen in similar luminosity radio galaxies at lower redshift, such as Cygnus A and 3C 295. This difference is consistent with the suggestion by Rees (1989) that the environments of $`z>2`$ radio galaxies consist of a multiphase medium of cool ($`10^4`$ K) clouds or filaments, embedded in a virialized medium at 106 to 107 K. The passage of the bow shock driven by the expanding radio source would have a number of effects in such a multiphase medium (Kaiser & Alexander 1999). Rees (1989) points out that the densest clouds are relatively unaffected by the passage of the shock, and may be induced to collapse and form stars by the high pressure medium in which they suddenly find themselves (cf. Icke 1999). The intermediate density clouds are shocked, but the cooling time is shorter than the jet lifetime, perhaps giving rise to some of the Ly$`\alpha `$ emission directly associated with bright radio features. The lower density, higher filling factor gas is shock heated to $`10^{78}`$, giving rise to the extended X-ray emission.

Barthel and Arnaud (1996) have proposed a selection effect in radio galaxies in which jets propagating through regions of higher gas density have a higher conversion efficiency of jet kinetic energy into radio luminosity (see also Eales 1992). If the ambient medium is anisotropically distributed, then radio jets of a given kinetic luminosity propagating along the direction of highest ambient density will be brighter than ones propagating through low density regions. In this case, preferential alignment between the jet axis and the axis of highest density of the ambient material will occur naturally in radio flux limited samples (West 1999).

A final comparison we make is with the recent observations of the extended X-ray emission around the $`z=1.79`$ radio galaxy 3C 294 by Fabian et al. (2001). The observations of 3C 294 show a similar alignment between the X-ray and radio axes as that seen in 1138–262, and a similar lack of one-to-one correspondence between high surface brightness radio and X-ray features. In particular, for 3C 294 the X-ray emission extends significantly beyond the radio source maximum radius. Fabian et al. (2001) interpret the extended X-ray emission as being thermal emission from a cooling flow cluster atmosphere. Using the temperature-luminosity relation for X-ray emitting cluster atmospheres, they show that the existence of such a system at this high redshift argues for a low $`\mathrm{\Omega }_M`$ universe. Inherent in this argument is that the gas temperature is a gross measure of the gravitational mass of system. For 1138–262 we have argued that the gas heating is due to the expanding radio source, and hence unrelated to the gravitational potential of the system.

The existence of a hot, high filling factor, thermal gas in high $`z`$ radio galaxies has long been hypothesized in order to confine both the radio source and the optical line emitting clouds (Chambers et al. 1990, Bremer et al. 1992). It may be that in 1138–262 the Chandra observatory has revealed this pervasive hot medium.

The National Radio Astronomy (NRAO) is a facility of the National Science Foundation, operated under cooperative agreement by Associated Universities, Inc.. DEH and CLC acknowledge support from NASA grant GO0-1137B. DEH acknowledges support from NASA contract NAS8-39073. We acknowledge support from a programme subsidy provided by the Dutch Organization for Scientific Research (NWO).

References

Arnaud, K.A., Fabian, A.C., Eales, S.A, et al. 1987, MNRAS, 227, 241

Barger, A.J., Cowie, L.L., Mushotzky, R.F., & Richards, E.A. 2001, AJ, 121, 622

Barthel, P.D. & Arnaud, K. 1996, M.N.R.A.S. (letters), 283, 45

Begelman, M.C., Blandford, R.D., & Rees, M.J. 1984, Rev. Mod. Phys., 56, 255.

Bicknell, G.V., Cameron, R.A., & Gingold, R.A. 1990, ApJ, 357, 373

Blundell, K.M. & Rawlings, S. 1999, Nature, 399, 330

Bremer, M. N., Fabian, A. C., Sargent, W. L. W., Steidel, C. C., Boksenberg, A., Johnstone, R. M. 1992, MNRAS, 258, 23L

Brunetti, G., Cappi, M., Setti, G., Feretti, L., & Harris, D.E. 2001, A&A, 372, 755

Carilli, C.L., Perley, R.A., Dreher, J.W., & Leahy, J.P. 1991, ApJ, 383, 554

Carilli, C.L., Röttgering, H.J.A., van Ojik, R., Miley, G.K., & van Breugel, W.J.M. 1997, Ap.J. (Supp.), 109, 1

Carilli, C.L., Harris, D.E., Pentericci, L., Roetterging, H.J.A., Miley, G.K., & Bremer, M.N. 1998, ApJ, 496, 57L

Carilli, C.L., Miley, G.K., Röttgering, H.J.A., et al. 2001, in Gas & Galaxy Evolution, (ASP: San Francisco), eds. J. Hibbard, M. Rupen, & J. van Gorkom.

Chambers, K. C., Miley, G. K., van Breugel, W. J. M. 1990, ApJ, 363, 21

Eales, S.A. 1992, ApJ, 397, 49

Fabian, A.C., Crawford, C.S., Ettori, S., & Sanders, J.S. 2001, MNRAS, 322, L11

Garmire, G.P. 1997, BAAS, 190, 3404

Harris, D.E. & Grindlay, J.E. 1979, MNRAS, 188, 25

Harris, D.E., Carilli, C.L., & Perley, R.A. 1994, Nature, 367, 713

Harris, D.E., Nulsen, P.E., Ponman, T.J. et al. 2000, ApJ, 530, L81

Icke, V. 1999, The Most Distant Radio Galaxies, (Royal Netherlands Acadamy: Amsterdam), eds. H. Roettgering, P. Best, & M. Lehnert, 217

Kaiser, C.R. & Alexander, P. 1999, MNRAS, 305, 707

Kurk, J.D., Roettgering, H.J.A., Pentericci, L., Miley, G.K., & Pentericci, L. 2001, Rev Mex. A& A, in press

McCarthy, P.J. 1993, A.R.A.A., 31, 639

Monet, D., Bird, A., Canzian, B., et al. 1996, USNO-SA1.0, (U.S. Naval Observatory, Washington DC)

Miley, G.K. 1980, A.R.A.A., 18, 165

Pentericci, L., Röttgering, H.J.A., Miley, G.K., Carilli, C.L., & McCarthy, P. 1997, A&A, 326, 580

Pentericci, L., Röttgering, H.J.A., Miley, G.K., et al. 1998, ApJ, 504, 139

Pentericci, L., Kurk, J.D., Röttgering, H.J.A. et al. 2000, A& A, 361, L25

Rees, M.J. 1989, MNRAS, 239, 1P

Reynolds, C.S. & Fabian, A.C 1996, M.N.R.A.S., 278, 479

Roettgering, H.J.A., Best, P., & Lehnert, M. 1999, The Most Distant Radio Galaxies, (Royal Netherlands Acadamy: Amsterdam)

Stanford, S.A., Holden, B., Rosati, P., Tozzi, P., Borgani, S., Eisenhardt, P.R., & Spinrad, H. 2001, ApJ, 552, 504

Stark, A.A., Gammie, C.F., Wilson, R.W. et al. 1992, Ap.J. (Supp.), 79, 77

Steidel, C., Adelberger, K., Dickinson, M., Giavalisco, M., Pettin, M., & Kellogg, M. 1998, ApJ, 492, 428

Taylor, G.B. & Carilli, C.L. 2002, ARAA, in press for Vol 40

Ueno, Shiro, Koyama, Katsuji, Nishida, Minoru, Shigeo, Yamauchi, & Ward, Martin J. 1994, Ap.J. (letters), 431, L1

Ward, M.J. 1996, in Cygnus A, eds. C. Carilli & D. Harris, (Cambridge Univ. Press: Cambridge), p.43

Weisskopf, M.C., O’Dell, S.L., & Van Speybroeck, L.P. 1996, Proc. SPIE, 2805, 2

West, M.J. 1999, in The Most Distant Radio Galaxies, (Royal Netherlands Acadamy: Amsterdam), eds. H. Roettgering, P. Best, and M. Lehnert p. 365

Figure Captions

Figure 1 – The contours show the total (0.3 to 10 keV) X-ray emission from the $`z=2.156`$ radio galaxy PKS 1138–262 as observed with the ACIS-S detector on the Chandra Observatory. The Chandra image has been convolved with a Gaussian of FWHM = 2′′. In this and subsequent images the contours are a geometric progression in square root two. The first contour level is 2.3 counts per Gaussian beam. The greyscale shows the radio continuum image at 1.4 GHz, 2′′ resolution made with the Very Large Array (Carilli et al. 1997). The greyscale range is 0 to 150 mJy beam-1.

Figure 2 – A comparison of the X-ray emission and emission at other wavelengths in PKS 1138–262. The contours in all cases are of total X-ray emission convolved with a Gaussian of FWHM = 1′′. The first contour level is 0.9 counts per Gaussian beam (note: the highest four contour levels in these images are by a factor two, not root two). The crosses mark the position of the radio galaxy nucleus, and of an X-ray point source located 5′′ northwest of the nucleus.
a. The greyscale shows the 5 GHz radio continuum emission at 0.5′′ resolution. The greyscale range is 0 to 2 mJy beam-1. (from Carilli et al. 1997).
b. The greyscale shows the HST F606W image (from Pentericci et al. 1998). Note that 606nm corresponds to 192nm at $`z=2.2`$.
c. The greyscale shows the VLT image of the Ly$`\alpha `$ emission (from Kurk et al. 2001).

Figure 3 – Three-band images of X-ray emission from PKS 1138–262 convolved with a Gaussian of FWHM = 2′′. The crosses are the same as in Figure 2.
a. The soft X-ray (0.3 to 1.2 keV) emission. The greyscale shows the radio continuum emission at 5 GHz. The first contour level is 1.1 counts per Gaussian beam.
b. The mid X-ray (1.2 to 2.4 keV) emission. The first contour level is 0.9 counts per Gaussian beam.
c. The hard X-ray ($`>`$2.4 keV) emission. The first contour level is 1.8 counts per Gaussian beam.

Figure 4 – The surface brightness profiles averaged azimuthally in rings centered on the nucleus of 1138–262. The solid circles are for the two quadrants along the radio axis, while the open squares are for the two quadrants perpendicular to the radio axis. Both profiles have been normalized at the first point, corresponding to a radius of 3.5′′, and the points at smaller radii are omitted due to confusion by AGN emission. The solid line is a $`\beta `$ model fit to the data along the radio axis, with $`\beta =2.5\pm 1`$ and $`\mathrm{r}_\mathrm{c}=12\pm 4^{\prime \prime }`$. The dotted line is the background level.

Figure 5 – The X-ray spectrum of the AGN emission from PKS 1138–262. The model is a power-law spectrum with $`\alpha =0.8`$ and N(HI) = $`3.5\times 10^{22}`$ cm-2 at $`z=2.156`$.


text:
# The binary components of GRS 1915+105

## 1. Introduction

GRS 1915+105 (Castro-Tirado et al. 1994) is the prototypical microquasar, a galactic X-ray binary ejecting plasma clouds at v$``$0.92 c (Mirabel & Rodriguez 1994). It exhibits unique X-ray variability patterns (Greiner et al. 1996) which have been interpreted as accretion disk instabilities leading to an infall of parts of the inner accretion disk (Belloni et al. 1997). Based on its X-ray properties GRS 1915+105 was suspected to be the most massive stellar black hole candidate in the Galaxy (Morgan et al. 1997). Besides GRO J1655-40, GRS 1915+105 is one of only two galactic sources which are thought to contain a maximally spinning black hole (Zhang et al. 1997). It is therefore of great importance to know some details about the system components in the GRS 1915+105 binary in order to understand the conditions which lead to the unique X-ray, radio, and infrared characteristics.

## 2. Observations and Results

### 2.1. The donor

GRS 1915+105 is located in the galactic plane at a distance of $``$11 kpc (Fender et al. 1999) and suffers extreme extinction of 25–30 mag in the visual band. We therefore obtained spectroscopic observations in the near-infrared H and K bands using the ISAAC instrument on the VLT-Antu telescope. Based on spectra taken in the 5 wavelengths bands 1.56–1.64 $`\mu `$m, 1.63–1.72 $`\mu `$m, 2.05–2.17 $`\mu `$m, 2.17–2.29 $`\mu `$m and 2.29–2.41 $`\mu `$m we successfully identified absorption features from the atmosphere of the companion (mass donating star) in the GRS 1915+105 binary (Greiner et al. 2001a). The detection of 12CO and 13CO band heads (Fig. 1) plus a few metallic absorption lines suggests a K-M spectral type and luminosity class III (giant). Comparing the equivalent widths of the 12CO(3,1) 12CO(4,2) and 12CO(5,3) bandheads between the GRS 1915+105 spectrum and that of a K2 III standard star (observed with the same settings) implies that the donor of GRS 1915+105 has an unveiled magnitude of 14.5–15.0 mag (prior to extinction correction). This in turn implies an absolute magnitude of $`M_\mathrm{K}`$ = –2…–3, consistent with the giant classification. This K-M III classification also implies a constraint on the mass of the donor of $`M_\mathrm{d}`$ = 1.2$`\pm `$0.2 M (Greiner et al. 2001a). This proves that GRS 1915+105 is a low-mass X-ray binary (LMXB), as has been suggested earlier based on the IR spectral variability and on its position on the IR H-R diagram (Castro-Tirado et al. 1996), though the contrary claim has been made by Mirabel et al. (1997) and Martí et al. (2000) that based on the detection of He I, the donor in GRS 1915+105 should be a high-mass O or B star, and that accretion occurs predominantly from the wind of the donor.

### 2.2. The radial velocity curve

Subsequently, a series of medium-resolution spectra in the 2.39–2.41 $`\mu `$m range using the VLT-Antu equipped with ISAAC was obtained between April and August 2000 (Greiner et al. 2001b). Radial velocities were measured for the individual spectra by cross-correlation of the major CO band heads, using as template a spectrum of the K2 III star HD 202135 taken with the same setting. The results of this cross-correlation (after heliocentric correction) are shown in Fig. 2: The radial velocity periodogram shows a clear peak at a period of 33.5 days (top panel) which is interpreted as the orbital period $`P_{\mathrm{orb}}`$ of the binary system (Greiner et al. 2001b). All measurements, folded over this period, are shown in the lower panel of Fig. 2. The velocity amplitude is measured to be $`K_\mathrm{d}`$ = 140$`\pm `$15 km/s. With these values, the mass function, i.e. the observational lower limit to the mass of the compact object is (Greiner et al. 2001b)

$$f(M)\frac{(M_\mathrm{c}\mathrm{sin}i)^3}{(M_\mathrm{c}+M_\mathrm{d})^2}=\frac{P_{\mathrm{orb}}K_\mathrm{d}^3}{2\pi G}=9.5\pm 3.0M_{}.$$
(1)
GRS 1915+105 has shown several major jet ejection events for which the brightness and the velocities of both the approaching and the receding blobs could be measured (Mirabel & Rodriguez 1994, Fender et al. 1999). From this, the angle of the jet axis relative to the line of sight has been determined to be $``$ 70 $`\pm `$ 2 and constant over several years. With the assumption that the jet is perpendicular to the accretion disk and orbital plane, this provides a reasonably accurate measure of the inclination angle $`i`$ of the binary system. Knowing the inclination $`i`$, Eq. 1 can be solved for the mass of the accreting compact object (Fig. 3), giving $`M_\mathrm{c}=14\pm 4`$ M. Table 1 summarizes the orbital parameters of GRS 1915+105. Even when accounting for the relatively large error (which is dominated by the error in the determination of the velocity amplitude $`K_\mathrm{d}`$), GRS 1915+105 is the galactic binary with the largest mass function and the largest mass of its compact object. Previous record holders were V404 Cyg ($`f(M)=6.07\pm 0.05`$, $`M_\mathrm{c}=710`$ M; Shahbaz et al. 1996) and XTE J1118+480 ($`f(M)=6.00\pm 0.36`$, $`M_\mathrm{c}=6.510`$ M; McClintock et al. 2001).

We note that given the high mass-loss of the donor (which is needed to explain the large X-ray luminosity), it is most certainly less luminous than a non-interacting star of the same spectral type (Langer et al. 2000). This implies a larger mass for the donor, thus making the above black hole mass a lower limit. While it is difficult to assess the exact impact of the mass loss on the donor mass, the very extreme and conservative case of doubling the mass is plotted in Fig. 3. This indicates that the black hole mass could be higher by 1–3 M at maximum.

Distortions of the radial velocity curve due to X-ray heating (e.g. Phillips et al. 1999) are likely to be unimportant because of the long orbital period in GRS 1915+105. Using their eq. 2, one finds that the maximum possible deviation for GRS 1915+105 would be 8%, as compared to their 14% for GRO J1655-40. In analysing their data, Phillips et al. (1999) finally find that the true velocity amplitude in GRO J1655-40 is smaller than the measured one by 4.5%. Thus, one could expect a 2–3% effect in GRS 1915+105, compared to our presently given error which amounts to 11%. However, since the infrared flux contribution from the secondary star in GRS 1915+105 is so small, there is the possibility that phase-dependent changes in the continuum near the absorption features may result in an additional source of systematic error in the measured value of the velocity amplitude $`K_\mathrm{d}`$.

The systemic velocity of the GRS 1915+105 binary system is $`\gamma `$ = –3$`\pm `$10 km/s which implies that based on the galactic rotation curve (Fich et al. 1989) the kinematic distance of GRS 1915+105 is $`d=12.1\pm 0.8`$ kpc, intermediate between earlier estimates (Mirabel & Rodriguez 1994, Fender et al. 1999).

## 3. Discussion

### 3.1. Formation of a 14 M black hole

The knowledge of the mass of the black hole in GRS 1915+105 has several implications for both the understanding of the physics in microquasars (for a review see Mirabel & Rodriguez 1998; Greiner 2000) as well as some broader astrophysical concepts. Most importantly, the formation of a 14 M black hole in a low-mass binary poses an interesting challenge for binary evolution. Stellar evolution of stars in a binary system proceeds differently from single star evolution due to primarily the mass transfer between the system components and/or common-envelope phases. There are, in general, two different paths for the black hole formation in a binary system. First, the progenitor system could be wide, and during the common envelope phase the low-mass (main sequence) star of $``$1 M will spiral into the envelope of the massive giant (progenitor of the black hole), causing the orbit to shrink (Kalogera 1999, 2001). Based on the measured system parameters (Tab. 1), the deduced orbital separation of the binary components in GRS 1915+105 is $`a`$ = 108$`\pm `$4 R. Thus, orbital contraction through a common-envelope phase caused by the expansion of the massive progenitor to typically $`>`$1000 R, is conceivable for GRS 1915+105. Second, the evolution could start with a progenitor system smaller than today, provided that the binary component interaction is delayed until after He burning has ceased (Brown et al. 1999). In this case, the time between the wind phase and the core-collapse is short, and black hole masses in the 5–10 M range have been reproduced for cases where the initial He-star progenitor is in the 10–25 M range, corresponding to initial primaries in the 25–45 M range (Brown et al. 2001, Kalogera 2001). How much mass is finally lost depends on the radii evolution of the two progenitor stars, and it remains to be shown whether black hole masses above 10 M can be achieved.

In order to produce BH masses $`>`$10 M, the progenitor might have been a massive Wolf-Rayet star. However, Wolf-Rayet stars are known to have a much larger wind-loss rate, and it is therefore unclear whether higher progenitor masses indeed will lead to higher final black hole masses. Also, the effects of rapid rotation on the mass loss become more serious. Even if one assumes no mass-loss at all from the helium star, the predicted masses for helium cores of massive stars can not account for BH masses $`>`$12–14 M (Wellstein & Langer 1999, Hurley et al. 2000). An alternative possibility to produce high-mass ($`>`$10 M) black holes may be to invoke hierarchical triples as progenitors (Eggleton & Verbunt 1986).

Whatever the formation path may be, it seems that the mass of the black hole in the binary GRS 1915+105 may provide the impetus to study the formation of black holes with masses above 10 M in more detail.

### 3.2. On the luminosity and X-ray variability of GRS 1915+105

With the values of Table 1, the implied Roche lobe size of the donor star is 21$`\pm `$4 R, in good agreement with the size of a K-M giant which therefore can be expected to fill its Roche lobe. Thus, accretion in GRS 1915+105 is likely to occur via Roche lobe overflow.

The knowledge of the mass of the accretor in GRS 1915+105 also allows us to possibly understand the rapid and large-amplitude X-ray variability which is seen only in GRS 1915+105 (Greiner et al. 1996). Combining the observed X-ray intensity with the measured distance immediately shows that this large-amplitude X-ray variability occurs near or even above the Eddington limit $`\dot{M}_{\mathrm{Edd}}`$. Such high accretion rates are seemingly never reached by other canonical black hole transients (e.g. GRO J1655-40) which usually operate in the 0.1-0.2 $`\dot{L}_{\mathrm{Edd}}`$ range. At these lower rates the accretion disks are likely gas pressure dominated, and thus viscously and thermally stable. In contrast, the uniquely high $`\dot{M}`$/$`\dot{M}_{\mathrm{Edd}}`$ ratio in GRS 1915+105 suggests that its inner part of the accretion disk is radiation pressure dominated, which in turn makes the inner disk quasi-spherical and thermally unstable. This property potentially provides the clue for the spectacular and unique X-ray variability in GRS 1915+105 (Greiner et al. 1996) which has not been found in any other X-ray binary. Indeed, under the assumptions of a corona dissipating 50% of the energy and a jet carrying a luminosity-dependent fraction of energy, Janiuk et al. (2000) show that the time-dependent behaviour of GRS 1915+105 can be well reproduced with the $`\alpha `$-viscosity prescription in a radiation pressure dominated region.

Also, while it is tempting to conclude that jet ejection occurs because the black hole can not accept this copious supply of matter, it is important to remember that jet ejection occurs also in these other sources, at e.g. 0.2 $`\dot{M}_{\mathrm{Edd}}`$ in GRO J1655-40, and thus near/super-Eddington accretion cannot be the determining factor for relativistic jets.

### 3.3. On the black hole spin in GRS 1915+105 and GRO J1655-40

The knowledge of the black hole mass also allows us to place constraints on the question of the black hole spin in GRS 1915+105 and GRO J1655-40. Information on the black hole spin in both sources has been deduced from two completely different sources: (1) Accretion disks around a (prograd) spinning black hole extend farther down towards the black hole, and thus allow the temperature of the disk to be higher. Since both GRS 1915+105 and GRO J1655-40 exhibit a thermal component in their X-ray spectra which has an unprecedently high temperature as compared to all other black hole transients (during outbursts), it has been argued that this is due to the spin of the black holes in these two sources, while the majority of black hole transients has non-rotating black holes (Zhang et al. 1997). (2) Several black hole binaries, including GRS 1915+105 and GRO J1655-40, show nearly-stable quasi-periodic oscillations (QPO) in their X-ray emission. The frequencies $`f`$ for these QPOs are 300 Hz in GRO J1655-40 (Remillard et al. 1999) and 67 Hz in GRS 1915+105 (Morgan et al. 1997). Very recently, a related second stable QPO has been discovered for each of these (Strohmayer 2001a,b). Most of the models proposed to explain these QPO either rely or depend on the spin of the accreting black hole.

If the black hole mass in GRS 1915+105 is indeed no larger than 18 M (Fig. 3), the deduction of the black hole spin by Zhang et al. (1997) is inconsistent to any of the four models on the origin of QPOs:

1. If associated with the Keplerian motion at the last stable orbit around a (non-rotating) black hole according to the simple relation $`f`$ (kHz)$`=2.2/M_{\mathrm{BH}}`$ (M), it gives a surprising agreement with the optically determined mass for GRO J1655-40 of 7 M (Orosz & Bailyn 1997), but is off by a factor of 2 for GRS 1915+105.
2. If associated with the trapped g-mode (diskoseismic) oscillations near the inner edge of the accretion disk (Okazaki et al. 1987, Perez et al. 1997, Nowak et al. 1997) would give consistency for a nearly maximally spinning black hole in GRO J1655-40, and a non-spinning black hole in GRS 1915+105, contrary to the argument that the spin in both systems should be very similar because of the extremely high temperature of the accretion disk (Zhang et al. 1997).
3. If associated with the relativistic dragging of inertial frames around a spinning black hole (Cui et al. 1998) which would cause the accretion disk to precess, the implied specific angular momentum (spin) of the black hole in GRS 1915+105 would be $`a`$ 0.8, thus considerably lower than the $`a`$ 0.95 deduced for GRO J1655-40. This is in conflict with the nearly identical accretion disk temperatures for both sources as deduced from the X-ray spectra, which in turn requires a larger spin for GRS 1915+105.
4. If associated with oscillations related to a centrifugal barrier in the inner part of the accretion disk (Titarchuk et al. 1998), the product of QPO frequency and black hole mass is predicted to be proportional to the accretion rate, implying that the accretion rate in GRO J1655-40 should be a factor $``$10 larger than in GRS 1915+105. This is certainly not the case, since the black hole mass in GRO J1655-40 is about 7 M (Orosz & Bailyn 1997).

A different phrasing of the above item (1) is that one cannot argue that based on the QPO frequency difference the mass of the black hole in GRS 1915+105 should be a factor of 4–5 larger than that in GRO J1655-40. This relation would only hold for the assumption that the QPO frequencies are correlated to the Kepler frequency. Our mass determination shows that the QPO frequency does not scale linearly with the mass of the black hole. As a consequence, their origin is very likely not related to the Keplerian frequency.

Thus, none of these four models provides a satisfactory solution if one adopts the interpretation that the high accretion disk temperatures are a measure of the black hole spin (Zhang et al. 1997). If, on the contrary, this latter interpretation is dropped, and the spin becomes a free parameter, the first three models could be applicable.

It should be noted that the applicability of the applied standard Shakura-Sunyaev disk model to deduce accurate accretion disk temperatures has been also questioned on other grounds (Sobczak et al. 1999, Merloni et al. 2000). In addition, there also exist alternative models, so-called slim disk models, which can reproduce high-temperature disks also for non-rotating black holes (see review by Balbus & Hawley, or a recent application by Watarai et al. 2000).

### 3.4. Other implications

Besides the above issues there are a number of other implications of the high value of the black hole mass in GRS 1915+105 which are summarized below: (1) It proves beyond doubt (and beyond any details in nuclear matter physics which occasionally argue even about 7 M neutron stars; see e.g. Negi & Durgapal 2001) the existence of stellar black holes. (2) It has been argued that the compact objects in most black-hole transients cluster around 7 M (Bailyn et al. 1998), contrary to the expectation (assuming that the initial mass function is weighted toward lower mass stars) of a monotonic distribution of black holes biased toward lower masses. The mass of GRS 1915+105 adds to the hints (based on V404 Cyg and Cyg X-1) that there is a range of black hole masses in binaries, though the bias towards low masses (3-5 M) still remains to be clarified. (3) It invalidates the counter-argument against these $``$7 M black holes being neutron stars plus a self-gravitating, 5 M accretion disk (Kundt 1999). An accretion disk with a mass of $``$13 M seems impossible, both because of the higher mass-transfer which would have to last many Myrs, as well as the unstable location of the neutron star within the disk.

Whatever the solutions to the individual implications may be, the microquasar GRS 1915+105 continues to surprise scientists with unforeseen challenges, and to be among the most exciting X-ray sources on the sky.

Acknowledgement: Though I haven’t had, unfortunately, any direct collaboration with Jan van Paradijs, I’m grateful for having known him and having had the chance to interact with and benefit from him. Thanks to the Organizing Committee and the opportunity to present this work.

## References

Bailyn C.D., Jain R.K., Coppi P., Orosz J.A. 1998, ApJ 499, 367

Balbus S. A., Hawley J.F. 1998, Rev. Mod. Phys. 70, 1

Belloni T., Mendez M., King A.R., van der Klis M., van Paradijs J., 1997, ApJ 488, L109

Brown G.E., Lee C.-H., Bethe H.A. 1999, New Astr. 4, 313

Brown G.E., Lee C.-H., Tauris T.M. 2001, New Astr. 6, 331

Castro-Tirado A.J., Brandt S., Lund N., et al. 1994, ApJS 92, 469

Castro-Tirado A.J., Geballe T.R., Lund N., 1996, ApJ 461, L99

Cui W., Zhang S.N., Chen W. 1998, ApJ 492, L53

Eggleton P.P., Verbunt F. 1986, MNRAS 220, 13p

Fender R.P., Garrington S.T., McKay D.J., et al. 1999, MNRAS 304, 865

Fich M., Blitz L., Stark A.A. 1989, ApJ 342, 272

Greiner J., Morgan E.H., Remillard R.A. 1996, ApJ 473, L107

Greiner J. 2000, in Cosmic Explosions, Proc. 10th Ann. Astrophys. Conf. in Maryland, eds. S. Holt, W.W. Zhang, AIP 522, p. 307

Greiner J., Cuby J.G., McCaughrean M.J., Castro-Tirado A.J., Mennickent R.E. 2001a, A&A 373, L37

Greiner J., Cuby J.G., McCaughrean M.J. 2001b, Nat. (in press)

Hurley J.R., Pols O.R., Tout C.A. 2000, MNRAS 315, 543

Janiuk A., Czerny B., Siemiginowska A., 2000, ApJ 542, L33

Kalogera V. 1999, ApJ 521, 723

Kalogera V. 2001, in Evolution of Binary and Multipe Star Systems, eds. P. Podsiadlowski et al., ASP Conf. Ser. (in press; astroph/0012064)

Kundt W., 1999, F. Giovannelli & L. Sabau-Graziati (eds.), Mem. Societa Astron. Italiana 70, Nos. 3/4, p. 1105

Langer N., Deutschmann A., Wellstein S., Höflich P., 2000, A&A 362, 1046

Martí J., Mirabel I.F., Chaty S., Rodriguez L.F. 2000, A&A 356, 943

McClintock J.E., Garcia M.R., Caldwell N., Falco E.E., Garnavich P.M. et al. 2001, ApJ 551, L147

Merloni A., Fabian A.C., Ross R.R. 2000, MNRAS 313, 193

Mirabel I.F., Rodriguez L.F. 1994, Nat. 371, 46

Mirabel I.F., Bandyopadhyay R., Charles P.A., Shahbaz T., Rodriguez L.F., 1997, ApJ 477, L45

Mirabel I.F., Rodriguez L.F. 1998, Nat. 392, 673

Morgan E.H., Remillard R.A., Greiner J. 1997, ApJ 482, 993

Negi P.S., Durgapal M.C. 2001, A&SS 275, 299

Nowak M.A., Wagoner R.V., Begelman M.C., Lehr D.E. 1997, ApJ 477, L91

Okazaki A.T., Kato S., Fukue J. 1987, PASJ 39, 457

Orosz J.A., Bailyn C.D., 1997, ApJ 477, 876

Perez C.A., Silbergleit A.S., Wagoner R.V., Lehr D.E., 1997, ApJ 476, 589

Phillips S.N., Shahbaz T., Podsiadlowski Ph. 1999, MNRAS 304, 839

Remillard R.A., Morgan E.H., McClintock J.E., Bailyn C.D., Orosz J.A. 1999, ApJ 522, 397

Shahbaz T., Bandyopadhyay R., Charles P.A., Naylor T. 1996, MNRAS 282, 977

Sobczak G.J., McClintock J.E., Remillard R.E., Bailyn C.D., Orosz J.A. 1999, ApJ 520, 776

Strohmayer, T.E. 2001a, ApJ 552, L49

Strohmayer, T.E. 2001b, ApJ 554, L169

Titarchuk L., Lapidus I., Muslimov A. 1998, ApJ 499, 315

Watarai K., Fukue J., Takeuchi M., Mineshige S. 2000, Publ. Astron. Soc. Japan 52, 133

Wellstein S., Langer N. 1999, A&A 350, 148

Zhang N.S., Cui W., Chen W. 1997, ApJ 482, L155