text:
# Constructing Phantom Cosmologies from Standard Scalar Field Universes

## Abstract

We illustrate how form-invariance transformations can be used for constructing phantom cosmologies from standard scalar field universes. First, we discuss how to relate two flat Friedmann-Robertson-Walker cosmologies with different barotropic indexes $`\gamma `$ and $`\overline{\gamma }`$. Then, we consider the particular case $`\overline{\gamma }=\gamma `$, and we show that if the matter content is interpreted in terms of self-interacting scalar fields, then the corresponding transformation provides the link between a standard and a phantom cosmology. After that, we illustrate the method by considering models with exponential potentials. Finally, we also show that the mentioned duality persists even if the typical braneworld modifications to the Friedmann equation are considered.

Current observations can accommodate matter which violates the weak energy condition $`\rho >0`$, $`\rho +p0`$ obs . One possible candidate for that sort of matter is a fluid with negative enough pressure. If the fluid has a barotropic equation of state of the form $`p=(\gamma 1)\rho `$, then $`\gamma <1`$. Alternatively, one can invoke matter of scalar nature so that the square of the time derivative of the scalar nature enters the energy-momentum with an opposite sign compared to the standard situation. The cosmological models filled with that kind of matter are referred to as phantom or ghost cosmologies phan . Even though the theoretical understanding of the subject is limited, we can rely on the motivation for phantom matter which is provided by string theory str . Moreover, we can rest assured that it does not go against the observations when considering such exotic matter. Those are the main reasons why phantom models have begun to attract the attention of an increasing number of theoreticians, including those concerned, with exact solutions DabStaSzy03 ; FeiJhi03 ; NojOdi03a  or cosmological dynamics SinSamDad03 . There is also an interesting debate on what would be the fate of a universe with large negative pressure fate .

One of the ways to obtain new exact solutions of the Einstein field equations is to exploit their form-invariance symmetries. In cosmology, that method has been used Chi02 ; AguChiJak03  for generating inflationary cosmologies from noninflationary ones. The Einstein equations for a flat Friedmann-Robertson-Walker (FRW) cosmological model with scale factor $`a`$, and filled with a perfect fluid with energy density $`\rho `$ and pressure $`p`$, are

$`3H^2=\rho ,`$ (1)
$`\dot{\rho }+3H(\rho +p)=0,`$ (2)
where a dot denotes differentiation with respect to the cosmic time, and $`H=\dot{a}/a`$. The units used in the paper are such that $`c=8\pi G=1`$.

A form-invariance transformation is a prescription to relate the quantities $`a`$, $`H`$, $`\rho `$ and $`p`$ of our original model, with the quantities $`\overline{a}`$, $`\overline{H}`$, $`\overline{\rho }`$ and $`\overline{p}`$ of what we will call the transformed model, so that it satisfies

$`3\overline{H}^2=\overline{\rho }`$ (3)
$`\dot{\overline{\rho }}+3\overline{H}(\overline{\rho }+\overline{p})=0.`$ (4)
In Chi02  it was shown that one such transformation is given by

$`\overline{\rho }=\overline{\rho }(\rho ),`$ (5)
$`\overline{H}=\left({\displaystyle \frac{\overline{\rho }}{\rho }}\right)^{1/2}H`$ (6)
$`\overline{p}=\overline{\rho }+\left({\displaystyle \frac{\rho }{\overline{\rho }}}\right)^{1/2}(\rho +p){\displaystyle \frac{d\overline{\rho }}{d\rho }}`$ (7)
where $`\overline{\rho }=\overline{\rho }(\rho )`$ is an invertible arbitrary function.

Let us consider now perfect fluids with a barotropic equation of state $`p=(\gamma 1)\rho `$, where $`\gamma `$ is the so-called barotropic index. The barotropic indexes of the original and transformed fluid are related by

$$\overline{\gamma }=\left(\frac{\rho }{\overline{\rho }}\right)^{3/2}\frac{d\overline{\rho }}{d\rho }\gamma $$
(8)
A scalar field $`\varphi `$ self-interacting through a potential $`V(\varphi )`$ can also be interpreted in the language of perfect fluids. It is well know that one can make the following identifications:

$`\rho ={\displaystyle \frac{1}{2}}\dot{\varphi }^2+V(\varphi ),`$ (9)
$`p={\displaystyle \frac{1}{2}}\dot{\varphi }^2V(\varphi ).`$ (10)
Under the transformation

$$\overline{\rho }=n^2\rho ,$$
(11)
$`n`$ being a constant, the scalar field and potential transform, respectively, as

$`\dot{\overline{\varphi }}^2=n\dot{\varphi }^2`$ (12)
$`\overline{V}=n^2\left({\displaystyle \frac{1}{2}}\dot{\varphi }^2+V\right){\displaystyle \frac{n}{2}}\dot{\varphi }^2,`$ (13)
Let us now investigate the simple transformations which preserve the energy density of the fluid; that is, $`\overline{\rho }=\rho `$. From Eqs. (6)-(8), we get

$`\overline{H}=\pm H,`$ (16)
$`\overline{p}=\rho \pm (\rho +p)`$
$`\overline{\gamma }=\pm \gamma .`$
Clearly, the $`(+)`$ branch is the identical transformation. However, as will be shown immediately, the $`()`$ branch leads to a phantom symmetry of the FRW equations (1) and (2), because it corresponds to

$$\overline{H}=H,\overline{\rho }+\overline{p}=(\rho +p),\overline{\gamma }=\gamma .$$
(17)
From these relations we see that the two quantities $`\rho `$ and $`H\gamma `$ are invariant under the phantom symmetry transformation.

On the other hand, the kinetic energy and the potential of the field transform as

$`\dot{\overline{\varphi }}^2=\dot{\varphi }^2`$ (18)
$`\overline{V}(\overline{\varphi })=\dot{\varphi }^2+V(\varphi ),`$ (19)
so in fact, up to a constant, we see that $`\overline{\varphi }=i\varphi `$. That is, the transformed scalar field is related to the original one by a Wick rotation, and we finally have

$$\overline{\rho }=\frac{1}{2}\dot{\varphi }^2+\overline{V}(\overline{\varphi }).$$
(20)
The sign in the kinetic part of the energy density indicates that what we actually have now is a phantom cosmology, with a phantom field $`\varphi `$, together with a phantom potential $`\overline{V}(\overline{\varphi })`$, which is a real function of $`\varphi `$. Note that the self-interaction potential in the phantom and standard cases are different, and are linked by the transformation rule (19).

Finally, in order to establish completely the link between flat FRW scalar field cosmologies with a barotropic equation of state and their phantom counterparts, we have to get $`\overline{a}`$ as a function of $`a`$. This follows directly from the relationship $`\overline{H}=H`$, which gives in turn $`\overline{a}a^1`$. As we have seen, this duality, which was noticed independently in DabStaSzy03 , is, in fact, a particular case of a more general form-invariance symmetry. Noticeably, this duality in the scale factor, which here allows one to establish a link between phantom and standard scalar field cosmologies, is identical to the one appearing in the pre-big bang models.

An straightforward extension would be the application of the method to the case of universes filled with several fluid components, or with multiple scalar fields, for instance in the fashion of o(n) , where an $`O(\mathrm{N})`$ field configuration was considered. The details of how this generalization could be done can also be found in Chi02 . Nevertheless, for the sake of simplicity, here we will restrict the discussion to single field configurations.

For illustration purposes, let us take an exponential potential

$$V=V_0\text{e}^{A\varphi },$$
(21)
where $`A`$ is a constant. As is well known, power-law solutions will exist if the scalar field has a logarithmic dependence on time. Such solutions to the Einstein-Klein-Gordon equation set will correspond to

$`a=t^{2/A^2},\varphi ={\displaystyle \frac{2}{A}}\mathrm{ln}t,V_0={\displaystyle \frac{2(6A^2)}{A^4}}.`$ (22)
If we apply now the phantom symmetry given by Eq. (19), the transformed scalar field and potential is

$$\overline{\varphi }=i\frac{2}{A}\mathrm{ln}t,\overline{V}=\overline{V}_0\text{e}^{iA\overline{\varphi }},$$
(23)
with

$$\overline{V}_0=\frac{2(6+A^2)}{A^4},$$
(24)
so that the transformed scale factor reads

$$\overline{a}=t^{2/A^2}.$$
(25)
Therefore, it is clear that solutions (22) y (25) satisfy the symmetry condition $`a\overline{a}=1`$.

A question that immediately comes to mind is whether the form-invariance symmetry can be extended. To this end, the next line of attack is the extension to brane scenarios (see NojOdi03b  for an study on brane cosmologies induced by phantom fields). The Friedmann equation for a perfect fluid obeying the standard conservation equation (2) on a flat FRW brane becomes

$$3H^2=\rho \left(1+\frac{\rho }{2\lambda }\right)$$
(26)
where bulk effects have been switched off, $`\lambda `$ being the brane tension.

The conservation equation (2) and the modified Friedmann equation (26) are form invariant under the symmetry transformation given by

$`\overline{H}=\sqrt{{\displaystyle \frac{\overline{\rho }(1+\overline{\rho }/(2\lambda ))}{\rho (1+\rho /(2\lambda ))}}}H`$ (27)
$`\overline{p}=\overline{\rho }+\sqrt{{\displaystyle \frac{\overline{\rho }(1+\overline{\rho }/(2\lambda ))}{\rho (1+\rho /(2\lambda ))}}}(\rho +p){\displaystyle \frac{d\overline{\rho }}{d\rho }}`$ (28)
It follows that, if we apply a symmetry transformation such that $`\overline{\rho }=\rho `$, we will obtain exactly the two same results as in the relativistic case, that is, $`\overline{H}=\pm H`$ and $`\overline{a}a^{\pm 1}`$, although the $`()`$ branch would correspond now to a brane phantom cosmology.

Summarizing, we have shown that there exists a duality between phantom and standard (or nonexotic) flat FRW cosmologies, both in the relativistic and braneworld frameworks. As a matter of fact, this duality is a particular case of the form-invariance symmetry of the corresponding gravitational field equations. It may be worth looking further into this problem, to see whether explicit relations could be obtained in less straightforward cases.

L.P.C. is supported by the University of Buenos Aires under project X223, and by the Consejo Nacional de Investigaciones Científicas y Técnicas. R.L. is supported by the Basque Government through fellowship BFI01.412, the Spanish Ministry of Science and Technology jointly with FEDER funds through research grant BFM2001-0988, and the University of the Basque Country through research grant UPV00172.310-14456/2002.


format: text

text:
# XMM-Newton observation of the 5.25 ms transient millisecond pulsar XTE J1807–294 in outburst

## 1 Introduction

In the past few years conclusive evidence that Low Mass X–ray Binaries (LMXRBs) contain fast spinning, weakly magnetic neutron stars has been gathered. Coherent periodicities have been discovered in four sources SAX J1808.4–3654 (401 Hz, Wijnands & van der Klis 1998), XTE J1751–305 (435 Hz, Markwardt et al. 2002), XTE J0929–314 (185 Hz, Galloway et al. 2002) and recently XTE J1807–294 (191 Hz, Markwardt, Smith & Swank 2003).

All these four sources are transients, i.e. exhibit only sporadic activity and for most of the time remain in a state of low level activity (for a review see e.g. Campana et al. 1998a). These transient pulsars share many similarities among them, but are rather peculiar among transient sources containing a neutron star (usually named Soft X–ray Transients, SXRT). Their peak X–ray luminosities are rather faint ($`0.01L_{\mathrm{Edd}}`$), whereas SXRTs have peak luminosities of $`0.11L_{\mathrm{Edd}}`$. They show X–ray pulsations, whereas for SXRTs and LMXRBs tight upper limits exist (e.g. Vaughan et al. 1994). The discovery of X–ray pulsation during the outburst phase has also allowed to reveal an orbital periodicity, which in all cases is shorted than $`2`$ hr (2.01 hr for SAX J1808.4–3654, Chakrabarty & Morgan 1998; 42.4 min for XTE J1751–305 Markwardt et al. 2002; 43.6 min for XTE J0929–314 Galloway et al. 2002; 40.1 min for XTE J1807–294, Markwardt, Juda & Swank 2003). These periods are much shorter than what usually found in SXRTs ($`\stackrel{>}{}4`$ hr). These differences led some authors to suggest that transient pulsars might form a separate (sub-)class among SXRTs, called faint transients (Heise et al. 1998; King 2000; in’t Zand 2001).

Faint transients are poorly known. The best studied source is SAX J1808.4–3658. It was studied across three outbursts (Wijnands & van der Klis 1998; Wijnands et al. 2001), in quiescence (Campana et al. 2002), in the optical (Homer et al. 2001) and radio (Gaensler et al. 1999) bands.

Concerning XTE J1807–294 very little is known. The source was discovered during a Galactic plane scan by RXTE/PCA (Markwardt, Smith & Swank 2003). The source was observed at a level of 33, 38, 58, 41 and 20 mCrab (2–10 keV) on 2003 Feb 16.7, 19.8, 21.6, 22.6 and Mar 13.9, respectively. A follow-up Chandra observation on 2003 Mar 10.8 allowed to obtain a much better position ($`1^{\prime \prime }`$ uncertainty). Pulsations were still detectable and allowed a refined orbital solution to be obtained from RXTE data (Markwardt, Juda & Swank 2003).

Here we report on an XMM-Newton (discretionary time) TOO observation of XTE J1807–294 on 2003 Mar 22 (see also Campana et al. 2003). In Section 2 we describe the observational data and the reduction strategy. In Section 3 we describe our results on spectral and timing analysis. Conclusions are reported in Section 4.

## 2 XMM-Newton observation and data analysis

The XMM-Newton Observatory (Jansen et al. 2001) includes three 1500 cm2 X–ray telescopes each with a European Photon Imaging Camera (EPIC, 0.1–15 keV) at the focus. Two of the EPIC imaging spectrometers use MOS CCDs (Turner et al. 2001) and one uses pn CCDs (Strüder et al. 2001). Reflection Grating Spectrometers (RGS, 0.35–2.5 keV, den Herder et al. 2001) are located behind two of the telescopes.

XTE J1807–294 was observed on March 22, 2003 from 13:57 to 18:40. MOS1 was operating in Small Window mode, MOS2 in Prime Full Window mode and pn in timing mode, all with thick filters. RGS1 and RGS2 were operating in the Spectro + Q standard mode. Due to a problem in the field acquisition the source was originally pointed 1 arcmin apart. This problem was then corrected, resulting in an overall efficiency decrease of the pn, MOS1 and RGS instruments. We generated final product using sas v. 5.4.1. A small soft proton flare event occurred at the beginning of the observation, we excluded the first 2.5 ks from the analysis. We obtained net exposure times of 14 ks and 9 ks for MOS2 and the other instruments, respectively.

### 2.1 MOS cameras

The source is extremely bright such that the MOS cameras are saturated and the radial support structure of the mirrors is visible. Piled-up rate for the MOS2 camera is $`8`$ counts s-1 (0.5–10 keV). Source position is consistent with that determined with Chandra (Markwardt, Juda & Swank 2003).

A simple way round to pile-up is to extract photons for spectral analysis in an annulus around the source, thus excluding the inner (piled-up) core. We followed this strategy for analysing the MOS2 data, whereas for MOS1 data the window is too small (see also the analysis carried out on Mkn 421, presenting similar problems, by Ravasio et al. 2003, in preparation). We extracted photons from an annulus of $`30^{\prime \prime }`$ and $`60^{\prime \prime }`$ inner and outer radii, respectively. Background was extracted from two different outer circular regions $`2^{}`$ each. We selected only single-pixel events and make sure that no pile-up is present (through the sas tool epatplot). This is also in line with prescriptions by Molendi & Sembay (2003)111http://xmm.vilspa.esa.es/docs/documents/CAL-TN-0036-1-0.ps.gz. The count rate within the 0.5–10 keV energy band is $`1.01\pm 0.01`$ counts s-1. Response matrix and ancillary files were then generated with rmfgen and arfgen. Spectral data were rebinned to have 80 counts per bin.

### 2.2 pn camera

Timing mode was selected to have the opportunity to reveal X–ray pulsations. Given the high count rate of the source, this mode is particularly useful to avoid any pile-up. We extracted the source spectrum from the raw image, taking a column 10-pixel wide and the background from two separate boxes 5-pixels wide. We consider only “single” and “double” events. No X–ray bursts were detected. The source count rate in the 0.5–10 keV band is $`34.0\pm 0.1`$ counts s-1. Background is at a level of $`1.4\%`$ in the full band. Pile-up is negligible and dead time is at a level of $`1.5\%`$. Response matrix and ancillary file for spectral analysis were generated using sas tasks. Spectral data were rebinned to have 80 counts per bin.

### 2.3 RGS

The standard sas procedure rgsproc was used to derive the first order RGS spectra. In the 7–36 Å band rates are $`0.39\pm 0.01`$ and $`0.56\pm 0.01`$ count s-1 for RGS1 and RGS2, respectively. For spectral analysis together with the other instruments we rebinned the data to 100 counts per bin. We further inspect the RGS spectra for emission or absorption lines. In this case we rebin the spectra to 20 counts per bin (see below).

## 3 Spectral analysis

The spectral analysis has been carried out in the 0.5–10 keV energy range for MOS2 and pn and in 7–36 Å for the two RGSs, using XSPEC (v11.2). All spectral uncertainties are given at $`90\%`$ confidence level for one degree of freedom ($`\mathrm{\Delta }\chi ^2=2.71`$). We first tried the standard model for SXRTs, i.e. an absorbed soft component (black body) plus a hard component (power law; e.g. Campana et al. 1998b). The fit is good with a $`\chi _{\mathrm{red}}^2=1.10`$ (1571 degrees of freedom, d.o.f.). Fit parameters can be found in Table 1. The column density is $`(6.3\pm 0.1)\times 10^{21}\mathrm{cm}^2`$ which is slightly higher than the Galactic value of $`3\times 10^{21}\mathrm{cm}^2`$. A comparable (but more physical) model is obtained by replacing the power law with a Comptonization model like COMPTT (Titarchuk 1994) as observed in SAX J1808.4–3658 (Gierlinski, Done & Barret 2002 but not in XTE 1751–305, Miller et al. 2003). In this case we obtain $`\chi _{\mathrm{red}}^2=1.10`$ (1569 d.o.f., see Fig. 1). Systematic errors at $`1.6\%`$ level provide an acceptable fit. The column density is $`4.6_{0.2}^{+0.3}\times 10^{21}\mathrm{cm}^2`$, consistent with the Galactic value. The black body temperature is $`kT=0.80\pm 0.03`$ keV and its radius (assuming spherical symmetry) is $`R=(2.2\pm 0.2)d_8^2f^2`$ km (with $`d_8`$ the unknown distance to XTE J1807–294 in units of 8 kpc222The source lies just in the direction of the Galactic center. and $`f`$ the spectral hardening factor obtained as the ratio of the color temperature to the effective temperature, e.g. Merloni, Fabian & Ross 2000). The temperature of the soft input photon is $`0.24\pm 0.02`$ keV. In the Comptonization model the temperature of the corona and the plasma optical depth are tightly related. The best fit is for a relatively cool corona ($`kT=9.6_{4.0}^{+102}`$ keV, i.e. $`kT>5.6`$ keV) and an optical depth ($`\tau <4.2`$). A very hot corona ($`kT>370`$ keV) with a negligible optical depth ($`\tau 0.02`$, $`\tau <2.68`$ at $`90\%`$ confidence level, c.l.) is also (formally) acceptable. However, with a detector sensitive up to 10 keV, it is impossible to test the presence of such a high coronal temperature. The unabsorbed 0.5–10 keV flux for our best fit spectrum amounts to $`2.6\times 10^{10}\mathrm{erg}\mathrm{s}^1\mathrm{cm}^2`$. Source luminosity is $`2.0\times 10^{36}d_8^2\mathrm{erg}\mathrm{s}^1`$. The extrapolated luminosity to the 0.01–100 keV band is a factor $`1.7`$ larger.

The hard component contributes almost the entire (unabsorbed) flux ($`87\%`$). The black body component is detectable only at low energies. Other soft components have been tried with similarly good results. We cannot discriminate between an emission from a black body or from an accretion disk ($`T_d0.2`$ keV, $`R_d10`$ km). A neutron star atmosphere (Gänsicke, Braje & Romani 2002) and a thermal spectrum however provided a worse fit to the data.

### 3.1 Line search

We searched for an iron lines in the pn and MOS data with negative results. Taking the line energy in 6.40–6.97 keV and 0.1 keV (best fit) or null values for the line width, we derive an upper limit on the equivalent width of 18 and 25 eV, respectively.

Following Miller et al. (2002) we searched the RGS spectra for emission or absorption lines. We divided the spectra (rebinned at 20 counts per bin) in slices of 3 Å and inspect them for lines in the energy band 7–20 Å. We fitted the sliced spectra with an absorbed (VARTBABS, Wilms, Allen & McCray 2000) power law model. All fits are acceptable. No emission or absorption lines are found. We directly searched for Ne K edge line at 14.25 Å following the evidence for enhanced absorption in ultra-compact binaries (Schulz et al. 2001; Juett, Psaltis & Chakrabarty 2001). No Ne edge is visible and the fit procedure is not sensible to its presence. We also searched for Ne Ly$`\alpha `$ line at 12.13 Å (in the range 11–13 Å). No line is present and an upper limit of 41 eV can be set ($`95\%`$ c.l.). In case of a narrow line the limit is 9 eV. Detected Ne lines in 4U 1626–67 had equivalent witdhs of 20 Å (Schulz et al. 2001).

## 4 Timing analysis

We used the sas task barycen to correct the pn data. The time resolution is 29.56 $`\mu `$s (Kuster et al. 2003). A significant signal in the power spectrum is readily visible in the data (Campana et al. 2003). However, the power spectrum shows a large number of peaks likely resulting from aliasing of the short orbital period (see Fig. 2). The correct pulse period can only be obtained after correcting for the orbital modulation. We assumed the orbital period derived by Markwardt, Juda & Swank (2003) using RXTE data, i.e. $`P_{\mathrm{orb}}=40.0741`$ min, and searched for the projected orbital radius and time of $`90^{}`$ mean longitude assuming the hypothesis of circular orbit. Stepping on those two parameters and, for each case, folding the events into a 10-bin light curve. The presence of a modulation in the folded light curve was then estimated with a $`\chi ^2`$ test. This is similar to the analysis carried out by Kirsch & Kendziorra (2003). We successfully recover the pulsation with a projected orbital radius $`a_x\mathrm{sin}i=4.7_{0.04}^{+1.2}`$ and an epoch of $`90^{}`$ mean longitude $`T_0=52720.68644(13)`$ ($`90\%`$ c.l., see Fig. 3). Pulsations were detected at 5.2459426(2) ms ($`90\%`$ c.l.).

After correction for the orbital motion and in the full band used for timing analysis (0.3–10 keV) the (background subtracted) pulsed fraction amounts to $`5.8\pm 0.3\%`$ ($`90\%`$ c.l., without correction the pulsed fraction is only $`1.9\%`$). The fit is not completely satisfactory ($`\chi _{\mathrm{red}}^2=4.8`$). An additional sinusoidal component at half of the spin period is highly significant (F-test probability $`3.5\sigma `$ and $`\chi _{\mathrm{red}}^2=1.2`$). The pulsed fraction of XTE J1807–289 is comparable to what observed in SAX J1808.4–3658 ($`4\%`$, Wijnands & van der Klis 1998), XTE 1751–305 ($`4\%`$, Markwardt et al. 2002) and XTE J0929–314 ($`5\%`$, Galloway et al. 2002). The pulsed fraction increases with energy, going from $`4.3\pm 0.7`$ in 0.3–2 keV to $`7.4\pm 0.7`$ in 2–5 keV and to $`7.5\pm 1.4`$ in 5–10 keV. Also for these curves two sinusoidal components provide a much better fit.

The orbital modulation is not directly visible in the data. Folding the barycentered pn light curve at the orbital period reported by Markwardt, Juda & Swank (2003), we can put a $`95\%`$ upper limit on the modulated component of $`1\%`$.

## 5 Conclusions

We report on XMM-Newton observation of the fourth millisecond low mass X–ray pulsar XTE J1807–294. We reveal the source still in a bright state $`30`$ d after its discovery. The source luminosity is $`2\times 10^{36}\mathrm{erg}\mathrm{s}^1`$ at a fiducial distance of 8 kpc.

Our results are very similar to what observed in XTE J1751–305 (Miller et al 2003). The hard spectral component contributes most of the observed flux ($`87\%`$), even though a soft component (a black body) is needed by the data. The equivalent radius of the black body component is small, $`R=(2.2\pm 0.2)d_8^2f^2`$ km, suggesting an emitting hot spot onto the neutron star surface.

No emission or absorption lines are observed in the data. We study in more detail Ne lines following the suggestion by Yungelson, Nelemans & van den Heuvel (2002) and Bildsten (2002) that the companions of faint transient sources are Ne-rich white dwarves. We did not find any evidence of emission lines (41 eV and 9 eV $`95\%`$ upper limits for a broad or narrow line) or absorption edges. Iron lines are not detected either with upper limits of $`25`$ eV.

The hard component (the main contributor to the source flux) can be better described with a Comptonization model rather than a simple power law model. This is at variance with the case of XTE J1751–305. The optical depth and the plasma temperature are highly covariant. In any case a firm upper limit on the optical depth of 4.3 ($`90\%`$ c.l.) can be set. This is in line with the predictions of Titarchuk, Cui & Wood (2002). They suggested that X–ray pulsations can be observed in sources with an optical depth lower than $`4`$ at a plasma temperature of $`20`$ keV (as in our case). For sources with a higher optical depth, e.g. $`\tau =11.7\pm 0.4`$ in the Z source GX 349+2 (Di Salvo et al. 2001) or $`\tau =5.3\pm 0.6`$ in the atoll source 4U 1728–34 (Di Salvo et al. 2000), pulsations would be more easily suppressed by scattering.

We are grateful to the XMM-Newton project scientist, F. Jansen, for granting time to this unsolicited target of opportunity observation. We thank the referee for her/his helpful comments. This work is partially supported through ASI and Co-fin grants.


text:
# Electromagnetic radiation from collisions at almost the speed of light: an extremely relativistic charged particle falling into a Schwarzschild black hole

## I Introduction

In a previous paper we have investigated the gravitational radiation emitted by a relativistic particle infalling into a four-dimensional Schwarzschild black hole cardosolemos0 , and afterwards extended it to the infall in a Kerr metric spacetime cardosolemos1 . The interest of these studies lies in the fact they can describe the gravitational radiation emitted in several collision phenomena, such as the collision between small and massive black holes, between stars and massive black holes, or between cosmic rays and small black holes, to name a few. The possibility of detecting gravitational radiation by the several operating antennae is now real, and a closer understanding of these collisional processes is accordingly important. Of course, to fully understand collisional processes, one has in principle to go beyond the infall of a particle and consider the collision between two black holes with comparable masses, in which case one has to take into account the full non-linearity and strong field regime of Einstein’s equations. Without resorting to numerical computations, but instead using topological arguments, Hawking was able to put upper limits on the gravitational radiation emitted in a black hole - black hole collsion in four dimensions. This calculation was then refined by D’Eath and Payne payne . Their calculation made use of the fact that, if one boosts the Schwarzschild metric to high velocities, then it approaches the Aichelburg-Sexl aichelburg  metric, which is a shock wave spacetime describing the gravitational field of a massless particle. It was found that the total efficiency of the shock would be $`16\%`$ (see Yoshino and Nambu yoshino  for the analysis in higher dimensions). In payne , the inclusion of a second term in the news function brings about a decrease in the efficiency from $`25\%`$ (using only the first term) to $`16\%`$, so one is entitled to ask what degree of confidence do these results and techniques offer. For instance, one might suspect the third term to lower even more the efficiency. One needs to have other means of computing this process, since D’Eath and Payne’s formalism is hard to pursue, due to its complexity. We have recently cardosolemos0 ; cardosolemos1  proposed the “point particle paradigm”, which may well be a good candidate. Handling high energy collisions of point particles with black holes is rather simple, since the point particle may be looked at as a perturbation in the black hole spacetime, and the full machinery to handle black hole perturbations is well developed bhperturbation . The interesting thing about the high energy collision of point particles with black holes is that in the limit that the mass $`\mu `$ of the particle goes to the mass $`M`$ of the black hole, the result for the total radiated gravitational energy agrees well with other predictions and never violates the area theorem cardosolemos0 ; cardosolemos1 . Moreover, one can use this method to treat the collision of rotating holes, which does not seem feasible using previous techniques.

Up to now, most studies have concentrated on the emission of gravitational radiation in the collision process, but if at least one of the particles is charged, electromagnetic radiation is also expected. Astrophysically this is perhaps less relevant for the time being, since the particles one can consider, such as stars or small black holes, are in principle discharged, although it is possible the latter ones can carry some charge (see, e.g., malheiro ). On the other hand, from the point of view of elementary black holes such a study is certainly of interest: elementary black holes and particles can easily carry electromagnetic charge, and the collision between such objects can be put to operate, either by a machine in some far future that can collide particles at center of mass energies of the order of the four-dimensional Planck mass, $`10^{16}`$TeV, or if there is any fundament in the TeV-scale gravity, by the near future generation accelerators such as the LHC. The latter case has been given a lot of consideration. The TeV-scale gravity hamed  requires extra large dimensions (of the order of sub-millimeters or smaller) in order to lower the higher-dimensional Planck mass to energies of the order of TeVs. In these TeV-scale gravity models, the gravitons are free to propagate in the higher dimensional spacetime, whereas, due to experimental constraints, the standard model fields live on a 3-brane, our universe. In this scenario, a particle collider with a center of mass energy of the order of TeVs can copiously produce higher dimensional small spherical type black holes, with radius of the order of fermis or less, in a space with large extra dimensions bhprod , or even black branes if some of the extra dimensions are large and other small (see cav  for a review). Moreover, in the collision process the gravitons escape to the extra dimensions and are therefore much harder, harder than usual, to detect. On the other hand, all the electromagnetic radiation emitted in such a collision can be detected.

Thus it is important to know the electromagnetic spectrum and the quantity of electromagnetic energy radiated in such an encounter. In this paper we extend the previous calculations into to the electromagnetic window and find the electromagnetic radiation emitted by a highly relativistic electrically charged particle infalling into a four-dimensional Schwarzschild black hole. This is a perturbation calculation: we suppose a small charged particle falling into a Schwarzschild black hole. Some former works that have dealt with the phenomenon of electromagnetic radiation from a charged particle falling from infinity into a Schwarzschild black hole are zerilli1 ; ruffini72 ; ruffini74 ; chandra . For example Ruffini ruffini72  first calculated the electromagnetic energy spectra and total radiated energy ruffini72 , for particles with low Lorentz factors, $`\gamma <<1`$. Here we go into high $`\gamma `$s. We use the point particle approximation and consider a charged point particle colliding head-on at high velocity with a Schwarzschild black hole. We show that the spectra is flat, and we also compare the total electromagnetic and gravitational energies emitted, and find that the former is supressed in relation to the latter for very high energies. The numerical results extracted by us are in very good agreement with Ruffini’s results ruffini72 . We shall also see that there is a classical calculation for this process that agrees extremely well with our numerical results.

Two comments are in order: (i) Our calculation is for the collision of a small particle with a black hole in a four-dimensional world. Therefore, in principle, it could apply to the astrophysical world in the case charged stars and charged black holes are out there, and to a very high energy collision experiment in a four-dimensional world (without extra dimensions) and with a fundamental Planck mass of $`10^{16}`$TeV. In this latter scenario the calculation is to be used for the moments just after the black hole formation, when the collision of charged debris with the newly formed black hole is certainly expected. Furthermore, the calculation can be extended to a black hole - black hole collision as we have already argued cardosolemos0 ; cardosolemos1 , and, in addition, it can give clues, although it does not aplly directly to the usual collision process in a collider, i.e., the collision between one non-charged particle with strong gravitational field (not necessarily a black hole) and one charged particle. These results may, however, serve as a model to the electromagnetic radiation emitted in the initial phase of a newly formed black hole, the stage in which the black hole sheds its hair by emitting gauge radiation, such as electromagnetic radiation. (ii) Our calculation is four-dimensional. It does not apply to Tev-scale gravity black holes, since these inhabit a world of six to eleven dimensions, with some of the dimensions being large, others perhaps being small. However, although there are certainly some differences, qualitatively our results should hold when extrapolated with some care to higher dimensions.

## II Basic Formalism

In this study, we assume the charged particle and the emitted radiation to be a small perturbation on the Schwarzschild spacetime, whose line element is given by

$`ds^2=e^\nu dt^2+e^\lambda dr^2+r^2(d\theta ^2+\mathrm{sin}^2\theta d\phi ^2),`$ (1)
where $`e^\nu =e^\lambda =12M/r`$, and $`M`$ stands for the mass of the background spacetime. Because a complete prescription to deal with this problem was given by Zerilli zerilli1 , here we briefly show the essential parts of the formalism. In virtue of the static property and spherical symmetry of the background spacetime, perturbations can be decomposed by using vector harmonics for the angular variables $`\theta `$ and $`\phi `$, and the Fourier components for the time variable $`t`$. After accomplishing a separation of variables, perturbations are therefore characterized by the harmonic indices $`l`$ and $`m`$, parity, and frequency $`\omega `$. For the present case, where a particle falls straight into a black hole along the $`z`$-axis or the $`\theta =0`$ line, no axial parity perturbations are excited due to the symmetry of the motion. Thus, only the polar parity perturbations are considered in this study. According to Zerilli zerilli1 , our master equation for determining the electromagnetic radiation is given by a single wave equation,

$$\frac{d^2\stackrel{~}{f}_{lm}(\omega ,r)}{dr_{}^2}+\left[\omega ^2e^\nu \frac{l(l+1)}{r^2}\right]\stackrel{~}{f}_{lm}(\omega ,r)=e^\nu \stackrel{~}{S}_{lm},$$
(2)
where $`r_{}`$ is the so-called tortoise coordinate, defined by $`dr/dr_{}=e^{(\nu \lambda )/2}`$, or $`r_{}=r+2M\mathrm{log}(r/2M1)`$ in our case, and $`\stackrel{~}{S}_{lm}`$ is the source term determined by the charge and motion of the particle. In the case of the radially falling particle, $`\stackrel{~}{S}_{lm}`$ is generally given by

$$\stackrel{~}{S}_{lm}=2q\sqrt{l+\frac{1}{2}}\frac{e^{i\omega T(r)}}{r^2}\delta _{m0},$$
(3)
where $`\delta _{ij}`$ means the Kronecker delta, and $`q`$ denotes the charge of the particle. Here, the function $`T(r)`$ is the coordinate time of the particle parameterized by the radial position $`r`$ of the particle. Since the charged particle is treated as a perturbation of the background spacetime, the particle traces a radial geodesic of the Schwarzschild geometry. Thus, the function $`T(r)`$ is determined by (see e.g., Chandrasekhar chandra )

$$\frac{dT(r)}{dr}=\frac{e^{(\lambda \nu )/2}}{\sqrt{1\gamma ^2e^\nu }},$$
(4)
where $`\gamma `$ is an integral of motion, given by $`\gamma =(1v_{\mathrm{}}^2)^{1/2}`$, where $`v_{\mathrm{}}`$ is the radial velocity of the particle at spatial infinity.

Once solutions $`\stackrel{~}{f}_{lm}(\omega ,r)`$ of equation (2) are calculated, time dependent functions $`f_{lm}(t,r)`$ can be obtained through the inverse Fourier transformation,

$$f_{lm}(t,r)=\frac{1}{\sqrt{2\pi }}_{\mathrm{}}^{\mathrm{}}e^{i\omega t}\stackrel{~}{f}_{lm}(\omega ,r)𝑑\omega .$$
(5)
From the functions $`f_{lm}(t,r)`$, all the components of the field strength $`F_{\alpha \beta }`$ can be derived except at the position of the charged particle. For example, the $`r`$$`\theta `$ component $`F_{r\theta }`$ is given by

$$F_{r\theta }(t,r,\theta ,\phi )=e^\nu \underset{l,m}{}f_{lm}(t,r)\frac{Y_{lm}}{\theta }(\theta ,\phi ).$$
(6)
In order to obtain a unique solution of equation (2), two boundary conditions must be specified, and physically acceptable conditions are purely outgoing wave at spatial infinity and purely incoming wave at the horizon, which are due to equation (5) given by

$$\stackrel{~}{f}_{lm}(r)\{\begin{array}{cc}A_{lm}^{\mathrm{in}}(\omega )e^{i\omega r_{}}\hfill & \mathrm{as}r_{}\mathrm{},\hfill \\ A_{lm}^{\mathrm{out}}(\omega )e^{i\omega r_{}}\hfill & \mathrm{as}r_{}\mathrm{}.\hfill \end{array}$$
(7)
where $`A_{lm}^{\mathrm{out}}(\omega )`$ and $`A_{lm}^{\mathrm{in}}(\omega )`$ are not dependent on $`r`$. Since we are interested only in solutions at spatial infinity, what we have to do is to obtain $`A_{lm}^{\mathrm{out}}(\omega )`$’s as functions of $`\omega `$ for our purpose of this study. Adapting a standard Green’s function technique, we can write $`A_{lm}^{\mathrm{out}}`$ as the integral, given by

$$A_{lm}^{\mathrm{out}}(\omega )=\frac{1}{2i\omega C_l(\omega )}_{2M}^{\mathrm{}}f_{\mathrm{L},l}\stackrel{~}{S}_{lm}𝑑r,$$
(8)
where $`f_{\mathrm{L},l}`$ is a homogeneous solution of equation (2) satisfying a boundary condition given by

$$f_{\mathrm{L},l}\{\begin{array}{cc}e^{i\omega r_{}}\hfill & \mathrm{as}r_{}\mathrm{},\hfill \\ B_l(\omega )e^{i\omega r_{}}+C_l(\omega )e^{i\omega r_{}}\hfill & \mathrm{as}r_{}\mathrm{}.\hfill \end{array}$$
(9)
The energy spectrum at spatial infinity is given by

$$\frac{dE}{d\omega }=\underset{l}{}\frac{dE_l}{d\omega }=\underset{l}{}\frac{l(l+1)}{2\pi }|A_{l0}^{\mathrm{out}}(\omega )|^2\mathrm{for}\omega 0,$$
(10)
In order to obtain numerical values of $`A_{lm}^{\mathrm{out}}(\omega )`$, we start the integration of $`f_{L,l}`$ and $`_{2M}^rf_{\mathrm{L},l}\stackrel{~}{S}_{lm}𝑑r^{}`$ at $`r=2M(1+10^6)`$ by using a Runge-Kutta method. Those functions are then integrated out to large values of $`r`$. The integration is stopped if the absolute value of $`_{2M}^rf_{\mathrm{L},l}\stackrel{~}{S}_{lm}𝑑r^{}`$ converges within the required numerical accuracy, and we simultaneously match $`f_{L,l}`$ with an asymptotic solution satisfying condition (9) at spatial infinity, which is given in chandra , to obtain a value of $`C_l(\omega )`$.

## III Numerical Results

Following the numerical procedure just outlined, we have computed the energy spectrum for several values of $`\gamma `$. In Figure 1 we show a typical result of the energy spectrum, here for $`\gamma =2`$ and for the four lowest radiatable multipoles. In Fig. 2 we show similar results but for a particle with $`\gamma 10`$. A general feature of the energy spectrum for high energy collisions is that it is flat up to some critical frequency, after which it rapidly (exponentially) decreases to zero. This was also verified for the gravitational energy spectrum resulting from the high-energy collision of a black hole with a point particle cardosolemos0 ; cardosolemos1 . This critical frequency is given, in a very good approximation, by the frequency of the fundamental quasinormal mode associated with the corresponding harmonic index $`l`$. The corresponding quasinormal frequencies are also indicated by a filled circle in Figures 1 and 2. The total energy $`\mathrm{\Delta }E`$ radiated away is given by

$$\mathrm{\Delta }E=\underset{l=1}{\overset{\mathrm{}}{}}\mathrm{\Delta }E_l,$$
(11)
where $`\mathrm{\Delta }E_l`$ is the total radiated energy associated with $`l`$-pole electromagnetic radiation, which is given by the area under the curve of the energy spectrum $`dE_l/d\omega `$ (see Figure 1). For the low $`\gamma `$ cases, the convergence of the infinite series (11) is relatively good in the sense that a finite summation of several $`\mathrm{\Delta }E_l`$’s is sufficient to obtain a value of $`\mathrm{\Delta }E`$ with sufficient accuracy. However, the convergence of the total radiated energy gets worse with increasing $`\gamma `$. This contrasts with the analogous calculation of the gravitational energy radiated in high energy collisions cardosolemos1  where a rescaling was possible to show that $`\mathrm{\Delta }E\gamma ^2`$. Here we shall see in the next section that the difficulty for the electromagnetic case lies in the fact that $`\mathrm{\Delta }E\mathrm{log}\gamma `$, for high $`\gamma `$’s. For a very fast particle, $`\gamma \mathrm{}`$, it is quite difficult to numerically obtain the total radiated energy, since all the $`l`$-pole energy spectra $`dE_l/d\omega `$ make a substantial contribution to the total radiated energy. In this study, we calculated the total energy for several values of $`\gamma `$ with $`\gamma 10`$. Some of these values are shown in Table 1.

Previous works concerned with the electromagnetic radiation from a charged particle falling from infinity into a Schwarzschild black hole include for example Ruffini ruffini72 . Ruffini first calculated the electromagnetic energy spectra and total radiated energy ruffini72 , for low $`\gamma `$ particles. The numerical results extracted by us are in very good agreement with Ruffini’s results ruffini72 .

## IV The Classical Calculation

There is a classical calculation in electromagnetism that suggests itself as a model to use in the high energy collision of charged particles: the radiation emitted when a charge, with constant velocity $`v_{\mathrm{}}`$ is suddenly decelerated to $`v=0`$. The deceleration is idealized as taking zero seconds. This model calculation has been applied with great success to, for example, beta decay chang49 ; jackson . The result for the energy spectrum per solid angle is chang49 ; jackson

$$\left(\frac{d^2E}{d\omega d\mathrm{\Omega }}\right)_{\mathrm{class}}=\frac{q^2v_{\mathrm{}}^2}{4\pi ^2}\frac{\mathrm{sin}^2\theta }{(1v_{\mathrm{}}\mathrm{cos}\theta )^2},$$
(12)
or, integrating over solid angle,

$$\left(\frac{dE}{d\omega }\right)_{\mathrm{class}}=\frac{q^2}{\pi }\left[\frac{1}{v_{\mathrm{}}}\mathrm{log}\left(\frac{1+v_{\mathrm{}}}{1v_{\mathrm{}}}\right)2\right].$$
(13)
This formula has been tested with great accuracy experimentally. To get the total energy, one has to integrate (13) over frequencies, and a naive procedure would lead to infinities in the total energy. To obtain reasonable results, one has to impose a cutoff $`\omega _c`$ in the frequency depending on the particular problem under consideration, and one obtains

$$\mathrm{\Delta }E_{\mathrm{class}}=\frac{q^2}{\pi }\left[\frac{1}{v_{\mathrm{}}}\mathrm{log}\left(\frac{1+v_{\mathrm{}}}{1v_{\mathrm{}}}\right)2\right]\omega _c.$$
(14)
How well do these classical formulas fit into our numerical results? Very well indeed. To allow a more direct comparison, we shall first decompose the energy spectrum (13) into spherical harmonics, i.e., into the energy spectrum associated with each $`l`$-pole radiation as follows:

$`\left({\displaystyle \frac{dE}{d\omega }}\right)_{\mathrm{class}}={\displaystyle \underset{l}{}}\left({\displaystyle \frac{dE_l}{d\omega }}\right)_{\mathrm{class}}`$
$`={\displaystyle \underset{l}{}}{\displaystyle \frac{q^2}{4\pi ^2l(l+1)}}\left|{\displaystyle \frac{v_{\mathrm{}}\mathrm{sin}\theta }{1v_{\mathrm{}}\mathrm{cos}\theta }\frac{Y_{l0}}{\theta }𝑑\mathrm{\Omega }}\right|^2.`$ (15)
When we consider the limit of $`v_{\mathrm{}}1`$, as shown in equation (13), the total energy spectrum $`\left(\frac{dE}{d\omega }\right)_{\mathrm{class}}`$ diverges. However, the energy spectrum due to $`l`$-pole radiation $`\left(\frac{dE_l}{d\omega }\right)_{\mathrm{class}}`$ converges to a finite value, given by

$$\left(\frac{dE_l}{d\omega }\right)_{\mathrm{class}}=q^2\frac{2l+1}{\pi l(l+1)}.$$
(16)
This yields for example $`\frac{dE}{d\omega }_{l=1}=0.477q^2`$ and $`\frac{dE}{d\omega }_{l=2}=0.265q^2`$ as $`\gamma \mathrm{}`$. This may be compared with the numerical values in Fig. 2. The numerical agreement is excellent not only for these lower multipoles, but also for higher ones, and one has therefore established numerically that for high energy collisions the classical result (13) for the energy spectrum is valid. Probably more interesting is the total energy radiated. Since one already knows that the classical result for $`\frac{dE}{d\omega }`$ is correct, let us now try to predict the cutoff frequency by identifying $`\mathrm{\Delta }E`$ in equation (11) with (14). Making use of this definition of $`\omega _c`$ and numerical values of $`\mathrm{\Delta }E`$, we evaluate the cutoff frequencies $`\omega _c`$ in (14) for several values of $`\gamma `$, and list the values of $`\omega _c`$ in Table II. In Figure 3, we also show the cutoff frequencies as a function of $`\gamma `$. It is found in this figure that for relatively high vales of $`\gamma `$, $`\omega _c`$ increase almost linearly with $`\gamma `$, and a good fit to our numerical values is given by

$$M\omega _c=0.224+0.0598\gamma .$$
(17)
In Figure 3, we have shown this linear function, and we can confirm that this linear function is in good agreement with our numerical values for high $`\gamma `$. This means that a good approximation to the total electromagnetic energy radiated is

$$\mathrm{\Delta }E=\frac{q^2}{M\pi }\left[\frac{1}{v_{\mathrm{}}}\mathrm{log}\left(\frac{1+v_{\mathrm{}}}{1v_{\mathrm{}}}\right)2\right](0.224+0.0598\gamma ),$$
(18)
or considering that $`v_{\mathrm{}}1`$ we have also

$$\mathrm{\Delta }E=\frac{2q^2}{M\pi }\left[\mathrm{log}2\gamma 1\right](0.224+0.0598\gamma ),$$
(19)
## V Comparing the Amount of Gravitational and Electromagnetic Energy radiated

Previous studies on the high energy collision of point particles with black holes led cardosolemos0 ; cardosolemos1  to the following conclusion: if a high energy point particle of mass $`\mu `$ collides head on with a mass $`M`$ non-rotating black hole, the total amount of gravitational energy radiated is

$$\mathrm{\Delta }E_g=0.26\frac{\mu ^2\gamma ^2}{M},$$
(20)
where $`\gamma `$ is again the Lorentz factor for the point particle. Furthermore, it was also found that the spectra is flat and also well described by a classical calculation (here classical means again working on a flat background). Now, for high $`\gamma `$’s (19) is

$$\mathrm{\Delta }E\frac{0.038q^2}{M}\gamma \mathrm{log}2\gamma ,\gamma \mathrm{}.$$
(21)
So, we get for the ratio electromagnetic energy over gravitational energy,

$$\frac{\mathrm{\Delta }E}{\mathrm{\Delta }E_g}0.146(\frac{q}{\mu })^2\frac{\mathrm{log}2\gamma }{\gamma },\gamma \mathrm{}.$$
(22)
This means that electromagnetic energy is supressed in relation to gravitational energy, for very high energies.

## VI Conclusions

We have computed the electromagnetic spectrum and total energy radiated during the high energy collision of a charged point particle with a Schwarzschild black hole. Our results show that the classical “instantaneous collision” calculation gives very good results, in accordance with the full numerical ones. We have dealt only with zero impact (i.e., head-on) collisions, but these results are easily generalized to non head-on collisions. For example, one expects that the classical results chang49 ; jackson  hold for such cases, and therefore that the total energy decreases as one increases the impact parameter. We stress that the results presented here are valid for any particle charge $`q`$, as long as the total effective stress energy of the particle is small compared to the total energy content of the black hole. Previous works cardosolemos0 ; cardosolemos1  have shown that the high energy collision between two black holes of equal mass may be well studied through the collision of a point particle with a black hole, and then taking the limit of equal mass, although this is formally not allowed. It is our belief that the same may be done here. One needs however some other method to attack this problem, to confirm or disprove this claim. The investigation carried here can also be carried over to higher dimensions, a case which is of more direct interest to TeV-scale gravity scenarios. Since Standard Model fields inhabit a four-dimensional brane, what one would need to generalize this construction would be to modify the induced four-dimensional metric (1), as was done for example in kanti . The generalization is straightforward.

Acknowledgements - We thank George Smoot for asking, during the talk of one of us (VC) in the Fourth International Workshop on New Worlds in Astroparticle Physics (University of Algarve, Faro, Portugal), how the gravitational radiation emitted in the process compares with the electromagnetic radiation, from which this work resulted. This work was partially funded by Fundação para a Ciência e Tecnologia (FCT) – Portugal through project CERN POCTI/FNU/49521/2002. VC acknowledges finantial support from FCT through PRAXIS XXI programme. JPSL acknowledges finantial support from ICCTI/FCT and thanks Observatório Nacional do Rio de Janeiro for hospitality. SY acknowledges finantial support from FCT through project SAPIENS 36280/99.


text:
# Observational Critical QCD

## Abstract

A detailed study of correlated scalars, produced in collisions of nuclei and associated with the $`\sigma `$-field fluctuations, $`(\delta \sigma )^2=<\sigma ^2>`$, at the QCD critical point (critical fluctuations), is performed on the basis of a critical event generator (Critical Monte-Carlo) developed in our previous work. The aim of this analysis is to reveal suitable observables of critical QCD in the multiparticle environment of simulated events and select appropriate signatures of the critical point, associated with new and strong effects in nuclear collisions.

preprint: UA-NPPS/04/2003

The existence of a critical point in the phase diagram of QCD, for nonzero baryonic density, is of fundamental significance for our understanding of strong interactions and so its experimental verification has become an issue of high priority WRB . For this purpose an extensive programme of event-by-event searches for critical fluctuations in the pion sector has already been started in experiments with heavy ions from SPS to RHIC energies NA49 . In ACDKK  we have emphasized, however, that in order to reveal critical density fluctuations in multiparticle environment, one has to reconstruct the $`\sigma `$-sector from pion pairs ($`\sigma \pi ^+\pi ^{}`$) and study correlations of sigmas as a function of the invariant mass near the two -pion threshold. In fact, the QCD critical point, if it exists, communicates with a zero mass scalar field ($`\sigma `$-field) which at lower temperatures ($`T<T_c`$) may reach the two-pion threshold and decay in very short time scales owing to the fact that its coupling to the two-pion system is very strong. Obviously, the fundamental, underlying pattern of $`\sigma `$-field fluctuations, built-up near the critical point by the universal critical exponents of QCD ACDKK , is phenomenologically within reach if and only if the study of correlated sigmas, reconstructed near the two-pion threshold, becomes feasible. In this Letter we perform a detailed feasibility study of the observables related to the detection of the QCD critical point in nuclear collisions. In order to proceed we summarize, first, the principles on which the behaviour of a critical system of sigmas is based ACDKK ; WRB .

(a) The geometrical structure of the critical system in transverse space (after integrating in rapidity) consists of $`\sigma `$-clusters with a fractal dimension $`d_F=\frac{2(\delta 1)}{\delta +1}`$ leading to a power law, $`<\sigma ^2>|\stackrel{}{x}|^{\frac{4}{\delta +1}}`$, for the $`\sigma `$-field fluctuations, within each cluster ($`\delta `$ : isotherm critical exponent)

(b) In transverse momentum space the $`\sigma `$-fluctuations obey a power law $`<\sigma ^2>|\stackrel{}{p}_{}|^{\frac{2(\delta 1)}{\delta +1}}`$ leading to observable intermittent behaviour of factorial moments: $`F_2(M)(M^2)^{\frac{\delta 1}{\delta +1}}`$ where $`M^2`$ is the number of $`2D`$ cells BialPes .

(c) The density fluctuations of pion pairs $`(\pi ^+\pi ^{})`$ with invariant mass close to two-pion threshold $`(2m_\pi )`$ simulate to a good approximation the sigma-field fluctuations, $`(\delta \sigma )^2<\sigma ^2>`$, at the critical point, under the assumption that the sigma mass reaches the two-pion threshold $`(m_\sigma \stackrel{>}{}2m_\pi )`$ at a freeze-out temperature close to the critical value.

(d) The QCD critical point belongs to the universality class of the $`3D`$ Ising system in which $`\delta 5`$.

On the basis of these principles and the fact that critical clusters in the above universality class interact weakly ACDPRE , one may construct a generator of critical events (Critical Monte-Carlo: CMC) containing only sigmas, correlated according to the above prescription ACDKK . The numerical experiment is completed by a self-consistent treatment of $`\sigma `$-reconstruction in the simulated events (after the process of decay of sigmas to pions) and a comparative study between the associated critical correlations-fluctuations and the corresponding behaviour of a conventional Monte-Carlo (HIJING).

The input parameters of the simulation are the size of the system in rapidity ($`\mathrm{\Delta }`$) and transverse space ($`R_{}`$) as well as the proper time scale ($`\tau _c`$) characteristic for the formation of the critical system. The CMC generator can be used to obtain a large set of critical events. The analysis of these events could serve as a guide for the search of the critical point in the data of real $`A+A`$ processes in SPS and RHIC experiments. This will be shown in a more transparent way in what follows.

We will investigate two sets of critical events. One set corresponding to a rather small system with mean charged pion multiplicity per event $`<n_{ch}>40`$ and another set describing a larger system ($`<n_{ch}>120`$). Each data set consists of $`10^5`$ CMC generated events. First we will consider the properties of the small system. For the values of the input parameters in this case we take $`\mathrm{\Delta }=6`$ adapted to the SPS energies and $`R_{}=15fm`$. In order to reveal underlying critical fluctuations, at the level of observation, one has first to perform factorial moment analysis in small cells of the phase space BialPes . We have chosen transverse momenta for this analysis, the reason being that density fluctuations in this plane are practically independent of the geometry of the critical system ACDKK . Applying factorial moment analysis to the transverse momenta $`(p_x,p_y)`$ of the negative pions in the sample of the $`10^5`$ critical events we obtain for the second moment the behaviour shown in Fig. 1. We observe a weak intermittency effect ($`F_2M^{2s_2}`$ with $`s_20.07`$) much smaller than the expected to occur in the critical system ($`s_2=2/3`$). For comparison we show in the same plot the corresponding factorial moment for the sigmas, before their decay, which follows closely the theoretical prediction: $`s_20.63`$. The pions analysed in this figure originate from the decay of the CMC sigmas using pure kinematics and taking into account the isospin properties of the decay amplitude. The mass of the decaying sigmas is treated as an almost uniformly distributed random variable ACDKK .

The reason for the suppression of fluctuations in the pionic sector is the kinematical distortion of the self-similar pattern formed in the sigma sector due to the sigma-decay. The strength of this distortion increases with momentum transfer: $`Q=\sqrt{m_\sigma ^24m_\pi ^2}`$ and becomes negligible near the two-pion threshold ($`m_\sigma 2m_\pi `$). Here $`m_\sigma ^2=(p_{\pi ^+}+p_\pi ^{})^2`$ where $`p_{\pi ^\pm }`$ are the four momenta of the produced charged pions.

In this treatment we consider only the visible mode of the sigma-decay to charged pions and the reconstruction of the sigmas, based on this sector, is simply a search for $`(\pi ^+,\pi ^{})`$ pairs with invariant mass close to $`2m_\pi `$. The aim of this procedure is to establish a self-consistent mechanism for producing critical sigmas either directly from the event generator (CMC) or through reconstruction from the final states $`(\pi ^+,\pi ^{})`$ in the same data set of simulated events. The significance of this self-consistent reconstruction approach in a search for critical $`\sigma `$-fluctuations in real events becomes self-evident. In practice this can be done by looking, in each recorded event, for $`(\pi ^+,\pi ^{})`$ pairs fullfiling the criterion A:

$$A=\{(\pi ^+,\pi ^{})|4m_\pi ^2(p_{\pi ^+}+p_\pi ^{})^2(2m_\pi +ϵ)^2\}$$
(1)
with $`\frac{ϵ}{2m_\pi }1`$. The momentum of the sigma-particle is then obtained as: $`\stackrel{}{p}_\sigma =\stackrel{}{p}_{\pi ^+}+\stackrel{}{p}_\pi ^{}`$. In general the pairing will not be one-to-one, allowing for different possibilities to pair the charged pions in order to get sigmas. It turns out that the full pairing forming all possible pairs $`(\pi ^+,\pi ^{})`$ fullfiling (1) for a given $`\pi ^+`$, produces more accurately the geometric properties of the critical sigma sector in the limit $`ϵ0`$. Due to finite statistics this limit can only be taken using data sets consisting of events with reconstructed sigmas, obtained for a sequence of $`ϵ`$-values approaching to zero, and then extrapolate smoothly to $`ϵ=0`$. Thus we have performed factorial moment analysis in transverse momentum space using values of $`ϵ`$ in the region $`815MeV`$. We have calculated the second moment $`F_2`$ for different sets of events with reconstructed sigmas, each obtained using a fixed $`ϵ`$. It exists a wide window in the resolution scale $`M`$ for which $`F_2(M^2)^{s_2(ϵ)}`$. Using linear extrapolation the full pairing leads to the value $`s_2(0)=0.70\pm 0.03`$ very close to the theoretically predicted $`s_2(0)=\frac{2}{3}`$. In an alternative rule for pairing, choosing for each $`\pi ^+`$ randomly the corresponding $`\pi ^{}`$ partner, from all possible ones fullfiling (1), and repeating the same moment analysis we find $`s_2(0)0.44\pm 0.01`$ deviating by $`35\%`$ from the theoretically expected value. Thus the full pairing reproduces succesfully the geometrical properties of the critical system due to the fact that in this way one reconstructs all critical sigmas even if it is impossible to identify them.

The above algorithm besides the reconstruction of the real critical sigmas also creates a number of fake sigmas. Some of these are generated by combining pions belonging to different critical sigmas and contribute significantly to the strong intermittency effect occuring in the reconstructed isoscalar sector. Therefore the intermittency effect in the factorial moment analysis is a necessary but not sufficient condition to identify critical, strongly correlated isoscalars. There is a need to substract the effect of the fake sigmas in order to reveal the correlated critical isoscalar sector. To achieve this goal one has to calculate the pair-correlation density of the sigmas defined by:

$$C(1,2)=<\rho _2(1,2)><\rho _2^{(m)}(1,2)>$$
(2)
where $`<\rho _2>`$ is the two-particle density averaged over the events. The arguments $`1,2`$ refer to the momentum coordinates of two sigma-particles. The second term is the two-particle density estimated using mixed events WolfDr . In fact for a given set of events one can construct a set of mixed events by shuffling the momentum coordinates of the pions in the initial data set. As a consequence the (possible) correlations in the original data, including the formation of sigmas, are destroyed in the mixed set. In general the two-particle density is given by the number of particle pairs per event and per phase space cell of volume $`d\mathrm{\Omega }`$. Therefore eq.(2) can be written as:

$$C(1,2)=\frac{1}{N_{ev}}\underset{i=1}{\overset{N_{ev}}{}}\frac{dn_{12,i}}{d\mathrm{\Omega }}W\frac{1}{N_{ev}^{(m)}}\underset{i=1}{\overset{N_{ev}^{(m)}}{}}\frac{dn_{12,i}^{(m)}}{d\mathrm{\Omega }}$$
(3)
where $`i`$ labels the events, $`N_{ev}`$ is the total number of CMC events fullfiling the condition (1) while $`N_{ev}^{(m)}`$ is the total number of mixed events, obtained from the shuffling of the tracks in the CMC events, fullfiling (1). The coefficient $`W`$ is the ratio of numbers of events: $`W=\frac{N_{ev}^{(m)}}{N_{ev}}`$. The correlator $`G(M)`$ is defined as the integral of $`C(1,2)`$ over a finite phase space cell with volume $`\delta \mathrm{\Omega }=\frac{\mathrm{\Omega }}{M^D}`$, where $`\mathrm{\Omega }`$ is the total available phase space in $`D`$-dimensions (here $`D=2`$). Performing the integration and normalizing properly we can express the correlator in terms of the second order factorial moments WolfDr . Taking into account the fact that the resulting distributions of the sigmas depend on the parameter $`ϵ`$ through (1) we finally get:

$$G(M,ϵ)=F_2(M,ϵ)W(ϵ)\left(\frac{<n_\sigma ^{(m)}>}{<n_\sigma >}\right)^2F_2^{(m)}(M,ϵ)$$
(4)
where $`<n_\sigma >`$ is the mean multiplicity of sigmas within a phase space cell for the CMC events while $`<n_\sigma ^{(m)}>`$ is the corresponding quantity for the mixed events. These two multiplicities are in fact almost equal and therefore their ratio appearing in eq.(4) is very close to one for all $`ϵ`$. The weight factor $`W`$ has a transparent physical interpretation: it describes the relative weight of the fake sigmas in the reconstructed set. If $`W(ϵ)1`$ the fake sigmas dominate while for $`W(ϵ)0`$ they are suppressed relative to the real sigmas. It turns out that the dependence of $`W`$ on $`ϵ`$ is crucial for the determination of the critical isoscalar sector. Accordingly, if $`W0`$ the correlator follows the properties of the second moment $`F_2`$ and the effect of the mixed events is marginal, while for $`W1`$ the correlator becomes negligible and no critical effect can be identified. At threshold ($`ϵ=0`$) we expect to fully recover the critical isoscalar sector therefore $`W(0)`$ must vanish if real critical sigmas exist in the data set of the reconstructed sigmas. This should be the case when analyzing the CMC events. In fact, it turns out that in order to avoid fake sigmas at threshold, in the reconstruction procedure, events with sufficient multiplicity of pion pairs $`(n_\sigma )`$ must be considered in the analysis $`(n_\sigma >2)`$. In addition the correlator $`G(M,ϵ)`$ should behave like $`M^{2\varphi _2(ϵ)}`$ if self-similar patterns of correlated sigmas are present in the corresponding phase space. The exponent $`\varphi _2(ϵ)`$ is a new intermittency exponent having the property $`\varphi _2(0)=s_2(0)`$ as $`G(M,0)F_2(M,0)`$. Finally, the strength of the correlations in the sigma sector is given by the magnitude of $`G(M,ϵ)`$. Therefore in terms of the correlator we have three very strict requirements in order to reveal the critical isoscalars in a data set of charged pions:

* Dominance of the real sigmas with respect to the fake sigmas when approaching the two-pion threshold ($`lim_{ϵ0}W(ϵ)=0`$)
* Large positive values of the correlator at threshold ($`G(M,0)1`$) at high resolution scales.
* Correct power-law dependence of the correlator on the resolution scale $`M`$, at threshold: $`G(M,0)M^{1.34}`$

In a conventional hadronic system with no critical isoscalar sector at least one of the above properties will not be met.

To test the above ideas we applied the sigma-reconstruction algorithm to $`10^5`$ CMC events describing a small system as well as to $`33000`$ HIJING events Hijing  for the $`C+C`$ system at the SPS. The mean charged pion multiplicity per event in these two systems is practically the same. The CMC data set is expected to possess a recoverable critical isoscalar sector. On the contrary the HIJING simulator describes, in the low $`p_{}`$ region, conventional hadronic dynamics and therefore no critical isoscalars are expected to occur in the corresponding data set Hijing .

We have first calculated the correlator for these two systems choosing $`n_\sigma 5`$ for the multiplicity of pion pairs. The results are given in Fig. 2 for $`ϵ`$-values in the region of $`815MeV`$, suggested by finite statistics.

As we can see the CMC-correlator behaves very differently from the HIJING-correlator for all values of $`ϵ`$. In Fig. 3 we also show the weight factor $`W(ϵ)`$ for both systems. It is obvious that in the CMC case $`W`$ approaches zero for $`ϵ0`$ while in HIJING $`W(ϵ)`$ is of order of one independently of the $`ϵ`$-value. The interpretation of this very important result is that in the case of HIJING, i.e. in a conventional hadronic system, the reconstructed sigmas are fake for every value of $`ϵ`$.

Finally in Fig. 4 we plot $`\varphi _2(ϵ)`$ for the CMC case. The linear extrapolation leads to $`\varphi _2(0)=0.69\pm 0.04`$ in very good agreement with the theoretically expected value of $`2/3`$. According to the analysis above we have recovered in a self-consistent manner the critical isoscalar sector in the CMC data set as a result of $`\sigma `$-reconstruction. On the contrary the HIJING events do not contain real critical sigmas, as expected.

We also examined the dependence of our results on the charged pion multiplicity. We have produced CMC events for a larger system with size in rapidity $`\mathrm{\Delta }=11`$ and $`R_{}=25fm`$. In this case the total charged pion multiplicity is $`120`$ per event coming from the decay of sigmas with mean multiplicity $`96`$. We found that $`W(0)`$ is very close to zero (linear extrapolation: $`W(0)=0.03`$) and $`\varphi _2(0)=0.664\pm 0.02`$ in very good agreement with the estimations in the smaller critical system (Fig. 4). Thus the results concerning the description of the critical isoscalar sector are not sensitive to the size of the system, as expected. What is changing is the correlation length which for a second order phase transition is proportional to the size of the system. In fact for the two systems considered here the critical correlation length is $`6fm`$ (small system) and $`10fm`$ (large system) respectively ACDPRE .

In order to support further our claim that in the CMC description the sector of $`\pi ^+\pi ^{}`$ pairs near two-pion threshold, has a dinstict dynamical behaviour, associated with $`\sigma `$-field fluctuations at the critical point, we have examined the behaviour of $`\pi ^+\pi ^+`$ pairs in the same region $`(ϵ\stackrel{<}{}4MeV)`$. As expected no effect is found, the pairing $`\pi ^+\pi ^+`$ has no dynamical origin $`(W1)`$ and the correlator $`G(M,ϵ)`$ of these pairs is strongly suppressed. We illustrate this in Fig. 5 where we compare, for $`ϵ=4MeV`$, the correlators of sigmas ($`\pi ^+\pi ^{}`$ pairs) and $`\pi ^+\pi ^+`$ pairs for the small system $`(<n_{ch}>40)`$. It must be noted that for this value of $`ϵ`$ the Coulomb correlations are negligible ($`Q^2=(p_{\pi ^+}p_\pi ^{})^250MeV`$ Seyboth ). Similar results are found if the system is larger $`(<n_{ch}>120)`$.

In conclusion, we have shown that a set of well prescribed observables (factorial moments, correlators, statistical weight of fake sigmas, intermittency exponents) associated with the existence of a critical point in quark matter, can be established in nuclear collisions. These observables belong to the reconstructed $`\sigma `$-sector describing massive scalars ($`\pi ^+\pi ^{}`$ pairs) near the two-pion threshold and their behaviour reveals strong critical effects suggested by $`\sigma `$-field fluctuations near the critical point. We claim that a search for such a critical behaviour in heavy ion experiments is feasible within the framework of a reconstruction procedure of massive scalars, discussed in this work. Since the critical effects in this sector are strong and their pattern remains robust for systems of different size, our proposal is to study, using the above observables, different processes at the SPS and RHIC with the aim to scan a substantial area of the phase diagram, in a systematic search for the QCD critical point in collisions of nuclei.

###### Acknowledgements.

We thank the NA49 Collaboration and especially K. Perl and R. Korus for their help in the treatment of the HIJING data for the $`C+C`$ system. One of us (N.G. A.) wishes to thank Larry McLerran, Dimtri Kharzeev and Jacek Wosiek for their valuable comments. This work was partly supported by the Research Committee of the University of Athens.


text:
# Effect of Semicore Orbitals on the Electronic Band Gaps of Si, Ge, and GaAs within the GW Approximation

## I Introduction

Since the early applications of the GW method to real materials (see Refs. 1-3 and References therein), the pseudopotential plane-wave approach has been the method of choice due to its accuracy and technical simplicity. Recent advances in LAPW and LMTO methodologies have allowed the implementation of “all-electron” applications of the GW method.blochl ; hamada ; kotani ; ku ; lebegue  One common feature of such calculations, using standard level of approximation for the self-energy, is an underestimation of the electronic energy gap compared to experimental measurements, whereas pseudopotential-based calculations show very good agreement with experiment.h&l ; wilkins ; szl  To explain this inconsistency, it was proposed that the pseudopotential approach does not correctly describe the effect of core orbitals in the self-energy corrections to the energy gaps, resulting in overestimated corrections.kotani ; ku ; lebegue

It is thus desirable to elucidate the effect of core orbitals in the quasiparticle band structure, and the preferred procedure is to perform a well converged all-electron calculation and compare its results with similarly converged pseudopotential-based calculations. Obviously, numerical precision should not be neglected. In this work, we explicitly include semicore orbitals in the pseudopotential plane-wave approach and calculate the quasiparticle energy gap for three semiconductors of technical importance: Si, Ge and GaAs. The underlying description of the ground-state electronic structure is based on Density Functional Theory in the Local Density Approximation (DFT/LDA).ihm ; p-z  Throughout this work, we are careful to converge all results systematically, and the final results are compared to previous pseudopotential results and to recent all-electron calculations. The paper is organized as follows: we outline the theoretical method in section II. Results are presented in section III and discussed in section IV.

## II Theoretical Method

The pseudopotential formalism has two advantages that make it convenient for practical ab initio calculations. First, degrees of freedom due to core electrons are removed from the system, resulting in a description that contains only valence electrons. Second, valence electronic pseudowavefunctions are smooth in the vicinity of atomic sites as well as in the interstitial region. Such valence wavefunctions can be expanded to convergence easily in a relatively small basis set such as a plane-wave basis. On the other hand, one possible problem with this formalism is that using pseudopotentials, instead of the true electron-ion potential, may not fully describe all effects produced by interactions between the valence and core electrons. In order to address this issue in GW calculations, we include in the present study electrons from the outermost valence shell as well as those from the second outermost (“semicore”) atomic shell as active “valence” electrons in the pseudopotential formalism. Only interactions between these “valence” electrons and electrons from the deeper core shells are described by pseudopotentials. Taking Si as an example, its “core” now contains only $`1s`$ electrons. Electrons from the $`2s`$,$`2p`$,$`3s`$,$`3p`$, and $`3d`$ shells are all treated on equal footing in the subsequent calculations.

Apart from the atomic configuration, we follow the standard prescription for generating ab initio pseudopotentials. tm2 ; kerker  Since the deeper core electrons have extremely large binding energy (e.g., 130 Ry for the $`1s`$ electrons in atomic Si), their interaction with valence electrons is expected to be much weaker than the already small interaction between the outermost valence and the semicore electrons.rohlfing

The Kohn-Sham DFT formalism within the LDA is used to solve for the ground state electronic structure and to provide a starting point for the calculation of the electron self-energy. We follow closely the GW method as developed by Hybertsen and Louie.h&l  In this method, the self-energy is given by the standard GW approximation,

$$\mathrm{\Sigma }(𝐫,𝐫^{};E)=i\frac{\mathrm{d}E^{}}{2\pi }e^{iE^{}0^+}G_0(𝐫,𝐫^{};EE^{})W_0(𝐫,𝐫^{};E^{}),$$
(1)
where $`G_0`$ is the one-electron Green’s function, calculated from LDA energy eigenvalues and eigenstates ($`\delta =0^+`$ for occupied states, $`\delta =0^{}`$ for unoccupied states),

$$G_0(𝐫,𝐫^{};E)=\underset{n𝐤}{}\frac{\phi _{n𝐤}(𝐫)\phi _{n𝐤}^{}(𝐫^{})}{E\epsilon _{n𝐤}^{DFT}i\delta }.$$
(2)
The screened Coulomb interaction $`W_0`$ is calculated within the Random Phase Approximation (RPA) as $`W_0=\left[1VP_0\right]^1V`$, with $`V`$ being the bare Coulomb interaction and the polarizability $`P_0=iG_0G_0`$.hedin  This is the commonly employed level of approximation for GW calculations (i.e. neglecting self-consistency and vertex corrections) wilkins ; aryasetiawan\_rev , and we specifically compare results from both pseudopotential and all-electron calculations at this particular level of approximation to help untangle the effects of core states.

The convolution integral in Eq. (1) is performed using the generalized plasmon-pole (GPP) model,h&l  which enables one to distinguish two contributions to the self-energy: a screened exchange part ($`\mathrm{\Sigma }_{sx}`$) arising from the poles of $`G_0`$, and a dynamical Coulomb interaction between an electron and the hole-like charge distribution around it ($`\mathrm{\Sigma }_{ch}`$) stemming from the poles of $`W_0`$. The former tends to increase the quasiparticle energy, after the bare exchange is excluded. The latter tends to decrease the quasiparticle energy. In particular, the dynamical contribution $`\mathrm{\Sigma }_{ch}`$ is highly sensitive to the number of bands $`n`$ included in the calculation of the Green’s function in Eq. (2). The final quasiparticle energy of a state $`\phi _{n𝐤}`$ is given by ($`V_{xc}`$ is the LDA exchange-correlation potential)

$$E_{n𝐤}^{qp}=\epsilon _{n𝐤}^{DFT}+\phi _{n𝐤}|\mathrm{\Sigma }(𝐫,𝐫^{};E_{n𝐤}^{qp})V_{xc}(𝐫)|\phi _{n𝐤}.$$
(3)
## III Results

A semicore, non-relativistic pseudopotential was generated for Si, using the Troullier-Martins scheme.tm2  For Ga, Ge and As, we constructed semi-relativistic pseudopotentials using the Kerker scheme.kerker  These choices resulted in stable, transferable pseudopotentials in the Kleinman-Bylander form, without ghost states. k-b  A summary of atomic parameters is presented in Table 1. A good expansion of electronic wavefunctions in a plane-wave basis was obtained using a cutoff energy of 700 Ry (600 Ry for Si), and the first Brillouin zone was sampled using a 4x4x4 Monkhorst-Pack grid. m-pack  These numerical parameters ensure convergence in LDA energy eigenvalues to 0.01 eV or better. The Ceperley-Alder exchange-correlation potential is used.p-z  For the lattice parameter, we used the experimental values: 5.43 $`\AA `$, 5.65 $`\AA `$, and 5.66 $`\AA `$ for Si, Ge, and GaAs, respectively. The polarizability was expanded in a plane-wave basis with an energy cutoff of 45 Ry (50 Ry for Si) and numerically inverted for the calculation of the screened Coulomb interaction $`W_0`$. Numerical precision in the calculation of the self-energy is 0.05 eV or better.

Table 2 shows some of the energy gaps obtained in the present approach for Si, compared with previous valence-only pseudopotential and all-electron calculations. Overall agreement between the present results (which explicitly include the effect of the semicore states) and experimental measurements is at the level of 0.1 eV, and discrepancies between our results and previous pseudopotential-based calculations of Ref. 9 are equally small. Recent all-electron calculations carried out at the same level of the GW approximation, however, systematically underestimate the minimum gap and the direct $`\mathrm{\Gamma }\mathrm{\Gamma }`$ gaps.kotani ; ku  We find that the convergence of the self-energy with respect to the number of unoccupied bands included in $`G_0`$ in Eq. (2) is an important factor. In Figure 1, we show the behavior of the calculated energy gap as function of the number of unoccupied bands, $`n_c`$, included in $`G_0`$. Convergence is typically very slow, and well-converged results require $`n_c120`$. Other energy transitions ($`\mathrm{\Gamma }X`$ and $`\mathrm{\Gamma }\mathrm{\Gamma }`$) show similar behavior and also approach the converged value from below. In contrast, the results of Ref. 7 were obtained with only $`n_c=24`$ and are closer to our results at approximately the same value of $`n_c`$ than to the converged results, as shown in Figure 1. This fact points to a lack of numerical convergence in the evaluation of the self-energy corrections in Ref. 7.

The energy gaps obtained for Ge are presented in Table 3. Spin-orbit interactions are included as first-order perturbations.h&l-so  A common feature of LDA-based calculations is the overlapping of the valence and conduction bands at the $`\mathrm{\Gamma }`$ point. The inclusion of $`3d`$ electrons moves the conduction bands further down resulting in a sizable negative direct gap.kotani ; ku ; szl ; rohlfing ; aryasetiawan\_prb  This feature is verified in our LDA calculation. As shown in Table 3, self-energy corrections are responsible for an opening of the gap and the correct positioning of the minimum, indirect gap between points $`\mathrm{\Gamma }`$ and $`L`$ in the Brillouin zone. Our GW results compare well with experiment, although the $`\mathrm{\Gamma }X`$ gap still shows a large discrepancy. On the other hand, there are significant differences between the pseudopotential and all-electron gaps at this level of the GW approximation, in particular regarding the $`\mathrm{\Gamma }X`$ gap. As found in previous studies, the self-energy corrections are needed to give the correct band topology.

Figure 2 shows direct and indirect energy gaps in Ge as function of the number of unoccupied bands included in the calculation of $`G_0`$. Whereas the indirect gaps $`\mathrm{\Gamma }L`$ and $`\mathrm{\Gamma }X`$ approach the converged value from below, the direct gap $`\mathrm{\Gamma }\mathrm{\Gamma }`$ approaches it from above. This particular convergence behavior arises from the fact that we are plotting differences of quasi-particle energies: taken individually, all quasi-particle energies $`E_{n,𝐤}^{qp}`$ converge monotonically from above, reflecting the attractive nature of the Coulomb-hole $`\mathrm{\Sigma }_{ch}`$ term.h&l ; wilkins  Additionaly, we note that in Ref. 7 the direct gap at $`\mathrm{\Gamma }`$ is strongly overestimated and the indirect gap $`\mathrm{\Gamma }X`$ is underestimated by $`0.8`$ eV. As expected, the same pattern of overestimation/underestimation is evident in Figure 2 when we reduce the number of unoccupied bands from 170 to 24, which was the value used in Ref. 7.

GaAs shows behavior similar to Ge. Table 4 summarizes our results. Agreement with experimental data is now within 0.15 eV, and all-electron calculations again underestimate the energy gaps. Regarding the convergence with $`n_c`$, Figure 3 shows that the indirect gaps converges slowly from below. On the other hand, the direct $`\mathrm{\Gamma }\mathrm{\Gamma }`$ gap converges more quickly.

In all systems studied, we applied the GPP model h&l  in two ways: (1) using only the valence charge density in the f-sum rule, and (2) using a total charge density from both the valence and semicore bands. Method (1) is physically more realistic, and predicts an energy dependence of the inverse dielectric function that is consistent with the RPA. Method (2) implies that semicore electrons are able to screen electric fields as efficiently as valence electrons, which is not physical. Nevertheless, we find agreement between the two schemes to better than 0.1 eV in the converged energy gap for all systems studied. This is a consequence of cancellation of errors: in method (2), the plasma frequency is overestimated, and therefore the contributions $`\mathrm{\Sigma }_{sx}`$ and $`\mathrm{\Sigma }_{ch}`$ are overestimated in absolute value. Since they have opposite sign, the final quasiparticle energies are weakly affected. We also observe that method (2) shows slower convergence of energy gaps with respect to $`n_c`$, which is expected since the $`\mathrm{\Sigma }_{ch}`$ term is enhanced. All results presented in this article were obtained using the physically more appropriate method (1).

## IV Discussion

As shown in Figures 1-3, energy gaps in the GW approximation have a significant dependence on the number of unoccupied bands that is slow to converge in many cases. This fact was already reported in calculations using standard pseudopotential techniques.wilkins ; fleszar  We can analyze the physics of this convergence by examining the GW approximation in the static limit, the so-called COHSEX approximation.h&l ; wilkins ; hedin  Bands gaps can be calculated more easily within the COHSEX approximation and show convergence behavior similar to the full dynamic calculation. Under this approximation, the Coulomb-hole term has a simple form in terms of the polarization potential, $`W_{pol}(𝐫,𝐫^{})W(𝐫,𝐫^{},\omega =0)V(𝐫𝐫^{})`$,

$$\mathrm{\Sigma }_{COH}(𝐫,𝐫^{})=\frac{1}{2}\delta (𝐫𝐫^{})W_{pol}(𝐫,𝐫^{})=\frac{1}{2}\underset{n𝐤}{}\phi _{n𝐤}(𝐫)\phi _{n𝐤}^{}(𝐫^{})W_{pol}(𝐫,𝐫^{}),$$
(4)
where the second equality follows from completeness of the basis of eigenvectors. In actual calculations, this sum over bands $`n`$ is always truncated, and the equality is violated. The Coulomb-hole energy evaluated at a given electronic state is calculated according to the expression

$$m𝐤|\mathrm{\Sigma }_{COH}|m𝐤=\frac{1}{2}\underset{\mathrm{𝐆𝐆}^{}𝐪}{}W_{\mathrm{𝐆𝐆}^{}}(𝐪)\underset{n}{}\left[_𝐆^{nm}(𝐤,𝐪)\right]^{}_𝐆^{}^{nm}(𝐤,𝐪),$$
(5)
where $`W_{\mathrm{𝐆𝐆}^{}}(𝐪)`$ are the coefficients in the plane-wave expansion of the polarization potential, h&l ; wilkins  and we define

$$_𝐆^{nm}(𝐤,𝐪)=n𝐤𝐪|\mathrm{e}^{i(𝐪+𝐆)𝐫}|m𝐤.$$
(6)
Physically, the summation in Eq. (4) describes virtual transitions produced when the quasiparticle induces a charge fluctuation around itself. $`\mathrm{\Sigma }_{COH}`$ is the energy associated to the interaction between the quasiparticle and the induced charge fluctuation. The matrix elements, Eq. (6), decay slowly as the energy difference between bands $`m`$ and $`n`$ increases. Monitoring the convergence of these matrix elements provides a good estimate of the relative error in the $`\mathrm{\Sigma }_{COH}`$, but gauging the absolute convergence of the self-energy requires knowing $`W_{pol}`$, which depends on the physical system.

Although we have not directly investigated the importance of imposing self-consistency in the calculation of the self-energy, this is an unsetled issue and we address it briefly below. Von Barth and Holm have investigated the effect of self-consistency in the electron gas,holm  and concluded that restricted self-consistency has small but significant effect on the full bandwidth and in the satellite structure of the electron gas. On the other hand, full self-consistency gives a poor description of the satellite structure and the bandwidth is drastically increased. Inclusion of vertex corrections are expected to recover the good, non-self-consistent results, but calculation of vertex corrections is not a simple task even for the electron gas system.wilkins ; holm  Self-consistency has been recently applied to real materials,ku ; eguiluz  and the valence bandwidth is also shown to increase when self-consistency is imposed. It appears that one must therefore include self-consistency and vertex corrections together in order to obtain a meaningful picture. Inclusion of vertex corrections and self-consistency is tangential for the purpose of our work, which is to compare pseudopotential-based and all-electron GW calculations and understand the role of semicore electrons.

## V Conclusion

We conclude that for the systems considered, the valence-only pseudopotential method does not suffer from large errors from the neglect of core states, as claimed in some all-electron calculations. kotani ; ku  While semicore effects are negligible in Si, they are important if one aims at good quantitative agreement with experiment in Ge and GaAs.wilkins ; aryasetiawan\_rev ; aryasetiawan\_prb ; rohlfing  However, discrepancies between pseudopotential and all-electron based GW calculations reported in recent works kotani ; ku  may be explained by a lack of numerical convergence in the latter. Specifically, the self-energy calculated within the GW method has slow convergence with respect to the number of energy bands included in the calculation of the Green’s function, as is demonstrated in this work and has been pointed out in the past. wilkins ; fleszar  This convergence behavior is present in the static limit to GW, the COHSEX approximation, and can be analyzed by comparing the COHSEX Coulomb-hole energy, Eq. (4), obtained with and without explicit summation over energy bands.

###### Acknowledgements.

We thank Andrew Canning and David Clatterbuck for fruitful discussions regarding the FLAPW method. This work was supported by National Science Foundation Grant No. DMR00-87088 and by the Director, Office of Science, Office of Basic Energy Sciences, Division of Materials Sciences and Engineering, U.S. Department of Energy under Contract No. DE-AC03-76SF00098. Computational resources have been provided by NSF at the National Partnership for Advanced Computational Infrastructure, and DOE at the National Energy Research Scientific Computing Center.