text:
# Evolution of the Inner Circumstellar Envelope of V838 Monocerotis

## 1 Introduction

The eruptive variable V838 Monocerotis, first reported as a possible nova on 2002 January 6.6 (Brown et al., 2002), experienced three significant photometric outbursts in early 2002 (Munari et al., 2002a; Kimeswenger et al., 2002). From pre-outburst to its maximum brightness during the second outburst, V838 Mon brightened by over nine magnitudes in V, from 15.85 to 6.66 (Goranskii et al., 2002). V838 Mon also exhibited significant spectroscopic variability in early 2002: neutral metal and s-process lines were observed following the first outburst (Zwitter & Munari, 2002), ionized metal lines were noted following the second outburst (Iijima & Della Valle, 2002; Morrison et al., 2002), and neutral metal and molecular lines were observed following the third outburst (Rauch et al., 2002; Banerjee & Ashok, 2002). V838 Mon exhibited intrinsic polarization on 2002 February 8 (Wisniewski et al., 2003). Polarimetric observations after 2002 February 13 only detected the presence of an interstellar polarization (ISP) component (Munari et al., 2002b; Wisniewski et al., 2003).

V838 Mon developed a light echo by 2002 February 17 (Henden et al., 2002). Bond et al. (2003) obtained imaging polarimetry of this light echo with HST ACS, and suggested a distance of 3 to 7 kpc to V838 Mon, although a wide range of other distances have been estimated (e.g. see discussion in Wisniewski et al. 2003 and Tylenda 2003). In late September and early October of 2002, V838 Mon’s spectral type had evolved to “later than M10-III” (Desidera & Munari, 2002). Furthermore, a weak blue continuum was detected by several groups, suggesting that a binary component of spectral type B3V might be present (Desidera & Munari, 2002; Wagner & Starrfield, 2002; Munari et al., 2002c). Recent MMT spectra confirm this secondary component (Starrfield, private communication).

Given the magnitude and complexity of the variability exhibited by V838 Mon, it is not surprising that there is still great uncertainty regarding the exact nature of this object. In this paper, we present imaging polarimetry observations of V838 Mon and its surrounding field obtained in 2002 October. These observations allow us to further interpret the evolution of the circumstellar matter near V838 Mon since its major outbursts.

## 2 Observations

We obtained polarimetric measurements of V838 Mon and its field in 2002 October at the CTIO 1.5 m telescope. The standard Cassegrain focus imaging camera with the f/13.5 (0.24 arcsec pixel-1) configuration was modified by the addition of a rotatable half-wave plate followed by a fixed analyzer placed before the second filter wheel. The analyzer was a double calcite block whose optical axes had been crossed to minimize astigmatism and color effects, which produced two orthogonally polarized images of each object in the field. One polarization modulation is covered for each 90 rotation of the waveplate. For our V, R, and I observations of V838 Mon, we took images with the waveplate rotated through 8 positions 22.5 apart. Further details about this polarimeter can be found in Magalhães et al. (1996), Melgarejo et al. (2001), & Pereyra & Magalhães (2002).

After basic image processing in IRAF111IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under contract with the National Science Foundation., we performed aperture photometry on the fields. The reduction package PCCDPACK (Pereyra, 2000), which calculates linear polarization from a least squares solution of the photometry in the 8 waveplate positions ($`\psi _i`$), was then applied. The residuals at each waveplate position with respect to the expected $`\mathrm{cos}4\psi _i`$ curve constitute the uncertainties in our data; these are consistent with the theoretically expected photon noise errors (Magalhães, Benedetti, & Roland, 1984). Each object whose polarization is reported has been carefully checked for contamination effects due to close neighbors; objects which contained such contamination were excluded from the present study.

Polarized and unpolarized standard stars were monitored nightly throughout the eleven nights of our observing run, and the stability of the polarimeter was reflected in the consistency of these standard star observations. Comparing our observations of polarized standard stars to available literature data, we transformed our polarization measurements to a standard equatorial system. This transformation is accurate to a position angle (PA) of better than 1. Instrumental polarization was consistently less than 0.05$`\%`$. Note that the simultaneous observation of the dual orthogonally polarized images of each target in each waveplate position allows for accurate polarimetry to be performed even under non-photometric skies, as all sky polarization is practically canceled (Magalhães et al., 1996). We summarize our observations of V838 Mon, along with polarimetry of this object from the literature, in columns 1-5 of Table 1.

## 3 Results

### 3.1 Polarization of Field Stars

Over small spatial regions, where variations in the distribution of dust and magnetic fields can be assumed to be small, field stars located at similar distances exhibit similar levels of interstellar polarization (ISP). To study the intrinsic polarization of an astrophysical object, it is often beneficial to examine the polarization of spatially nearby field stars, located at similar distances, to get a statistical measure of the ISP along a particular line of sight.

We measured the polarization in V, R, and I of most of the field stars surrounding V838 Monocerotis; a sample of these results is given in Table 2 and a complete listing of these data is available online. Figure 1a shows a V band image obtained on 2002 October 16 at the CTIO 0.9m, and includes numerical references to stars shown in Figure 1b and tabulated in Table 2. A plot of the spatial distribution of the polarization of these field stars in the V band is shown in Figure 1b. We note that three stars near V838 Mon show very similar polarizations to each other, star #2, star #3, and star #5. Furthermore, these stars exhibit very similar polarizations to the ISP previously estimated for V838 Mon (Wisniewski et al., 2003). We plot the wavelength dependence of the polarization of these three objects in Figure 2, and overlay the modified Serkowski law (Serkowski et al., 1975; Wilking et al., 1982) description of V838 Mon’s interstellar polarization with $`P_{max}=2.746\pm 0.011\%`$, $`\lambda _{max}=5790\pm 37\AA `$, $`PA=153.43\pm 0.12^{}`$, $`\delta PA(\lambda )=0`$, and $`K=0.971`$ (Wisniewski et al., 2003). The good agreement shown in this figure supports the previously determined ISP, derived entirely from spectropolarimetry. Furthermore, these results imply that stars #2, #3, and #5 are located at similar distances as V838 Mon. Future efforts to determine the distance to V838 Mon should ensure consistency with these three neighboring objects.

### 3.2 Intrinsic Polarization of V838 Mon

We subtracted the ISP along the line of sight (Wisniewski et al., 2003) from our observations of V838 Mon to determine its intrinsic polarization. The non-zero results (columns 6-8, Table 1) show that on 2002 Oct 22-24 V838 Mon once again exhibited a significant intrinsic polarization component. Interestingly, the magnitude of intrinsic polarization present on 2002 October 22-24 is similar to that observed on 2002 February 8; however, the position angle of this scattered light differs by roughly 90, e.g. a PA “flip” has occurred.

The presence of an intrinsic polarization component can be interpreted as a signature of light scattering off circumstellar material which is distributed in an asymmetrical geometry. The intrinsic polarization reported here and by Wisniewski et al. (2003) is produced in a region immediately surrounding the unresolved central star, and as such is different from the reported polarization of V838 Mon’s light echo (Bond et al., 2003). In our observations of the unresolved source, the observed polarization is roughly the ratio of the scattered to total (i.e. direct plus scattered) light received. Thus the renewed intrinsic polarization component suggests that either a) a new source of asymmetrical circumstellar material, i.e. a new source of scatterers, formed around V838 Mon in 2002 October; b) the opacity of the circumstellar envelope surrounding V838 Mon evolved significantly between February 5 to February 13 to October 22-24, thereby changing the dominant source of scattered light over time; c) the illuminating source had changed; or d) a combination of the aforementioned scenarios.

## 4 Discussion

Schulte-Ladbeck et al. (1992) interpreted an observed wavelength-dependent PA flip in polarimetric observations of an unresolved B\[e\] star as evidence of a bipolar nebula. These authors argued that at short wavelengths an optically thick circumstellar disk blocked out most starlight and thus little scattered light (with a PA aligned along the polar axis) originated from the equatorial disk. Most of the polarized light at short wavelengths originated from the polar regions, where the resulting PA would be aligned with the equatorial disk. At long wavelengths, where the disk was optically thin, the dominant source of scattered light was the equatorial disk region, not the polar regions. The net effect of such a wavelength dependent opacity effect was the production of a 90 PA flip in the polarization signal at the wavelength where the equatorial and polar regions contributed equal amounts of scattered light.

We suggest that a similar type of scenario might explain the renewed intrinsic polarization component and PA flip in V838 Mon. On 2002 February 8, V838 Mon had an intrinsic polarization component, indicating the presence of an asymmetrical distribution of circumstellar material, initially suggested to be a flattened circumstellar envelope (Wisniewski et al., 2003). We postulate that at this initial stage, the observed scattered light originated primarily from specific physical locations, e.g. the polar region. As the opacity of the circumstellar material evolved over time, specifically as the opacity of the equatorial material decreased, we suggest that the contribution of scattered light from the polar and equatorial regions were nearly equivalent. Thus, viewed as an unresolved object, the scattered light from the circumstellar envelope would appear to be unpolarized, consistent with observations in late 2002 February and March (Wisniewski et al., 2003; Munari et al., 2002b). As the envelope continued to evolve and the equatorial region experienced a further decline in opacity, one would expect the equatorial region to slowly emerge as the dominant source of scattered light. This re-emergence of a dominant (and secondary) scattering region would produce an intrinsic polarization component oriented 90 from the original position angle. The projection of the intrinsic PA of this envelope onto the sky (2 $`\theta 75^{}`$, measured N to E) is not inconsistent with the overall morphology seen in the HST ACS images (Bond et al., 2003). The intrinsic PA of the polar region is also not inconsistent with the direction of the hole in the HST ACS light echo images.

We note that the opposite relative opacity changes would produce a similar PA flip. The scattered light in the circumstellar envelope could have initially been dominated by an optically thin disk. As this equatorial material dispersed, the scattered light from the equatorial and polar regions could have balanced, producing zero intrinsic polarization. Over time, expansion of the disk could have evacuated the equatorial region, causing the polar region to emerge as the dominant source of scattered light.

Followup infrared (IR) imaging polarimetry would be valuable in providing additional constraints on the nature of V838 Mon’s circumstellar envelope. Based upon infrared (IR) light and color curves, Crause et al. (2003) suggest that a dust shell had formed around the central star by 2002 April. If in early 2002 this dust formed in the equatorial regions as described in the aforementioned first scenario, i.e. an initial thick equatorial and thin polar region, then the IR polarization position angle would be the same as that quoted in this paper.

As described in the introduction of this paper, V838 Mon appears to have a binary companion. Polarimetric observations of mass transfer binary systems show some degree of periodicity and scattering at a constant PA, defined by the orbital plane of the system (Hoffman et al., 1998). The observed PA flip in V838 Mon’s intrinsic polarization, along with the lack of any apparent periodicity, suggest that binarity is not the dominant source of scatterers responsible for the observed intrinsic polarization component.

In summary, the polarization of field stars confirms the previously suggested interstellar polarization along the line of sight to V838 Mon. We find strong evidence of a renewed intrinsic polarization component in V838 Mon, at a position angle nearly 90 from that present on 2002 February 8. We suggest that these observations indicate an evolution of the circumstellar environment of V838 Mon, possibly characterized by: 1) an initially optically thick equatorial disk which contributed a limited amount of scattered light; 2) a subsequent decline in disk opacity resulted in a balance of scattered light produced by the polar and equatorial regions; and 3) the further decline in the opacity of the equatorial disk finally allowed the equatorial region to become the dominant source of scattered light. Continued spectropolarimetric monitoring of V838 Mon is strongly encouraged, as such observations would enable detailed modeling of the circumstellar environment of this unique object to be performed.

We thank the referee, S. Starrfield, for his suggestions to improve this paper. We also thank Antonio Pereyra and Brian Babler for assisting with the data reduction. This work has been supported in part by NASA LTSA grant NAG5-8054 to the University of Toledo. KSB is a Cottrell Scholar of the Research Corporation and gratefully acknowledges their support. JPW is supported by a NASA GSRP fellowship and thanks NOAO for supporting his travel to CTIO. AMM acknowledges travel support by FAPESP; he is also partially supported by CNPq. Polarimetry at the University of São Paulo (USP) is supported by FAPESP. AMM and KSB acknowledge partial travel support by USP.


format: text

text:
# Phonon and plasmon excitation in inelastic electron tunneling spectroscopy of graphite

## Abstract

The inelastic electron tunneling spectrum (IETS)of highly oriented pyrolitic graphite (HOPG) has been measured with scanning tunneling spectroscopy (STS) at 6K. The observed spectral features are in very good agreement with the vibrational density of states (vDOS) of graphite calculated from first principles. We discuss the enhancement of certain phonon modes by phonon-assisted tunneling in STS based on the restrictions imposed by the electronic structure of graphite. We also demonstrate for the first time the local excitation of surface-plasmons in IETS which are detected at an energy of 40 meV.

Initially applied to measure inelastic tunneling currents across metal-insulator-adsorbate-metal junctions wolf , Inelastic Electron Tunneling Spectroscopy (IETS) has in the recent past been revived and its concept extended to Scanning Tunneling Microscopy (STM-IETS) stipe . Superimposed on the elastic tunneling current, the inelastic tunneling process opens up alternative tunneling channels increasing the total conductance at the onset of its excitation. Very recently, elegant experiments (e.g. stipe ; kim ; lee ; pascnat ) have been performed in detecting local inelastic electron energy-loss spectra of individual molecule adsorbed on metallic substrates. The observed features in the tunneling spectra are found to correspond to the vibrational modes of the adsorbed species. Vibrations of single molecules as large as C60 have have been detected pascual . Despite the increasing number of publications of STM-IETS experiments on single molecules adsorbed on metals, to our knowledge, only two studies were reported on collective vibrational excitations of surfaces smith ; olejn . In both these studies graphite was used as a substrate and several spectral features remained unexplained or were attributed to vibrational excitation in the tungsten tip. However, these tip effects have not been detected in many other STM-IETS experiments. In order to resolve the pure inelastic response of a graphite surface, we have performed STM-IETS experiments on a highly oriented pyrolitic graphite crystal (HOPG) at cryogenic temperature. In our experiments the whole vibrational spectrum of HOPG from the rigid layer shear mode at $``$ 5meV to the optical modes at $``$ 200 meV is clearly observed without being perturbed by inelastic effects from the tip. The phonon modes of graphite can be clearly identified by the direct comparison of experimental data with first principles density function theory (DFT) calculations. We report for the first time that also surface-plasmon losses can also be observed directly with STM-IETS.

Freshly cleaved HOPG has been introduced into the vacuum system and shortly annealed at 500K. The experiments were performed at 6K with a home built UHV-STM using a tungsten tip. The inelastic tunneling spectrum of HOPG was measured over a topographically clean and flat terrace with a lock-in technique. In this technique a sinusoidal signal with amplitude $`\mathrm{\Delta }V`$ is superimposed to the applied sample bias. In IETS the second harmonic component of the modulated signal, which is proportional to the second derivative vs. bias of the tunneling current, $`d^2I/dV^2`$, is recorded. In this way the inelastic loss features in the conductance are transformed into peaks at the onset of the loss modes for positive sample bias and into minima for negative sample bias, respectively. A typical example for HOPG is shown in fig.1.

We obseve phonon features that are symmetrically located with respect to zero bias as expected for the inelastic tunneling process. Possibly due to a particular tip-end configuration teh phonon on the negative bias side are more pronounced. In Fig. 2a we present a high-resolution spectrum of the $`d^2I/dV^2`$ signal of HOPG taken at negative bias. The data was obtained with a bias modulation of $`\mathrm{\Delta }V=10`$ meV RMS at a frequency of 2.7kHz and inverted to render the phonon excitations as peaks to facilitate the comparison to theory. Since the electron density of states of HOPG is flat in the energy range close to the Fermi level, the features in the spectrum in Fig. 1 are exclusively due to inelastic processes.

Along with the measured inelastic spectrum, we present in Fig. 2b and c the calculated vibrational density of states (vDOS) and the phonon dispersion relation of graphite along the A-$`\mathrm{\Gamma }`$-M-K directions. Graphite is a layered structure with strong carbon bonds in the plane and weak Van der Waals-like forces between the planes. As a result of this anisotropy, the phonon spectrum of graphite covers a very wide range. We have performed first principles calculations using density-functional theory (DFT) in the local density approximation (LDA) calcdetails . Aim of this calculation is to resolve some of the deviations in the vDOS calculated by the different recent force-constant parametrizations for graphite jishi ; gruneis ; aizawa ; sieben  and to have a solid foundation for the assignment of the measured peaks. The first principles calculation has been shown kresse ; pavone ; dubay  to yield very close agreement with phonon dispersion relations measured by high resolution electron energy loss spectroscopy (HREELS) along the $`\mathrm{\Gamma }`$-M direction aizawa  and along the $`\mathrm{\Gamma }`$-K direction sieben . Only for the acoustic shear-mode (SH) it was observed dubay  that at the M-point, the HREELS data of Ref.aizawa,  yields an energy which is higher by 25 meV than the ab-initio value. A similar underestimation of of the shear-mode frequency might occur as well at the K-point where a mode with energy 134 meV was observed by double-resonance Raman scattering kawashima . Our calculated vDOS is also at variance with the HREELS data at M aizawa  but is in very good agreement with the vDOS measured by neutron scattering on a powdered sample of graphite rols  where a separate SH(M)-feature was not observed. However, since the SH K-point mode contributes only weakly to the calculated DOS the comparison to the experimental vDOS in Ref. rols,  cannot answer the question whether the mode comes out too soft in the calculations.

The observed vibrational spectral features, their comparison with values obtained with other experimental techniques and their assignment to phonon modes at the various points of the surface Brillouin zone (SBZ) are summarized in Table 1. Most features coincide with peaks in the calculated vDOS which arise from the flat dispersion relation around the high-symmetry points. In particular, the higher peaks due to phonon modes around the M-point are clearly mapped. The first peak at 16 meV is in good agreement with the calculated energy of the out of plane acoustic (ZA) mode at $`\mathrm{\Gamma }`$ where two neighboring planes are oscillating out-of-phase. While in the vDOS and in neutron scattering rols  this peak is only weak, in the IETS experiment it is stronlgy pronounced which can be explained by the direction of oscillation which is orthogonal to the surface plane. The weak shoulder at 6 meV corresponds to the low frequency horizontal mode at $`\mathrm{\Gamma }`$ where two neighboring planes are oscillating out-of-phase parallel to each other. The prominent feature at 40 meV is not due to a vibrational feature but is assigned to a plasmon excitation (see below). In Ref.smith, , a phonon feature at 39 meV was observed which was associated with a double-phonon excitation in the surface (16 meV) and in the tungsten tip (23 meV). However, in our experiment the tip-phonon peak is clearly not present and all the observed features can therefore be assigned to surface excitation. The vibrational feature observed at 58 meV can be assigned to the ZA at the M-point. The optical out-of-plane (ZO) mode and the shear acoustic mode SH (also called transverse acoustic mode) display a saddle in the dispersion relation around the M-point adding up to a very high peak in the vDOS at 78 meV. This explains the peak at 81 meV in the measurement. However, as noted above, HREELS data aizawa  suggests the energy of the SH mode at M to be about 100 meV. This gives rise to the alternative explanation that the peak at 81 meV is only due to the ZO mode while the shoulder at 100 meV may be a signal of the SH mode at the M-point. The ZO mode at $`\mathrm{\Gamma }`$ which gives rise to a peak at 108 meV in the theoretical vDOS seems to be missing in the spectrum. Instead, the spectrum reveals actually a dip ($`dI^2/dV^2<0`$) at this energy which will be discussed below. At 124 and 150 meV two low-intensity phonon modes, the SH, and the longitudinal acoustic (LA) modes at the K-point are predicted, respectively. Experimentally we have observed a large asymmetric inelastic feature at about 137 meV. This disagreement between experiment and theory might point to an underestimation of the SH mode also at the K point. If, as could be the case for the SH mode at M dubay , the real energy is higher by 10 to 20 meV, the peak at 137 meV can be assigned to the SH mode at K which corresponds to the above mentioned observation of a K-point phonon at 134 meV kawashima  by double resonance Raman spectroscopy. In the highest energy range of the spectrum two loss-features protrude the background signal at 180 and 200 meV, respectively. The first is related to the optical shear mode (SH, also called transverse optical mode) around the M-point. The neighboring feature at 200 meV corresponds to the longitudinal optical (LO) branch which displays a strong overbending in the direction $`\mathrm{\Gamma }`$-M and $`\mathrm{\Gamma }`$-K giving rise to a peak in the vDOS.

In general, the measured IET spectrum follows very closely the vibrational density of states with two notable exceptions which will be discussed in the following. These exceptions are related to the question if all phonons can be excited or if the geometry of the tunneling system leads to a restriction to phonons in a certain region of the Brillouin zone, as it is the case in phonon-assisted electronic tunneling through narrow p-n junctions of semiconductors. In e.g. Ge, the valence-band maximum is centered at $`\mathrm{\Gamma }`$ while the conduction-band minimum is located at the $`L`$ point. This allows only vibrational modes at L to be excited by a phonon-assisted tunneling process payne . The electronic structure of graphite is also quite unique being a semimetal with a band crossing only at the K-point of the SBZ. Therefore the elastic current is carried by electrons that tunnel with a relative large $`k_{}=`$ K (1.7 Å-1). In an STM geometry the expansion of the initial state, i.e., the tip wave function in terms of surface plane-waves tersoff , leads to a parallel momemtum distribution centered at $`k_{}=0`$. Consequently the electrons of the elastic current have a low transmission probability. Electrons with any other $`k_{}`$ can only tunnel by exchange of finite crystal momentum with a phonon. The phonon-assisted tunneling is relevant to all the dominant features observed besides the $`\mathrm{\Gamma }`$ phonons at 16 and 200 meV. An enhanced tunneling contribution should especially occur for electrons with initial parallel momentum $`k_{}=0`$ that tunnel to the Fermi level at the K-point by exciting (or absorbing) a K-point phonon. This explains why the peak at 137 meV which is probably due to a K-point phonon excitation shows up so dominantly in the spectrum even though the vDOS of states only displays weak features for K-point phonons.

The second deviation of the measured spectrum from the vDOS of graphite is the dip at 111 meV where the vDOS displays a peak due to the ZO phonon at $`\mathrm{\Gamma }`$. The theory for inelastic tunneling spectroscopy persson ; lorente  distinguishes between an inelastic contribution (real phonon emission or absorption) and an elastic contribution (virtual phonon emission with subsequent reabsorption or phonon absorption with subsequent re-emission). In some cases, the elastic contribution leads to a decrease of the tunneling conductance (dI/dV) due to back-scattering of the electron persson . This effect has been measured recently as dips in the $`dI^2/dV^2`$ curve of O2 on a silver surface hahn . We believe that the dip at 111 meV is due to a similar elastic contribution at the onset of vDOS singularity corresponding to the ZO($`\mathrm{\Gamma }`$) phonon branch. In the light of the above discussion of phonon-assisted tunneling on graphite, it is also not surprising that this effect is observed for phonons at $`\mathrm{\Gamma }`$. The only matching pair for a virtual phonon scattering process are electrons initially at $`k_{}=K`$ and phonons at $`\mathrm{\Gamma }`$. All other combinations are forbidden due to the unavailability of intermediate electronic states. Furthermore, the steady reduction of $`d^2I/dV^2`$ from the ZO (M) peak at 81 meV, in contrast to the vDOS, shows that symmetry of the phonon displacement (the ZO branch) should play an important role in the supression of the IET current lorente2 . We note that one of the degenerate modes of the ZO branch at $`\mathrm{\Gamma }`$ is the only $`\mathrm{\Gamma }`$ phonon mode that has $`A_{2u}`$ symmetry.

In addition to the predicted phonon modes, the measured spectrum contains an energy loss feature at 40 meV, which cannot be explained by phonon excitations. Other than the previous features, this loss feature is of purely electronic origin and corresponds to an out-of-plane oscillation of the electrons. A tunneling electron traversing a vacuum-metal interface has a significant probability of interacting with collective electron oscillations. The overlapping of the $`\pi `$ and $`\pi ^{}`$ band at the K point of the SBZ yields a low density of free carriers (i.e. electrons and holes) near the Fermi level and allows low-energy plasmon excitation. A strong temperature dependence of the excitation energy of this plasmon mode was observed in HREELS jensen . The excitation energy shifts from 63 meV to 45 meV when going from room temperature to 150 K, respectively. This effect has been attributed to the low density of states at the Fermi level in conjunction with the thermal excitations of carriers. Thermal excitations from the valence into the conduction band can radically alter the electron population at the Fermi level and yields the strong dependence of surface plasmon energy and line width on temperature. The temperature dependence of this excitation energy was also calculated in Ref. jensen, . Within the random phase approximation the plasmon energy excitation was predicted to be lowered to 36 meV at 0K in rather good agreement with our experimental value. The low-energy plasmon was also observed in transmission electron energy loss spectroscopy, which also detected a second plasmon loss, assigned to an excitation polarized parallel to the layer at 128 meV geiger . The shoulder observed around 125 meV might be assigned to this plasmon excitation.

In conclusion, the IET spectrum of HOPG has been measured by low temperature scanning tunneling spectroscopy and is compared to an ab-initio calculation of the vibrational density of states in graphite. The observed vibrational loss features can be identified with peaks in the vDOS which are due to the flat dispersion of the phonon branches around the high-symmetry points of the Brillouin zone. Additionally, a prominent plasmon feature at 40meV has been identified in the inelastic tunneling current measured with STM. In the past, the local capability of STM-IETS has been proven to be a valuable tool to access the vibrational properties of adsorbates at surfaces. In this work we have shown that STM-IETS resolves the vibrational and collective electronic excitations of graphite with high accuracy. This is an important step for the application of the method to confined geometries and individual carbon nanostructures.

###### Acknowledgements.

L.W. and A.R. thank support from the EC research training network COMELCAN (HPRN-CT-2000-00128), Spanish MCyT(MAT2001-0946) and UPV/EHU (9/UPV 00206.215-13639/2001).


text:
# X-ray Spectroscopy of Thermally Emitting Neutron Stars

## 1. Introduction

ROSAT discovered a number of nearby neutron stars whose emission appears to be entirely thermal, uncontaminated by accretion or magnetospheric processes. At present, six (possibly seven) sources are known (Treves et al. 2000; Haberl 2003). For four sources, optical counterparts have been identified. The high X-ray to optical flux ratios leave no model but an isolated neutron star.

As a class, the sources are interesting because of the implied existence of a fair number of neutron stars different from the usual radio pulsars and X-ray binaries. At present, their nature remains unclear. The original idea, of old neutron stars accreting slowly from the interstellar medium, has become unlikely because of the very low accretion rates implied by high proper motions. Instead, they might be radio pulsars beamed away from us, although in this case the long periods are surprising. Perhaps they have very strong magnetic fields, and are descendants of anomalous X-ray pulsars and soft gamma-ray repeaters.

As individuals, the sources are of particular interest because of the opportunity to study uncontaminated emission from a neutron-star atmosphere. The hope is that this will allow one to infer precise values of the temperature, surface gravity, gravitational redshift, and magnetic field strength. In turn, these could be used to constrain the interior, and hence learn about the equation of state of cold, ultradense matter, an unexplored region in QCD parameter space.

## 2. X-ray Spectra

Given the interest, long spectroscopic observations were taken with the Chandra X-ray Observatory and XMM-Newton. The first results were discouraging: no lines were found in XMM spectra of RX J0720.4$``$3125 (Paerels et al. 2001), nor in Chandra spectra of the prototype of the class, RX J1856.5$``$3754 (Burwitz et al. 2001), not even after 500 ks (Drake et al. 2002; Braje & Romani 2002). Instead, for both sources the spectra were found to be remarkably well described by mildly extincted black bodies. This was unexpected, since for light-element atmospheres one would expect a hard tail (because the opacities decrease towards higher energies), while for anything else one would expect to see lines (for reviews, see Pavlov, Zavlin, & Sanwal 2002; Zavlin & Pavlov 2002).

The situation changed this year. For RX J1308.6+2127, the fourth-brightest object in the class, Haberl et al. (2003) discovered a broad absorption feature in XMM data. It extends from $`0.5`$keV to lower energies, and can be described by a Gaussian centred at $`<0.3`$keV. Furthermore, the third-brightest object, RX J1605.3+3249, was found to show a weaker, less wide absorption feature centred at a somewhat higher energy of $`0.45`$keV (Van Kerkwijk et al. 2003), as well as a narrower, marginally significant absorption line at 0.55 keV.

In Fig. 1, the RGS spectra obtained with XMM for both sources are shown. In addition, the spectrum of RX J0720.4$``$3125 is shown. This source has a very similar temperature, $`kT90`$eV, but a featureless spectrum.

## 3. Proton Cyclotron Absorption and Vacuum Resonance

For RX J1308.3+2127, Haberl et al. (2003) suggested that the feature was due to proton cyclotron absorption. Given the energy of $`<0.3`$keV, the implied magnetic field strength is $`<5(1+z)\times 10^{13}`$G (where $`1+z=(12GM/Rc^2)^{1/2}1.3`$ is the gravitational redshift factor). This is not unreasonable, given the observed slow spin period, $`P=10.3`$s. From the temperature, the neutron star should be about half a million years old. Equating this to the characteristic age, one infers a current spin-down rate $`\dot{P}P/2t3\times 10^{13}`$, and a magnetic field of $`3.2\times 10^{19}(P\dot{P})^{1/2}6\times 10^{13}`$G, consistent with what is required.

The large width of the feature, $`\sigma _E/E>1/3`$, is as expected for proton cyclotron absorption, since the cyclotron energy scales linearly with the magnetic field strength, which will vary over the surface (by a factor two for a centred dipole). Furthermore, as mentioned by Haberl et al. (2003), the equivalent width is consistent with model calculations of Zane et al. (2001).

If we assume the absorption in RX J1605.3+3249 is due to proton cyclotron absorption as well, we infer $`B7(1+z)\times 10^{13}`$G. Unfortunately, no pulsations have been found – to a limit of 3% in the frequency range 0.001–800 Hz (Van Kerkwijk et al. 2003) – so we cannot verify this. Puzzling in this case, however, is the difference in strength and width of the features in the two sources. As mentioned, the temperatures are very similar, so that cannot be the explanation.

The fields strengths inferred above are in excess of the critical quantum electrodynamics field $`B_{\mathrm{QED}}=4.4\times 10^{13}`$G, at which the electron cyclotron energy equals the electron rest mass. One effect that may become important is that photons propagating down the density gradient in the atmosphere can change polarisation mode at “vacuum resonance,” where the plasma contribution to the dielectric properties is compensated by the QED effect of vacuum polarisation.

Vacuum resonance has recently been studied in detail (e.g., Lai & Ho 2003; see also Lai, these proceedings), and it was found to lead to a reduction in the contrast of spectral features when it occurs between the deeper photosphere for the extraordinary mode photons and the shallower one for the ordinary mode photons. For the relevant energies of $`0.6`$keV at the surface, this will be the case for magnetic fields in the range 0.7–$`50\times 10^{14}`$ G.

If this is correct, it will not be very important for RX J1308.6+2127, and it makes sense that that source’s feature could be reproduced by Zane et al.’s 2001 models, which do not take vacuum polarisation into account. It should affect RX J1605.3+3249, however. And indeed, it may explain why the observed feature is so narrow: we might be seeing only absorption from regions with relatively low field, $`B<9\times 10^{13}`$G, the contrast of the absorption in regions with higher field being reduced due to the vacuum resonance.

## 4. Future Prospects

The discovery of absorption features should help further theoretical work on the vacuum resonance (in progress; Ho & Lai 2003b). In particular, in both sources, the absorption extends up to $`0.5`$keV. Might it be that the maximum energy out to which absorption is seen is set by vacuum resonance? If so, the observations provide a direct measurement of the critical magnetic field strength: about $`7(1+z)\times 10^{13}`$G.

Theoretical work is also needed on neutral hydrogen, which is strongly bound in high magnetic fields (for a review, Lai 2001). Hence, it should be present, and could lead to observable features. So far, it has only been studied at lower fields (e.g., Pavlov & Meszaros 1993; Zavlin & Pavlin 2002). It might possibly also help understand the puzzling overall spectral energy distributions.

Observationally, the obvious steps are a survey of the remaining sources, and using Chandra/LETGS to obtain spectra over a larger wavelength range and with phase resolution. Timing and astrometric studies could help constrain magnetic field strengths, distances, space velocities, places of origin, and ages.

All in all, the future looks bright, even if we may learn more about QED than about QCD.

### Acknowledgments.

I thank Dong Lai and George Pavlov for the good discussions during the conference, during which the ideas described here arose.

## References

Burwitz, V., Zavlin, V. E., Neuhäuser, R., Predehl, P., Trümper, J., & Brinkman, A. C. 2001, A&A, 379, L35

Braje, T. M. & Romani, R. W. 2002, ApJ, 580, 1043

Drake, J. J. et al. 2002, ApJ, 572, 996

Haberl, F. 2003, in High Energy Studies of Supernova Remnants and Neutron Stars, ed. W. Becker & W. Hermsen, in press \[astro-ph/0302540\]

Haberl, F., Schwope, A. D., Hambaryan, V., Hasinger, G., & Motch, C. 2003, A&A, 403, L19

Ho, W. C. G. & Lai, D. 2003, ApJ, submitted

Lai, D. 2001, Reviews of Modern Physics, 73, 629

Lai, D. & Ho, W. C. G. 2003, ApJ, 588, 962

Paerels, F. et al. 2001, A&A, 365, L298

Pavlov, G. G. & Meszaros, P. 1993, ApJ, 416, 752

Pavlov, G. G., Zavlin, V. E., & Sanwal, D. 2002, in Neutron Stars, Pulsars, and Supernova Remnants (Garching: MPE), 273 \[astro-ph/0206024\]

Treves, A., Turolla, R., Zane, S., & Colpi, M. 2000, PASP, 112, 297

Van Kerkwijk, M. H., Kaplan, D. L., Durant, M., Kulkarni, S. R., Paerels, F. 2003, ApJ, submitted

Zane, S., Turolla, R., Stella, L., & Treves, A. 2001, ApJ, 560, 384

Zavlin, V. E. & Pavlov, G. G. 2002, in Neutron Stars, Pulsars, and Supernova Remnants (Garching: MPE), 263 \[astro-ph/0206025\]


text:
# Nonlocal potentials for short-range electronic correlation in atoms, molecules, and solids

## I Introduction

Without invoking the electronic density as an intermediate, variational equations for an independent-electron model can be derived as an orbital functional theory (OFT), in agreement with Kohn-Sham equations in the local density approximation (LDA). In contrast to density functional theory (DFT), an exact energy functional for OFT can be derived from many-body theory. The orbital Euler-Lagrange equations (OEL) of OFT are determined by orbital functional derivatives, which take the general form of linear operators acting on occupied orbital functions of a model state. Because the exclusion principle requires independent normalization of the orbital partial densities, derivation of the OEL equations requires orbital wave functions or densities to be varied independently. For N electrons, there are N normalization constraints. Independent Lagrange multipliers (the orbital energy eigenvalues) are determined by these normalization constraints, as they are for Hartree-Fock and Kohn-Sham equations.

In a density-based theory that incorporates the exclusion principle, Euler-Lagrange equations for the orbital densities are determined by partial (Gâteaux) density functional derivatives, indexed by the orbital densities. Unless they are independent of the orbital index, these functional derivatives do not determine a unique total (Fréchet) density functional derivative, equivalent to a multiplicative local potential function. The required Gâteaux functional derivatives can be constructed from the corresponding orbital functional derivatives of OFT. The implied generalized Thomas-Fermi equations for orbital partial densities are operationally equivalent to the OEL equations.

Unless all orbital energies are equal, the OEL equations cannot be derived by a theory which considers only the total electronic density. Equal orbital energies conflict with the exclusion principle for the lowest state of any compact system with more than one electron of each spin. This is easily verified in the example of the $`1s2s^3S`$ state of an atom with two noninteracting electrons of the same spin. Independent variation of orbital densities, required by the exclusion principle, defines Gâteaux functional derivatives. It follows from this analysis that extension of Kohn-Sham theory beyond the LDA cannot be done in general without introducing nonlocal potentials. The required nonlocal potentials are well-defined in OFT, formulated as an application of the many-body theory of electrons.

The Schrödinger kinetic energy operator $`\widehat{t}=\frac{1}{2}^2`$ and the Fock exchange operator are well-known nonlocal ”potentials” in standard theory. The present work includes nonlocal potentials for short-range correlation based on the correlation energy Ansatz of Colle and Salvetti. The plan being followed is to implement these nonlocal potentials for calculations within atomic cells, then to use energy-linearized variational multiple scattering theory to match such local calculations together in large molecules and solids. This methodology is outlined here. Unrestricted Hartree-Fock (UHF) theory can be implemented for atoms or atomic cells by modifying any standard local-potential atomic Slater or Kohn-Sham program, to replace local exchange by the nonlocal Fock operator. A similar modification allows inclusion of the correlation energy functional considered here. Total and orbital energies computed with this nonlocal exchange-correlation model are reported for ground states of atoms He through Ne.

## II Nonlocal potentials in orbital functional theory

The simplest example of a nonlocal functional derivative is provided by the Schrödinger kinetic energy orbital functional

$`T={\displaystyle \underset{i}{}}n_i(i|\widehat{t}|i).`$ (1)
The notation here sums over occupied orbitals $`\varphi _i`$ of a Slater-determinant reference state $`\mathrm{\Phi }`$, with occupation numbers $`n_i=1`$ for $`iN`$. Spin indices and sums are assumed, but suppressed in the notation. The total electronic density of the model state $`\rho (𝐫)=_in_i\rho _i(𝐫)`$ is a sum of orbital densities $`\rho _i(𝐫)=\varphi _i^{}(𝐫)\varphi _i(𝐫)`$. The partial (Gâteaux) density functional derivative

$`{\displaystyle \frac{\delta T}{n_i\delta \rho _i}}=v_{ti}(𝐫)={\displaystyle \frac{\varphi _i^{}\widehat{t}\varphi _i}{\varphi _i^{}\varphi _i}}.`$ (2)
is determined as an orbital-indexed local potential function by the orbital functional derivative

$`{\displaystyle \frac{\delta T}{n_i\delta \varphi _i^{}}}=\widehat{t}\varphi _i.`$ (3)
Applied to the functional $`V=_in_i(i|v|i)`$, for an external local potential $`v(𝐫)`$, this analysis determines the Gâteaux derivative $`v_{vi}(𝐫)=v(𝐫)`$. Since this is independent of the orbital index, it reduces to a Fréchet derivative, a multiplicative local potential. The implied generalized Thomas-Fermi equations for orbital densities $`iN`$ of noninteracting electrons are

$`v_{ti}(𝐫)=ϵ_iv(𝐫).`$ (4)
The Lagrange multipliers $`ϵ_i`$ are to be determined such that $`d^3𝐫\rho _i=1`$ for each orbital density. These equations are operationally equivalent to the noninteracting Kohn-Sham or OEL equations

$`\widehat{t}\varphi _i=\{ϵ_iv(𝐫)\}\varphi _i.`$ (5)
Two-electron functionals are defined by $`U=E_h+E_x`$, where

$`E_h={\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j}{}}n_in_j(ij|u|ij);E_x={\displaystyle \frac{1}{2}}{\displaystyle \underset{i,j}{}}n_in_j(ij|u|ji),`$ (6)
and $`u=1/r_{12}`$. The Hartree functional $`E_h`$ is an integral of an explicit function of total density, which determines the Coulomb potential $`v_h(𝐫)`$ as a Fréchet derivative. Similarly, if $`E_{xc}`$ is defined as in the LDA, the present analysis determines a Fréchet functional derivative, and verifies the LDA Kohn-Sham equations. The orbital functional derivative of the exchange functional $`E_x`$ defines the Fock exchange operator $`\widehat{v}_x`$ such that $`\widehat{v}_x\varphi _i=_jn_j(j|u|i)\varphi _j`$. The functional derivative of $`U=E_h+E_x`$ defines $`\widehat{u}=v_h(𝐫)+\widehat{v}_x`$. Explicitly,

$`{\displaystyle \frac{\delta U}{n_i\delta \varphi _i^{}}}=\widehat{u}\varphi _i={\displaystyle \underset{j}{}}n_j(j|\overline{u}|j)\varphi _i,`$ (7)
where $`\overline{u}=u(1𝒫_{12})`$, and $`𝒫`$ is the exchange operator for equal spin indices.

Given $`(HE)\mathrm{\Psi }=0`$ for an N-electron eigenstate and any rule $`\mathrm{\Psi }\mathrm{\Phi }`$ that determines a model state $`\mathrm{\Phi }`$, unsymmetric normalization $`(\mathrm{\Phi }|\mathrm{\Psi })=(\mathrm{\Phi }|\mathrm{\Phi })=1`$ implies $`E=(\mathrm{\Phi }|H|\mathrm{\Psi })=E_0+E_c`$. Here $`E_0=(\mathrm{\Phi }|H|\mathrm{\Phi })=T+U+V`$ is an explicit orbital functional, and $`E_c=(\mathrm{\Phi }|H|\mathrm{\Psi }\mathrm{\Phi })`$ defines the correlation energy. If $`𝒬=I\mathrm{\Phi }\mathrm{\Phi }^{}`$, $`E_c=(\mathrm{\Phi }|H|\mathrm{\Psi }\mathrm{\Phi })=(\mathrm{\Phi }|H|𝒬\mathrm{\Psi })`$. This implies an exact but implicit orbital functional

$`E_c=(\mathrm{\Phi }|H[𝒬(HE_0E_ci\eta )𝒬]^1H|\mathrm{\Phi }),`$ (8)
for $`\eta 0+`$. In practice, some parametrized approximate $`E_c`$ must be used, defining a correlation potential operator such that

$`{\displaystyle \frac{\delta E_c}{n_i\delta \varphi _i^{}}}=\widehat{v}_c\varphi _i.`$ (9)
This defines an indexed local potential, the Gâteaux derivative

$`{\displaystyle \frac{\delta E_c}{n_i\delta \rho _i}}=v_{ci}(𝐫)={\displaystyle \frac{\varphi _i^{}\widehat{v}_c\varphi _i}{\varphi _i^{}\varphi _i}}.`$ (10)
Defining universal functional $`F=EV=T+U+E_c`$, and an operator $`=\widehat{t}+\widehat{u}+\widehat{v}_c`$, the general OEL equations are

$`{\displaystyle \frac{\delta F}{n_i\delta \varphi _i^{}}}=\varphi _i=\{ϵ_iv(𝐫)\}\varphi _i,`$ (11)
reducing in the LDA to Kohn-Sham equations. The OEL equations imply generalized Thomas-Fermi equations,

$`{\displaystyle \frac{\varphi _i^{}\varphi _i}{\varphi _i^{}\varphi _i}}=v_{fi}(𝐫)=ϵ_iv(𝐫).`$ (12)
If $``$ is hermitian, the indexed potential $`v_{fi}`$ is the Gâteaux derivative $`\delta F/n_i\delta \rho _i`$. There is no implication in general that $`E_x+E_c`$ defines a Fréchet derivative.

## III Plan for large molecules and solids

The general case of nonlocal potentials determined explicitly by the idempotent Dirac density matrix

$`\widehat{\rho }(1,2)={\displaystyle \underset{i}{}}\varphi _i(𝐫_1)n_i\varphi _i^{}(𝐫_2)`$ (13)
is considered. In full-potential multiple scattering theory (MST), local basis functions are constructed by integrating the Schrödinger or semirelativistic equation for specified orbital angular momentum $`\mathrm{}`$ and energy $`ϵ`$ within the enclosing sphere ($`r=r_1`$) of each atomic cell. Following ideas of canonical energy-band theory, reviewed by Skriver, the energy-dependent radial wave function $`u_{\mathrm{}}(ϵ;r)`$ is characterized by its logarithmic derivative $`D_{\mathrm{}}(ϵ)=r_Su_{\mathrm{}}^{}(r_S)/u_{\mathrm{}}(r_S)`$, evaluated at the radius of a sphere whose volume equals that of a polyhedral atomic cell. As discussed originally by Wigner and Seitz, this cellular wave function can continue smoothly across a cell interface if $`D_{\mathrm{}}(ϵ)`$ is negative, implying that $`D_{\mathrm{}}(ϵ_B)=0`$ and $`D_{\mathrm{}}(ϵ_A)=\mathrm{}`$ define the lower and upper energy limits of an energy band. The band center is estimated to occur at $`D_{\mathrm{}}(ϵ_C)=\mathrm{}1`$. An equivalent parameter is $`p_{\mathrm{}}(ϵ)=\frac{D_{\mathrm{}}(ϵ)+\mathrm{}+1}{D_{\mathrm{}}(ϵ)\mathrm{}}`$, which varies nearly linearly over the width of a band. It is found that orbital wave functions are well-approximated by linear interpolation over energies $`ϵ_Bϵϵ_A`$.

This behavior indicates that an initial self-consistent calculation should be carried out for spherically averaged potentials, subject to the orbital boundary condition $`D_{\mathrm{}}(ϵ)=\mathrm{}1`$ at $`r_S`$, in each inequivalent atomic cell. Because they satisfy a fixed boundary condition, the set of eigenfunctions $`u_{\mathrm{}}(ϵ_C;r)`$ can be extended to completeness within $`r_S`$, defining a local basis set for the linearized variational cellular method (LVCM). The computed nonlocal potential must be extended out to the enclosing sphere $`r=r_1`$, using the density matrix $`\widehat{\rho }`$ constructed from the self-consistent cell orbitals, normalized to unity within $`r_S`$. Orbital functions should be computed, using this fixed potential, for energies $`ϵ_C`$. Energy-derivative functions $`\dot{u}_{\mathrm{}}(ϵ_C;r)`$ must be computed as a basis for expanding the LVCM global matching function.

At the UHF level of calculation, in each self-consistent iteration an average local exchange potential (Slater exchange) is defined by

$`v_x(𝐫)={\displaystyle \underset{i}{}}n_iv_{xi}(𝐫)\rho _i(𝐫)/{\displaystyle \underset{i}{}}n_i\rho _i(𝐫),`$ (14)
an exact formula if $`v_{xi}`$ reduces to a Fréchet derivative. This approximation is corrected by incremental inhomogeneous terms

$`\{\widehat{v}_xv_x\}\varphi _i`$ (15)
evaluated from the previous iteration. This procedure converts a local-potential algorithm into UHF, and is valid for the indexed correlation potential considered below.

## IV Short-range correlation energy

Electronic correlation energy arises from two quite different sources. At distances larger than an atomic radius, multipolar response produces correlation effects evident in polarization potentials and dispersion forces. Long-range correlation, not considered here, requires computation of first-order multipole response pseudostates for each basis function $`u_{\mathrm{}}(ϵ_C;r)`$. A different approach is required for the short-range correlation due to the singularity of the interelectronic Coulomb potential $`u=1/r_{12}`$. Expansion in $`r_{12}`$ about such a singularity shows that an N-electron wave function must have specific cusp behavior in order to cancel the singularity. The wave function must vary as $`1+\frac{1}{2}r_{12}+\mathrm{}`$ . Colle and Salvetti (CS) impose this cusp condition through a symmetrical factor

$`\mathrm{\Pi }_{i<j}[1\xi (𝐫_i,𝐫_j)],`$ (16)
multiplying an antisymmetric model wave function. Using coordinates $`𝐪=𝐫_i𝐫_j`$ and $`𝐫=\frac{1}{2}(𝐫_i+𝐫_j)`$ the CS Ansatz is

$`\xi (𝐫,𝐪)=\mathrm{exp}(\beta ^2q^2)[1\mathrm{\Gamma }(𝐫)(1+{\displaystyle \frac{1}{2}}q)].`$ (17)
This ensures the limiting forms

$`(1\xi )|_{q0}`$ $`=`$ $`\mathrm{\Gamma }(𝐫)(1+{\displaystyle \frac{1}{2}}q+\mathrm{})`$ (18)
$`\xi |_q\mathrm{}`$ $`=`$ $`0.`$ (19)
In adapting this Ansatz to orbital functional theory (OFT), it is desirable to retain the unsymmetric normalization condition $`(\mathrm{\Phi }|\mathrm{\Psi })=(\mathrm{\Phi }|\mathrm{\Phi })=1`$, where $`\mathrm{\Psi }`$ is the correlated state, and $`\mathrm{\Phi }`$ is a reference state Slater determinant. Then the correlation energy is $`E_c=(\mathrm{\Phi }|H|\mathrm{\Psi }\mathrm{\Phi })=_{i<j}n_in_jE_{ij}`$, a sum of electron-pair correlation energies. This suggests parametrization for each pair of occupied orbital functions $`\varphi _i,\varphi _j`$, such that

$`E_{ij}=(ij|\overline{u}\xi _{ij}(q)|ij).`$ (20)
A parametrized form similar to CS is

$`\xi _{ij}(q)=\mathrm{exp}(\beta _{ij}^2q^2)[1\gamma _{ij}(1+{\displaystyle \frac{1}{2}}q)].`$ (21)
Parameter $`\gamma _{ij}`$ is determined by the normalization condition $`(\mathrm{\Phi }|\mathrm{\Psi }\mathrm{\Phi })=0`$, or $`(ij|\xi _{ij}|ij)=0`$ for each pair:

$`\gamma _{ij}={\displaystyle \frac{(ij|\mathrm{exp}(\beta _{ij}^2q^2)|ij)}{(ij|\mathrm{exp}(\beta _{ij}^2q^2)(1+\frac{1}{2}q)|ij)}}.`$ (22)
The free parameter $`\beta _{ij}`$ can be chosen to minimize Bethe-Goldstone (BG or IEPA) energy for specified $`i,j`$, or can be treated as a semiempirical parameter and fitted to known correlation energies of atoms.

This Ansatz can be incorporated into OFT by using Eq.(20) to define the model correlation energy. At any stage of the self-consistency iteration, parameter $`\gamma _{ij}`$ is determined by the consistency condition given above, and $`\beta _{ij}`$ is either a fixed empirical parameter, or is to be updated by solving a 2-electron BG equation indexed by orthogonal occupied orbitals $`ij`$. The indexed local correlation potential (Gâteaux functional derivative) is

$`v_{ci}(𝐫)={\displaystyle \frac{\varphi _i^{}(𝐫)\widehat{v}_c\varphi _i(𝐫)}{\varphi _i^{}(𝐫)\varphi _i(𝐫)}},`$ (23)
where

$`\widehat{v}_c\varphi _i={\displaystyle \frac{\delta E_c}{n_i\delta \varphi _i^{}}}={\displaystyle \underset{j}{}}n_j(j|\overline{u}\xi _{ij}|j)\varphi _i,`$ (24)
a direct generalization of the operator $`\widehat{u}`$. Both $`v_h`$ and $`\widehat{v}_x`$ are modified by short-range correlation, and antisymmetry is built in. Because terms $`j=i`$ vanish, this Ansatz for Coulomb-cusp correlation does not produce self-interaction. The indexed local potential $`v_{ci}(𝐫)`$ is singular at nodes of $`\varphi _i(𝐫)`$. As in Hartree-Fock methodology, $`\widehat{v}_c\varphi _i`$ should be treated as an inhomogeneous term in numerical solution of the orbital Euler-Lagrange (OEL) equations.

Coupling of electron pairs is significant in standard CI methods for electronic correlation. Although such coupling is inherent in the self-consistency of the proposed mean-field model, a more accurate extension of the method may be needed. A possible procedure is to implement the coupled electron-pair approximation (CEPA) within each atomic cell, using the modified CS Ansatz as a closure formula for the CI expansion.

### A Notes on integrals

The normalized orbital basis functions are of the form

$$\varphi _a(𝐫)=N_a\chi _a(r)Y_{\mathrm{}_am_a}(\theta ,\varphi ),$$
(25)
where $`r\chi _a(r)`$ is a numerical solution of the radial OEL equation. The normalization constants are $`N_a=[r^2𝑑r\chi _a^2(r)]^{\frac{1}{2}}`$. Definite and indefinite integrals are required for two-electron generalized potential functions $`F(q)`$, where $`q^2=r_1^2+r_2^22r_1r_2\mathrm{cos}\theta `$. Integration over angles follows the standard derivation of Condon and Shortley for $`F(q)=1/q`$. For normalized radial functions,

$`(ab|F|cd)={\displaystyle \underset{k}{}}c^k(ac)c^k(db)F^k(ac;db),`$ (26)
where $`c^k`$ denotes a Gaunt coefficient and

$`F^k(ac;db)={\displaystyle _0^{\mathrm{}}}r_1^2𝑑r_1{\displaystyle _0^{\mathrm{}}}r_2^2𝑑r_2`$ (27)
$`f^k(r_1,r_2)\chi _a^{}(r_1)\chi _c(r_1)\chi _d(r_2)\chi _b^{}(r_2).`$ (28)
The factor $`f^k`$ is

$`f^k(r_1,r_2)={\displaystyle \frac{2k+1}{2}}{\displaystyle _0^\pi }\mathrm{sin}\theta d\theta P_k(\mathrm{cos}\theta )F(q).`$ (29)
Given $`r_1`$ and $`r_2`$, $`\mathrm{cos}\theta =\frac{r_1^2+r_2^2q^2}{2r_1r_2}`$ and $`\mathrm{sin}\theta d\theta =\frac{qdq}{r_1r_2}`$. In internal coordinates, $`f^k(r_1,r_2)`$ is an integral of the form

$`{\displaystyle \frac{2k+1}{2r_1r_2}}{\displaystyle _{|r_1r_2|}^{r_1+r_2}}P_k({\displaystyle \frac{r_1^2+r_2^2q^2}{2r_1r_2}})qF(q)𝑑q.`$ (30)
These integrals are needed for $`qF(q)`$ equal to $`1,q,q^2`$ times the factor $`\mathrm{exp}(\beta ^2q^2)`$, so that the integrand of $`f^k`$ is this Gaussian factor times a polynomial in $`q`$. It can easily be verified that $`f^k`$ reduces to $`r_<^k/r_>^{k+1}`$ for $`\beta 0`$ if $`qF(q)=1`$. The Legendre polynomial factors of the integrand are determined by the recurrence formula

$`P_{k+1}(x)={\displaystyle \frac{2k+1}{k+1}}xP_k(x){\displaystyle \frac{k}{k+1}}P_{k1}(x),`$ (31)
with $`P_1(x)=0,P_0(x)=1`$. The elementary integrals required are of the form $`I_n(\beta ;\kappa )=_0^\kappa e^{\beta ^2q^2}q^n𝑑q`$. By a change of variables such that $`t=\beta ^2q^2`$, this reduces to

$`I_n(\beta ;\kappa )`$ $`=`$ $`{\displaystyle \frac{1}{2\beta ^{n+1}}}{\displaystyle _0^{\beta ^2\kappa ^2}}e^tt^{\frac{1}{2}(n1)}𝑑t`$ (32)
$`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }(\frac{1}{2}(n+1))}{2\beta ^{n+1}}}p({\displaystyle \frac{1}{2}}(n+1),\beta ^2\kappa ^2),`$ (33)
in terms of the incomplete gamma function $`p(a,x)=1\mathrm{\Gamma }(a,x)/\mathrm{\Gamma }(a)`$, Sect.6.2. $`p(a,x)`$ can be computed effficiently using a power series for small $`x`$ and a continued fraction for large $`x`$. The continued fraction terminates if $`n`$ is odd.

## V Calculations on light atoms

For applications to molecules and solids, using multiple scattering theory and an atomic cell model, values of the parameters $`\beta _{ij}`$ can be obtained by calculations on atoms. Results of such calculations, for light atoms He through Ne, are reported here. A numerical Dirac-Slater program was modified as described above for UHF (exchange-only) and OFT calculations, the latter incorporating the modified CS correlation energy functional described above. The program was used in its nonrelativistic mode. Angular coefficients in the total energy functional were computed such that effective potentials are spherically averaged, but retain a spin index. This equivalence restriction implies that radial orbital functions with indices $`\mathrm{},m_s`$ are well-defined.

For each atom considered, three sets of self-consistent calculations were carried out: for the He-like ion, for the Be-like ion (for $`N4`$), and for the neutral atom. To verify the computational method, computed UHF energies are compared in Table(I) with established RHF energies and with total energies including correlation. The OFT calculations were used to determine $`\beta _{ij}`$ parameters such that the computed total energies agreed with the ”experimental” values shown in Table(I). Pair-indexed parameters $`\beta _{1s1s}`$ were determined for He-like ions, parameters $`\beta _{2s2s}`$ for Be-like ions, and parameters $`\beta _{2p2p}`$ for the neutral atoms. In each case, inner parameters were frozen and intershell parameters such as $`\beta _{1s2s}`$ were scaled to the geometric mean of the corresponding intrashell parameters. Values of the latter that fit total energies of the ions and atoms considered are listed in Table(II). These parameters are the principal result of the present calculations, intended to define parametrized correlation functionals for extended systems.

Self-consistent UHF orbital energies are tabulated in Table(III) and OFT orbital energies in Table(IV). Because Janak’s theorem is valid in OFT, these energies have a physical meaning in the context of a theory in which orbital occupation numbers are allowed to change continuously and to have fractional values. They are derivatives of the total energy with respect to infinitesimal changes of the occupation numbers. Physical energy differences, by implication, correspond to integrals of these derivatives over finite increments of occupation numbers.

## VI Conclusions

This paper considers an independent-electron model that incorporates a theoretically motivated Ansatz for correlation energy, expressed as a parametrized orbital functional. Computed results extend exchange-only theory (UHF) to a formalism parametrized by exact atomic ground-state energies. Parameters are obtained that make it possible to apply this formalism to calculations of the electronic structure of molecules and solids.

## Acknowledgments

This work was initiated at the University of Connecticut, supported by a grant from the University of Connecticut Research Foundation. The author is grateful to H. H. Michels for discussions of the new methodology proposed here, and to Prof. Michels and Prof. Wm. C. Stwalley for encouraging and implementing a visiting appointment at Storrs.


text:
# Space-Time Properties of Hadronization from Nuclear Deep Inelastic Scattering

## 1 Introduction

Hadronization has been a ubiquitous feature of high energy physics experiments for more than two decades. The subculture that has evolved around the study of this process has succeeded in characterizing the topological features of high energy jets, such as the fragmentation functions that describe the probabilities of quarks of a given flavor evolving into particular hadrons. The string picture is an example of a successful phenomenology that captures the important features over a wide range of energies. Few aspects of hadronization, however, are presently understood at the most fundamental level, since it is a non-perturbative process.

One important type of experimental data that is generally not easily accessible for hadronization is its space-time development. While some phenomenological models can be used to predict this development, the process takes place on a microscopic distance scale that is not directly observable. However, this type of information can be inferred by implanting the hadronization process into a nucleus, which acts as a spatial filter. Fully formed hadrons in the nuclear medium interact via ordinary hadronic cross sections; by contrast, propagating quarks do not appear to interact in a way that modifies hadron yields significantly. Therefore, much can be learned from studying hadron production in deep inelastic scattering kinematics (DIS) from nuclei of varying sizes. The degree of modification of fragmentation functions, for instance, can be studied as a function of nuclear size A, energy transfer $`\nu `$, hadron energy fraction $`z=E_h/\nu `$, and the momentum component perpendicular to $`\stackrel{}{q}`$, which is labelled $`p_T`$. These studies can yield important insights into the fundamental processes driving hadronization, such as gluon emission and the temporal evolution of the truncated color field following hard interactions.

Several previous measurements of semi-inclusive leptoproduction of hadrons from nuclei have been performed. The earliest generation of these were carried out with electron beams at SLAC and muon beams at CERN and FNAL. Prominent features displayed by these efforts are an attenuation of hadrons which increases with A and vanishes with increasing $`\nu `$. Following these early efforts, recent data from the HERMES experiment at DESY have inspired a new wave of theoretical interest. Unlike the historical efforts, these new data offer excellent hadron identification and high precision, allowing the flavor dependence of hadron formation to be studied.

Finally, existing and future measurements at Thomas Jefferson National Accelerator Facility (Jefferson Lab) will offer unique capabilities to study hadron production in DIS. Measurements at the CEBAF Large Acceptance Spectrometer (CLAS) at Jefferson Lab will provide data on the widest possible range of nuclear target masses at high luminosity. The first data at 5 GeV were taken in 2003 with sustained luminosities approaching $`10^{34}cm^2s^1`$, while 11 GeV data are planned for the future, accompanied by an order-of-magnitude luminosity increase and improved particle identification using the upgraded CLAS (“CLAS++”).

The scientific topics connected with these studies are also important themes in other physics communities. Quark propagation through strongly interacting media is a foundational element in predicting signatures of the quark gluon plasma at RHIC through the phenomenon of jet quenching; in this connection, it has been explored in ultrarelativistic d-A collisions. It has also been investigated in the Drell-Yan process at FNAL in studies of quark energy loss in the nuclear medium.

## 2 Experimental Tools

A primary experimental tool for nuclear hadron formation is the *hadronic multiplicity ratio:*

$`R_M^h(z,\nu ,p_T^2,Q^2,\varphi )={\displaystyle \frac{\left\{\frac{N_h^{DIS}(z,\nu ,p_T^2,Q^2,\varphi )}{N_e^{DIS}(\nu ,Q^2)}\right\}_A}{\left\{\frac{N_h^{DIS}(z,\nu ,p_T^2,Q^2,\varphi )}{N_e^{DIS}(\nu ,Q^2)}\right\}_D}}`$ (1)
where $`N_h^{DIS}`$ and $`N_e^{DIS}`$ denote the number of hadrons and electrons measured in DIS kinematics, $`Q^2`$ is the four-momentum transfer, and $`\varphi `$ is the angle between the lepton scattering plane and the photon-hadron plane. The primary connection to theory has traditionally been through $`R_M^h`$ interpreted in a partonic framework, for particular flavor combinations. Additional tools include *energy and angular correlations* between outgoing hadrons involved in the reaction. There are at least two classes that may be of interest: correlations between two energetic hadrons formed directly from the hadronization of the struck quark, and correlations between the energetic hadrons and very low energy protons. The $`p_T`$ distribution is expected to become broader for larger nuclei; this *transverse momentum broadening* is expected to be large enough to be directly measurable. In addition to these quantities, *polarization observables* such as single spin asymmetries may also play a role.

A variety of physical properties can be derived from the above observables. Hadron *formation lengths*111In some models, referred to as “production length,” the distance to produce a “pre-hadron” in a two-step picture of hadron formation., the characteristic distances over which hadrons form, can be extracted from $`R_M^h`$. *Quark energy loss* in passing through the nuclear medium may also be accessible, particularly via transverse momentum broadening, which is also closely connected with a *quark-gluon correlation function*. Correlations between final-state particles are expected to shed light on the reaction mechanisms, and the nature of the interaction with the nuclear medium, whether partonic or hadronic.

A variety of experimental capabilities are required in order to fully exploit the observables discussed above. Access to the detailed multivariate dependences of $`R_M^h`$ requires a large statistical sample; this can be most economically accessed at high luminosity. Isolating the energy and angular correlations requires large acceptance, since the kinematic variables of the particles of interest are not strongly correlated. Good particle identification for several hadron species is required to analyze the flavor dependence of observables. A large energy range is needed to understand the $`\nu `$ dependence as well as test the validity of parton model assumptions. Polarization is needed to study such quantities as single spin asymmetries. Accessing all of these experimental ingredients will require combining information from several different facilities.

## 3 Model Descriptions and Contact with Observables

Lattice computations can address nonperturbative processes in QCD. However, there are several difficulties in applying this technique to hadronization. A natural choice to describe the time development of the process is Minkowski space, however, this introduces some significant technical obstacles. In addition, the relevant space-time volume needed requires a very large lattice compared to calculating static properties of hadrons. These calculations may become feasible in the 5-10 year time frame.

Therefore, to make progress on this topic, modeling is needed. Most relevant are models possessing close contact to QCD, and which can predict a wide variety of observables. A synopsis of model approaches is given in the following section.

### 3.1 Models

Historically, a number of models have been developed to describe hadron production in nuclear DIS. While it is not practical to review the older models here, a representative approach by Bialas and Chmaj is instructive. In this model, the struck quark evolves into a hadron with a characteristic time $`\tau `$. The probability that the propagating object is a quark at a given time $`t`$ after it is struck is $`P_q=e^{t/\tau }`$ and the corresponding probability that it has become a hadron is simply $`P_h=1P_q`$. Phenomenological cross sections for the quark and hadron interacting with the nuclear medium are postulated: $`\sigma _q`$ and $`\sigma _h`$. The probability of the quark or hadron interacting with the nuclear medium is taken as proportional to the density of the medium for a spherical nucleus, using a standard density parameterization. $`R_M^h`$ can then be calculated by numerical integration. This model approach was qualitatively successful in describing the $`\nu `$ dependence of $`R_M^h`$ for the HERMES nitrogen data, but failed to describe the $`z`$ dependence even qualitatively.

In a more sophisticated approach, the gluon bremsstrahlung model divides the hadronization process into several stages. As the struck quark propagates, an initial stage of uninhibited gluon emission occurs, followed by a second stage of reduced gluon emission. The suppression of gluon emission in the second stage is imposed to conserve energy, since the propagating quark begins with only a finite amount of energy, $`\nu `$. At some point in the second stage, a color dipole is formed from the struck quark and the last emitted gluon. The color dipole subsequently evolves into a hadron. The role of the nuclear medium is to modulate the rate at which the color dipole turns into a hadron, including the effect of color transparency; in later versions of the model, medium-stimulated gluon emission is included, which has a small effect on the prediction for $`R_M^h`$. The assumptions of the model are valid for pions with $`z>0.5`$. Its predictions agree well with the HERMES nitrogen data in both the $`z`$ dependence and the $`\nu `$ dependence.

A more recent model approach to describing hadron production in nuclear DIS is based on a continuous series of developments since the late 1980’s addressing jet quenching in relativistic heavy ion reactions. In this approach, the prediction for the behavior of $`R_M^h`$ is based solely on medium-stimulated gluon emission; hadronization is assumed to occur outside the nuclear medium. The calculation of $`R_M^h`$ is performed in pQCD; in addition to the leading order terms, the dominant twist-four term is also included. The twist-four term represents multiple scattering in the nuclear medium, and the contribution due to this term is a free parameter that can be fixed by comparison to data. This model was able to describe the HERMES krypton data by fitting the HERMES nitrogen data. The modification of the fragmentation functions represented by $`R_M^h(z,\nu )`$ emerges naturally from this description. In this approach, a novel quadratic dependence of the energy loss on path length is predicted, originating in the non-Abelian analog of the LPM effect in QED. However, the assumption that hadronization occurs outside the nucleus will break down at low energies or for more massive hadrons.

Another recent model invokes deconfinement of nucleons in the nuclear medium as a mechanism for the modification of fragmentation functions in the nuclear medium. The basic method follows the rescaling models developed in the 1980’s. This approach combines partial deconfinement with absorption of the formed hadron in the nuclear medium. The nucleonic deconfinement effectively extends the range in $`Q^2`$ for which gluon radiation takes place, modifying the fragmentation functions accordingly. The model successfully describes the EMC data and the HERMES data for pions and $`K^+`$ in nitrogen and krypton, although the prediction for $`K^{}`$ systematically deviates from the data.

A fourth recent model adapts and improves stationary string model (SSM) techniques developed in the 1980’s and 1990’s. This model, applied to valence quarks, successfully fits the HERMES pion data for nitrogen and krypton, and predictions are given for kaons and heavier nuclei. As with the previous model, the prediction for $`K^+`$ and $`K^{}`$ attenuations are very similar, while the HERMES data show a systematic difference.

Finally, an approach that emphasizes accurate modeling of the final state interactions (FSI) is able to explain the HERMES data without invoking modification of fragmentation functions. The struck quark propagates, evolves into a hadron, and its interaction with the nuclear medium is given by a semi-classical coupled channel transport calculation that has been tested in other reactions. These calculations, although in preliminary form, seem to have enough flexibility to reproduce the HERMES nitrogen and krypton data for charged pions and kaons, and protons and anti-protons. There is good qualitative agreement throughout, and quantitative agreement for some kinematics.

### 3.2 Connections between Models and Observables

The models highlighted above are generally able to describe the published data, in spite of the fact that they rely on quite different physical pictures. In order to make good progress in understanding hadron formation in nuclear DIS, it is necessary to discriminate among these models using more data that span a wider range of kinematics or that introduce new observables. Similarly, more complete calculations are needed. For instance, model predictions for the flavor and mass dependence of $`R_M^h(z,\nu )`$ can be tested by studying pions, kaons, and protons. In fact, a wide range of particles can be considered, as seen in Table 1. In this table, hadrons with $`c\tau `$ much greater than nuclear dimensions are listed. To the extent that hadronization is the production mechanism, inter-comparison of $`R_M^h(z,\nu )`$ for these particles will shed light on such topics as whether semi-classical calculations of FSI suffice, or whether a quantum mechanical description is required. Likewise, study of moments of transverse momentum may distinguish between models with different assumptions about the number of multiple scatters or the shape of the emitted gluon momentum distribution. Measurements of $`R_M^h`$ at low $`z`$ may be useful, particularly at low energies, to distinguish between pion flux entering other reaction channels, as in FSI, from a change in the intrinsic fragmentation function, i.e., due to gluon emission. A clean evaluation of the $`Q^2`$ dependence of $`R_M^h`$ will help distinguish between different reaction mechanisms. Studies of the angular correlations between two high-energy pions may also clarify whether the reaction mechanism is dominated by FSI in the nuclear medium or by the fragmentation process. Spin observables, such as single-spin asymmetries, may offer indications of whether partonic multiple scattering is an important component of the process. Finally, the twist-4 pQCD model makes the assumption that hadronization takes place outside the nucleus altogether. The point at which this assumption is no longer valid is an important unknown that can be addressed by the CLAS low-energy data.

While these are rich and exciting prospects for learning much new physics, there are a few experimental and theoretical complications that require further study and clarification. For instance, hadron formation lengths can be derived from the data, however, the definition of the formation length depends on which model is considered. As an example, hadron formation is conceived as a multi-step process with multiple time constants in some models, but a single-step process in others. Reference frame considerations also enter the interpretation of the data; a lower limit on x is required to preserve the simple picture that a single quark propagates with initial energy $`\nu `$. Resonances in the residual system may need to be suppressed, which reduces the kinematic coverage available. There may be an issue of isolating current fragmentation from target fragmentation, although the nuclear medium may perhaps suppress the latter. Finally, the degree of validity of QCD factorization needs to be evaluated. The last few issues are common to all semi-inclusive deep inelastic scattering, and as such are being investigated by several physics communities; further progress is therefore to be expected.

## 4 Future Prospects

Over the next few years it is expected that new datasets will become available which will provide further stringent tests of these models. Some HERMES data were taken with a 12 GeV beam on both nitrogen and krypton targets; 27 GeV data on helium and neon are also under analysis. A new program at Jefferson Lab will provide even lower energy data at high luminosity on targets ranging from carbon to lead. This 5 GeV data, from the CLAS EG2 running period, will have sufficient statistics to map out the full multivariate character of $`R_M^h`$ for several hadron species, and to directly measure transverse momentum broadening.

Other future prospects include the possibility of additional targets in HERMES such as xenon, to complete the study of nuclear mass dependence at high energy transfer, and a large program which would be possible at Jefferson Lab following the energy upgrade. The ultimate moderate-energy program will be feasible with CLAS++, due to the planned upgrade in luminosity (to $`10^{35}cm^2s^1`$) and improved particle identification capabilities. These improvements will permit mapping out multivariate $`R_M^h`$ for all the hadrons in Table 1 as well as $`p_T`$ broadening for many of these, spanning a kinematic range of $`\nu `$ from $`29GeV`$ and $`29GeV^2`$ in $`Q^2`$.

## 5 Conclusions

The study of hadronization from semi-inclusive nuclear DIS is a new topic that is drawing increasing interest from several different physics communities. The recent wave of interest stimulated by the HERMES measurements has resulted in a number of theoretical efforts to describe the data. These efforts, based on a variety of physical pictures, demonstrate the exciting possibility of gaining new insight into hadronization as well as other aspects of quark propagation through nuclei such as quark energy loss. However, more data is needed to definitively discriminate among these model approaches. In addition, much more theoretical effort is needed to address all of the recent data, to reduce the number of assumptions in the models, and to make greater contact with the fundamental theory, QCD.