text:
# Magnetic Collective Mode Dispersion in High Temperature Superconductors

## Abstract

Recent neutron scattering experiments in the superconducting state of YBCO have been interpreted in terms of a magnetic collective mode whose dispersion relative to the commensurate wavevector has a curvature opposite in sign to a conventional magnon dispersion. The purpose of this article is to demonstrate that simple linear response calculations are in support of a collective mode interpretation, and to explain why the dispersion has the curvature it does.

Inelastic neutron scattering, which measures the momentum and energy dependence of the dynamic spin susceptibility, has emerged as a powerful probe of high temperature cuprate superconductors. One of the most striking features of the data is the appearance below the superconducting transition temperature of a commensurate resonance in the case of the bilayer cuprates, YBCO and Bi2212. In addition, an incommensurate magnetic response has been seen for energies below the resonance energy in both cases. Recently, it has been argued that these two structures are part of the same collective mode dispersion. This interpretation is unusual in that the curvature of the purported collective mode dispersion relative to the commensurate wavevector is opposite in sign to that of a conventional magnon dispersion.

One of the motivations of this interpretation is the similarity of the measured dispersion to theoretical results based on linear response theory (RPA) . In these works, though, it was unclear whether the incommensurate response below resonance was actually due to a pole in the response function as would be expected for a collective mode interpretation. The purpose of this article is to investigate this issue in detail.

The methodology is standard, and details can be found in Ref. . An effective quasiparticle dispersion is used to calculate the bare bubble in the superconducting state, $`\chi _0`$. The RPA susceptibility is then obtained as $`\chi _0/(1J\chi _0)`$. The most crucial input is the form of the dispersion. In Ref. , two dispersions were analyzed, one based on angle resolved photoemission (ARPES) data in the normal state (tb1), another on ARPES data in the superconducting state where the Fermi surface was flattened around the node to enhance the effects of incommensurability (tb2). In addition, in this paper, a slight modification of this second dispersion was made to further enhance the incommensurability (tb3). All three dispersions are listed in Table I for completeness sake. For the superconducting gap, the standard $`\mathrm{cos}(k_x)\mathrm{cos}(k_y)`$ form was assumed, with a maximum value, $`\mathrm{\Delta }_{max}`$, taken from recent tunneling measurements in YBCO. $`J`$ was adjusted to obtain a resonance at 41 meV, the resulting $`J`$ being similar in value to estimates based on neutron scattering in underdoped samples .

To begin, it is helpful to discuss the frequency dependence of the bare bubble, $`\chi _0`$. In Figure 1, $`Re\chi _0`$ and $`Im\chi _0`$ are shown at a typical incommensurate wavevector for the three dispersions of Table 1. In $`Im\chi _0`$, there are two steps, a low energy one at threshold, and another at slightly higher energy. The two energies correspond to the spin gap energies $`\omega _g(\pm q)=\mathrm{min}_k(E_k+E_{k\pm q})`$, where $`k`$ is confined to the first quadrant of the Brillouin zone and $`E_k`$ are the quasiparticle energies in the superconducting state. The steps are a consequence of the coherence factors associated with a d-wave order parameter. For the commensurate wavevector, $`Q=(\pi ,\pi )`$, these two energies are degenerate, so there is only one step. By the Kramers-Kronig relations, these steps lead to logarithmic divergences in $`Re\chi _0`$ which become simple peaks when damping is taken into account. (All plotted results are obtained by replacing $`\omega `$ by $`\omega +i\mathrm{\Gamma }`$ in the energy denominators when calculating $`\chi _0`$).

In Figure 2, $`Re\chi _0`$ and $`Im\chi _0`$ are shown as a function of $`q`$ at an energy of 35 meV for the three dispersions of Table 1. In the first case (tb1), one finds the typical result that $`Re\chi _0`$ is quadratic relative to the commensurate wavevector $`Q`$ with a negative curvature (i.e., a maximum at $`Q`$). Since $`Re\chi _0`$ increases with $`\omega `$ (up to the spin gap energy $`\omega _g(Q)`$), one expects that the RPA pole condition, $`1JRe\chi _0=0`$, will be satisfied first at $`Q`$, then at higher energies at wavevectors increasingly displaced from $`Q`$. This would then give a conventional magnon like dispersion. The actual situation is more complicated, though. Note in Figure 2a that there is a second maximum in $`Re\chi _0`$ at an incommensurate wavevector, associated with $`\omega _g(q)`$ as discussed in the context of Figure 1. As the energy increases, this second maximum becomes the global maximum.

In Figure 3a, the $`q`$ vector where the RPA $`Im\chi `$ is maximal at a given $`\omega `$ ($`\omega (q)`$) is plotted for this dispersion. In agreement with experiment, one finds incommensurability below the resonance energy, but the RPA pole condition is not satisfied below resonance. Instead the incommensurability simply tracks $`\omega _g(q)`$. As a consequence, the intensity plummets below resonance, as can be seen in Figure 4a. Above resonance, the RPA pole condition is satisfied, but as can be seen in Figure 3a, the dispersion crosses into the continuum (i.e., beyond $`\omega _g(q)`$) just above resonance, so the mode becomes damped, as obvious from the intensity plot of Figure 4a.

The above behavior is not what is observed, in that the experimental intensity decays more slowly below resonance. This is because the incommensurability effects associated with dispersion tb1 are too weak. This can be contrasted with the dispersion analyzed by Brinkmann and Lee . In that case, the constant energy contours of $`E_k`$ are quite flat for energies below $`\mathrm{\Delta }_{max}`$, leading to an enhancement in the incommensurate effects. This can be seen in Figures 2b and 2c, where $`Re\chi _0`$ is plotted for the other two dispersions considered here (tb2 and tb3). In the first case (tb2), the dispersion used in earlier work , one finds a rather flat behavior around $`Q`$ due to the flat quasiparticle dispersion around $`k=(\pi ,0)`$. In addition, a global maximum now occurs at an incommensurate wavevector due to the flattening of the Fermi surface around the d-wave nodal $`(\pi ,\pi )`$ directions relative to dispersion tb1.

The resulting dispersion, $`\omega (q)`$, of $`Im\chi `$ from tb2 is plotted in Figure 3b. This is in better agreement with experimental results in YBCO . In addition, it can be seen from this plot that the RPA pole condition is now satisfied below the resonance energy for a few meV until the spin gap energy, $`\omega _g(q)`$, is encountered. For energies below this, the mode is strongly diminished in intensity as it traces out the edge of the continuum. This result is consistent with earlier work .

Although Figure 3b is in good agreement with experiment as far as $`\omega (q)`$ is concerned, it is still deficient in that the intensity of $`Im\chi `$ below resonance still drops off faster than experiment, as plotted in Figure 4b and also noted earlier . To analyze this further, dispersion tb2 was modified by making the dispersion less flat in one direction around $`(\pi ,0)`$, and flattening the Fermi surface around the d-wave node even more. The resulting $`Re\chi _0`$ is shown in Figure 2c, and is similar to that obtained earlier by Brinkmann and Lee . One now finds a truly quadratic behavior around $`Q`$ with positive curvature, a direct consequence of the stronger incommensurate peaks. As a consequence, the RPA pole condition is satisfied for an even greater energy range below resonance (Figure 3c), and moreover, the intensity below resonance falls off much more slowly (Figure 4c), as observed experimentally.

The results in this paper concentrated on the $`(\pi ,q_y)`$ direction, since this is where the maximum in the incommensurability was observed. Two dimensional $`q`$ plots from dispersions tb2 and tb3 are consistent with experiment  in showing a baseball diamond shaped incommensurability pattern below resonance with global maxima along $`(\pi ,q_y)`$ and $`(q_x,\pi )`$, as reported earlier . Recently, an anisotropy was observed in this intensity pattern for partially detwinned YBCO samples . If this is a true magnetic anisotropy not related to the phonon background, then within the RPA context, it would have to be due to the influence of the chain bands wiping out two of the incommensurate spots, not surprising given the sensitivity of the incommensurability to the electronic structure. On the other hand, this is conjecture at this stage, as a proper calculation taking into account both plane and chain bands has yet to be performed. Of course, the most natural explanation of this anisotropy would be due to stripe formation, but it is unclear whether such a picture can reproduce the results demonstrated here (in particular, the energy dependence of the incommensurability).

In conclusion, RPA calculations of the dynamic susceptibility are in support of a collective mode interpretation of the magnetic dispersion relation observed in bilayer cuprate superconductors. The observed incommensurability has strong implications for the effective quasiparticle dispersion, and it would be of interest to verify whether the dispersions analyzed here, in particular in regards to the flattening of the Fermi surface around the d-wave node, are consistent with high resolution angle resolved photoemission data, particularly in the case of YBCO.

The author thanks Andrei Chubukov and Oleg Tchernyshyov for discussions. This work was supported by the U.S. Dept. of Energy, Basic Energy Sciences, under Contract No. W-31-109-ENG-38.


format: text

text:
# Critical Depinning Force and Vortex Lattice Order in Disordered Superconductors

## I Introduction

A large class of condensed matter systems can be represented as an elastic lattice interacting with a quenched random substrate. Vortices in type II superconductors, which can be pinned by random defect sites in the material, are a particularly ideal example of such a system since the substrate and lattice interactions can be readily tuned experimentally. In this case, a rich variety of static  and dynamical phases  can occur due to the competition between the repulsive vortex-vortex interactions, which tend to order the vortex lattice, and the attractive vortex-pin interactions, which tend to disorder the vortex lattice. Unlike the early studies of collective pinning of a defect-free elastic lattice, recent experimental and theoretical work point to the importance of topological defects in pinned vortex matter. If the pinning in the sample is sufficiently strong, the vortex lattice will be highly defected; conversely, in a sample with weak pinning the vortex lattice can be relatively ordered and contain only a few defects. In contrast with the purely elastic model, however, the relationship between the critical force and the density of flux line defects is not known. Moreover, applications of superconducting materials require that the critical current be high, and so a key question is how the number of defects in the vortex lattice can affect the critical current.

In a given sample, i.e., with a given realization of quenched disorder, the relative strength of the pinning can be increased by lowering the vortex-vortex interaction strength, i.e., by softening the lattice. Vortex-vortex interactions become weak at low applied fields, when vortices are dilute, as well as very near $`H_{c2}`$ and $`T_c`$, where the induction of the vortices overlap significantly. In the latter regime one encounters the peak effect: a rapid increase in the critical current before it collapses to zero at the superconducting - normal phase boundary. At these two extreme field values, the pinning interaction dominates and the vortex lattice is expected to be highly disordered. Indeed, direct evidence for vortex lattice disordering through the peak effect regime has been observed in neutron scattering measurements . A reentrant disordering has been found in the low field regime where the vortex lattice again softens . Bitter decoration experiments have also provided evidence for the disordering of the vortex lattice at low fields . In each of these cases the critical current $`J_c`$ is high due to the dominance of pinning.

Experiments on superconducting samples with different amounts of pinning have revealed some systematic trends. The peak effect is much sharper and more pronounced in cleaner samples , while in samples with stronger pinning, the width of the peak region increases and the relative valley-to-peak height decreases, although the overall critical current is higher . In samples with the strongest pinning, $`J_c`$ is high through most of the applied field range, and the peak is completely absent. A possible interpretation of this behavior connects vortex lattice defects with high critical currents. A detailed understanding of the observed effects as a possible transition/crossover between regimes with low and high critical currents, i.e., between, say, a Bragg glass-like ordered phase with few defects and a Vortex glass/ pinned liquid-like disordered vortex phase is still lacking. A key reason for this uncertainty is the absence of a systematic connection between the number of defects in the vortex lattice, the softness of the lattice, and the critical current.

In this work, we directly examine the effect of vortex lattice softness on the critical current through a series of numerical simulations of samples with different pinning strengths. For strong pinning the vortex lattice is highly defective for the entire range of vortex-vortex interactions investigated. The critical depinning force, which is proportional to the critical current, only weakly increases with decreasing vortex-vortex interaction strength, reaching a maximum for non-interacting vortices. For the weakest pinning strengths, on the other hand, the lattice is almost defect-free over a large range of vortex-vortex interaction strengths. When defects begin to appear in the softest vortex lattices, the critical depinning force increases sharply from its low value to a peak at zero vortex-vortex interaction. The behavior for intermediate pinning is intermediate between these extremes. We find that the relative height of the critical current peak increases for decreasing pinning while the width of the peak decreases. Thus, the most pronounced peaks occur for the weakest pinning. These results are in excellent agreement with the experimental trends described above. We discuss how these results can be connected to the behavior of the peak effect in superconductors with varying pinning strength.

## II Simulation

We consider a 2D slice of a system of superconducting vortices interacting with a random pinning background. The applied magnetic field $`𝐇=H\widehat{𝐳}`$ is perpendicular to our sample, and we use periodic boundary conditions in $`x`$ and $`y`$. The $`T=0`$ overdamped equation of motion for a vortex $`i`$ is

$$𝐟_i=\eta \frac{d𝐫_i}{dt}=𝐟_i^{vv}+𝐟_i^{vp}+𝐟_d=\eta 𝐯_i,$$
(1)
where $`𝐯_i`$ is the velocity of vortex $`i`$ and $`\eta =1`$ is the damping coefficient. The total force on vortex $`i`$ from the other vortices is $`f_i^{vv}=_{j=1}^{N_v}A_vf_0K_1(|𝐫_i𝐫_j|/\lambda )\widehat{𝐫_{\mathrm{𝐢𝐣}}}`$, where $`𝐫_i`$ is the position of vortex $`i`$, $`\lambda `$ is the penetration depth, $`f_0=\mathrm{\Phi }_0^2/8\pi ^2\lambda ^3`$, the prefactor $`A_v`$ is used to vary the vortex lattice softness, and $`K_1`$ is the modified Bessel function. The pinning $`𝐟_i^{vp}`$ is modeled as randomly placed attractive parabolic wells of radius $`r_p=0.15\lambda `$ with $`f_i^{vp}=(f_p/r_p)(|𝐫_i𝐫_k^{(p)}|)\mathrm{\Theta }(r_p|𝐫_i𝐫_k^{(p)}|)\widehat{𝐫}_{ik}^{(p)}`$, where $`𝐫_k^{(p)}`$ is the location of pin $`k`$, $`f_p`$ is the maximum pinning force, which is varied from $`0.1f_0`$ to $`3.0f_0`$, $`\mathrm{\Theta }`$ is the Heavside step function and $`\widehat{𝐫}_{ik}^{(p)}=(𝐫_i𝐫_k^{(p)})/|𝐫_i𝐫_k^{(p)}|`$. The pin density is $`n_p=3.0/\lambda ^2`$ and the vortex density is $`n_v=0.75/\lambda ^2`$. We simulate a $`36\lambda \times 36\lambda `$ system containing $`N_v=864`$ vortices and $`N_p=3887`$ pins. We initialize the vortex positions by performing simulated annealing, starting from a high temperature and slowly cooling to $`T=0`$. This method of preparing the lattice is similar to field cooled experiments. To identify the depinning force $`f_c`$ we apply a slowly increasing uniform driving force $`𝐟_d`$ on the vortices in the $`x`$ direction, which would correspond to a Lorentz force from an applied current $`𝐉=J\widehat{𝐲}`$. For each drive increment we measure the average vortex

velocity in the direction of drive, $`V_x=(1/N_v)_1^{N_v}v_x`$. The $`f_d`$ versus $`V_x`$ curve corresponds experimentally to a $`V(I)`$ curve. The depinning force $`f_c`$ is defined as the drive at which $`V_x>0.03`$.

## III Vortex Order and Pinning for Varied Vortex-Vortex Interaction Strength

We first consider the effect of the vortex lattice softness on the stationary vortex lattice. In Fig. 1 we show the Delaunay triangulation for a system with $`f_p=0.25f_0`$ for decreasing vortex-vortex interaction $`A_v`$, after the

lattice has been annealed and with no driving force applied. Defect sites in the vortex lattice are indicated by circles. In Fig. 1(a) for $`A_v=4.0`$, the vortex lattice contains no defects. In Fig. 1(b) for $`A_v=3.0`$ a small number of 5-7 defect pairs appear. In Fig. 1(c), at $`A_v=2.0`$ the vortex lattice is considerably disordered but some domains of order are still present. For softer vortex lattices in Fig. 1(d-f), the vortex lattice becomes progressively more disordered by the underlying random pinning. Since the simulation is for a two-dimensional system, we expect the system to be defective for any pinning strength. The defect-free situation likely results from the finite size of the sample.

By applying a transport current to the annealed lattices, we determine the critical depinning force $`f_c`$ that must be applied before the vortices begin to move. In Fig. 2(a) we plot $`f_c`$ versus $`A_v`$ for a sample with $`f_p=0.25f_0`$. For comparison, we measure the amount of order in the lattice before depinning using the Delaunay triangulation, and in Fig. 2(b) we plot $`P_6`$, the fraction of six-fold coordinated vortices, versus $`A_v`$. Here $`P_6=1.0`$ indicates a perfect triangular lattice. As $`A_v`$ is lowered $`f_c`$ increases and simultaneously the order in the lattice, $`P_6`$, decreases. The maximum value of $`f_c`$ occurs at $`A_v0`$ which coincides with the minimum value of $`P_6`$. Thus the softer lattices with low values of $`A_v`$ and large amounts of disorder are more strongly pinned than stiffer, more ordered lattices.

## IV Height and Width of Critical Depinning Force Peak

The shape and magnitude of the critical current peak at $`A_v=0`$ are affected by the strength of the pinning in the sample. To demonstrate this, in Fig. 3(a) we plot $`f_c/f_p`$ versus $`A_v`$. Scaling the curves with $`f_p`$ in this way causes all of the curves to approach $`f_c/f_p=1`$ at $`A_v=0`$. We find that the relative height of the peak increases as the pinning becomes weaker. For example, comparing $`f_p=0.25`$ (bottom curve) with $`f_p=3.0`$ (top curve), the critical current $`f_c`$ increases by a factor of 8 from $`f_c=0.1`$ at $`A_v=6.0`$ to $`f_c=0.8`$ at $`A_v=0.01`$ in the weakly pinned sample, whereas in the strongly pinned sample, $`f_c`$ increases only by a factor of 2 from $`f_c=0.38`$ at $`A_v=6.0`$ to $`f_c=0.8`$ at $`A_v=0.01`$. The peak not only becomes higher for weaker pinning, but also it becomes much narrower as can be seen by comparing the widths of the peaks in Fig. 3(a). As shown in Fig. 3(b), in the strongly pinned samples such as $`f_p=3.0f_0`$ with broad low peaks, $`P_6`$ is low over the entire range of $`A_v`$, indicating that the vortex lattice is always strongly disordered by the pinning. In contrast, for weakly pinned samples such as $`f_p=0.25f_0`$, $`P_6`$ drops a large amount from $`P_61`$ at large $`A_v`$, indicating a nearly perfect

lattice, to $`P_60.4`$ at the lowest $`A_v`$, indicating a large amount of disorder in the lattice.

In Fig. 4 we show explicitly how the relative peak height, $`H=f_c(A_v=0.01)/f_c(A_v=6.0)`$, and the width at half maximum, $`dW`$, vary with pinning strength $`f_p`$. Here we see that the relative strength and sharpness of the peak increase with weaker pinning. The width $`dW`$ depends linearly on $`f_p`$, while the relative height $`H`$ shows a nonlinear increase that can not be fitted by a simple functional form.

The results in Fig. 3 and Fig. 4 suggest that in clean superconductors, where the vortex lattice softness has a much more significant impact on the effectiveness of the pinning, the critical current enhancements should be both sharper and of larger relative height than in strongly pinned superconductors. This behavior agrees with the results of recent experiments  on the peak effect. Our results also suggest that the peak can be made arbitrarily sharp simply by adjusting the disorder strength, and that such large peaks in the critical current must be accompanied by the proliferation of defects.

## V Critical Depinning Force Dependence on Defect Density

Using the data obtained by varying $`A_v`$, we can determine the relationship between the critical current and the density of defects. In Fig. 5(a) we plot $`f_c/f_p`$ versus $`P_6`$, showing that the curves for samples with different

pinning strengths roughly collapse. At lower values of $`P_6`$ when the vortex lattice becomes defected, the critical current increases rapidly with the onset of defects. In Fig. 5(b) we plot $`f_c/f_p`$ versus the distance between vortex lattice defects, $`d_d=1/\sqrt{1.0P_6}`$. In samples with weak pinning there is a range of high $`A_v`$ values that produce well ordered vortex lattices. In this regime, when there are very few defects in the vortex lattice which are far apart, we find that the $`f_c/f_p`$ values do not collapse, and that the rate at which $`f_c/f_p`$ changes with $`d_d`$ is slow. For $`d_d<3.0`$, or when the lattice contains approximately 10% or more defects, all the curves collapse and $`f_c/f_p`$ increases rapidly with decreasing $`d_d`$. The behavior of $`f_c`$ versus $`d_d`$ suggests that there is a critical distance between defects below which $`J_c`$ begins to rapidly increase.

We can compare the dependence of the critical depinning force on the defect density with the collective pinning results of the Larkin-Ovchinnikov (LO) theory . In Fig. 6(a) we plot $`f_c`$ as predicted by LO, using $`f_cf_p^2/A_v`$. Fig. 6(b) shows $`f_c`$ computed from the distance between dislocations, using $`f_cf_p/d_d`$. The value of $`f_c`$ measured in the simulations is shown in Fig. 6(c). We compare $`f_c`$ obtained by the LO and dislocation distance methods with the actual $`f_c`$ in Fig. 6(c). The agreement with the LO prediction is poor, but that with the dislocation distance is good out to values $`A_v/f_p5`$, beyond which $`d_d`$ becomes of order the system size.

## VI Conclusion

We have investigated the dependence of the critical depinning force $`f_c`$ and vortex lattice topological order $`P_6`$ on the vortex lattice rigidity for different values of pinning strength. The study is restricted to the two dimensional case which is appropriate for vortex lattices where the longitudinal correlation length is large enough to be effectively given by the sample thickness. In all cases, $`f_c`$ increases with decreasing vortex-vortex interaction strength, reaching a maximum for non-interacting vortices, as expected. Although for very weak pinning and strong vortex-vortex interaction, the lattice is relatively defect-free, for strong pinning the lattice is defective for the entire vortex interaction range investigated. These results are in accord with previous studies. We have, in addition, studied the systematic variation of $`f_c`$ with varying pinning and inter-vortex interaction. The results show that the pinning force variation departs significantly from the expectations of an elastic picture. The increase in $`f_c`$ towards the maximum occurs very rapidly with decreasing interaction or increasing pinning. The relative height increases while the width decreases with decreasing pinning strength, strongly remniscent of the peak effect phenomenon. Furthermore, a pseudo-Larkin picture where a plastic length, the distance between topological defects, replaces the elastic correlation length, provides a good account of the variation of $`f_c`$. We also find evidence that there is a critical value of the plastic length, typically spanning 3-5 vortices on each side, at which the rapid crossover occurs. These results should be relevant to experiments in very weak pinning (quasi-2D) flux lattices, as well as to the general systems of two-dimensional lattices with quenched disorder where the lattice interaction can be tuned.


text:
# First principles study of the adsorption of C60 on Si(111)

\[

## Abstract

The adsorption of C60 on Si(111) has been studied by means of first-principles density functional calculations. A 2x2 adatom surface reconstruction was used to simulate the terraces of the 7x7 reconstruction. The structure of several possible adsorption configurations was optimized using the ab initio atomic forces, finding good candidates for two different adsorption states observed experimentally. While the C60 molecule remains closely spherical, the silicon substrate appears quite soft, especially the adatoms, which move substantially to form extra C-Si bonds, at the expense of breaking Si-Si bonds. The structural relaxation has a much larger effect on the adsorption energies, which strongly depend on the adsorption configuration, than on the charge transfer.

\]

Recent studies using scanning tunneling microscopy (STM) have reported images that reveal intra-molecular features of C60 molecules deposited, at very low coverages, on Si(100)-(2$`\times `$1) , and Si(111)-(7$`\times `$7) . These observations are very interesting as they might be used to infer the orientation of the molecules , and other details of the bonding configuration. They also confirm the strong interaction of the fullerenes with these surfaces. Furthermore, in the images reported by Pascual et al.  it is clear that there exist two types of molecules: “large” molecules with an apparent height of $``$0.6 nm and width of $``$2.0 nm, which appear more round and fuzzy in the constant current STM images, and “small” molecules with a height of $``$0.5 nm and width of $``$1.5 nm, which present a more clearly discernible internal structure and a larger variety of shapes. A similar observation of two different adsorption states has also been done by Yao et al  on the Si(100)-(2$`\times `$1) substrate. Interestingly enough, after annealing at 870 K the “large” molecules evolve to “small” ones, indicating that they were probably in a weaker adsorption state.

In this paper we apply first-principles electronic structure methods to study the adsorption of fullerenes on the terraces of the Si(111)-(7$`\times `$7) reconstruction, identifying candidates for the two adsorption states observed. There are some previous electronic structure calculations of fullerenes over silicon substrates . However, in the sole case where a structural optimization was performed , a semi-empirical force field was used. Therefore, to the best of our knowledge this is the first systematic study of the adsorption geometry of C60 on a silicon surface using ab initio atomic forces.

The calculations have been performed with the SIESTA program , which allows standard calculations within density functional theory  (DFT) for systems with hundreds of atoms. It uses norm-conserving pseudopotentials , and a basis set of numerical atomic orbitals, obtained from the solution of the atomic pseudopotentials at a slightly excited energy .*** We have used here an ‘energy shift’ of 0.5 eV. The corresponding radii are 4.65 and 5.73 a.u. for the $`s`$ and $`p`$ orbitals of Si, 3.79 and 4.41 in the case of C, and 4.26 for the $`s`$ states of H.  In this work we have used the local density approximation to DFT  and a minimal $`sp^3`$ basis set for both C and Si. With this basis, the Si-Si and C-C bond lengths differ only 1% and 2% from experiment, respectively. Our method also uses a real-space grid to compute the Hartree and the exchange correlation contributions to the self-consistent potential and the Hamiltonian matrix. Real and reciprocal-space integrations were performed with a 70 Ry-cutoff grid and with 2 inequivalent k-points. For a given basis set, these meshes guarantee a convergence better than $``$ 1 meV/atom for the total energies and of $``$0.1 eV for the C60 binding energies  The binding energies have been corrected for the basis set superposition error as described, for example, in Ref.

The terraces of the Si(111)-(7$`\times `$7) surface have been modeled by slabs of two double layers, covered with a 2$`\times `$2 adatom reconstruction on one face, and saturated with hydrogen on the other face. The C60 molecules are arranged in a $`2\sqrt{3}\times 2\sqrt{3}`$ periodic supercell, big enough to avoid strong interactions between them. First, we relax the structure of the clean surface with the lattice constant parallel to the surface fixed at the calculated bulk value, to avoid artificial stresses. The relaxed structure is similar to that of other authors . The distance between the pedestal atoms (restatoms bound to the adatoms) is reduced, and a strong downward relaxation for the atoms beneath the adatom is observed. The rest atom lies $``$0.4 Å higher than the ideal bulk height. The main discrepancy with the results of Ref.  is that the position of the adatoms is approximately 9% higher in our case. This difference may be due to the use of a minimal basis. However, this is not crucial, since we are interested in the changes induced by the adsorption of the C60, rather than in the surface itself. Our results also reproduce the main features of the electronic structure of the surface.

First we place the C60 molecule on top of one of the rest atoms, with different orientations, suggested by the STM images of Pascual et al. . In one of them (labelled S3) the molecule is oriented with an hexagon on top, to preserve the threefold symmetry of the substrate. In the other orientation, a C60 double-bond (DB) is placed directly over the rest atom. Fig 1(a), shows the binding energy as a function of the height of the center of mass of the fullerene relative to the initial position of the rest atom. If no internal relaxation of the molecule or the surface are allowed, the binding energies are 0.64 eV and 1.47 eV for the S3 and DB configurations, respectively. The DB molecule is more bound because it is lower and closer to the rest atom.

We then investigate the structural changes, shown in Fig 2, in both the surface and the molecule after the adsorption, relaxing the atomic positions according to the ab initio forces. The carbon atoms move very little from their ideal positions, the molecule remains almost spherical, and its height over the substrate decreases from $``$5.9 to 5.35 Å. In contrast, the three nearest adatoms to the S3 molecule move substantially, approaching the closest carbon atoms from $``$3.0 Å to 2.02 Å. The rest atom below the molecule relaxes downward 0.61 Å, while the pedestal atoms closer to the molecule relax upwards 0.08 Å, and their bond distance to the adatom decreases from 2.57 Å to 2.51 Å. The pedestal atoms more distant to the molecule relax upwards 0.24 Å, and their bond length to the adatoms becomes 2.75 Å. The relaxation energy is large, increasing the binding energy to 3.96 eV.

The case of the DB molecule in rather different. It rotates during the relaxation, binding to the rest atom and to two adatoms. The bond length with the adatoms is 2.01 Å, while the bond with the rest atom is a little shorter (1.98 Å) and stronger (larger bond charges). The relaxation of these atoms, and those bound to them, is similar to that for the S3 molecule, but the displacements are somewhat smaller. The deformation of the molecule is again very small. The third adatom is largely displaced towards the C60, but the large distance to the molecule (2.63 Å), and the very low bond charge do not indicate the formation of covalent bonding. Its displacement is probably driven by the movement of the pedestal atoms, and the ionic interaction with the negatively charged molecule. The equilibrium height is reduced from $``$5.5 to 5.39 Å. The DB molecule has a binding energy of 5.96 eV, quite larger than the S3 molecule. This difference is due to the stronger bonds formed and to the somewhat less strained surface.

Fig. 1 shows data for two other relaxed structures. The structure labelled as SB (not shown) has initially a single-bond pointing towards the rest atom. Like the DB structure, it forms, after the relaxation, a strong convalent bond with the rest atom. However, in this case the adatom more distant to the molecule suffers a larger displacement, breaking its bond with one of the pedestal atoms, and forming another covalent bond with the fullerene. This fourth bond is longer (2.31 Å) and weaker than the other three. In fact, in spite of the formation of this extra-bond, the binding energy (5.93 eV) and its height (5.41 Å) are almost identical to those of the DB molecule. This indicates that the adatom most loosely bound to the molecule can probably oscillate between two positions with the same energy: one where it is covalently bound to the fullerene, and another where the interaction with the molecule is mainly ionic.

The S3 configuration might be a good candidate for the “large” (less stable) molecules observed by STM , whose resolved internal structure present a triangular shape , consistent with the threefold symmetry. The DB and SB configurations are more bound to the substrate, but their height is almost identical to that of S3. Although STM does not directly measure the atomic positions, our estimations, using the Tersoff-Hamann  theory, also indicate a similar “electronic height” for the DB and S3 configurations ($``$6.8Å). A reasonable candidate for the “small” molecules should have a higher binding energy and a lower heigth over the substrate than the structures studied to this point. The relaxed structure labelled L, which fulfills these conditions, is shown in Fig. 3. This is just an example of the several configurations, with different molecular orientations, that can be obtained by allowing a larger rearragement of the substrate than the structures previously studied. These new structures are characterized by a lower height, and the formation of more bonds between the molecule and the surface atoms. The L configuration has a binding energy of 6.36 eV (0.4 eV larger than the DB), its height is 4.56 Å ($``$0.8 Å lower than the other configurations). The molecule presents an hexagon on the top and occupies almost a “bridge” position between two adatoms. Each of these adatoms has a broken bond with the substrate, allowing for the formation of an extra-bond with the molecule (the new bonds have been highlighted in Fig. 3). The molecule has five bonds with the substrate, two with the adatoms (2.13 Å), two with the freed pedestal atoms (1.98 and 2.01 Å), and one with the rest atom (1.96 Å).

Our findings indicate that adsorption states like S3, DB, or SB, can be accessed whitout overcoming any energetic barrier, in other words we have not found any physisorbed precursor for the adsorption of C60 over Si(111). However, it must exist an energetic barrier between the “large” and “small” molecules observed experimentally, so both can be simultaneously visualized at room temperature. The nature of the barrier has been clarified by our calculations: some of the bonds of the adatoms with the substrate have to be broken to allow the formation of extra-bonds with the molecule. Only after annealing the sample the majority of the molecules are found in the lower configuration. It is also worth noting here that some recent evidence supports a strong rearrangement of the Si(111) surface upon the adsorption of C60 . In these experiments, after annealing one monolayer of C60 at 670 K, the (7$`\times `$7) diffraction pattern is lost and a (1$`\times `$1) pattern appears in its place. This transition is accompanied by a change to a more strongly bound adsorption state of the molecules.

The values of the charge transfer, from the Si(111)-(7$`\times `$7) substrate to the C60, deduced from different experimental techniques, oscillate between 3$`\pm `$1 electrons , and $``$0.21 electrons . The theoretical calculations might be important to understand this discrepancy. The charge transfers reported here are obtained by Mulliken analysis . Fig. 1(b) shows the evolution of the charge transferred to the molecule versus the distance to the rest atom. The charging of the molecule increases rapidly as the molecule approaches the surface. For the relaxed structures the charge transfers are: 1.00 electrons for the S3, 0.94 and 0.97 for the DB and SB molecules, and 1.13 for the L configuration. Most of the transferred charge is in the atoms directly bound to the surface: 0.16-0.17 electrons in the carbon atoms bound to the rest atom, 0.12-0.14 electrons in those bound to the adatoms, and 0.16 in those bound to the pedestal atoms in the L molecule. The rest of the charge, $``$0.5 electrons in all the cases, is mainly in the atoms closer to the surface, rather than uniformly distributed. This is an indication of the nature of the charge transfer. The charging of the molecule does not come through the occupation of the initially unoccupied states of the fullerene, but through the hybridization with the surface states. This result agrees with recent experimental reports of the transport through C60 on Si(111). .

A previous calculation with the DV-Xα-LCAO method , has reported a considerably larger charge transfer of 3.35 electrons. We believe that this discrepancy stems from the inherent arbitraryness of the Mulliken analysis to split the charge to covalently bound atoms, which makes the results strongly dependent on the basis set used. However, we stress that both calculations agree on the general pattern of the charge transfer.

In conclusion, the adsorption of C60 on Si(111)-(7$`\times `$7) has been studied by ab initio density functional calculations. Different adsorption configurations have been explored, finding good candidates for the two adsorption states experimentally observed. The adsorption energies range between 4 and 6.5 eV, while the charge transfer is always very close to one electron and mainly localized on the carbon atoms bound to the substrate. Work is in progress to simulate the STM images of the different adsorption structure described in this paper. We hope this will allow a detailed comparison with the experiments. Molecular dynamics simulations are also in progress to explore the mechanism, and energetic barriers for the transition between different adsorption configurations.

This work was supported by Grants No. DOE 8371494, and No. DEFG 02/96/ER 45439.


text:
# Spinodal Decomposition in Finite Temperature 𝑆𝑈(2) and 𝑆𝑈(3)

## 1 Introduction

Much of the experimental relevance of finite temperature QCD comes from non-equilibrium situations, as in the case of heavy ion collisions. Current lattice techniques are restricted to equilibrium simulations, and cannot address directly non-equilibrium time evolution. However, there are aspects of non-equilibrium behavior that are determined by the equilibrium effective potential. After a rapid change in temperature, induced by some external agency, it is the effective potential that determines the process of re-equilibration, in particular whether relaxation, nucleation or spinodal decomposition is the dominant equilibration mechanism. We have shown that a large, rapid increase of temperature across the deconfinement temperature results in spinodal decomposition in pure $`SU(2)`$ and $`SU(3)`$ gauge theories. In particular, we have found the existence of exponentially growing long wavelength modes in the approach to equilibrium, characteristic of spindodal decomposition. This behavior depends only on the features of the equilibrium effective action in an unstable region.

We consider first the generic case of a pure $`SU(N)`$ gauge theory in which the temperature is raised rapidly from a temperature less than $`T_d`$, the deconfinement temperature, to a temperature above $`T_d`$. We refer to such a rapid increase in system temperature as a quench, a borrowing of nomenclature from condensed matter physics. The equilibrium order parameter for the deconfinement transition is the Polyakov loop $`L_F`$. At temperatures $`T<T_d`$, the $`Z(N)`$ global symmetry associated with confinement is unbroken, and $`L_F=0`$. When the temperature is rapidly increased to $`T>T_d`$, $`L_F=0`$ is no longer the stable state of the system, and must evolve to a new equilibrium state with $`L_F0`$, reflecting the transition to the gluon plasma phase.

## 2 Effective Potential

The instability of $`L_F=0`$ at high temperatures is seen directly from the one-loop finite temperature effective potential for the Polyakov loop, as derived by Gross, Pisarski and Yaffe and Weiss. For simplicity, consider the case of $`SU(2)`$. It is convenient to parametrize the Polyakov loop as

$$L_F=\mathrm{cos}\left[\pi \left(1\psi \right)/2\right].$$
(1)
where $`1\psi 1`$. The effective potential at one loop for gauge bosons in a constant Polyakov loop background is

$$V(\psi )=\frac{\pi ^2T^4}{15}+\frac{\pi ^2T^4}{12}\left(1\psi ^2\right)^2.$$
(2)
The one loop result dominates the effective potential for $`T\mathrm{\Lambda }_{QCD}`$ due to asymptotic freedom. The first term is the standard black body result, obtained when $`\psi =1`$. The use of the variable $`\psi `$ makes the $`Z(2)`$ symmetry of the potential under $`\psi \psi `$manifest. Note that the equilibrium value of $`\psi `$ is $`\pm 1`$, corresponding to $`L_F=\pm 1`$; $`\psi =0`$, corresponding to $`L_F=0`$, is a maximum of $`V(\psi )`$.

Our picture of the quenching process is that the system is initially in a state where $`\psi `$is equal to zero at some temperature below $`T_d`$. When the system is quickly raised to a new temperature $`T>T_d`$, the system is still in the state with $`\psi =0`$. However, the system is unstable, and must eventually find its way to either $`\psi =+1`$or $`\psi =1`$. Because we quench into a region of the phase diagram where $`V^{\prime \prime }(\psi )<0`$, the system will decay to the equilibrium state via spinodal decomposition.

Because pure $`SU(2)`$ gauge theory has a second-order deconfining transition, spinodal decomposition will occur after quenching to any temperature $`T>T_d`$. The situation is more subtle for $`SU(3)`$, where the deconfinement phase transition is first order. This implies the existence of a metastable confined phase for some range of temperatures above $`T_d`$, which in turn implies that nucleation is the dominant equilibration mechanism just above $`T_d`$. Perturbation theory gives no hint of this behavior, and the one-loop effective potential is unstable at $`L_F=0`$. For temperatures sufficiently large that the one-loop effective potential for $`L`$ is reliable, spinodal decomposition will occur.

## 3 Langevin Model

In order to study the dynamics of this transition, we use the effective action of Bhattacharya combined with Langevin dynamics. For $`SU(2)`$, the effective action takes the form

$$S_{eff}\left[\psi \right]=d^3x\left[\frac{\pi ^2T}{2g^2}\left(\psi \right)^2+V(\psi )/T\right].$$
(3)
The equilibrium distribution will be $`\mathrm{exp}\left[S_{eff}\left[\psi \right]\right]`$. We postulate Langevin dynamics of the form

$$\frac{\psi (x,\tau )}{\tau }=\mathrm{\Gamma }\frac{\delta S_{eff}\left[\psi \right]}{\delta \psi (x,\tau )}+\eta (x,t)$$
(4)
where the white noise $`\eta `$ is normalized to

$$\eta (x,\tau )\eta (x^{},\tau ^{})=2\mathrm{\Gamma }\delta ^3(xx^{})\delta (\tau \tau ^{}).$$
(5)
The late-time relaxational behavior of $`\psi `$ is controlled by the Debye screening mass $`m_D`$, given by $`m_D^2=2g^2T^2/3`$. Near equilibrium, any effects of initial conditions decay exponentially as $`\mathrm{exp}\left[\frac{2\pi ^2\mathrm{\Gamma }T}{g^2}(k^2+m_D^2)\tau \right]`$. However, for early times, the initial conditions contribute to $`\stackrel{~}{\psi }(k,\tau )\stackrel{~}{\psi }(k,\tau )`$ a term of the form

$$\stackrel{~}{\psi }(k,0)\stackrel{~}{\psi }(k,0)\mathrm{exp}\left[\frac{2\pi ^2\mathrm{\Gamma }T}{g^2}(k^2k_c^2)\tau \right]$$
(6)
where $`k_c^2`$$`=g^2T^2/3=m_D^2/2`$. Modes with $`k<k_c`$ are initially not damped but grow exponentially, with the $`k=0`$ mode growing the fastest. This is a characteristic feature of spinodal decomposition with a non-conserved order parameter.

## 4 Monte Carlo Results

We have verified the existence of spinodal decomposition with simulations of rapidly quenched pure $`SU(2)`$ gauge theory. Lattices of size $`32^3\times 4`$ and $`64^3\times 4`$ were equilibrated at $`\beta =2.0`$, below the deconfinement transition at $`\beta _d=2.2986\pm 0.0006`$. The coupling constant was increased instantaneously to $`\beta =3.0`$, and the approach to equilibrium monitored via the Polyakov loop and other observables. The heat bath algorithm was used both for equilibration at low temperature before the quench and for subsequent dynamical evolution after the quench. While this time evolution is not the true time evolution of the non-equilibrium quantum field theory, features such as spinodal decomposition which depend only on the equilibrium effective action will occur with any local updating algorithm which converges to the equilibrium distribution. The abrupt change in $`\beta `$ is a potential cause of concern with this procedure, since the lattice spacing, and hence the physical volume, changes in all directions when $`\beta `$ is changed. However, the large spatial sizes used should mitigate this effect.

Figure 1 shows the Fourier transform of the connected Polyakov two point function $`S(k,\tau )`$for low values of the wave number as a function of Monte Carlo time for the same simulation. Note the early exponential rise in these modes, followed by a sharp disappearance as the Polyakov loop reaches its equilibrium value, characteristic of spinodal decomposition. Only the low momentum modes exhibit this growth; above $`k_c`$ no such growth occurs. Although the general behavior of $`S(k,\tau )`$ is the same for each run, many details are run dependent. In this particular run, the $`k/T=0.68`$ mode achieves a larger amplitude than $`k/T=0.56`$, which is atypical. In some runs, there is clear evidence for mode-mode coupling, reflecting the nonlinearity of the system. For each run, we have estimated the rate of growth of each low-momentum mode by fitting $`log(S(k,\tau ))`$ to a straight line in $`\tau `$ for early times. We can extract $`k_c`$ as the value where the growth rate is zero. From equation $`(6)`$, the growth rate of each line is proportional to $`k_c^2k^2`$ for the linearized theory, but this may not represent the true time evolution of the simulation. In any case, the growth rates measured in each run are highly dependent on initial conditions. In figure 2, we plot these growth rates versus $`k^2/T^2`$ for the same run used in figure 1. The error bars are naive estimates of the error for each growth rate for this particular run. The x intercept provides an estimate of $`k_c^2`$ for each run. Using multiple $`64^3\times 4`$ runs at $`\beta =3.0`$, we have estimated $`k_c/T`$ to be $`1.14\pm 0.02`$. The principle errors in this estimate come from the sensitivity of individual modes to initial conditions and the discrete character of $`k`$ on the lattice. We can compare this with lattice measurements of the Debye screening length, assuming the one-loop relation

$$\frac{k_c}{T}=\frac{m_D(T)}{\sqrt{2}T}$$
(7)
holds in general. Using the results of Heller et al. for $`m_D(T)`$, we obtain $`k_c/T=1.35(5)`$. We consider this to be reasonable agreement, given the many uncertainties involved.

It is natural to ask what happens when rapid cooling takes place, as might occur in the late stages of the expansion of a quark-gluon plasma, or in the early universe. At low temperatures, we expect that $`V(\psi )`$ has a single minimum and $`V^{\prime \prime }(\psi )>0`$ everywhere. Thus theory predicts the absence of unstable modes, and that relaxational processes should dominate the approach to equilibrium. Simulations of such a cooling process, in which $`32^3\times 4`$ lattice configurations equilibrated at $`\beta =3.0`$ are suddenly cooled to $`\beta =2.0`$, shows no sign of spinodal decomposition.

We have begun studies of the more physically relevant case of $`SU(3)`$. Since nucleation dominates at temperatures just above $`T_d`$, it is necessary to go to higher temperatures to observe spinodal decomposition. Working with $`32^3\times 4`$ lattice configurations equilibrated at $`\beta =5.5`$, we performed quenches to various values of $`\beta >\beta _d5.69`$. For values of $`\beta `$ close to $`\beta _d`$, we believe we see evidence for metastability. We have observed exponentially growing low wavelength modes for quenches to $`\beta 5.80`$.

In figure 3, we show the evolution of the three lowest non-zero modes as a function of simulation time for a quench from $`\beta =5.5`$ to $`\beta =5.92`$. Note the similarity to the behavior seen in $`SU(2)`$. A systematic investigation is underway, emphasizing $`k_c`$ as well as the limit of metastability.

## 5 Conclusions

Analytical and simulation results both indicate the relevance of spinodal decomposition in the equilibration of a gluon plasma after a rapid quench to high temperature. In the case of $`SU(N)`$ gauge theories, we have shown that the physical parameter $`k_c`$ which controls domain growth can be determined from simulations. In the case of $`SU(3)`$, determination of $`k_c`$ as a function of $`\beta `$ will allow us to define the limit of metastability as the point $`k_c(\beta _m)=0`$. Note that small values of $`k_c`$ require correspondingly large lattice sizes to see non-zero exponentially growing modes.

Just as the formation of a DCC may lead to enhanced production of low-momentum pions, spinodal decomposition could lead to enhanced production of low-momentum gluons in the early stages of plasma formation. The characteristic scale for such a phenomenon would be $`gT`$.

Full QCD has quarks as well as gluons, and the phase structure is different. It would be very interesting to study the behavior of both the Polyakov loop and the chiral condensate after abrupt changes in temperature. The Polyakov loop can be studied fairly easily if large equilibrated unquenched lattice field configurations are available. However, the chiral condensate requires a reliable estimator for the local condensate as a function of space and simulation time.


text:
# A Search for Jovian Planets around Hot White Dwarfs

## 1 Introduction

Searches for extrasolar planets have employed four techniques: direct imaging of the star and planet, photometric monitoring for transits, and astrometric and Doppler measurements of the reflex motion of the star due to planets. The imaging technique is difficult because a Jupiter-like planet is $``$109 and $``$104 times fainter than a solar-type host star at the visible and thermal IR wavelengths, respectively. The photometric method searches for occultations, but can only detect systems with the orbital plane seen edge-on, which is an exceedingly rare circumstance. For example, a Jupiter-sized planet orbiting a solar-type star at 1 AU has only $``$0.5% chance of being in the proper alignment for photometric detection. The reflex motion of a star due to planets is typically small. To date, only the Doppler measurements have been successful in discovering planets around a pulsar (Wolszczan & Frail, 1992) and solar-type stars (e.g., Mayor & Queloz, 1995; Marcy & Butler, 1996; Butler & Marcy, 1996; Cochran et al., 1997). With the exception of the planets around the pulsar PSR1257+12, all other known extrasolar planets are Jupiter-sized, 0.2–11 Jupiter masses (MJ), with orbital semi-major axis ranging from 0.04 AU to 3 AU.

Current searches for extrasolar planets do not explore whether intermediate-mass (2 – 7 M) stars also possess planetary systems, because these stars, of MK spectral types A and B, have fewer stellar absorption lines available for accurate velocity measurement. Furthermore, because they are more massive than solar-type stars, the stellar reflex motion is even smaller. These difficulties can be circumvented, if one searches for planets around the descendants of these intermediate-mass stars: hot white dwarfs with $`T_{eff}>50,000`$ K. A hot white dwarf’s UV radiation will photoionize the planetary atmosphere and recombination lines will be emitted by the planet. The contrast between the planet’s recombination line emission and the stellar spectrum will be the limiting factor for detecting the planet. The strength of the planet’s recombination line emission is proportional to the UV flux where the stellar emission peaks, but the recombination line is emitted at a longer wavelength where the stellar emission is lower; therefore, the contrast between the planet’s recombination line emission and the stellar spectrum is not as extreme as the contrast between the planet and star both in continuum.

Thus, it may be possible to find planets around hot white dwarfs by searching for a periodically Doppler-shifted recombination line emission superposed on the stellar spectrum. We have carried out calculations to demonstrate the feasibility of this search, taking into account the stellar evolutionary effects. In this letter we report the results of our calculations, discuss practical considerations, and present observations of a white dwarf to illustrate observational considerations.

## 2 Detectability of Planets around Hot White Dwarfs

Our proposed method to search for planets around hot white dwarfs utilizes recombination lines emitted by the ionized atmosphere of orbiting planets. In order to detect the planetary emission, the planet’s atmosphere needs to contain an abundant element that has recombination lines in easily observable wavelength ranges. Given currently available facilities, hydrogen provides the most observable recombination lines. Therefore, we will consider searching for only Jupiter-like gas giants around hot white dwarfs.

The progenitor of a hot white dwarf with $`T_{eff}>50,000`$ K has gone through various evolutionary stages, among which the asymptotic giant branch (AGB) phase is the stage at which the star has the most pronounced effects on its planetary system, because the star swells up to a few hundred R, and goes through episodes of copious mass loss (Iben & Renzini, 1983). At least four competing factors may affect the planetary system: (1) Planets engulfed in the stellar envelope will experience a drag force and migrate inward. (2) The decreasing central star mass will cause the planets to migrate outward. (3) Tidal interaction between a planet and the stellar envelope will slow down the planet’s orbital motion and cause it to migrate inward. (4) Planets orbiting an AGB star can experience accretion and/or evaporation. These factors have been taken into account by Livio & Soker (1984) in their calculation of the evolution of a star-planet system. They found that Jovian planets engulfed in the AGB envelope will either evaporate or accrete mass to become a close companion, depending on the mass of the stellar envelope, the initial mass of the planet, and the initial separation between the planet and the star. The calculations of Livio & Soker (1984) are admittedly qualitative, and there are no theoretical predictions which preclude the existence of Jovian planets within a few AU of a post-AGB star, or a hot white dwarf. Furthermore, a Jovian planet has been suggested to be a companion to the hot white dwarf PG 0950$`+`$139 (Dopita & Liebert, 1989).

Hot white dwarfs with $`T_{eff}>50,000`$ K are powerful sources for H-ionizing UV radiation. In this environment, the hydrogen-rich outer envelopes of any Jovian planets present will be photoionized. The resulting recombinations will cause the planet to emit hydrogen recombination lines. To calculate the recombination line strengths we have made a number of assumptions and approximations. For the stellar radiation, we use a blackbody radiation model to approximate the stellar spectrum and use 1 R as the stellar radius to determine the ionizing photon flux emitted by the star. We approximate the planetary atmosphere as a plane-parallel slab with physical dimension equivalent to the cross section of the planet. Within the slab, the UV flux from the white dwarf will ionize the outer layer of the atmosphere, to the point where the number of incident ionizing photons balances the number of recombinations. In other words, we treat the atmosphere as a “Strömgren slab.” Assuming the ionized atmosphere is at a temperature of 10,000 K, the expected H$`\alpha `$ luminosity is then given by:

$$L_{\mathrm{H}\alpha }=7.7\times 10^{20}Q(\frac{r_\mathrm{P}}{a})^2\mathrm{ergs}\mathrm{s}^1,$$
(1)
where $`Q`$ is the ionizing flux of the white dwarf in number of ionizing photons per second, $`r_\mathrm{P}`$ is the radius of the planet in Jupiter radii (RJ), and $`a`$ is the planet’s orbital semi-major axis in AU.

We have estimated the expected brightness contrast between a Jovian planet and the host star at the H$`\alpha `$ line using the model described above. The calculations were performed for planets with radii of 0.5, 1.0, 2.0, and 5.0 RJ and orbital semi-major axes of 0.5, 1.0, 2.0, and 5.0 AU over a range of white dwarf temperatures, 20,000–200,000 K. The ratio of the Jovian planet’s H$`\alpha `$ line flux to the underlying stellar flux (within a comparable bandwidth, $`0.5`$ Å, the expected thermal width of the H$`\alpha `$ line) is shown in Figure 1.

From this simplistic model, it can be seen in Figure 1 that the ratio between the planet’s H$`\alpha `$ line emission and the stellar continuum emission can be as high as $``$1/20 for a Jupiter-sized planet at 1 AU from a $``$200,000 K star. This model, however, ignores many effects which will lower the contrast between the planet and star. For instance, in a dense, molecule-rich planetary atmosphere, H+ can react with atoms or molecules rather than recombine with a free electron; thus, the planetary H$`\alpha `$ emission will be reduced. Furthermore, many hot white dwarfs show narrow H$`\alpha `$ and He II $`\lambda `$6560 line emission superposed on broad H$`\alpha `$ and He II absorption (Reid & Wegner, 1989) because of non-LTE (non-Local Thermal Equilibrium) effects in the stellar atmosphere (Lanz & Hubeny, 1995; Hubeny et al., 1999). This will complicate the search for planetary H$`\alpha `$ emission. Therefore, the results in Figure 1 portray an overly optimistic view of the detectability of Jovian planets around a hot white dwarf. On the other hand, the photoionization-driven ablation of a planetary atmosphere may produce an extended low-density ($`10^6`$ H-atom cm-3) nebula, and the forbidden lines emitted form this nebula can be easily detected (Dopita & Liebert, 1989).

## 3 Observational Considerations

To detect the H$`\alpha `$ emission from a Jovian planet around a hot white dwarf, multi-epoch, high-dispersion, high-S/N spectroscopic observations of a large sample of hot white dwarfs are needed. Below are some of the observational considerations for such a search.

* Velocity resolution: The orbital velocity of a planet within 5 AU of the host star will be $``$10 km s-1. To be able to detect this motion, the precision of the velocity measurements should be $``$1–2 km s-1. This can be easily achieved with a high-dispersion spectrograph with R = 30,000.
* Large telescope: A sensitive search for recombination line emission from a planet requires a high signal to noise ratio (S/N). To achieve a S/N ratio of 100 for a 13.5 mag star with a reasonable exposure time ($``$1.5 hours), the telescope must be at least $``$4 m in diameter.
* Long slit: Due to the possible faintness of the planetary recombination line emission, telluric and interstellar line emission can interfere with the detection of the planetary lines. To facilitate the removal of telluric and interstellar line emission, a long slit is required for the observations.
* Multi-epoch observations: Given the range of semi-major axes over which this method is likely to detect H$`\alpha `$ emission from a planet (a few AU), multiple epoch observations separated by 4–6 months are required to detect any change in the H$`\alpha `$ emission line due to the orbital motion of the planet. Planetary H$`\alpha `$ emission may be distinguished from stellar H$`\alpha `$+He II emission by the periodic velocity change of the planetary emission. Monitoring observations over a period of 2–3 years will be needed to determine the planetary orbital elements.
* Planetary mass measurements: Unlike reflex motion planet searches, the velocity shift observed is that of the planet rather than that of the host star. The mass of the planet cannot be accurately determined from this motion. Therefore, astrometric measurements of the reflex motions of the host white dwarfs are needed to determine the masses of discovered planets.

## 4 Example Observation of a Hot White Dwarf

The hot white dwarf Feige 34 was observed on 2000 April 23 (UT) with the echelle spectrograph on the 4-m telescope at Kitt Peak National Observatory. The high-dispersion spectra were obtained using the 79 l mm-1 echelle grating in combination with a 226 l mm-1 cross disperser and the long focus red camera. A reciprocal dispersion of 3.5 Å mm-1 at H$`\alpha `$ was achieved with this setup. The spectra were imaged with the T2KB CCD detector. The pixel size is 24 $`\mu `$m, corresponding to 0$`\stackrel{}{\mathrm{.}}`$24 pixel-1 along the slit and $``$3.7 km s-1 pixel-1 along the dispersion axis. Two 1800 s exposures were taken of Feige 34 with an east–west oriented 20″-long slit of width 1$`\stackrel{}{\mathrm{.}}`$5. The H$`\alpha `$ region of the echellogram is presented in Figure 2. All of the data were reduced using standard routines in IRAF111Image Analysis and Reduction Facility – IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation..

The echellogram of Feige 34 shows broad H$`\alpha `$+He II absorption from the white dwarf superposed with a narrow H$`\alpha `$ emission component and a faint He II emission component (marked in Figure 2). The H$`\alpha `$ emission component has a FWHM of $``$70 km s-1; it appears offset blueward from the center of the absorption because the absorption is a blend of a strong H$`\alpha `$ line and a weak He II line. Using the observations of Massey et al. (1988) for a flux calibration, we determine the H$`\alpha `$ flux of the emission component to be $`1.4\times 10^{14}`$ erg cm-2 s-1. This large flux and the lack of interstellar emission features on the Palomar Sky Survey plates rules out an interstellar origin for this emission component. Furthermore, the lack of forbidden line emission excludes the possibility that the H$`\alpha `$ emission arises from a typical circumstellar nebula. It is possible that this H$`\alpha `$ emission is dominated by the stellar H$`\alpha `$ emission due to a non-LTE atmosphere (Werner, 2000). This does not, however, exclude the possibility that the non-LTE stellar emission is superposed with weak H$`\alpha `$ emission from an orbiting planet or low-mass dwarf companion.

A low-mass red dwarf companion to Feige 34 has been suggested by Thejll, MacDonald, & Saffer (1991) based on an infrared excess in JHK bands. Alternatively, a large planet close to the star, heated by the UV flux, may also produce the observed infrared excess. To distinguish between these two possibilities, further theoretical and observational work is required. Detailed modeling of a non-LTE white dwarf atmosphere for Feige 34 is needed to quantitatively assess whether all or part of the H$`\alpha `$ emission feature in the spectrum of Feige 34 arises from the star. Follow-up spectroscopic observations are necessary to determine whether the H$`\alpha `$ emission line profile shows temporal variations suggesting the existence of a companion. If temporal variations are detected, further monitoring observations would allow the orbital parameters and nature of the companion to Feige 34 to be determined.

## 5 Summary

We suggest that high-dispersion spectroscopic observations of hot white dwarfs can be used to search for H$`\alpha `$ emission from Jovian planets within a few AU of the star. H$`\alpha `$ emission from the ionized outer envelope of a Jovian planet as weak as 1/100 the stellar emission can be detected with modern large ($``$4-m) telescopes. Non-LTE effects in the hot white dwarf atmosphere may produce stellar H$`\alpha `$ emission. Stellar and planetary H$`\alpha `$ emission can be disentangled by multi-epoch observations if changes in the H$`\alpha `$ emission line profile due to the orbital motion of the planet are detected. Monitoring observations to establish the orbital parameters are required to distinguish between Jovian planets and a low-mass dwarf companion. High-dispersion spectroscopic observations of Feige 34 demonstrate the complexity and feasibility of this method to search for these planets. Such a search will provide targets for future high-precision astrometric observations (e.g., SIM, adaptive optics, etc.), and valuable insight can be obtained on the development and evolution of planetary systems around stars more massive than the Sun.

We would like to thank Klaus Werner for enlightening discussions on non-LTE effects in a hot white dwarf atmosphere and sharing his model calculations of Feige 34’s H$`\alpha `$+He II lines. We also thank Ron Webbink, Charles Gammie, and Wing Ip for helpful conversations, the anonymous referee for insightful comments, and Caty Pilachowski for sharing unpublished high-dispersion spectra of white dwarfs.