text:
# Atomic Scale Imaging and Spectroscopy of a CuO2 Plane at the Surface of Bi2Sr2CaCu2O8+δ

\[

## Abstract

We have used a scanning tunneling microscope to demonstrate that a single CuO2 plane can form a stable and atomically ordered layer at the surface of Bi2Sr2CaCu2O8+δ. In contrast to previous studies on high-Tc surfaces, the CuO2-terminated surface exhibits a strongly suppressed tunneling conductance at low voltages. We consider a number of different explanations for this phenomena and propose that it may be caused by how the orbital symmetry of the CuO2 plane’s electronic states affects the tunneling process.

PACS numbers: 74.50.+r,73.20.Hb,74.72Hs

\]

High-temperature cuprate superconductors are layered compounds containing one or more CuO2 planes and other layers that act as charge reservoirs for these planes. It is believed that the strong electronic interactions in the CuO2 planes are responsible for the cuprates’ unusual electronic properties, including unconventional superconductivity with a $`d`$-wave order parameter and the non-Fermi liquid behavior found at low carrier concentrations . Consistent with this view, recent efforts have demonstrated that the phenomena observed in the bulk cuprates can be induced in a single CuO2 plane using charge injection devices fabricated at the surface of the insulating infinite-layer compound CaCuO2 . Therefore, it is both timely and of fundamental importance to study the electronic properties of a single CuO2 plane at the surface. Surfaces of most materials from simple metals to compound semiconductors exhibit a variety of novel electronic phenomena; however systematic studies of layered cuprate surfaces have been limited. In particular, a study of surfaces terminated with different layers has not been accomplished.

Here we report the first imaging and electron spectroscopy measurements of a single CuO2 plane at the surface of a cuprate superconductor using a scanning tunneling microscope (STM). Similar to most previous experiments that have probed the surfaces of cuprates , we carry out our studies on Bi2Sr2CaCu2O8+δ. However, previous efforts have focused almost exclusively on surfaces terminated with the BiO layer, 4.5Å above the pair of CuO2 planes in this compound. Although there have been some reports of STM experiments on surfaces terminated with layers other than BiO, these efforts have failed to definitively identify the CuO2 plane as the surface layer . Such an identification must rely on a successful correlation of STM topographic images with crystallographic data, which depends, for example, both on preparing the samples’ surface and on making measurements under ultra-high vacuum (UHV) conditions .

Our measurements demonstrate that a CuO2 plane in Bi2Sr2CaCu2O8+δ can form a stable and well-ordered surface layer with a lattice similar to that of the same plane in the bulk. In contrast to typical spectroscopic measurements on BiO-terminated surfaces, which exhibit an inchoate gap , we find that tunneling into a CuO2 plane at the surface is strongly suppressed at energies close to the Fermi level. This behavior, which is unexpected from a $`d`$-wave superconductor, raises the fundamental question of how a single CuO2 plane behaves when it is at the surface and asymmetrically doped by layers beneath the surface. We consider a number of possible explanations for this phenomena including a model based on the effects of CuO2 plane’s orbitals on the tunneling process, which can produce suppression of tunneling at low energies.

We performed our measurements using a home-built UHV STM operating at 4.2K on Bi2Sr2CaCu2O8+δ thin films (1000Å) which were grown using molecular beam epitaxy (MBE). Resistivity measurements show the samples to have a Tc of 84K and a transition width of 4K. The samples were introduced into the UHV chamber and cleaved mechanically at room temperature prior to making STM measurements. The majority of the areas of the nearly optimally doped samples show STM data similar to those previously obtained on the BiO-terminated surface of Bi2Sr2CaCu2O8+δ . This surface is typically exposed as a result of cleaving, a process in which the Bi2Sr2CaCu2O8+δ samples break at the weakest bond in the crystal structure, namely between the BiO planes of two adjacent half-unit cells . However, in addition to the BiO plane, we occasionally find smaller terraces terminated with crystallographic layers other than the BiO plane.

We focus on situations such as that shown in Fig. 1, where we find a BiO terrace ending in a sub-unit cell height step. In Fig. 2, we show the spectra taken along a line perpendicular to the step edge in Fig. 1a. Tunneling spectra measured on the top terrace (on the right, Fig. 1a) are characteristic of those measured throughout the majority of the surface (Fig. 3c), and hence the top terrace is identified as the BiO plane. The spectra measured on the bottom terrace (on the left, Fig. 1a) are remarkably different from those measured on the top terrace, even when taking into account the extremes in variation of the BiO spectra (Fig. 3c). To identify the atomic terrace on the left side of Fig. 1a, we compare the average of the topographic line scans measured perpendicular to the step edge with the known separation between the different crystal planes in Fig. 1b . We have estimated various sources of inaccuracy in relating STM topography to structural information, such as those due to piezo calibration or to the difference in the electronic structure of the two planes. All such factors are far smaller (0.1-0.2Å) than the expected distances between adjacent layers in the crystal structure ($``$2Å). Hence, we can correlate the STM data and plane spacing as shown in Fig. 1b and conclude that the terrace on the left is the lower CuO2 plane of Bi2Sr2CaCu2O8+δ. Using this procedure, we have identified the CuO2 plane as the terminating layer of distinct regions of several different samples grown under the same conditions. The coupling between CuO2 planes in adjacent half-unit cells along the crystallographic c-axis in Bi2Sr2CaCu2O8+δ is extremely weak. Our measurements thus represent an examination of the electronic structure of a single CuO2 plane that is asymmetrically doped via the layers beneath the surface.

Having identified a surface terminated with a CuO2 plane, we examine its electronic and structural properties far from step edges. The high-resolution images of this surface in Fig. 3b show a well-ordered lattice with an atomic spacing consistent with that expected for a CuO2 plane in the bulk crystal structure. The atomic lattice also shows a periodic distortion due to the incommensurate b-axis supermodulation known to run throughout the bulk of the Bi2Sr2CaCu2O8+δ crystal.  Spectroscopic measurements (Fig. 3d) made on the surface CuO2 plane show a wider energy gap as compared with those typically measured on the BiO surface, as well as a strong suppression of the tunneling conductance within 10 meV of the Fermi level EF. The unusual tunneling characteristics measured on this plane are found to be homogeneous on length scales from interatomic distances to several hundreds of angstroms (gap value 60$`\pm `$10 meV). The strong contrast between the spectra taken on the BiO and CuO2 terminated surfaces is the main result reported in this paper, which demonstrates the dramatic influence of the layered structure on their surface electronic properties.

Surface sensitive measurements on the BiO-terminated surfaces is commonly interpreted to reflect the density of state (DOS) of the ”bulk” CuO2 layers beneath surface. However, in the case of tunneling measurements performed on the BiO plane, the large separation between the tip and the pair of CuO2 planes (greater than 10Å) precludes the possibility of tunneling directly to these planes. In this situation, the electronic states of planes other than CuO2 planes play a role in the tunneling process and must be taken into account. Despite the lack of a clear model for this tunneling process, experimental support for interpretation of BiO spectra in terms of the DOS of the sub-surface CuO2 layers has come from strong resemblance of STM spectra with those obtained by angle-resolved photoemission (ARPES) at $`\stackrel{}{k}=(0,\pi )`$ on the same surface.  However, in this region of momentum space, the DOS of Bi2Sr2CaCu2O8+δ could be remarkably different from that expected for a single CuO2 plane. Electronic structure calculations for Bi2Sr2CaCu2O8+δ have long predicted the presence of BiO-derived bands and a splitting of CuO2 bands due to bilayer interactions at the Fermi level near $`\stackrel{}{k}=(0,\pi )`$.  These two effects may have to be included to interpret both the STM and the ARPES data on BiO-terminated surfaces. To date, the contribution of the BiO layer to the DOS near EF has not yet been observed; however, various ARPES experiments have recently accumulated strong evidence for the predicted bilayer splitting of the CuO2 bands.  In general, the exact nature of ARPES measurements of the band structure near $`(0,\pi )`$ continues to be a subject of intense debate.  The resolution of this debate is critical to the interpretation of tunneling spectra on the BiO surface and how they relate to our spectroscopic measurements of a single CuO2 plane at the surface.

In contrast to tunneling on the BiO plane, direct tunneling into the surface CuO2 plane does not involve the electronic structure of layers other than a single CuO2. Despite this simplification, to interpret our measurements we must consider the possibility that this surface layer may have a dramatically different electronic character than either the same plane in the bulk or those underneath the BiO-terminated surface. Such a situation may arise as a result of structural modifications or changes in the doping level when the CuO2 plane is at the surface. The STM images of the surface-terminated CuO2 plane rule out the possibility of a structural reconstruction, as they show this plane has the atomic structure expected for a bulk CuO2 plane in Bi2Sr2CaCu2O8+δ. However, we can not rule out the possibility that the doping level is different at the surface and hence consider various electronic states that can give rise to the measured spectra.

One possibility would be that the surface CuO2 plane is an insulator and the observed energy gap is associated with making single electron excitation in such a system. While we can not rule out this possible explanation, based on the measurements presented here we can place stringent constraints on possible insulating phases. First, the insulator must be electronically homogenous and have an energy gap of about 60 meV, a value unusually close to previously reported values of energy gap for highly underdoped (0.1 holes per Cu-O plaquette) but still superconducting samples. Second, to explain our data, the insulator must have a particle-hole symmetric electronic structure with the chemical potential exactly in between the occupied and unoccupied bands. For bulk cuprates, at very low doping, there is an antiferromagnetic insulating ground state that has a 700 meV Mott wide gap , which is inconsistent with our measurements reported here. In addition to this phase, there have been proposals for a variety of interesting insulating bulk phases, such as spin or charge density waves states at low doping Many of these models not only involve inhomogeneous electronic states but also are expected to give rise to an inchoate pseudogap.

An alternative possibility would be that the surface CuO2 plane is conducting, but the measured gap is due to Coulomb blockade, such as those observed when tunneling into disordered or granular thin films. Such a possibility can be diagnosed, as the size of the charging gap should be a sensitive function of capacitance. Our observation that the spectra is independent of tip-sample separation rules out charging effects as a potential explanation.

Finally, we consider the possibility that the surface CuO2 plane is in the superconducting state at low temperatures, although with a doping level that may be different than that of the bulk of our sample. Assuming a similar superconducting state as in the bulk, we would expect the density of states (DOS) of the surface layer at low energies to have a form $`\rho _s(E,\stackrel{}{k})=\rho _0E/\sqrt{E^2\mathrm{\Delta }^2(\stackrel{}{k})}`$, with an energy gap $`\mathrm{\Delta }(\stackrel{}{k})=\mathrm{\Delta }_0[cos(k_xa)cos(k_ya)]`$ that has the symmetry of a $`d`$-wave order parameter (a is the lattice constant). If we assume that the STM tunneling averages the DOS equally around the Fermi surface, the differential conductance $`dI/dV`$ is expected to be a linear function of energy near the EF–a behavior that is in strong contrast to that measured on the CuO2 plane. However, there have been theoretical efforts which suggest tunneling into a CuO2 plane along the c-axis may not be averaging equally over the entire Fermi surface.

These previous efforts have shown that tunneling in the perpendicular direction into a CuO2 plane is facilitated by the unoccupied 4s orbital of the Cu atoms, which extends further out-of-plane than the other orbitals in this plane. With the STM tip perpendicular to the CuO2 plane, we would thus expect the tip’s electronic states to overlap most strongly with the Cu 4s orbitals and probe $`\rho _s(E,\stackrel{}{k})`$ near EF through these orbitals (Fig. 4a). However, the CuO2 plane states near EF originate primarily from the Cu 3d band, which because of its d-orbital symmetry, couples to the Cu 4s orbitals with the matrix element $`|M(\stackrel{}{k})|^2[cos(k_xa)cos(k_ya)]^2`$ (Fig. 4a).  Including the consequences of the orbital symmetry of the CuO2 plane states on the tunneling process, we can relate the tunneling conductance to the DOS of a single CuO2 plane as $`dI/dV(V)|M(\stackrel{}{k})|^2\rho _s(E_F+V,\stackrel{}{k})e^{(E_F+VU)^2/2\sigma ^2}𝑑\stackrel{}{k}𝑑U`$,

where we have introduced Gaussian broadening (as characterized by $`\sigma `$) to account for energy broadening due to various sources. The effect of the matrix element is to filter the measurement of $`\rho _s(E,\stackrel{}{k})`$ preferentially for $`\stackrel{}{k}`$ directions in which it has a larger energy gap.

Figure 4b shows a comparison of the above model with the experimental data, using parameters to best match the tunneling conductance at low voltages, close to EF. The best fit to the experimental data is with $`\mathrm{\Delta }_0=60meV`$, which is larger than that expected for the slightly underdoped samples (Tc=84K) used in our experiments. It is, however, close to previously reported values of the energy gap for highly underdoped (0.1 holes per Cu-O plaquette) but still superconducting, samples from various measurements.  The underdoped behavior of a single CuO2 plane reported here may be the consequence of this layer being at the surface, where doping is provided to it only from below the surface–a situation different from that of a bulk CuO2 plane. The model described above has the obvious shortcoming in that it does not address features at higher energies in the data. However, it does provide a simple physical argument, based on the orbital symmetry of the electronic states of a CuO2 plane, for the apparent suppression of tunneling conductance at low biases.

In conclusion, we report the first STM imaging and spectroscopy measurements of a single CuO2 plane at the surface of a high-Tc cuprate. We find this layer to be atomically ordered and to exhibit an unusual suppression of tunneling DOS at low energies. We have considered a number of possible explanations and have proposed that the orbital symmetry of the electronic states in the CuO2 may be influencing the tunneling process and be responsible for the observed behavior at the surface.

We acknowledge discussions with P.W. Anderson, A.V. Balatsky, S.L. Cooper, J.C. Davis, M.E. Flatt$`\stackrel{´}{e}`$, E. Fradkin, L.H. Greene, M.V. Klein, R.B. Laughlin, A.J. Leggett, Z.-X. Shen, D. Pines, P.W. Phillips, and D.J. Van Harlingen. This work was supported under NSF-CAREER(DMR-98-75565), ONR(N000140110071), NSF-FRG(DMR-99-72087), and DOE through FSMRL(DEFG-02-96ER4539).


format: text

text:
# Measurement of electron wave functions and confining potentials via photoemission

## Abstract

Wave functions and electron potentials of laterally-confined surface states are determined experimentally by means of photoemission from stepped Au(111) surfaces. Using an iterative formalism borrowed from x-ray diffraction, we retrieve the real-space wave functions from the Fourier transform of their momentum representations, whose absolute values in turn are directly measured by angle-resolved photoemission. The effective confining potential is then obtained by introducing the wave functions into Schrödinger’s equation.

The electron wave functions and potentials in low-dimensional systems are of primary importance for tailoring electronic properties of nanostructures. The electron energy levels and the probability density are the physical observables. These can be obtained in real space from local conductance maps via scanning tunneling microscopy/spectroscopy (STM/STS) Avouris1 ; Buergi1 , or in reciprocal space by angle-resolved photoemission Kawakami1 ; Mugarza1 . The latter has been used to study thin films Kawakami1  or two-dimensional arrays of nano-objects Mugarza1 , whereas STM/STS resolves individual nanostructures on a surface.

The question arises whether it is possible to directly derive wave functions and/or electron potentials from experimental data. The standard procedure consists of assuming a model potential, whose parameters are obtained by fitting the experiment. For instance, in core-level photoelectron diffraction muffin-tin potentials are assumed for retrieving atomic positions that fit the experiment KFH78\_1 ; GVF01\_1 . The electron density of surface states confined between two steps has been measured by STM Avouris1 ; Buergi1  and has been modeled by semitransparent mirrors located at the step edges Buergi1 . It has been argued that the position of the mirror has to be moved slightly away from the steps to account for electron spill over Otero1 .

We introduce a direct method to determine the effective electron potential from angle-resolved photoemission data, without any a priori knowledge about its nature. All the necessary information is contained in the momentum distribution of the photoemission intensity. This allows us to unambiguously determine the potential and wave functions of one-dimensional quantum well states on stepped Au(111). The real-space wave functions are derived by Fourier transform of their momentum-space representations, the square of which is proportional to the photoelectron intensity under the conditions discussed below. However, the phase in momentum space is not measured and this leads to the well-known phase problem in optics, shared with many other techniques, such as x-ray and electron diffraction Saldin . Various iteration methods have been devised to retrieve real-space objects from the modulus of their momentum space representation, in particular the oversampling method Miao . These methods are valid for a confined wave function, because the phase is obtained by repeatedly diminishing the amplitude of the wave function outside the confinement region. The one-electron potential can be obtained after dividing the Schrödinger equation by the wave function.

As a test system we use vicinal noble metal surfaces characterized by equally-spaced, linear step arrays and one-dimensional surface states confined within individual terraces of nominal size $`L`$. That is the case of Au(111) vicinal surfaces with relatively wide terraces, such as Au(788) ($`L=38`$ Å) and Au(11,12,12) ($`L=56`$ Å) shown in Fig. 1 Mugarza1 ; Mugarza2 . On Au(11,12,12), the lowest three quantum-well levels ($`N=13`$) lie below the Fermi energy and they can be probed by photoemission, as shown in Figs. 2(a)-(b). The emission angle $`\theta `$ is given relative to the terrace normal in the plane perpendicular to the step array. The data have been taken at the PGM beam line in the Synchrotron Radiation Center (SRC) of the UW-Madison, using a hemispherical Scienta SES200 spectrometer. The energy and angular resolution were 20+7 meV (photons + electrons) and 0.3, respectively. The light is p-polarized and incident 60 off-normal. By tuning the photon energy to 60 eV we avoid overlap between surface umklapp replicas Mugarza2 . Fig. 2(c) shows the angular photoemission intensity scan for each quantum level as obtained from the areas of the peaks in Fig. 2(a). The scale for the component of the photoelectron wave vector parallel to the terrace $`q_x`$ has been determined from the emission angle and the kinetic energy $`E`$ as $`q_x=\mathrm{sin}\theta \sqrt{2mE}/\mathrm{}`$. The parallel momentum scale is shifted by the reciprocal step lattice vector $`g=2\pi /L=0.11`$ Å-1 to bring it into the first Brillouin zone Mugarza2  and corrected for the small parallel photon momentum of $`0.026`$ Å-1.

The photoemission intensity in reciprocal space $`𝐪`$ can be understood in the framework of standard photoemission theory, assuming a one-electron description. We start with the photoemission intensity, which is proportional to the matrix element

$`I_N(𝐪_{})|<𝐪|\mathrm{e}^{\mathrm{i}𝐤𝐫}𝐀|\mathrm{\Psi }_N>|^2.`$ (1)
$`\mathrm{\Psi }_N`$ is the initial surface state wave function, $`|𝐪>`$ is the final electron state of momentum $`𝐪=(𝐪_{},q_z)`$, $`𝐪_{}`$ is the momentum component parallel to the surface, $`𝐀`$ is the light polarization vector, $`𝐤`$ is the momentum of the photon, and $`N`$ refers to the quantum number of the wave function. To proceed further, it will be assumed that the surface state wave functions can be factorized into components that are either parallel or perpendicular to the surface, i.e. factorization :

$`\mathrm{\Psi }_N(𝐫)=\varphi _N(𝐑)\phi (z),`$
where $`𝐑`$ denotes the coordinates along the surface. If, in addition to the surface-state confinement in the perpendicular direction $`z`$, the electron is also bound to a 2D region of the surface, the parallel component of the wave function $`\varphi _N(𝐑)`$ must adopt a form which depends on the detailed shape of the confining region and on the boundary conditions at its border. Thus, $`I_N`$ can be written as

$`I_N(𝐪_{})`$ $`C|<𝐪_{}|\varphi _N>|^2,`$
where

$`C=|<q_z|\mathrm{e}^{\mathrm{i}k_zz}A_z{\displaystyle \frac{}{z}}|\phi >|^2`$
depends very weakly on $`𝐪_{}`$, so that it can be absorbed into a normalization factor that we set to 1. Therefore, there is a direct relation between photoemission intensity $`I_N(𝐪_{})`$ and the Fourier transform of the wave function $`\stackrel{~}{\varphi }_N(𝐪_{})`$:

$`\stackrel{~}{\varphi }_N(𝐪_{})=<𝐪_{}|\varphi _N>={\displaystyle 𝑑𝐑\mathrm{e}^{\mathrm{i}𝐪_{}𝐑}\varphi _N(𝐑)}.`$ (2)
Thereby we neglect multiple scattering in the final state, which corresponds to Fourier components $`𝐪_{}`$+$`𝐠_{}`$ with a reciprocal lattice vector $`𝐠_{}`$. These affect the intensity distribution between different Brillouin zones but not the momentum distribution within one Brillouin zone (compare this to the spot profiles in LEED which are not affected by multiple scattering Henzler ). This permits obtaining the wave function in real space $`\varphi _N(𝐑)`$ from the intensity in reciprocal space as

$`\varphi _N(𝐑)={\displaystyle \frac{d𝐪_{}}{(2\pi )^2}\mathrm{e}^{\mathrm{i}𝐪_{}𝐑}\sqrt{I_N(𝐪_{})}\mathrm{e}^{\mathrm{i}\delta _N(𝐪_{})}}.`$ (3)
However, this equation still contains an unknown phase $`\delta _N(𝐪_{})`$. We have explored two methods to obtain this phase: (1) an iterative procedure using oversampling Saldin ; Miao  and (2) an expansion of the wave function into a Fourier series combined with a least squares fit. In both cases, the strategy is minimizing $`|\varphi |`$ outside the confinement region.

For the iterative method, we begin with a constant phase $`\delta _N(𝐪_{})=0`$ to start the iteration with Eq. (3). The resulting $`\varphi _N(𝐑)`$ is corrected outside the confining region, and transformed back into $`𝐪_{}`$ space using Eq. (2). The phase of $`\stackrel{~}{\varphi }_N(𝐪_{})`$ is extracted and inserted again into Eq. (3) to start a new iteration. The confinement length is obtained directly from the photoemission data (see below), and the noted correction consists of subtracting from the newly calculated real-space wave function the wave function at the previous step multiplied by a factor 0.1. This ensures convergence to a wave function that vanishes outside the confinement region. As a strong test of convergence of this iterative technique, we have used the phase of the $`N=2`$ state to start the iteration of the $`N=1`$ state. The first iteration step leads to a wave function that resembles that of the $`N=2`$ state, a clear indication of the importance of the phase Saldin , and convergence to the true $`N=1`$ state is achieved after several hundred iterations.

This iterative method leads to results that are in agreement with those derived from an expansion of the wave function into sine functions that vanish at the edges of the confinement region. Moreover, the length of the region where the wave function takes non-negligible values is not sensitive to the input value for the length of the confinement region.

The size of the confining region is obtained from the Fourier transform of the experimental intensity distribution, which can be written as the self-convolution of the real-space wave function,

$`{\displaystyle \frac{d𝐪_{}}{(2\pi )^2}\mathrm{e}^{\mathrm{i}𝐪_{}𝐑}I_N(𝐪_{})}={\displaystyle 𝑑𝐑^{}\varphi _N(𝐑^{})\varphi _N^{}(𝐑^{}𝐑)}.`$ (4)
This convolution takes non-zero values in an area twice the size of the confining region. For the data of Fig. 2 in particular, one obtains curves that are basically confined within a region of $`60`$ Åin diameter (not shown), in agreement with the nominal terrace width in Au(11 12 12), $`L=56`$ Å.

The natural normalization for both the two-dimensional wave function $`\varphi _N`$ and the measured intensity is provided by Eq. (4) if one sets $`𝐑=0`$:

$`{\displaystyle \frac{d𝐪_{}}{(2\pi )^2}I_N(𝐪_{})}=<\varphi _N|\varphi _N>=1.`$
Applying this procedure to the angular scans of Fig. 2(c), we have obtained the surface state wave functions for quantum-well levels in Au(11 12 12) that are shown in Fig. 1. All three wave functions are confined to a region whose width matches the terrace width $`L`$. Thus the wave functions exhibit a clear terrace confinement that has not been assumed by our reconstruction procedure, but rather it has emerged from the information contained in the photoemission data.

The actual effective potential of the terrace can be retrieved from the Schrödinger equation as

$`V(𝐑)E_N={\displaystyle \frac{\mathrm{}^2}{2m^{}\varphi _N(𝐑)}}^2\varphi _N(𝐑),`$ (5)
where $`m^{}=0.26m`$ is the effective mass of the electron in the initial state of energy $`E_N`$, taken from the dispersion of the surface state on flat Au(111). The potential becomes independent of the quantum number $`N`$ if many-electron effects are absent, as we assumed initially.

The experimental wave functions of Fig. 1 have been introduced into Eq. (5) and the resulting electron potential has been represented in Fig. 3. In order to compare the potential derived from the quantum states with various $`N`$ we add on the left side of Eq. (5) the respective experimental values $`E_N^{\mathrm{exp}}`$ shown in Fig. 1 (horizontal bars). Excellent mutual agreement is obtained for the shape of the potential derived independently from each of the three wave functions $`N=13`$, proving the validity of our method. The potential exhibits a smooth central region and sharp boundaries that force electron confinement. It must be stressed that the validity of Eq. (5) is limited to regions where the wave function is not too small, and therefore, the asymptotic limit of the step barrier potential cannot be determined.

There is a slight uphill/downhill asymmetry of $`\varphi `$ and $`V`$, which is related to an asymmetry of $`|\stackrel{~}{\varphi }|`$ for $`\pm q`$ . This effect shows up for the highest level ($`N=3`$), which is the least confined. Such an asymmetry explains the asymmetric reflectivity observed in STM Buergi1  where the downhill step reflects more strongly than the uphill. The asymmetry increases for higher-lying empty states which start spilling over the weaker potential barrier on the uphill side.

The same iterative procedure has been applied to the one-dimensional quantum well of Au(788), measured in a different system Mugarza1 . In this case the terraces are 38 Åwide and only the first two quantum levels are occupied. We also obtain a good consistency in the electron potential between the two levels. The resulting average is displayed in Fig. 1, and compared with the average potential of Au(11,12,12). In contrast to the latter case, the width of the potential well in Au(788) is narrower than the terrace size, as also observed in thin Pb films Otero1 .

In summary, we have introduced a simple procedure for directly obtaining wave functions and effective potentials of confined electronic states from the momentum distribution of the photoemission intensity. This procedure, which relies on photoemission theory within the one-electron approach, has been successfully tested in the case of laterally confined surface states on Au(111) vicinal surfaces. The method is of general applicability for discrete states because they are spatially confined and allow an iterative determination of the phase by forcing $`\varphi `$=0 outside the confinement region. It can be applied to any type of surface nanostructure confined in at least one dimension, such as arrays of quantum wires and quantum dots. While STM provides direct images of the low spatial frequencies, photoemission determines the high spatial frequencies. It would be interesting to pursue the reverse transformation from the real space $`|\varphi |`$ in STM to $`\stackrel{~}{\varphi }`$ in $`q`$ space using the technique described here. A connection between $`|\varphi |`$ and the $`E(q)`$ band dispersion has been made already Buergi1 ; Horn

A.Mu. and J.E.O. are supported by the Universidad del País Vasco (9/UPV 00057.240-13668/2001) and the Max Planck Research Award Program. J.G.deA. acknowledges support by the Spanish Ministerio de Ciencia y Tecnología. F.J.H. acknowledges support by the NSF under Award Nos. DMR-9815416 and DMR-0084402.


text:
# The local structure of topological charge fluctuations in QCD

## Abstract

We introduce the Dirac eigenmode filtering of topological charge density associated with Ginsparg-Wilson fermions as a tool to investigate the local structure of topological charge fluctuations in QCD. The resulting framework is used to demonstrate that the bulk of topological charge in QCD does not appear in the form of unit quantized lumps. This means that the mixing of “would-be” zeromodes associated with such lumps is probably not the prevalent microscopic mechanism for spontaneous chiral symmetry breaking in QCD. To characterize the coherent local behavior in topological charge density at low energy, we compute the charges contained in maximal coherent spheres enclosing non-overlapping peaks. We find a continuous distribution essentially ending at $`0.5`$. Finally, we study, for the first time, the overlap-operator topological-charge-density correlators and find consistency with non-positivity at nonzero physical distance. This represents a non-trivial check on the locality (in gauge paths) of the overlap Dirac operator for realistic gauge backgrounds.

The aim of this work is to set up a physically meaningful and unbiased framework to probe the low-energy structure of topological charge fluctuations in QCD, and use it to study some related physics . As a starting point, we will consider the following question. Can one confirm or disprove a qualitative picture of vacuum fluctuations where most topological charge appears in the form of identifiable coherent unit quantized lumps? This question is interesting because if the answer is yes, then one could possibly argue that S$`\chi `$SB, the $`\theta `$-dependence, and the $`\eta ^{}`$ mass have common qualitative origin modeled in the instanton picture. However, if the answer is no, then this frequently-used intuitive reasoning has no verifiable basis. The main reason why a lumpy vacuum (vacuum dominated by unit quantized lumps) represents a dividing line is that such vacuum has a very specific effect on fermions. In particular, the chiral “would-be” zeromodes $`\chi ^i`$ implied by the index theorem and localized on the lumps $`_i`$ would mix and form a special subspace of states $`\psi ^i_ja^{ij}\chi ^j,i,j=\mathrm{\hspace{0.17em}1},\mathrm{}N_L`$ which will be referred to as the topological subspace.

The usual reasoning associated with topological mixing contains a subtlety rooted in the fact that the Euclidean $`q(x)`$ is odd with respect to reflections and consequently $`<q(x)q(0)>\mathrm{\hspace{0.17em}0}`$ for $`|x|>0`$,  (see also ). This means that any four-dimensional coherent structure in the long-wavelength components of $`q(x)`$ will be obscured by anticorrelated background fluctuations. This fact appears to be at the heart of all the problems associated with studying topology locally and, since we discuss this in the continuum, it has nothing to do with lattice as is sometimes incorrectly argued.

Since the regions $`_i`$ cannot be identified in $`q(x)`$, the topological mixing scenario is ill-defined (and thus can not be verified) until one defines the coherent lumps for realistic backgrounds in the first place. The issue here is that, naively, the topological mixing picture relates the local behavior of infrared modes, where the ultraviolet fluctuations of the underlying gauge field are filtered out , to the behavior of $`q(x)`$ containing all ultraviolet fluctuations. This suggests that a logically consistent framework can be built around the possiblity of defining an effective topological charge density, based on the infrared eigenmode expansion. This can be conveniently achieved starting from lattice topological charge density associated with GW fermions

$$q_x=\frac{1}{2}\text{tr}\gamma _5D_{x,x}=\text{tr}\gamma _5(1\frac{1}{2}D_{x,x})$$
(1)
and defining the effective densities through

$$q_x^{(k)}\underset{i=1}{\overset{N_0}{}}c_x^{0,i}\underset{j=1}{\overset{k}{}}(2\text{Re}\lambda _j)c_x^{\lambda _j}$$
(2)
where $`D\psi ^\lambda =\lambda \psi ^\lambda `$, and $`c_x^\lambda =\psi _x^{\lambda +}\gamma _5\psi _x^\lambda `$ is the local chirality. $`N_0`$ zeromodes and $`k`$ conjugate complex eigenpairs (ordered by increasing real part of the eigenvalue) contribute to $`q^{(k)}`$.

Like the full density $`q_x`$, the filtered densities $`q_x^{(k)}`$ satisfy the exact index theorem , and yield identical global fluctuations since $`_xq_x^{(k)}=Q`$, $`k`$. However, the singular short-distance behavior of $`q_x`$ is smoothed out , and the coherent structures in $`q_x^{(k)}`$ can potentially survive the continuum limit. The topological mixing picture can now be formulated consistently and checked. In particular, the picture asserts that for a typical gauge background, the lowest effective densities $`q_x^{(k)}`$ will be dominated by $`N_L`$ sign-coherent regions $`_i`$, and the topological charge contained in these regions will saturate at $`k_L`$ satisfying $`N_L=N_0+2k_L`$, when $`_{x_i}q_x^{(k_L)}\pm 1`$. Note that the existence of unit lumps in effective densities would then be consistent with the basic assumption of topological mixing, according to which the low-lying modes are built of “would-be” zeromodes $`\chi ^i`$ localized on $`_i`$ so that $`_{x_i}|\chi _x^i|^21`$. Indeed, since $`\lambda _i0`$ for modes in the topological subspace, it follows that

$$\underset{x_i}{}\underset{j=1}{\overset{N_L}{}}\psi _x^{j+}\gamma _5\psi _x^j\underset{x_i}{}q_x^{(k_L)}\pm 1$$
(3)
We now test this proposition on QCD backgrounds. If the low-energy structure in a configuration with topological charge $`Q`$ is dominated by $`|Q|+k_L`$ coherent unit lumps (antilumps) and $`k_L`$ antilumps (lumps), then we should observe that $`_x|q_x^{(k_L)}|Q^{k_L,abs}N_L=|Q|+2k_L`$. In fact, the quantity $`|Q|+2kQ^{k,abs}`$ would be close to zero for all $`kk_L`$ (each mode in the topological subspace contributes $`1`$ to $`Q^{k,abs}`$), and it would start to increase rapidly for $`k>k_L`$ (modes outside the topological subspace contribute $`1`$). If the topological mixing scenario was relevant, monitoring the $`k`$-dependence of $`|Q|+2kQ^{k,abs}`$ could thus serve as a procedure for determining the dimension of the topological subspace for a given configuration. We have computed this $`k`$-dependence using overlap Dirac operator for a few dozen of configurations at various lattice spacings, but have never observed the expected break at definite $`k_L`$. Instead, the typical robust behavior is shown below for an ensemble of six $`12^4`$ Wilson gauge configurations at $`\beta =5.91`$.

Given our volume, $`|Q|+2kQ^{k,abs}`$ should be close to zero for $`k1,2`$, as the consistency with the value of topological susceptibility would require, but we observe large values instead. This apparent inconsistency leads us to conclude that the effective densities for QCD configurations are not dominated by coherent unit lumps.

An important feature of effective densities is that they can be used to study the low-energy structure of topological charge fluctuations without prejudice to a specific picture. As a first step towards this we calculate the typical amounts of topological charge contained in the maximal coherent spheres enclosing the peaks observed in effective densities. We fix the resolution for identification of such peaks by requiring that most of the resulting structures do not overlap . Shown below is the distribution of charges obtained from six $`12^4`$ configurations at $`\beta =5.91`$ using $`q^{(9)}`$.

The charges contained in the peaks are reasonably stable with respect to inclusion of additional eigenmodes. Interestingly, the distribution essentially ends at $`0.5`$ and most of the charge is carried by structures containing smaller amounts. This is consistent with the scenario emerging from the assumption of dominance by center vortices .

Finally, we come back to the question of whether the overlap $`q_x`$ complies with non-positivity of the correlator at non-zero distances. For any valid reflection-positive lattice action $`S(U)`$ and well-defined local $`q_x`$, the consistency of the continuum limit requires that $`lim_{a0,a|x|=r_p}<q_xq_0>\mathrm{\hspace{0.17em}0}`$, where $`r_p`$ is a non-zero physical distance. Since the Wilson gauge action is reflection positive, in our case this property relies on the locality of overlap $`q_x`$, and hence on the locality of $`D_{x,x}`$. This is not guaranteed since, contrary to the locality in fermionic variables, the locality in gauge paths has not been studied for realistic backgrounds. In the plot above we show preliminary results for $`<q_xq_0>`$ (normalized at the origin) at two lattice spacings. Even though the number of configurations is small, the behavior is very robust and reveals that the effect of chiral smoothing  is about two lattice spacings. This appears to be scale-independent. Indeed, the size of the positive core drops accordingly (see figure below), and is consistent with being a contact term in the continuum limit. This justifies the use of overlap $`q_x`$ and provides a nontrivial check on the locality properties of overlap Dirac operator. A detailed account of these observations is in preparation.


text:
# On the Ejection Mechanism of Bullets in SS 433

## 1 Introduction

SS 433 remains one of the most enigmatic objects in the sky. Even twenty-five years after its first appearance in the catalogue of Stephanson & Sanduleak (1977), it is not clear whether the compact object is a black hole or a neutron star. However, there is ample evidence that the companion is an OB type star with an orbital period of 13.1 days, which is losing mass at the rate of about $`10^4M_{\mathrm{}}\mathrm{yr}^1`$ (van den Heuvel 1981), corresponding to extremely super-Eddington accretion regardless of the mass of the compact object.

One of the most curious properties of the jets of SS 433, which first made their presence distinctly felt through the emission of variable H$`\alpha `$ lines, is that they are apparently ejected as bullets (e.g. Borisov & Fabrika 1987; Vermeulen et al. 1993; Paragi et al. 1999; 2002; Gies et al. 2002), with a surprisingly nearly constant radial velocity of about $`0.26c`$. The absence of a significant intrinsic rotational velocity (i.e., $`v_\varphi `$) component is clear from the fact that the kinematic model (e.g. Abell and Margon 1979), which assumes only radial injection, quite accurately explains the time variation of the red- and blue-shifts of the H$`\alpha `$ emission from the jets with a period of 162 days, which is attributed to the precession of the accretion disk about the compact object. The radial velocity is less than the maximum allowed sound speed of $`c/\sqrt{3}`$ and thus hydrodynamic acceleration could, in principle, explain it. Therefore one may not require a magnetic or electrodynamic acceleration processes (e.g. Belcher & MacGregor 1976; Lovelace 1976). However, the rather good collimation (Margon 1984; Paragi et al. 1999) supports the hypothesis that a substantial degree of confinement produced by toroidal flux tubes may be present. Gies et al. (2002) showed that the ratios of the H$`\alpha `$ emission equivalent widths from the approaching and receding jets as a function of precessional phase only could be nicely fit if these emission components are bullet-like. Indeed, the recent Chandra X-ray Observatory discovery of X-rays at a distance of about $`10^{17}`$cm from the center may result from the collision of such bullets (Migliari, Fender, & Mendez 2002).

SS 433 poses another interesting problem: it was pointed out by Chakrabarti (1999) and Das & Chakrabarti (1999) that significant outflows are produced only when the accretion rate is such that the X-ray source is in a low/hard state, and all the observational indications in other micro-quasars also suggest that the jets are indeed produced in low/hard states (Corbel et al. 2001; Klein-Wolt et al. 2002). However, it is difficult to imagine how SS 433 manages to remain in the low/hard state with $`10^4M_{\mathrm{}}\mathrm{yr}^1`$ of wind matter ejected from its companion. The answer to this quandary probably lies in the recent results of Paragi et al. (1999) and Blundell et al. (2000), whose high resolution radio maps show that there is a large region of roughly $`50\mathrm{A}\mathrm{U}`$ in radius which is filled with enough gas and dust to obscure the accretion disk and the base of the jets. They also found an equatorial outflow. Gies et al. (2002) present additional evidence from observations of the “stationary” H$`\alpha `$ and He I lines for an extended “disk wind”. So it is distinctly possible that most of the matter from the donor is rejected either by centrifugal force (Chakrabarti 2002) or by radiation force far outside the central accretion disk, and thus the compact object receives only a few times the Eddington rate ($`\dot{M}_{\mathrm{Ed}}`$) of its companion’s wind matter to accrete. This consideration finds further support from the fact that the kinematic luminosity of the jet itself is around $`10^{39}`$ erg s-1 (Margon 1984), which corresponds to about one Eddington rate for a $`10M_{}`$ compact object.

In numerical simulations of supercritical winds by Eggum, Coroniti & Katz (1985) designed to model SS 433, it was shown that only a fraction of a percent of the infalling matter is ejected from a radiation pressure supported Keplerian disk, which indicates that the accretion rate must be at least $`100\dot{M}_{\mathrm{Ed}}`$ if the accretion takes place through a Keplerian disk. On the other hand, numerical simulations of a sub-Keplerian disk by Molteni, Lanzafame & Chakrabarti (1994) suggest that about $`1520`$ percent of matter is ejected as an outflow, indicating that the accretion rate onto the compact object in SS 433 need be at most a few $`\dot{M}_{\mathrm{Ed}}`$. Similar simulations with different parameters yield situations where no steady shocks can form, even though two saddle-type sonic points are present (Ryu, Chakrabarti, & Molteni 1997, hereafter RCM); under these conditions large scale shock oscillations produce intermittent outflows instead of continuous outflows. Since the compact object is a wind accretor, a low-angular momentum, sub-Keplerian flow is the most likely description of the accretion flow. Indeed, the presence of sub-Keplerian flows in several other high mass X-ray binaries has now been verified (Smith, Heindl & Swank 2002).

In this Letter, we present a few scenarios leading to ejection of matter as bullets in SS 433. We discuss four possible ways to create blobs of matter emerging from the disk and conclude that periodic ejection of the blobs by the large scale oscillation of an accretion shock (something like a piston) may be the fundamental production mechanism of the “normal” bullets. The irregularly observed rapid flaring (Vermeulen et al. 1993) could be understood in terms of the catastrophic collapse of toroidal magnetic flux tubes, very similar to what has been argued to be occurring in GRS 1915+105 (Vadawale et al. 2001; Nandi et al. 2001). In the next section we discuss these processes and their suitability or unsuitability for SS 433. In §3, we present concluding remarks.

## 2 Mechanisms to produce bullet-like ejecta from accretion flows

In both the works of Eggum et al. (1985) and Molteni et al. (1994) continuous ejection was reported when a radiation pressure dominated Keplerian disk, or a sub-Keplerian disk capable of producing a steady shock, were considered. However, in SS 433 the basic ejection is bullet-like and since the size of the X-ray emitting region is smaller than $`l_x10^{12}`$cm within which the material in the jets is already accelerated to $`v_{jet}0.26c`$ (Watson et al. 1986; Stewart et al. 1987), the bullets are not expected to be delayed by more than $`l_x/v_{jet}100`$s. Indeed, recent Rossi X-Ray Timing Experiment (RXTE) observations of hard X-rays from SS 433 indicated variability on time scales of 50–1000 s (Safi-Harb & Kotani 2002), roughly corroborating this picture. In fact, a simultaneous measurement of a flare at $`2`$ GHz in the radio (Kotani & Trushkin 2001) and in hard X-rays (Safi-Harb & Kotani 2002) indicated a strong anti-correlation of radio and X-ray fluxes, similar to what is observed in GRS 1915+105 (Mirabel & Rodriguez 1994). Moreover, the X-ray luminosity is very low ($`10^{36}`$erg s-1) and is believed to come from the base of the jets (Watson et al. 1986). It is believed to have a thermal origin and EXOSAT (Watson et al. 1986) and GINGA (Yuan et al. 1995) observations were adequately fitted with a thermal bremsstrahlung model with $`kT\stackrel{>}{}30`$ keV. The overall spectral shape suggests that the source has always been in a standard low/hard state and so far no quasi-thermal emission expected from a ‘Keplerian disk’ has been detected. From the interaction of the jet with the supernovae remnant W50, the lower limit of kinematic luminosity is found to be at least $`10^{39}`$erg s-1 (Biretta et al. 1983; Davidson & McCray 1980). This means that the mean mass outflow rate is around $`10^{18}`$g s-1, and if most of it is in the form of bullets ejected at 50–1000 s intervals, the mass accumulated in each bullet should be in the range of $`10^{19}`$$`10^{21}`$g.

The above data implies that the essential features that one must explain when attempting to produce bullets out of the accretion disks are: (a) the disk should be a sub-Keplerian flow; (b) the object (black hole or a neutron star) and its surroundings should be in a low/hard state; (c) bullets should be ejected in 50–1000 s time-scales under normal circumstances; (d) the mass of each bullet should be around $`10^{1921}`$g; and finally, (e) there should be occasional flaring with an anti-correlation of radio and X-ray emission. We now discuss several scenarios and present what we believe to be the most probable picture of what is going on in SS 433. The four processes are schematically shown in Fig. 1(a-d).

### 2.1 Cooling of the jet-base by Comptonization and separation of blobs

It was shown by several numerical simulations that significant outflows are produced from regions very close to the inner edge of the accretion flow, possibly from the centrifugal pressure dominated region (Molteni et al. 1994; Molteni et al. 2001). These jets are launched subsonically, but quickly pass through the inner sonic point to become supersonic. In the subsonic region while the matter moves slowly, the density is high and the optical depth could be large enough ($`\tau >1`$) to undergo Compton cooling (Fig. 1a) provided there is a Keplerian disk underneath to supply soft photons. A part of the outflow, which was subsonic previously becomes supersonic because of this rapid cooling and separates from the base of the jet. This separation of blobs is expected to occur at the sonic surface $`r_c`$ which is $`23r_s`$, where, $`r_s`$ is the size of the centrifugal barrier; see Chakrabarti (1999).

This possibility, though attractive, and in fact likely to be a major mechanism for rapid state change in objects like GRS 1915+105 (Chakrabarti & Manickam 2000), is untenable in SS 433 because the latter is a wind accretor: thus no significant Keplerian disk is expected in this system to supply the soft photons, and indeed none has been detected so far (Watson et al. 1986; Yuan et al. 1995).

### 2.2 Resonance oscillation of accretion shocks in the presence of bremsstrahlung cooling

Numerical simulations of accretion flow show that in cases where the cooling time-scale nearly matches the infall time-scale, a shock forms, but it then starts oscillating and ejects matter quasi-periodically (Langer, Chanmugam & Shaviv 1983; Molteni, Sponholz & Chakrabarti 1996, hereafter MSC; see Fig. 1b). In order to have an oscillation period of around $`50`$s, the shock must be located at the large distance of $`r_{s,\mathrm{MSC}}6400r_g`$ for a black hole of mass $`M=10M_{\mathrm{}}`$, where $`r_g=2GM/c^2`$. The mass of the post-shock region is computed by equating the bremsstrahlung (which we assume to be the major cooling mechanism) cooling time and the infall time in the post-shock region (MSC):

$$T_{\mathrm{MSC}}\frac{}{\dot{}}\frac{r_{s,\mathrm{MSC}}}{v_f}\left(\frac{Rr_{s,\mathrm{MSC}}}{r_g}\right)^{3/2}\frac{r_g}{c},$$
(1)
where, $``$ is the specific thermal energy, $`v_f`$ is the infall velocity and $`R=(\gamma +1)/(\gamma 1)47`$ (these limits are for a strong shock with $`\gamma =5/3`$ and $`\gamma =4/3`$, respectively) is the compression ratio at the shock. Assuming the gas density ($`n`$) and temperature ($`T`$) scale as $`nr^{3/2}`$ and $`Tr^1`$ respectively, the mass of the sub-Keplerian region of $`r<r_s`$ turns out to be $`7\times 10^{19}`$g (with $`M=10M_{\mathrm{}}`$, $`\gamma =5/3`$). This is indeed of the same order as the mass of the bullets observed in SS 433. However, one has to have both the angular momentum and energy of the injected material comparable to the marginally bound values determined by the central object in order to achieve such an oscillation. On the other hand, if the mass expulsion from the system takes place at the similar radius of $`r_{ex}10^4r_g=2\times 10^9M/M_{}`$cm due to the centrifugal force, the specific angular momentum of the flow is approximately $`[r_{ex}/(2r_g)]^{1/2}r_gc70r_gc`$, which is very large compared to the marginally bound value of $`2r_gc`$. So it is unlikely that this mechanism works in SS 433.

### 2.3 Non-steady and non-linear shock oscillation

A standing shock can form in a sub-Keplerian flow only if there are two saddle type sonic points and the Rankine–Hugoniot relation is satisfied at least at one point in between these two sonic points. However, Chakrabarti (1990) showed that there is a large region of the parameter space where there are two saddle type sonic points but the shock-conditions are not satisfied. Even an initially supersonic accretion (such as the wind from the companion) can fall into this category.

What will happen to such a realistic flow, especially when the specific entropy at the inner sonic point is greater than that at the outer sonic point? RCM discovered that a flow injected with these parameters exhibits yet another type of shock-oscillation (Fig. 1c). Here the shock searches for a stable location and oscillates without finding it. In the first half of the cycle, the shock recedes far away, the post-shock region fills up, but the accretion is essentially completely blocked. In the second half of the cycle, the shock pushes the matter into the black hole, thereby evacuating the post-shock region. In a realistic simulation, RCM find that while the ratio of actually accreted matter to the amount available from the companion, $`R_{\mathrm{ai}}\dot{M}_{\mathrm{acc}}/\dot{M}_{\mathrm{inj}}`$ would be around $`0.2`$ during the first half-cycle, $`R_{\mathrm{ai}}1.3`$ in the second half-cycle. The outflow was also found to be very large. The time scale of oscillation was found to be $`T_{\mathrm{RCM}}(40006000)r_g/c`$ for a $`r_s20r_g`$ whose infall time is only about $`T_{\mathrm{MSC}}(Rr_s/r_g)^{3/2}(r_g/c)(350400)r_g/c`$. Thus, this type of oscillation takes about a factor of $`R_T=T_{\mathrm{RCM}}/T_{\mathrm{MSC}}15`$ times longer than the resonance oscillation discussed in §2.2. For a $`50\mathrm{s}`$ oscillation, the location of the shock should be obtained from $`(r_{s,\mathrm{RCM}}/r_g)^{3/2}(1/R)(50\mathrm{s}/R_T)(c/r_g)10^4`$, which gives, $`r_{s,\mathrm{RCM}}450r_g`$ for a $`10M_{}`$ black hole, a more physically reasonable value. Even though the size of the oscillating region goes down by a factor of $`10`$ or so, compared to that involved in the resonance oscillation, the ejected mass need not go down (even for the same accretion rate as in the earlier case). This is because nearly all of the accretion flow is accumulated in half the cycle ($`25`$s in this case) before being ejected (see Fig. 2 of RCM).

Another advantage of this type of non-steady shock oscillation is that it is driven by centrifugal force and not by thermal cooling. Hence the result is generally independent of the accretion rate. Thus as long as the viscosity remains low, equivalent to having the Shakura-Sunyaev (1973) parameter, $`\alpha \alpha _c0.015`$ (Chakrabarti 1990), and $`\dot{M}_{\mathrm{inj}}`$ remains fairly constant, this oscillation, once established, could be sustained indefinitely.

### 2.4 Magnetic rubber-band effect

In the event of increase in magnetic activity of the disk, as could happen for instance when the accretion disk bends towards the binary companion during its precessional motion, it is not unlikely that a strong magnetic field will be first intercepted, and then advected, toward the inner edge of the disk. In this case the field will preferentially become toroidal due to shear in the rotating flow. Then, as has already been pointed out (Chakrabarti & D’Silva 1994; Nandi et al. 2001) the acceleration due to magnetic tension

$$a_T=\frac{B_\varphi ^2}{4\pi r(\rho _e+\rho _i)}\frac{B_\varphi ^2}{4\pi r\rho _e},$$
(2)
would be the dominant force in the post-shock region of the sub-Keplerian flow (Fig. 1d). Here, $`r`$ is the major radius of the toroidal flux tube and $`\rho _i`$ and $`\rho _e`$ are the densities of the medium internal and external to the flux tube respectively. The last step in Eqn. (2) is written because $`\rho _i<<\rho _e`$ for a strong flux tube. Since $`B_\varphi 1/r`$ and $`\rho _er^{3/2}`$, we get

$$a_Tr^{3/2},$$
(3)
thus increasing rapidly as the tube comes closer to the black hole, and even surpassing the magnetic buoyancy,

$$a_{MB}=\frac{1X}{1+X}\left[\frac{\lambda _{\mathrm{Kep}}^2\lambda ^2}{r^3}\right]\frac{\lambda _{\mathrm{Kep}}^2\lambda ^2}{r^3},$$
(4)
where $`X=\rho _i/\rho _e0`$, and $`\lambda _{\mathrm{Kep}}`$ and $`\lambda `$ are respectively the specific angular momentum of a Keplerian disk and the disk under consideration. the accelerations in Eqns. (3) and (4) do cross over, since when very close to a black hole, $`\lambda \lambda _{\mathrm{Kep}}`$ for a sub-Keplerian flow.

The effect of magnetic tension is dramatic, and the inner part of the disk is evacuated in the Alfvén time scale: $`r/v_A(r/a_T)^{1/2}\stackrel{<}{}0.1`$s, for a $`10M_{\mathrm{}}`$ black hole with a realistic Alfvén speed, $`v_A0.1c`$ (Nandi et al. 2001). The enhanced plasma ejection along the axis presumably causes sporadic magnetic flare events which would be observable as radio outbursts, at the same time reducing the X-ray emission from the disk which forms the base of the jet. Recently, such effects may have been seen (Safi-Harb & Kotani 2002) where simultaneous observations of 2 GHz radio and 2–20 keV X-ray fluxes from SS 433 have been made, and a clear dip in X-ray flux is seen at the same time a strong radio flare is observed. It is worth noting that similar anti-correlated variations are common during flares in GRS $`1915+105`$ (Feroci et al. 1999; Naik et al. 2001) and we suggest that the flares in SS 433 originate in the same way.

## 3 Concluding remarks

In this Letter, we have studied various competing processes for the creation of bullets which move ballistically in the jet of SS 433. We showed that blobs may be separated by: (1) Comptonization; (2) shock oscillations due to resonance; (3) oscillations due to inherent unsteady accretion solutions; (4) intense magnetic tension of the toroidal flux tubes. We reject the first possibility because it requires a large Keplerian disk, which is unlikely. We are unable to distinguish at this stage which type of shock oscillation is more capable of producing bullet formation in SS 433, but we prefer the third possibility due to its impulsive and generic nature and smaller involved region. We believe that the fourth possibility of the inner disk evacuation should produce flaring events, but will occur rather rarely, perhaps only once in a single precession period, when the magnetic field of the companion is preferentially tilted towards the accretion disk during precessional motion. This fourth mechanism gives rise to an anti-correlation between radio and X-rays, perhaps already observed in SS 433 (Safi-Harb & Kotani 2002).

SKC, AN, and SD acknowledge a grant from the Department of Science and Technology, India, and financial support from CEA/Saclay where part of this work was performed. PJW is grateful for hospitality at the Department of Astrophysical Sciences, Princeton University, and for support from the Research Program Enhancement program at Georgia State University.


text:
# Quantum Monte Carlo method using phase-free random walks with Slater determinants

## Abstract

We develop a quantum Monte Carlo method for many fermions that allows the use of any one-particle basis. It projects out the ground state by random walks in the space of Slater determinants. An approximate approach is formulated to control the phase problem with a trial wave function $`|\mathrm{\Psi }_T`$. Using plane-wave basis and non-local pseudopotentials, we apply the method to Si atom, dimer, and 2, 16, 54 atom (216 electrons) bulk supercells. Single Slater determinant wave functions from density functional theory calculations were used as $`|\mathrm{\Psi }_T`$ with no additional optimization. The calculated binding energy of Si2 and cohesive energy of bulk Si are in excellent agreement with experiments and are comparable to the best existing theoretical results.

Quantum Monte Carlo (QMC) methods based on auxiliary fields (AF) are used in areas spanning condensed matter physics, nuclear physics, and quantum chemistry. These methods BSS ; Koonin  allow essentially exact calculations of ground-state and finite-temperature equilibrium properties of interacting many fermion systems. The required CPU time scales in principle as a power law with system size, and the methods have been applied to study a variety of problems including the Hubbard model, nuclear shell models, and molecular electronic structure. The central idea of these methods is to write the imaginary-time propagator of a many-body system with two-body interactions in terms of propagators for independent particles interacting with external auxiliary fields. The independent particle problems are solved for configurations of the AF and averaging over different AF configurations is then performed by Monte Carlo (MC) techniques.

QMC methods with auxiliary fields have several appealing features. For example, they allow one to choose any one-particle basis suitable for the problem, and to fully take advantage of well-established techniques to treat independent particles. Given the remarkable development and success of the latter DFT , it is clearly very desirable to have a QMC method that can use exactly the same machinery and systematically include correlation effects by simply building stochastic ensembles of the independent particle solutions. Vigorous attempts have been made from several fields to explore this possibility Silvestrelli93 ; shiftcont\_rom97 ; Wilson95 ; AFQMC\_appl .

A significant hurdle exists, however: except for special cases (e.g., Hubbard), the two-body interactions will require auxiliary fields that are complex. As a result, the single-particle orbitals become complex, and the MC averaging over AF configurations becomes an integration over complex variables in many dimensions. A phase problem thus occurs which ultimately defeats the algebraic scaling of MC and makes the method scale exponentially. This is analogous to but more severe than the fermion sign problem with real AF loh90 ; Zhang  or in real-space methods schmidt84 . No satisfactory, general approach exists to control the phase problem. As a result, only small systems or special forms of interactions can be treated.

In this paper we address this problem. We develop a method for many-fermions that allows the use of any one-particle basis. It projects out the ground state by random walks in the space of Slater determinants. The phase problem is eliminated with an approximation that relies on a trial wave function $`|\mathrm{\Psi }_T`$. We demonstrate the method by applying it to electronic systems using a plane-wave basis and non-local pseudopotentials, which can be implemented straightforwardly in this method. We calculate the binding energy of Si2 and the cohesive energy of bulk Si using fcc supercells consisting of up to 54 atoms (216 electrons). These calculations represent the first application of AF-based QMC to solids. The results are in excellent agreement with experiments and are comparable to the best existing theoretical results. Particularly worth noting is that our results were obtained with a trial wave function which is a single Slater determinant formed by orbitals from density functional theory (DFT) calculations (with the local density approximation (LDA)), with no additional parameters or optimization.

The Hamiltonian for any many-fermion system with two-body interactions can be written in any one-particle basis in the general form

$$\widehat{H}=\widehat{H}_1+\widehat{H}_2=\underset{i,j}{\overset{N}{}}T_{ij}c_i^{}c_j+\frac{1}{2}\underset{i,j,k,l}{\overset{N}{}}V_{ijkl}c_i^{}c_j^{}c_kc_l,$$
(1)
where $`N`$ is the size of the chosen one-particle basis, and $`c_i^{}`$ and $`c_i`$ are the corresponding creation and annihilation operators. Both the one-body ($`T_{ij}`$) and two-body matrix elements ($`V_{ijkl}`$) are known.

To obtain the ground state $`|\mathrm{\Psi }_G`$ of $`\widehat{H}`$, QMC methods use the imaginary time propagator $`e^{\tau \widehat{H}}`$ acting on a trial wave function $`|\mathrm{\Psi }_T`$: $`lim_n\mathrm{}(e^{\tau \widehat{H}})^n|\mathrm{\Psi }_T|\mathrm{\Psi }_G`$. $`|\mathrm{\Psi }_T`$ must not be orthogonal to $`|\mathrm{\Psi }_G`$, and we will assume that it is of the form of a single Slater determinant or a linear combination of Slater determinants. The time step $`\tau `$ is chosen to be small enough so that $`\widehat{H}_1`$ and $`\widehat{H}_2`$ in the propagator can be accurately separated with the Trotter decomposition.

The propagator $`e^{\tau \widehat{H}_1}`$ is the exponential of a one-body operator. A propagator of this form acting on a Slater determinant is straightforward to calculate, and it simply yields another determinant. The two-body propagator $`e^{\tau \widehat{H}_2}`$ can be expressed as an integral of propagators of this form, as follows. Any two-body operator can be written as a quadratic form of one-body operators: $`\widehat{H}_2=\frac{1}{2}_\alpha \lambda _\alpha \widehat{v}_{\alpha }^{}{}_{}{}^{2}`$, where $`\lambda _\alpha `$ is a real number and $`\widehat{v}_\alpha `$ is a one-body operator. The Hubbard-Stratonovich (HS) transformation HS  then allows us to write

$$e^{\tau \widehat{H}_2}=\underset{\alpha }{}\left(\frac{1}{\sqrt{2\pi }}_{\mathrm{}}^{\mathrm{}}e^{\frac{1}{2}\sigma _\alpha ^2}e^{\sqrt{\tau }\sigma _\alpha \sqrt{\lambda _\alpha }\widehat{v}_\alpha }𝑑\sigma _\alpha \right).$$
(2)
Introducing vector representations $`\sigma \{\sigma _1,\sigma _2,\mathrm{}\}`$ and $`\widehat{𝐯}=\{\sqrt{\lambda _1}\widehat{v}_1,\sqrt{\lambda _2}\widehat{v}_2,\mathrm{}\}`$, we have the desired form

$$e^{\tau \widehat{H}}=P(\sigma )B(\sigma )𝑑\sigma ,$$
(3)
where $`P(\sigma )`$ is the normal distribution in Eq. (2) and

$$B(\sigma )e^{\tau \widehat{H}_1/2}e^{\sqrt{\tau }\sigma \widehat{𝐯}}e^{\tau \widehat{H}_1/2}$$
(4)
is a one-body propagator.

The imaginary-time propagation thus requires evaluating the multidimensional integral in Eq. (3) over time slices $`n`$ and the corresponding auxiliary fields. MC techniques are the only way to evaluate such integrals efficiently. We use a random walk approach Zhang . In each step, a walker $`|\varphi `$, which is a single Slater determinant, is propagated to a new position $`|\varphi ^{}`$: $`|\varphi ^{}(\sigma )=B(\sigma )|\varphi `$, where $`\sigma `$ is a random variable sampled from $`P(\sigma )`$. After a sufficient number of steps (iterations), the ensemble of random walkers is a MC representation of the ground-state wave function: $`|\mathrm{\Psi }_G_\varphi ^{}|\varphi ^{}`$.

In general $`\lambda _\alpha `$ cannot be made all positive in Eq. (2) note\_shift\_potential . The one-body operators $`\widehat{𝐯}`$ are therefore complex. As the projection proceeds, the orbitals in the random walkers will become complex. As a result, the statistical fluctuations in the MC representation of $`|\mathrm{\Psi }_G`$ increase exponentially with projection time $`\beta n\tau `$. This is the phase problem referred to earlier. It is of the same origin as the sign problem that occurs when $`B(\sigma )`$ is real. The phase problem is more severe, however, because for each $`|\varphi `$, instead of a $`+|\varphi `$ and $`|\varphi `$ symmetry Zhang , there is now an infinite set $`\{e^{i\theta }|\varphi \}`$ ($`\theta [0,2\pi )`$) from which the random walk cannot distinguish. At large $`\beta `$, the phase of each $`|\varphi `$ becomes random, and the MC representation of $`|\mathrm{\Psi }_G`$ becomes dominated by noise. This problem is generic, and the same analysis would apply if we had chosen, instead of the random walk, the standard AF QMC sampling approach Koonin . In Fig. 1, the curves labeled “free projection” illustrate the phase problem.

Existing fixed-node type approximations have often worked very well to control the sign/phase problem in real space Foulkes01 ; Ortiz93  or in Slater determinant space when the propagator is real Zhang . The phase problem here is unique because not only do the determinants acquire overall phases, but the internal structures of their orbitals become complex. The real-space analogy would be to have walkers whose coordinates become complex. This makes straightforward generalization of existing approaches ineffective. For example, similar to the constrained path approximation Zhang  we could impose the condition $`\mathrm{Re}\mathrm{\Psi }_T|\varphi >0`$. Or, in the spirit of the fixed-phase approximation in real space Ortiz93  we could project the walker by including a factor $`\mathrm{cos}(\mathrm{\Delta }\theta )`$ in the weight, where $`\mathrm{\Delta }\theta `$ is the phase of $`\mathrm{\Psi }_T|\varphi ^{}/\mathrm{\Psi }_T|\varphi `$. They give similar results and do not work well Wilson95 . The former is shown in Fig. 1 (“simple constraint”). Importance sampling with $`\mathrm{Re}\mathrm{\Psi }_T|\varphi `$ or $`|\mathrm{\Psi }_T|\varphi |`$ does not change the results ZKunpublished .

To formulate a new method that can better separate the overall phase from the determinant, we first borrow from the idea of importance sampling Kalos74 , although our choice of the so-called importance function, $`\mathrm{\Psi }_T|\varphi `$, is actually complex. We modify Eq. (3) to obtain the following new propagator for $`|\varphi `$:

$$\mathrm{\Psi }_T|\varphi ^{}(\sigma \overline{\sigma })P(\sigma \overline{\sigma })B(\sigma \overline{\sigma })\frac{1}{\mathrm{\Psi }_T|\varphi }𝑑\sigma ,$$
(5)
where we have included a constant shift shiftcont\_rom97  $`\overline{\sigma }`$ in the integral in Eq. (3), which does not affect the equality. Eq. (5) can be re-written as

$$P(\sigma )W(\sigma ,\varphi )B(\sigma \overline{\sigma })𝑑\sigma ,$$
(6)
where

$$W(\sigma ,\varphi )\frac{\mathrm{\Psi }_T|\varphi ^{}(\sigma \overline{\sigma })}{\mathrm{\Psi }_T|\varphi }e^{\sigma \overline{\sigma }\frac{\overline{\sigma }\overline{\sigma }}{2}}.$$
(7)
The new propagator in Eq. (6) defines a new random walk. In each step the walker $`|\varphi `$ is propagated to $`|\varphi ^{}`$ by $`B(\sigma \overline{\sigma })`$: $`|\varphi ^{}(\sigma \overline{\sigma })=B(\sigma \overline{\sigma })|\varphi `$, where $`\sigma `$ is again sampled from $`P(\sigma )`$. $`W(\sigma ,\varphi )`$ is a c-number which can be accounted for by having every walker carry an overall weight factor and updating them according to: $`w_\varphi ^{}=W(\sigma ,\varphi )w_\varphi `$. Formally the MC representation of $`|\mathrm{\Psi }_G`$ in the new random walk is:

$$|\mathrm{\Psi }_G\underset{\varphi ^{}}{}w_\varphi ^{}\frac{|\varphi ^{}}{\mathrm{\Psi }_T|\varphi ^{}}.$$
(8)
For any choice of the shift $`\overline{\sigma }`$, the new random walk is an exact procedure to realize the imaginary time propagation, in the sense of Eq. (8). The optimal choice of $`\overline{\sigma }`$ is determined by minimizing the fluctuation of $`W(\sigma ,\varphi )`$ with respect to $`\sigma `$. To $`𝒪(\sqrt{\tau })`$ this yields

$$\overline{\sigma }=\sqrt{\tau }\frac{\mathrm{\Psi }_T|\widehat{𝐯}|\varphi }{\mathrm{\Psi }_T|\varphi }.$$
(9)
With this choice the leading $`\sigma `$-dependent term in $`W`$ is reduced to $`𝒪(\tau )`$ and, by expanding $`B(\sigma \overline{\sigma })`$ in $`|\varphi ^{}`$ in Eq. (7), we can manipulate $`W`$ into the following form:

$$W(\sigma ,\varphi )\mathrm{exp}\left[\tau \frac{\mathrm{\Psi }_T|\widehat{H}|\varphi }{\mathrm{\Psi }_T|\varphi }\right]\mathrm{exp}[\tau E_L(\varphi )],$$
(10)
where the term $`E_L`$ parallels the local energy in real-space QMC methods. Both $`E_L`$ and the shift $`\overline{\sigma }`$ in Eq. (9) are independent of any overall phase factor of $`|\varphi `$.

The weight of the walker in the new random walk is determined by $`E_L`$. In the limit of an exact $`|\mathrm{\Psi }_T`$, $`E_L`$ is a real constant, and the weight of each walker remains real. The so-called mixed estimate for the energy is phaseless:

$$E_G=\frac{\mathrm{\Psi }_T|\widehat{H}|\mathrm{\Psi }_G}{\mathrm{\Psi }_T|\mathrm{\Psi }_G}\frac{\underset{\varphi ^{}}{}w_\varphi ^{}E_L(\varphi ^{})}{_\varphi ^{}w_\varphi ^{}}.$$
(11)
With a general $`|\mathrm{\Psi }_T`$ which is not exact, a natural approximation is to replace $`E_L`$ in Eq.’s (10) and (11) by its real part, $`\mathrm{Re}E_L`$. We have thus obtained a phaseless formalism for the random walk, with real and positive weights in Eq.’s (8) and (11).

Despite this, an additional constraint is still required. To illustrate the problem we consider the overlap $`\mathrm{\Psi }_T|\varphi ^{}`$ during the random walk. Let us denote the phase of $`\mathrm{\Psi }_T|\varphi ^{}(\sigma \overline{\sigma })/\mathrm{\Psi }_T|\varphi `$ by $`\mathrm{\Delta }\theta `$, which is in general non-zero (of order $`\sigma \mathrm{Im}\overline{\sigma }`$). This means that, the walkers will undergo a random walk in the complex plane defined by $`\mathrm{\Psi }_T|\varphi ^{}`$. At large $`\beta `$ they will therefore populate the complex plane symmetrically, independent of their initial positions. It is useful to contrast the situation with the special case of a real $`\widehat{𝐯}`$. For any $`\widehat{𝐯}`$ the shift $`\overline{\sigma }`$ diverges as a walker approaches the origin in the complex plane, i.e., as $`\mathrm{\Psi }_T|\varphi ^{}0`$. The effect of the divergence is to move the walker away from the origin. With a real $`\widehat{𝐯}`$, $`\mathrm{\Delta }\theta =0`$ and the random walkers move only on the real axis. If they are initialized to have positive overlaps with $`|\mathrm{\Psi }_T`$, $`\overline{\sigma }`$ will ensure that the overlaps remain positive throughout the random walk, much like in fixed-node diffusion Monte Carlo (DMC) in real space. Thus in this case the phaseless formalism reduces to the constrained path Monte Carlo method of Ref. Zhang , and it alone is sufficient to control the sign problem. For a complex $`\widehat{𝐯}`$, however, the random walk is “rotationally invariant” in the complex plane, and the divergence of $`\overline{\sigma }`$ is not enough to prevent the build-up of a finite density at the origin. Near the origin the local energy $`E_L`$ diverges, which causes diverging fluctuations in the weights of walkers. To address this we make an additional approximation. We project the random walk to “one-dimension” and multiply the weight of each walker in each step by $`\mathrm{max}\{0,\mathrm{cos}(\mathrm{\Delta }\theta )\}`$. Imposing instead $`\mathrm{Re}\mathrm{\Psi }_T|\varphi ^{}>0`$ gave similar results, but with somewhat larger variance.

We apply the new method to Si atom, molecule, and bulk. The Si4+ ions are represented by a norm-conserving LDA Kleinman-Bylander (KB) non-local pseudopotential rappe . We use periodic boundary conditions, and a plane-wave basis with a kinetic energy cut-off $`E_{\mathrm{cut}}=12.25`$ Ry. The error resulting from $`E_{\mathrm{cut}}`$ was estimated through LDA calculations and is smaller than the MC statistical errors. The pseudopotential can be applied in essentially the same way as in plane-wave-based LDA calculations ZKunpublished . Calculations involving $`\widehat{𝐯}`$ and the local part of the pseudopotential are efficiently handled using fast Fourier transforms. The separable KB form of the non-local pseudopotential makes its application as efficient as in LDA plane-wave codes. Our $`|\mathrm{\Psi }_T`$ is a single Slater determinant consisting of LDA orbitals.

In Table 1, we show results for the atom and molecule. Additional calculations with $`a=22a_B`$ supercells show that finite-size errors at $`a=19a_B`$ were smaller than the MC statistical errors. Our calculated Si2 binding energy is in excellent agreement with the experimental value Schmude .

In the bulk calculations, we use fcc supercells consisting of 2, 16, 54 atoms (5209 plane waves). As Fig. 1 shows, the new method leads to a large improvement. Results for 16 and 54 atoms are shown in Table 2. Our calculation for 54 atoms took several days on 20 Compaq Alpha 667 MHz processors. For the bulk cohesive energy, we first included a correction for the independent-particle finite-size error from the LDA results. We then corrected for the remaining Couloumb finite-size error note\_CFSE  using the results of Kent et al. Kent99 . Our result is again in excellent agreement with the experimental value (from Ref. Foulkes01 ). It also compares very well with the result of a recent fixed-node DMC calculation leung99 , which also used a 54-atom supercell and gave $`4.63(2)`$ eV per atom after similar finite-size and zero-point energy corrections.

Without an exact solution to the sign/phase problem, reducing the reliance on trial wave functions is clearly of key importance to increasing the predictive power of QMC. For continuum electronic systems such as our test cases above, fixed-node DMC has often been the most accurate theoretical method Foulkes01 . It is encouraging that the new method, using simple LDA trial wave functions, gave comparable results to DMC. For similar supercells DMC often uses trial wave functions with 30-100 additional parameters Foulkes01 . Obtaining a good enough $`|\mathrm{\Psi }_T`$ is instrumental to a successful DMC calculation, and often constitutes a substantial effort. The quality of $`|\mathrm{\Psi }_T`$ controls the systematic errors from the fixed-node approximation and the variance. It also affects errors due to the locality approximation Mitas91 , which has been employed by most DMC calculations with non-local pseudopotentials. In the new method the latter approximation is eliminated. It remains to be seen whether the present method could lead to more accurate results than fixed-node DMC for continuum systems. This has been possible in some cases Zhang99  with real AF in the Hubbard model.

We have presented a general framework. Various possibilities exist for further improvement of the method. An improved $`|\mathrm{\Psi }_T`$ will give improved results. The freedom to choose the one-particle basis and the form of HS transformation, both of which can impact the quality of the results, offers significant opportunities. For periodic systems it should be possible to generalize the formalism to allow $`𝐤`$-point sampling.

In conclusion, we have described a method for ground-state QMC calculations that allows the use of any one-particle basis. The method is general and applies to any Hamiltonian of the form in Eq. (1). It provides an approximate way to control the phase problem in all AF-based QMC methods, while allowing many of their advantages to be retained that lead to their applications spanning several areas. We have shown that the method gave accurate results for systems from an atom to a large supercell, using a simple trial wave function.

We are very grateful to E. J. Walter for help with LDA calculations and with several programming issues. We thank P. Kent for sending us his DMC data and J. Carlson for useful conversations. This work is supported by the NSF under Grant No. DMR-9734041 and the Research Corporation, and by ONR Grant DMR-N000149710049. We acknowledge computing support by the Center for Piezoelectrics by Design.