text:
# Flux-loss of buoyant ropes interacting with convective flows

## 1 Introduction

The concept of buoyant magnetic flux tubes is an essential part of the framework of current theories of dynamo action in stars, particularly in the case of cool dwarf stars such as the Sun. Results from studies of buoyant magnetic flux tubes carried out within the essentially 1-d thin flux tube approximation (e.g. Spruit Spruit1981 (1981) and Moreno-Insertis Moreno86 (1986)) are consistent with the observed latitudes of emergence and tilt angles of bipolar regions on the surface of the Sun (Fan et al. Fan+ea94 (1994) and Caligari et al. Caligari+ea95 (1995)). More general 2-d simulations of flux tube cross-sections have shown that cylindrical tubes are disrupted by a magnetic Rayleigh-Taylor instability (e.g. Schüssler Schussler1979 (1979), Tsinganos Tsinganos1980 (1980), Cattaneo et al. Cattaneo+ea90 (1990), Matthews et al. Matthews+ea95 (1995), Moreno-Insertis & Emonet Moreno+ea96 (1996)). This renders them unlikely to reach the surface unless the presence of fieldline twist introduces a sufficient amount of magnetic tension to suppress this effect (Emonet & Moreno-Insertis Emonet+ea96 (1996, 1998) and Dorch & Nordlund Dorch+ea98 (1998)). On the one hand, 3-d simulations of buoyant, twisted flux ropes have confirmed several of the results from 2-d simulations (Matsumoto et al. Matsumoto+ea98 (1998) and Dorch et al. Dorch+ea99 (1999)), and have further shown that the S-shaped structure of a twisted flux tube as it emerges through the upper computational boundary is qualitatively similar to the sigmoidal structures observed in EUV and soft X-ray by the Yohkoh and SoHO satellites (e.g. Canfield et al. Canfield+ea99 (1999) and Sterling et al. Sterling+ea00 (2000)). Moreover, tightly packed $`\delta `$-spots may be interpreted as the emergence of highly twisted, kinking flux ropes (e.g. Fan et al. Fan+ea99 (1999)). On the other hand, it has been suggested that the value of the critical degree of twist needed to prevent the Rayleigh-Taylor instability may be unrealistically high in the 2-d case, and a smaller twist may be sufficient in the case of sinusoidal 3-d magnetic flux loops (Abbett et al. Abbett+ea00 (2000)).

In the solar convection zone, buoyant flux structures are constantly interacting with the surrounding convective downdrafts and updrafts, and the question remains whether the quasi-steady behavior that the flux ropes reach in the later phase of their rise in 2-d simulations (Emonet & Moreno-Insertis Emonet+ea98 (1998) and Dorch & Nordlund Dorch+ea98 (1998)) is stable towards perturbations from the surroundings, and whether the results found for 3-d flux ropes moving in a essentially 1-d static stratification are valid in the more realistic case.

In this paper, we present our first results regarding the behavior of buoyant, twisted flux ropes embedded in a fully dynamic model of solar-like convection.

## 2 Numerical model

The set-up of the model is twofold, consisting of a snapshot of a time-dependent, but statistically relaxed “local box” convection zone model (sandwiched between two stable layers), and of an idealized twisted magnetic flux rope. We solved the full resistive and compressible MHD-equations on a staggered mesh of 150 vertical $`\times `$ $`105^2`$ horizontal grid points, using the method by Galsgaard and others (e.g. Galsgaard & Nordlund Galsgaard+ea97 (1997), Nordlund et al. Nordlund+92 (1992)):

$`{\displaystyle \frac{\rho }{t}}`$ $`=`$ $`\rho 𝐮,`$ (1)
$`{\displaystyle \frac{\rho 𝐮}{t}}`$ $`=`$ $`(\rho \mathrm{𝐮𝐮}\tau )P+𝐅_{\mathrm{grav}}+𝐅_{\mathrm{Lorentz}},`$ (2)
$`{\displaystyle \frac{𝐁}{t}}`$ $`=`$ $`\times (𝐮\times 𝐁)+\times (\eta \times 𝐁),`$ (3)
$`{\displaystyle \frac{e}{t}}`$ $`=`$ $`(e𝐮)+P(𝐮)+Q_{\mathrm{rad}}+Q_{\mathrm{visc}}+Q_{\mathrm{Joule}}.`$ (4)
Here $`\rho `$, $`𝐮`$, $`P`$, $`𝐁`$, and $`e`$ represent the density, velocity, pressure, magnetic field, and internal energy respectively. $`\eta `$ and $`\tau `$ denote the magnetic diffusivity and the viscous stress tensor; the source terms $`Q_{\mathrm{visc}}`$, $`Q_{\mathrm{Joule}}`$, and $`Q_{\mathrm{rad}}`$ refer to the viscous, Joule, and diffusive heating. In the upper part of the domain, $`Q_{\mathrm{rad}}`$ includes an additional term that provides for a simple, isothermal cooling layer.

The computational method employs a finite difference staggered mesh with 6th order derivative operators, 5th order centering operators, and a 3rd order time-stepping routine. The diffusive terms are quenched in regions with smooth variations, to reduce the diffusion of well-resolved structures. Magnetic Reynolds numbers in non-smooth regions are of the order a few times $`10^2`$, but can be much higher in smooth regions.

The computational box is horizontally periodic with sides of 250 Mm and a height of 313 Mm (of which 166 Mm is the convection zone, that covers 2.5 orders of magnitude in pressure). The flows are turbulent through-out the convection zone, and the kinetic energy spectrum displays a power law at intermediate wavenumbers ($`k`$ 3 – 10). As it is typical for over-turning stratified convection, a cellular granulation pattern is generated on the surface of the convection zone (Fig. 1). The typical length scale of this pattern is about 50 Mm, somewhat larger than the canonical size of 32 Mm of solar super-granules (e.g. Leighton et al. Leighton+ea62 (1962)). The typical velocity of the granulation is 200 m/s in the narrow downdrafts at the surface and slightly less in the upwelling regions, close to what is found for solar super-granulation (e.g. Worden & Simon Worden+Simon76 (1976)).

We choose an initially isentropic flux rope with a buoyancy of $`1/\gamma `$ (with $`\gamma =5/3`$) lower than the case of temperature balance, where the buoyancy is $`1/\beta `$ (with $`\beta `$ being the classical plasma beta). This is computationally advantageous, since we avoid the costly process of perturbing the flux rope from a state of mechanical equilibrium.

The initial twist of the flux rope is given by

$$B_z=B_0e^{(r/\mathrm{R})^2}\mathrm{and}B_\varphi =\alpha (r/\mathrm{R})B_z,$$
(5)
where $`B_z`$ is the parallel and $`B_\varphi `$ the transverse component of the magnetic field with respect to the rope’s horizontal main axis. The coordinate system is chosen so that $`z`$ corresponds to the initial axial direction of the rope.

The wavelength $`\lambda `$ of the flux rope is equal to the horizontal size of the domain so that $`\lambda =3.2\mathrm{H}_{\mathrm{P0}}`$ at the initial position of the rope. Thus, the flux rope is not undular Parker-unstable even though the stratification permits this instability for longer wavelengths (Spruit & van Ballegooijen Spruit+ea82 (1982)). The rope is initially twisted, with a pitch angle (at $`r=\mathrm{R}`$) of $`\psi _\mathrm{R}=\mathrm{arctan}(\alpha )`$ and a radius $`\mathrm{R}_0=0.177\mathrm{H}_{\mathrm{P0}}`$ which corresponds to a half-width at half-maximum of $`B_z`$ (henceforth HWHM) of $`0.1\mathrm{H}_{\mathrm{P0}}`$.

To avoid problems associated with the large ratio of thermal to dynamic time scales, our convection model has a much higher luminosity than the Sun, and thus, all variables must be scaled to compare with solar values. The choice of the magnetic field strength is somewhat problematic in this regard. The ratio of kinetic to thermal energy density $`e_\mathrm{K}/e`$ is much larger in our model than in the Sun (though the convective flows remain subsonic with an average Mach number of 0.01). This requires a choice of $`\beta `$ that is smaller than its presumed value at the base of the solar convection zone so that the ratio of magnetic to kinetic energy density $`e_\mathrm{M}/e_\mathrm{K}`$ is the proper order of magnitude. However, a small $`\beta `$ is what is needed so that the time it takes to complete a simulation is not prohibitively long. We choose $`\beta =100`$, which yields a solar-like $`e_\mathrm{M}/e_\mathrm{K}`$ of 100.

## 3 Results

We have performed several fully convective 3-d simulations, as well as a number of 2-d convection-less simulations. The results of the latter agree with previous 2-d findings, and are used for reference in the following. We discuss results from a 3-d simulation with $`\psi _\mathrm{R}=45^{}`$ ($`\alpha =1`$). In that case, the degree of twist is small enough to prevent the onset of the kink instability (the linear growth rate vanishes for $`\alpha =1`$, e.g. Fan et al. Fan+ea99 (1999)), yet it is large enough to prevent the onset of the Rayleigh-Taylor instability. Thus, the rope retains its cohesion without distorting its shape by any of these two instabilities, and we can focus our attention on the effects of the convective flows on the rope.

Fig. 2 compares our 3-d simulation to a corresponding 2-d convection-less reference simulation, and a simple analytic flux tube. As the 3-d rope rises, convective flows perturb its motion, preventing it from entering a well-defined terminal rise phase with a constant rise speed, as in the 2-d reference simulation (see Fig. 2a). The rope remains straight and the maximum excursion of its axis, at the end of the simulation, is $`0.04\lambda `$, where we define the rope’s axis as the set of positions along the rope, where the magnetic field strength is maximum. With the chosen super-equipartition axial field strength, the main action of the large-scale convective flows is to push the rope both left and right of the central plane (Fig. 2b; see also the mpeg-movie at Dorch Dorch-mpeg (2001)), while the effect of the small-scale downdrafts (of the order of the rope’s radius) is to locally deform its equipartition boundary.

Initially the rope is located in a general updraft region. This explains why the rise speed of the rope is slightly greater than that of the 2-d reference simulation, which reaches a terminal speed of $`0.1v_{\mathrm{A0}}`$ ($`v_{\mathrm{A0}}`$ being the Alfvén speed at the initial position of the rope). Nevertheless, the two ropes arrive at the same final height at the end of the runs (though in the 3-d case, we note that the rope follows a longer path). Our 3-d rope also expands more quickly than the rope in the 2-d simulation, and its rate of expansion is closer to what is expected from an adiabatically expanding, non-stretching tube with constant flux ($`B/\rho =`$const.):

$$\mathrm{R}(x)=\mathrm{R}_0\left(1_a\frac{xx_0}{\mathrm{H}_{\mathrm{P0}}}\right)^{1/2_a}.$$
(6)
Fig. 2c shows the rope’s characteristic size $`\mathrm{R}_{hwhm}`$ defined by the average between the vertical and horizontal HWHM along its axis (the short period oscillations due to differential buoyancy, that are not well-resolved, have been filtered out by smoothing over two grid points). As the rope rises and expands, its magnetic field strength $`\mathrm{B}_c`$, here defined as the average axial field strength along the rope, decreases at a rate close to that determined by Eq. (6). At later times, the field strength of the 3-d and 2-d ropes decrease at nearly the same rate (Fig. 2d). The deviation can be attributed to the fact that, during its ascent, a significant amount of the magnetic flux within the 3-d rope is lost to its surroundings.

This is illustrated by Fig. 3 (left), which shows the total normalized magnetic flux within the HWHM-boundary $`\mathrm{\Phi }_i`$ as a function of time for both the 2-d and 3-d ropes. We note that as the 3-d simulation progresses, the total flux-loss from the computational domain is only 0.3$`\%`$. The flux content of the rope, however, decreases much more quickly.

Also shown in Fig. 3 (right) is the magnetic flux external $`\mathrm{\Phi }_e`$ to the rope both above and below its center $`\mathrm{\Phi }_u`$ and $`\mathrm{\Phi }_l`$ respectively. Since the sum $`\mathrm{\Phi }_e+\mathrm{\Phi }_i`$ is nearly conserved, as $`\mathrm{\Phi }_i`$ decreases, $`\mathrm{\Phi }_e=\mathrm{\Phi }_u+\mathrm{\Phi }_l`$ must increase by an equal amount. However, the distribution of the flux-loss is not symmetric: more flux is lost to the surroundings below the rope than above. The e-folding time of the increase of flux $`\mathrm{\Phi }_l`$ in the lower domain is $`20\tau _\mathrm{A}`$ (with $`\tau _\mathrm{A}=\mathrm{R}_0/v_{\mathrm{A0}}`$).

This asymmetry also exists in the 2-d reference simulation, even though the total flux-loss is much smaller in that case. The asymmetry is a result of two factors. First, as the rope rises, the total volume above it decreases, while the volume below it increases. Second, there is an anti-symmetry of the relative velocity across the rope. When the rope ascends, there is a tendency for flux to be advected towards the rope near its apex, and transported away from the rope in its wake. The more pronounced asymmetry in the 3-d case can be attributed to the pumping effect that transports the weak field downwards (Dorch & Nordlund Dorch+Nordlund01 (2001) and Tobias et al. Tobias+ea01 (2001)).

We have defined the flux rope as the magnetic structure that lies within the HWHM-boundary. This boundary is not, however, a contour that moves with the fluid in the classical sense: the flux within the latter kind of contour is naturally conserved (neglecting resistive effects) and equal to the total flux $`\mathrm{\Phi }_0=\mathrm{\Phi }_e+\mathrm{\Phi }_i`$. The HWHM-boundary is a convenient way of defining the flux rope and a characteristic size $`\mathrm{R}_{hwhm}`$, that behaves more or less as it is expected from the analytical expression Eq. (6). The evolution of the flux within the rope’s core $`\mathrm{\Phi }_i`$ is determined by

$$\dot{\mathrm{\Phi }}_i=_{\mathrm{boundary}}\mathrm{\Delta }𝐯\times 𝐁d𝐥,$$
(7)
where $`\mathrm{\Delta }𝐯`$ is the difference between the fluid velocity and the motion of the HWHM-boundary. The average “slip” $`\overline{\mathrm{\Delta }v}`$ of the rope’s boundary in the simulation is only a small fraction of the rise speed, and varies between the range of plus or minus a few times $`10^4`$ and $`10^3v_{\mathrm{A0}}`$.

Making the rather crude assumptions that the boundary only moves radially relative to the fluid and that the circumference of the boundary is circular (which it is not), Eq. (7) reduces to

$$\dot{\mathrm{\Phi }}_i=\pi \mathrm{R}_{hwhm}\mathrm{\Delta }v\mathrm{B}_c,$$
(8)
where $`\mathrm{B}_c`$ is the field strength at the center of the flux rope. Integrating Eq. (8) numerically with $`\mathrm{R}_{hwhm}`$ and $`\mathrm{B}_c`$ determined from the simulation (see Fig. 2), the result is a perhaps surprisingly good fit to the actual flux-loss, see Fig. 3 (left), if $`\mathrm{\Delta }v`$ is set to 3 10$`{}_{}{}^{4}v_{\mathrm{A0}}^{}`$ throughout the time span of the simulation except for a short interval of $`3\tau _\mathrm{A}`$ around $`t=30\tau _\mathrm{A}`$, where $`\mathrm{\Delta }v=10^3v_{\mathrm{A0}}`$, when the rope passes from one updraft to another (see the discussion below).

## 4 Discussion and conclusions

The 3-d numerical simulations show that the interaction of a buoyant twisted flux rope with stratified convection leads to a considerable loss of magnetic flux from the core of the flux rope (as defined by the rope’s HWHM-boundary).

The initial position of a flux rope in the convection zone is significant for the subsequent detailed history of its rise: with the present convective flows and the initial location of the flux rope, most of the rope starts out located inside or close to a convective updraft. Thus, the ascent of the rope is likely to be influenced by this fact, and we are therefore not able to draw any conclusions on the detailed path of its rise. However, in the course of the simulation, the flux rope rises 96 Mm, and loses about 25% of its original flux content. This, ceteris paribus, leads to an increase in the amount of toroidal flux that must be stored at the bottom of the convection zone during the course of the solar cycle.

In the Sun, toroidal flux ropes rise about 200 Mm through the convection zone before emerging as bipolar active regions. One may thus expect them to lose even more of their initial flux, which would then be pumped back down toward the bottom of the convection zone. We can quantify this subsequent flux-loss by assuming that Eq. (8) is valid through-out the rise, that the ropes expand according to the simple analytical estimate of Eq. (6), and that the ratio of the slip $`\mathrm{\Delta }v`$ to the rise speed remains constant. Given these assumptions, the flux-loss at a height of 200 Mm is 26 % of the initial flux, i.e. not much more than in our simulation. However, the relative slip may not remain constant throughout the rope’s rise. For example, $`\mathrm{\Delta }v`$ and thus $`\dot{\mathrm{\Phi }}_i`$ changes at the time around $`t=30\tau _\mathrm{A}`$, which corresponds to the time when the rope is at its maximum (rightward) excursion from a vertical ascent (see Fig. 2b). At that time, the rope exits the convective updraft with which it was initially associated, and enters a different ascending “plume” to the left of the its original position. This leads to a transient compression of the rope ($`\mathrm{\Delta }v<0`$, in the simplified expression Eq. 8). After entering the new plume, the average slip returns to its previous positive value for the remainder of the rise.

Petrovay & Moreno-Insertis (Petrovay+Moreno97 (1997)) suggested that turbulent erosion of magnetic flux tubes may take place within the solar convection zone due to the “gnawing” of turbulent convection. They propose a mechanism whereby a flux tube is eroded by a thin current sheet that forms spontaneously within a diffusion time. That we do not see a loss of flux via this type of enhanced diffusion should not, however, be taken as a dis-proof of the feasibility of turbulent erosion: it requires the turbulence to be resolved down to much smaller scales $`\mathrm{}\lambda `$, than in our simulations. Instead, the flux-loss is completely due to the advection of flux away from the core of the flux rope by convective motions. Most of the flux that is “gnawed-off” ends up in the trailing wake and some of this flux is mixed back into the upper layers by ascending flows. We speculate that both types of flux-loss may take place simultaneously in the Sun, and as a result, the amount of toroidal flux stored near the bottom of the solar convection zone may currently be underestimated.

###### Acknowledgements.

SBFD and BVG was supported through an EC-TMR grant to the European Solar Magnetometry Network. WPA was supported through the NSF and NASA’s SECT and Solar Physics Research and Analysis programs. Computing time was provided by the Swedish National Allocations Committee. The authors thank Kristóf Petrovay for discussions on flux-loss mechanisms.


format: text

text:
# References

The wave packet propagation using wavelets

Andrei G. BORISOV a and Sergei V. SHABANOV b

aLaboratoire des Collisions Atomiques et Moléculaires,
Unité mixte de recherche CNRS-Université Paris-Sud UMR 8625,
Bâtiment 351, Université Paris-Sud, 91405 Orsay CEDEX, France
b Institute for Fundamental Theory, Departments of Physics
and Mathematics, University of Florida, Gainesville, FL 23611, USA

## Abstract

It is demonstrated that the wavelets can be used to considerably speed up simulations of the wave packet propagation in multiscale systems. Extremely high efficiency is obtained in the representation of both bound and continuum states. The new method is compared with the fast Fourier algorithm. Depending on ratios of typical scales of a quantum system in question, the wavelet method appears to be faster by a few orders of magnitude.

1. Owing to the fast development of the computational tools the direct solution of the time-dependent Schroedinger equation has become one of the basic approaches to study the evolution of quantum systems. Thus, the wave packet propagation (WPP) method is successfully applied to time dependent and time-independent problems in gas-surface interactions, molecular and atomic physics, and quantum chemistry . One of the main issues of the numerical approaches to the time-dependent Schroedinger equation is the representation of the wave function $`|\mathrm{\Psi }(t)`$ of the system. Development of the pseudospectral global grid representation approaches yielded a very efficient way to tackle this problem. The discrete variable representation (DVR)  and the Fourier grid Hamiltonian method (FGH)  have been widely used in time-dependent molecular dynamics , S-matrix  or eigenvalue  calculations. The FGH method based on the Fast Fourier Transform (FFT) algorithm is especially popular. This is because for the mesh of $`N`$ points the evaluation of the kinetic energy operator requires only $`NlogN`$ operations and can be easily implemented numerically. In the standard FGH method the wave function is represented on a grid of equidistant point in coordinate and momentum space. It is well suited when the the local de Broglie wavelength does not change much over the physical region spanned by the wave function of the system . There are, however, lot of examples where the system under consideration spans the physical regions with very different properties. Consider, for example, the scattering of a slow particle on a narrow and deep potential well. Despite the short wave lengths occur only in the well, in the FGH method they would determine the lattice mesh over entire physical space leading to high computational cost. In fact, pseudospectral global grid representation approaches are difficult to use in multiscale problems.

This is why much of work has been devoted recently to the development of the mapping procedures in order to enhance sampling efficiency in the regions of the rapid variation of the wave function . Though mapping procedure, based on the variable change $`x=f(\xi )`$ is very efficient in 1D, it is far from being universal. One obvious drawback is that it is difficult to implement in higher dimensions. In the case of several variables, there is no simple procedure to define the mapping functions $`f`$ so that the lattice would be fine only in some designated regions of space. Next, the topology of the new coordinate surfaces can be different from that of Cartesian planes. Therefore the Jacobian may vanish at some points (e.g., spherical coordinates). This leads to singularity in the kinetic energy and imposes small time step in simulations .

In this letter we propose a novel approach to the wave packet propagation in multiscale systems. It is based on the use of wavelets as basis functions in the projected Hilbert space. In contrast to the plane waves, the wavelets are localized in both real and momentum spaces. This characteristic property of wavelets allows us to accurately describe short wave length components of the wave function in designated spatial regions, while keeping only few basis elements in the bulk. This leads to a drastic reduction of the projected Hilbert space dimension without any loss of accuracy. In some cases it may even become possible to diagonalize the Hamiltonian matrix so that the time evolution can be elementary followed in the basis of eigenstates. We illustrate our approach by simulations of a simple multiscale one-dimensional scattering problem. A generalization to higher dimensions is straightforward and based on the standard mathematical construction of multidimensional wavelets . Wavelet bases have been successfully used for the systematic treatment of the multiple length scales in the electronic structure of matter  (for a comprehensive review see  and references therein). Here we apply, for the first time, wavelets to study the time evolution in multiscale quantum systems and demonstrate the advantages of the wavelet method over the conventional FGH method.

2. To illustrate the efficiency of the wavelet method in the wave packet propagation, we consider the 1D scattering of an electron on a narrow and deep potential well. Despite we use an electron as a projectile, procedure described below readily applies to molecular or atomic wave packets. The potential has the form $`\widehat{V}=`$ $`64\mathrm{exp}[(100x)^2]`$ (atomic units are used throughout the paper) and supports single bound state with the energy $`0.63200`$. We are interested in the transmission and reflection coefficients for the energies of the scattered electron within $`0.21.1`$ range. This corresponds to a typical wave length $`\lambda =2\pi `$ which is much larger than the potential well width. To calculate the scattering properties of the system we use the Gaussan wave packet impinging on a potential well. The calculation requires a typical size of the box $`x[60,60]`$, where 18 a.u. from each side are allocated for an optical potential $`\widehat{V}_{opt}`$ that absorbs reflected and transmitted waves .

We first introduce a uniform lattice and simulate the wave packet propagation by means of the FGH method. Split-operator technique  is used to calculate the action of the evolution operator $`(\widehat{W}=\widehat{V}+\widehat{V}_{opt})`$:

$`|\mathrm{\Psi }(t+\mathrm{\Delta }t)`$ $`=`$ $`\mathrm{exp}(i\mathrm{\Delta }t\left[\widehat{H}+\widehat{V}_{opt}\right])|\mathrm{\Psi }(t)`$ (1)
$`=`$ $`\mathrm{exp}(i\mathrm{\Delta }t\widehat{W}/2)\mathrm{exp}(i\mathrm{\Delta }t\widehat{T})\mathrm{exp}(i\mathrm{\Delta }t\widehat{W}/2)|\mathrm{\Psi }(t)+O(\mathrm{\Delta }t^3)`$
Numerical convergence is obtained with $`N=2^{14}`$ points at the grid and a time step $`\mathrm{\Delta }t=0.0002`$. Calculated transmission and reflection coefficients are shown in Fig. 1. The time evolution of the wave packet is presented in Fig.2. The colors represent the magnitude of the wave packet in the logarithmic scale: $`10\mathrm{ln}|\mathrm{\Psi }(t,x)|0.3`$. As the color changes from red to violet, $`\mathrm{ln}|\mathrm{\Psi }(t,x)|`$ decreases from $`0.3`$ to $`10`$. Without the potential present, we would see a colored ray (a trajectory of a free particle is a straight line) spreading for larger values of $`t`$.

To compare the above results with our new wavelet method, we take the Daubechies wavelets $`𝒟_{10}`$ which have ten first vanishing moments and the filter of length 20 . The basis is generated by the scaling function $`\varphi (x)`$ whose support lies in the interval $`x[0,19]`$. The lowest resolution level is set so that the initial wave packet is reproduced with high accuracy in the orthonormal basis of the functions $`\varphi _{1,j}(x)=\sqrt{2}\varphi (2xj)`$ where $`j`$ runs over integers inside the interval $`(120,120]`$. The wavelets $`\psi _{k,j}(x)=\sqrt{2^k}\psi (2^kxj)`$ are used to describe short wave length components of the wave function in the vicinity of the potential well. For a moment, they should be regarded as a special orthonormal basis with compact support and the following properties. If the basis of scaling functions $`\varphi _{1,j}(x)`$ spans well functions whose Fourier transforms are concentrated in a wave length band around $`\lambda _1`$, the corresponding wavelets from the $`k`$th resolution level form a good basis for a band centered at $`2^k\lambda _1`$. In our case we have used $`k=1,2,3,4`$. The number of wavelets needed on each resolution level is determined by the potential width (by its shape, in general). The technical details are explained in Section 3. It appears necessary to take 20 wavelets on each resolution level. Thus, all together we use only $`240+420=320`$ coefficients (or, equally basis functions) to describe the wave packet propagation instead of $`2^{14}`$ in the fast Fourier method. Due to such a tremendous reduction of the size of the problem, the Hamiltonian matrix can be directly diagonalized and the time dependent wave function can be easily obtained as:

$$|\mathrm{\Psi }(t+\mathrm{\Delta }t)=e^{i\widehat{V}_{opt}\mathrm{\Delta }t}_Ee^{iE\mathrm{\Delta }t}|EE|\mathrm{\Psi }(t),$$
(2)
where $`|E`$ and $`E`$ stands for the eigenvector and eigenvalue of the Hamiltonian, respectively. Convergent results are obtained with time step $`dt=0.025`$, i.e., much larger than in the fast Fourier method above. The results are given on Figs. 1 and 3. There is a perfect coincidence of the fast Fourier and our wavelet results, while the wavelet method is nearly 100 times faster.

For the sake of comparison, we also used Eq. (1) to simulate the time evolution. The same accuracy is achieved for a time step comparable to the time step in FGH method. However there is an alternative, much better splitting scheme thanks to the localization properties of wavelets. Consider the decomposition $`\widehat{H}=\widehat{H}_1+\widehat{H}_2`$ where $`\widehat{H}_1`$ contains matrix elements of the basis elements localized in the vicinity of the potential well, while $`\widehat{H}_2`$ corresponds to the free space and, therefore, contains matrix elements only of the scaling functions. The propagation can then be done accordingly:

$$|\mathrm{\Psi }(t+\mathrm{\Delta }t)=e^{i\widehat{V}_{opt}\mathrm{\Delta }t}e^{i\widehat{H}_2\mathrm{\Delta }t/2}\left\{_{E_1}e^{iE_1\mathrm{\Delta }t}|E_1E_1|\right\}e^{i\widehat{H}_2\mathrm{\Delta }t/2}|\mathrm{\Psi }(t),$$
(3)
where $`|E_1`$ and $`E_1`$ stands for the eigenvector and eigenvalue of the $`\widehat{H}_1`$. In the spit method (1) the error depends on the spectral range of operators involved. For larger spectral ranges, the errors get larger. The conventional way to cope with this problem is to reduce the time step. The approach (3), which we call the “wavelet tower diagonalization”, offers a better alternative. The convergence here is drastically improved because (a) $`[\widehat{H}_1,\widehat{H}_2]`$ is small and (b) the operator $`\widehat{H}_1`$ with a large spectral range can be diagonalized. Both the properties are hardly achievable without the wavelet basis. In our example the convergence is reached for a time step $`dt=0.01`$ (vs $`dt=0.0002`$ in (1)). Corresponding results are presented in Fig. 1. This approach applies even better to evolving heavy-particle wave packets because the Chebyshev  or Lanczos  schemes can be used. These schemes are, first, more efficient than the split method and, second, allow one to take full advantage of the sparse structure of the Hamiltonian in the wavelet basis (see below).

Note also that by making the potential depth greater so that the minimal wave length becomes twice shorter, we would have to increase the lattice size by factor two in the Fourier method, thus increasing the number of coefficients needed to describe the wave packet by $`2^{14}`$ (!), while in our wavelet approach one more resolution level is to be added, implying only 20 extra wavelet coefficients in the wave packet decomposition. Thus, the wavelet approach offers a systematic and easy way to improve the accuracy of simulations.

3. The description of technical details of our approach is limited to essential practical steps necessary to reproduce our results and to apply the algorithm to new systems. A general theoretical analysis of wavelets bases in multidimensional spaces can be found in . So we only discuss the Hilbert space $`L_2`$ of square integrable functions of a real variable $`x`$, $`𝑑x|\mathrm{\Psi }(x)|^2<\mathrm{}`$.

(i). The scaling function $`\varphi (x)`$ is defined by the equation $`\varphi (x)=2_lh_l\varphi (2xl)`$ (called the scaling relation), where the coefficients $`h_l`$ are called a filter. For a finite filter, the scaling function has compact support. Define $`\varphi _{n,j}(x)=2^{n/2}\varphi (2^nxj)`$ for all integers $`j`$ and $`n`$. The filter $`h_l`$ satisfies an equation obtained by combining the scaling relation with the required orthonormality condition of the scaling functions $`\varphi _{n,j}`$ for fixed $`n`$. Given a filter, numerical values of $`\varphi (x)`$ can be generated by an iteration procedure .

(ii). The subspace of $`L_2`$ spanned by $`\varphi _{n,j}`$ is denoted $`V_n`$. An important property $`V_n`$ is that $`V_{n1}V_n`$ and the projection of any function $`\mathrm{\Psi }`$ from $`L_2`$ onto $`V_n`$ converges to $`\mathrm{\Psi }`$ in the $`L_2`$ sense as $`n\mathrm{}`$. Consider an orthogonal decomposition: $`V_{n+1}=V_nW_n`$. There exists a function $`\psi W_0`$ called the mother wavelet such that $`\psi _{n,j}(x)=2^{n/2}\psi (2^nxj)`$ form an orthonormal basis in $`W_n`$. The numerical values of $`\psi `$ can be generated from the scaling relation $`\psi (x)=2_l(1)^lh_{1l}\varphi (2xl)`$. A finite dimensional subspace of $`V_n`$ is used to approximate $`L_2`$ in simulations. By construction, the functions $`\varphi _{l,j}`$, $`l<n`$, and $`\psi _{k,j}`$ for $`k=l,l+1,\mathrm{},n1`$ form an orthonormal basis in $`V_n`$: $`\psi _{k,j}|\psi _{k^{},i}=\delta _{kk^{}}\delta _{ij}`$, $`\varphi _{l,j}|\varphi _{l,i}=\delta _{ij}`$, and $`\psi _{k,j}|\varphi _{l,i}=0`$.

(iii). From the definition of $`\varphi _{l,j}`$ and $`\psi _{k,j}`$ it is clear that the index $`j`$ indicates the position of support of the scaling function or wavelet. The index $`k`$ of $`\psi _{k,j}`$ can be understood as follows. Let the Fourier transform of the mother wavelet $`\psi `$ be peaked at a momentum $`p_0`$. Then the Fourier transform of $`\psi _{k,j}`$ is peaked at the momentum $`p_k=2^kp_0`$. Since the wavelets with different $`k`$ are orthogonal, the coefficients $`d_{k,j}=\psi _{k,j}|\mathrm{\Psi }`$ determine relative amplitudes of successively shorter wave length components of $`\mathrm{\Psi }`$ in the vicinity of $`x=2^kj`$ as $`k`$ increases. For this reason, the index $`j`$ is called a position index, and $`k`$ is called a resolution level.

(iv). From the physical properties of a system one can estimate a wave length band, $`\lambda [\lambda _{min},\lambda _{max}]`$, required for simulations. The lowest resolution level $`l`$ is identified by the condition that $`\varphi _{l,j}`$ span functions with Fourier components in the vicinity of $`\lambda _{max}`$. The necessary maximal resolution level is determined by $`n\mathrm{log}_2(\lambda _{max}/\lambda _{min})`$. In our example $`l=1`$ and $`n=5`$. If $`N`$ scaling functions are needed (to cover a given physical volume), then, in general, on every wavelet resolution level there will be $`J_k=2^{kl}N`$ basis functions. The wave packet is decomposed in the corresponding basis of $`V_n`$

$$\mathrm{\Psi }(t)=\underset{j=1}{\overset{N}{}}s_{1,j}(t)\varphi _{1,j}+\underset{k=1}{\overset{n1}{}}\underset{j=1}{\overset{J_k}{}}d_{k,j}(t)\psi _{k,j}.$$
(4)
The decomposition coefficients as functions of time are to be found from the Schroedinger equation. Yet, the total number of coefficients in $`\mathrm{\Psi }(t)`$ is about the same as that in the uniform finite lattice approach.

The advantage of using wavelet bases becomes significant whenever the volume of regions where short wave lengths can appear during the time evolution is much smaller than the total physical volume of the system. This allows one to substantially reduce $`J_k`$ in Eq.(4) by taking higher resolution level wavelets only where they are needed. In our case we need only 20 wavelets at each resolution level. We shall refer to all the wavelets needed in the vicinity of one local minima of the potential as a wavelet tower. Higher tower floors correspond to higher resolution levels $`k`$. Each wavelet is thought of as a building block of width equal $`2^k`$.

(v). An important part of the algorithm is the projection of the Hamiltonian onto the wavelet towers. Matrix elements of the potential as well as the initial coefficients $`s_{l,j}=\varphi _{l,j}|\mathrm{\Psi }`$, $`d_{k,j}=\psi _{k,j}|\mathrm{\Psi }`$ are computed via Riemann sums. The matrix elements of the kinetic energy are computed using the the fast Fourier transform to obtain the second derivative of the function. The Hamiltonian matrix is sparse because of compact support of the basis functions. For instance, in our example $`\varphi _{n,j}|H|\varphi _{n,i}=0`$, if $`|ij|20`$. Finally, the time evolution can be computed by standard techniques such as split, Lanczos, or Chebychev methods . If after the projection on wavelet towers, the Hamiltonian matrix is small enough for a direct diagonalization, Eq. (2) can be used for the propagation.

In conclusion, we have demonstrated that wavelet bases can be extremely efficient for solving the time-dependent Schroedinger equation for multiscale systems. The bound and continuum states are accurately described with a much less number of basis functions as would be required in the standard Fourier-grid methods. In our example we were able to reduce the basis size by a factor of 50, which allowed us to scale down the computation time by a factor of 100. The method can easily be implemented in higher dimensions where a wavelet basis is built by taking the direct product of one-dimensional wavelets . The wavelet approach becomes more efficient as ratios of typical scales get larger. Finally we would like to stress that the wavelet method is ideally suitable for simulating wave packet propagation in multiscale problems with complicated form of the potential. This is because the wavelet towers can be custom designed for any topology of the potential minima (e.g., potential valleys) regardless of the dimensionality of the system.

Acknowledgments. S.V.S. thanks the LCAM of the University of Paris-Sud for warm hospitality. We are also grateful to Victor Sidis and Jean Pierre Gauyacq for support of this project, and to John Klauder for reading the manuscript and useful comments.

Figure captions

Fig. 1. Calculated transmission (black) and reflection (red) coefficients. Circles: The standard approach based on the split propagation and Fourier grid with uniform mesh. Solid curves: The wavelet approach with Hamiltonian matrix diagonaliztion. Triangles: The split propagation with the “wavelet tower diagonalization” as in Eq. (3).


text:
# On a theorem of Goussarov

## 1 Introduction

The heart of this paper is a simple proof of a theorem of Goussarov, using an elementary tool developed in my thesis. His theorem relates two possible notions of finite type invariants, the standard one which he co-invented with Vassiliev, and a more subtle one based on what he calls “interdependent” modifications of links. He proves that these invariants coincide up to a degree shift by a factor of two.

The reason I wrote this paper is Corollary 4, that a grope cobordism preserves Vassiliev invariants up to roughly half the degree. This corollary is needed to prove Theorem 2 of \[CT\], which states that grope cobordism of class $`c`$ is precisely the same as surgery on *simple* claspers (see \[H\] or \[CT\]) with *grope degree* $`c`$ or more. It is much easier to see that grope cobordism coincides with *rooted* clasper surgeries, for which Goussarov’s interdependent modifications are well-suited.

In my thesis I considered a problem similar to Corollary 4, the difference being that I considered knots bounding gropes, which is a much stronger notion than cobounding a grope with the unknot. In that case, I obtained a sharp answer, namely a knot bounding a grope of class $`k`$ is $`k/2`$-trivial. The current paper implies $`(k1)/2`$, which is not very far away. Achieving the sharp result is a lot of work for not much gain.

Goussarov’s main idea was to define filtrations using alternating sums over “disjoint variations.” For instance, one may take the variations to be crossing changes of a knot, as he did in \[G1\], defining the usual finite type filtration independently of Vassiliev. In \[G2\] he considered a more subtle notion of variation, an “interdependent modification,” which involves replacing an arc of a circle with an arbitrary other arc in the link complement with the same endpoints. Evidently he was trying to find more subtle invariants than those which are of standard finite type. But his discovery, which is the main subject of this paper, was that they are the same up to a degree change. Most recently, he invented and developed a notion of 3-manifold finite type invariants, where the variation is a surgery along a “Y.” (Equivalently, a Matveev move on a genus 3 handlebody.) This notion appears in \[GGP\] for example. In \[GR\] Garoufalidis and Rozansky use a variation of this notion to study pairs (Homology 3-sphere, Knot), where the leaves of the “Y” must link the knot trivially. If it weren’t for his untimely death, Goussarov surely would have produced more exotic and beautiful notions.

Acknowledgments: It is a pleasure to thank Thang Le and Peter Teichner for useful discussions. I also wish to dedicate this to the late Mikhail Goussarov, whose paper \[G1\], gave me an inspiring start in this field.

## 2 Two filtrations

Let $`X`$ be the set of isotopy classes of links or string links with a fixed number of components. The Vassiliev filtration of $`[X]`$ is a descending filtration $`[X]=_0^v_1^v_2^v\mathrm{}`$. Each piece $`_k^v`$ is spanned by alternating sums of the following form. Let $`xX`$ be a knot, link, or string link. Choose $`k+1`$ framed arcs from $`x`$ to itself, which guide homotopies of $`x`$ supported in neighborhoods of these arcs, which take the little piece of $`x`$ at one end of the arc and push this across the other end. Such a move is called a finger move. (Such finger moves can always be realized as crossing changes in some projection of $`x`$.) Let the set of these finger moves be called $`S`$. Then define $`[x;S]=_{\sigma S}(1)^{|\sigma |}x_\sigma `$ where $`x_\sigma `$ is $`x`$ modified by the finger moves in $`\sigma `$. Now, by definition $`_k^v`$ is additively generated by elements of the form $`[x,S]`$, where $`S`$ has $`k+1`$ finger moves in it.

One can replace single finger moves by groups222I use the term group to indicate that a collection of moves is to be viewed as a single variation. of finger moves in the above definition, and this is well-known to give the same filtration. See \[G1\] for a proof.

An alternative filtration can be defined in the same way using moves more general than finger moves. One way to think of these moves is that they are guided by circles which are attached to $`x`$ along a subarc. The move is to replace the part of $`x`$ running across the circle with the other arc of the circle. $`k+1`$ disjoint circles lead to generators of this alternative filtration, which I denote by $`_k^{alt}`$. Goussarov calls these interdependent moves. The reason is that such moves cannot neccesarily be realized by independent, i.e. disjoint, groups of finger moves.

As before, one can replace single moves along circles by groups of moves along circles and the same filtration is achieved, the proof being modified *mutatis mutandis*.

The main result of \[G2\], which is stated there in the dual setting, is that $`_k^v=_{2k}^{alt}=_{2k+1}^{alt}`$. It is proved in two steps.

###### Theorem 1

$`_k^v_{2k+1}^{alt}`$.

###### Theorem 2

$`_{2k}^{alt}_k^v`$.

Theorem 1 is the easier of the two. I provide a repackaged version of Goussarov’s proof for completeness. Theorem 2 is harder and this is the one for which I provide a new simplified proof.

## 3 Proofs of Theorems 1 and 2

*\[Proof of Theorem 1\]*

Each finger move has two interdependent moves associated to it as follows. A symmetric way to view a finger move along an arc $`\eta `$ is to push the strands of $`x`$ at both endpoints of $`\eta `$ along $`\eta `$ so that they crash through each other at some point in the middle of $`\eta `$. For each of the two modified strands of $`x`$, there is a disk cobounding that strand with the strand before the finger move, the disk having been swept out by the isotopy. These two disks intersect in a single clasp. The boundaries of these two disks are the circles guiding the two interdependent moves associated to the finger move. Now I claim the alternating sum over a set of finger moves is, up to a sign, the same as the alternating sum over the associated interdependent moves. This is essentially because if one only does one of a given pair associated to a finger move, this is just an isotopy. Let $`S=\{s_i\}`$ be the set of $`k+1`$ finger moves, and $`T=\{a_i,b_i\}`$ the associated set of $`2k+2`$ interdependent moves. Inductively I show that

$$\underset{\sigma T}{}(1)^{|\sigma |}x_\sigma =(1)^k\underset{\mu S}{}(1)^{|\mu |}x_\mu ,$$

which will be sufficient to prove the theorem. This left hand sum breaks up into

$$\underset{\tau \{a_2,\mathrm{},b_{k+1}\}}{}(1)^{|\tau |}(x_\tau x_{\tau \{a_1\}}x_{\tau \{b_1\}}+x_{\tau \{a_1,b_1\}}),$$

and since the first three terms are equal, we get

$$\underset{\tau \{a_2,\mathrm{},b_{k+1}\}}{}(1)^{|\tau |}x_\tau +\underset{\tau \{a_2,\mathrm{},b_{k+1}\}}{}(1)^{|\tau |}(x_{s_1})_\tau ,$$

and by induction this is just

$$(1)^k\underset{\mu \{s_2,\mathrm{},s_{k+1}\}}{}((1)^{|\mu |}x_\mu (1)^{|\mu |}x_{\{s_1\}\mu })$$

which equals $`(1)^k_{\mu S}(1)^{|\mu |}x_\mu `$ as desired. $`\mathrm{}`$

*\[Proof of Theorem 2\]*

Given a set of $`2k+1`$ interdependent moves, $`S`$ on $`x`$, we wish to show that $`[x;S]_k^v`$. The strategy is to show that $`[x;S]`$ is congruent modulo $`_k^v`$ to sums of simpler alternating sums $`[x;S^{}]`$, where the complexity of such alternating sums is recorded by a graph I will define in a moment. After iteration, we reduce the problem to alternating sums for which a direct argument is possible to show that they are congruent to $`0`$.

Fix a projection of $`x`$ and the circles in $`S`$. (That is, think of it as a planar picture with over and under crossing data.) The projection is assumed to be such that $`x`$ does not cross itself where circles are attached. In order to keep track of the complexity of the interdependence of the moves, I now define a graph associated to the picture. It has $`2k+1`$ vertices corresponding to the $`2k+1`$ circles that guide the interdependent moves. Fix an ordering on these vertices, say they are called $`v_1,\mathrm{},v_{2k+1}`$. Draw an edge between $`v_i`$ and $`v_j`$ if $`i<j`$, but the circle corresponding to $`v_i`$ crosses over the one corresponding to $`v_j`$. Draw an edge from a vertex to itself if the corresponding circle is knotted. Finally, label a vertex with a star if $`x`$ crosses over the corresponding circle.

For each edge of the graph, one can do an obvious group of crossing changes of $`x`$ union the circles, which eliminates that edge. That is, to eliminate an edge between $`v_i`$ and $`v_j`$, where $`i<j`$, do those crossing changes which make $`v_i`$ always pass under $`v_j`$ at each crossing. Similarly, to eliminate a star, there is a group of crossing changes that always make $`x`$ pass under a given circle. Finally there is a group of crossing changes eliminating a self-edge, i.e. by unknotting the corresponding circle.

Suppose a graph has at least $`k+1`$ edges plus stars. Let $`T`$ denote a set of $`k+1`$ groups of crossing changes on $`x`$ union the circles, each of which removes an edge or star from the graph. Each such group of crossing changes induces a group of crossing changes on each summand $`x_\sigma `$ of $`[x;S]`$. Thus each $`x_\sigma `$ is congruent modulo $`_k^v`$ to $`_{\mathrm{}\tau T}\pm x_{\sigma }^{}{}_{\tau }{}^{}`$. Then

$$\begin{array}{cc}\hfill [x;S]& =\underset{\sigma S}{}(1)^{|\sigma |}x_\sigma \hfill \\ & \underset{\sigma S}{}(1)^{|\sigma |}\underset{\mathrm{}\tau T}{}\pm (x_\sigma )_\tau \hfill \\ & =\underset{\mathrm{}\tau T}{}\pm \underset{\overline{\sigma }S_\tau }{}(1)^{|\overline{\sigma }|}x_{\overline{\sigma }}\hfill \\ & =\underset{\mathrm{}\tau T}{}\pm [x,S_\tau ],\hfill \end{array}$$

where $`S_\tau `$ is the set of interdependent moves guided by the circles modified by $`\tau `$. Thus it suffices to show that $`[x;S_\tau ]_k^v`$ for $`\tau \mathrm{}`$. Each of the resulting sets $`S_\tau ,\tau \mathrm{}`$, of interdependent moves has fewer edges plus stars in the resulting graph. We can always iterate this simplification unless the number of edges plus stars is less than or equal to $`k`$. It is hence sufficient to consider this case.

On such a graph, I claim that there are $`k+1`$ unstarred vertices, no pair of which has a connecting edge. This follows from the claim that there are at least $`k+1`$ connected components without any stars on them. To see this, let $`st`$ denote the number of stars. Then $`E+stk`$ implies $`Ekst`$. The Euler characteristic can be computed in two different ways: $`VE=b_0b_1`$, so that $`b_0b_0b_1=2K+1E(k+1)+st`$. Thus $`b_0stk+1`$.

The fact that the vertices are unstarred means they each bound a disk, and the lack of edges implies the disks are disjoint. There are now $`k+1`$ groups of finger moves of $`x`$ union the circles which push everything out of the $`k+1`$ disks. Thus, as above, modulo $`_k^v`$, we can assume that at least one of the circles guiding a move bounds an embedded disk with no intersections with anything in its interior. This means that the move is an isotopy. However $`[x;S]`$ is obviously $`0`$ if any of the moves in $`S`$ are isotopies. $`\mathrm{}`$

## 4 Application to rooted claspers and grope cobordism

For the reader’s convenience in following the proof of Theorem 3, I remind him that a rooted clasper is a clasper in the sense of Habiro(see \[H\]), embedded in the complement of a knot or link, such that there is one root leaf, which is a zero framed leaf linking the knot as a little meridian, whereas the other leaves can be embedded arbitrarily.

###### Theorem 3

Suppose two knots $`K_1`$ and $`K_2`$ are related by a rooted clasper surgery with $`2k+1`$ non-root leaves. Then $`K_1K_2_k^v=_{2k}^{alt}`$.

*\[Proof\]*

Surgering along the clasper, the original knot is modified inside a regular neighborhood of the clasper union the disk bounding its root leaf to get $`K_2`$. We will find $`2k+1`$ groups of interdependent moves as follows. For each non-root leaf of the clasper there is an interdependent move of a subarc of the leaf so that the result is a tiny leaf bounding a disk avoiding intersections with anything. Each of these moves on the clasper descend to a group of interdependent moves which take strands of the knot running through the neighborhood of the leaf and move them to a position corresponding to the modified clasper. Thus we have found a set $`S`$ of $`2k+1`$ groups of interdependent moves on $`K_2`$. Thus the surgered knot $`K_2`$, is congruent modulo $`_{2k}^{alt}`$ to the sum $`_{\mathrm{}\sigma S}(1)^{|\sigma |}(K_2)_\sigma `$, but for each nonempty $`\sigma `$, $`(K_2)_\sigma `$ is by construction just $`K_1`$ modified by a clasper with at least one trivial leaf. Such clasper surgeries do not change the knot, hence $`K_2`$ is congruent to $`_{\mathrm{}\sigma S}(1)^{|\sigma |}K_1=K_1`$.$`\mathrm{}`$

###### Corollary 4

Suppose two knots $`K_1`$ and $`K_2`$ are related by a grope cobordism of class $`2k+1`$. Then $`K_1K_2_k^v`$.

*\[Proof\]*

By Proposition 6 of \[CT\], we may assume that the grope cobordism has all of its stages of genus one. By Theorem 9 of the same paper, this is just a rooted clasper surgery with $`2k+1`$ non-root leaves. $`\mathrm{}`$

###### Remark:

Even though this Corollary is used in \[CT\] to prove Theorem 2 of that paper, Theorem 2 is not used to prove Theorem 9 or Proposition 6, so there is no logical circuity. Corollary 4 could have been proven directly in the same way as Theorem 3, using techniques of my thesis \[C\], but as there is no logical neccessity I have avoided the added complication.


text:
# Ten Low Mass Companions from the Keck Precision Velocity Survey¹

## 1 Introduction

All $``$70 known extrasolar planets have been discovered by the precision Doppler technique over the last six years (Marcy, Cochran, Mayor 2000). These include four published multiple planet systems (Butler et al. 1999; Marcy et al. 2001a, 2001b; Fischer et al. 2002) and one transiting planet (Henry et al. 2000; Charbonneau et al. 2000). The mass function of giant planets is emerging (Marcy & Butler 2000; Jorissen, Mayor & Udry 2001) and their occurrence correlates with metallicity (Butler et al. 2000; Santos, Israelian, & Mayor 2001).

Precision Doppler surveys are biased toward finding massive companions in small orbits, as these two characteristics enhance the Doppler amplitude and allow many orbits to be observed quickly. At the current epoch, all of the confirmed extrasolar planets orbit within 4 AU of their host stars. This is primarily due to the time baseline of the active surveys, but also because the Doppler amplitude of an orbiting companion decreases with increasing orbital radii. A 1 MJUP planet in a 0.05 AU orbit induces a maximum Doppler amplitude of 127 m s$`^1`$on its star, but at 3 AU this is reduced to 16 m s$`^1`$. While Doppler surveys with precision of 10 to 20 m s$`^1`$can easily detect 1 MJUP planets at 0.05 AU, they are hard pressed to detect similar planets at 3 AU.

Due to their larger mass, brown dwarf companions are much easier to detect than planets. A 20 MJUP brown dwarf companion at 0.05 AU induces a maximum Doppler amplitude of 2500 m s$`^1`$, a signal that has been in the detectable range for most of the last century. Yet no companion having mass between 10 and 80 MJUP has ever been found orbiting within 0.1 AU of a star.

This paper reports the discovery of five planet candidates, two brown dwarf candidates, and three low mass stellar companions from the Keck precision velocity survey. Hipparcos astrometric data is used to solve for, or set limits on, the inclination angle and maximum mass of these companions. Section 2 describes the Keck precision velocity program. The stellar properties and Keplerian orbital fits for the ten new low mass companions are presented in Section 3. A discussion follows, including an updated substellar companion mass function.

## 2 The Keck Planet Search Program

The Keck Planet Search Program makes use of the HIRES echelle spectrometer (Vogt et al. 1994) on the Keck I telescope. The resolution of these spectra is R $``$ 80000, spanning wavelengths from 3900–6200 Å. Wavelength calibration is carried out by means of an iodine absorption cell (Marcy & Butler 1992) which superimposes a reference iodine spectrum directly on the stellar spectra (Butler et al. 1996; Valenti et al. 1995). This system currently achieves photon–limited measurement precision of 3 m s$`^1`$(Vogt et al. 2000). Efforts are under way to improve the single–shot precision of this system to the 2 m s$`^1`$level.

Observations from the Lick Observatory precision velocity program are combined with the Keck observations for two of the stars presented here. Lick observations are made with either the 3–m Shane or 0.6–m Coudé Auxilliary (CAT) Telescopes, both of which feed the “Hamilton” echelle spectrometer (Vogt 1987). Wavelength calibration at Lick is also acomplished with an iodine absorption cell.

The Keck Planet Survey began in July 1996, and is currently surveying about 650 main sequence dwarfs ranging in spectral type from late F to mid–M. The spectrum of stars earlier than F7 do not contain enough Doppler information to achieve precision of 3 m s$`^1`$, while stars later than M5 are too faint even for Keck.

Nearly all suitable northern hemisphere G dwarfs within 50 pc, and K dwarfs within 30 pc, are included in either the Keck survey or the Lick 3–m survey (Fischer et al. 2002). Evolved stars have been removed from the observing list based on Hipparcos distances (Perryman et al. 1997, ESA 1997). The list has been further sieved to remove chromospherically active stars as these stars show velocity “jitter” of 10 to 50 m s$`^1`$, related to rapid rotation, spots, and magnetic fields (Saar et al. 1998). The Ca II H&K lines are used as a chromospheric diagnostic (Noyes et al. 1984). We measure the strength of the H&K line reversal directly from our Keck HIRES spectra. Keck H&K measurements are placed on the Mt Wilson “S” scale by calibration with previously published results (Duncan et al. 1991; Baliunas et al. 1995; Henry et al. 1996).

Stars with known stellar companions within 2 arcsec are removed from the observing list as it is operationally difficult to get an uncontaminated spectrum of a star with a nearby companion. Otherwise there is no bias against observing multiple stars. The Keck target stars also contain no bias against brown dwarf companions.

## 3 New Companions from the Keck Survey

Five planet–mass candidates, 2 brown dwarf candidates, and 3 stellar companions have emerged from the Keck survey. The stellar properties of the 10 host stars are given in Table 1. The first two columns provide the HD catalog number and the Hipparcos catalog number respectively. Spectral types are from a calibration of $`BV`$ and Hipparcos derived absolute magnitudes. The stellar masses are estimated by interpolation of evolutionary tracks (Fuhrmann 1998, Fuhrmann et al. 1997). The \[Fe/H\] values are drawn from a variety of sources (given in Section 3), including spectral synthesis matched directly to our Keck HIRES spectra.

The R’HK values, a chromospheric activity indicator (Noyes et al. 1984), are measured from the CaII H&K line cores in our Keck spectra. The level of Doppler “jitter” is correlated with R’HK (Saar et al. 1998). Slowly rotating, chromospherically inactive stars are intrinsically stable to at least the 3 m s$`^1`$level, while the Doppler “jitter” for young rapidly rotating stars ranges from 10 to 50 m s$`^1`$.

Figure 1 shows the H line for the 7 G dwarfs reported in this paper. The Sun is shown for comparison. Of these stars, HD 33636 is the most active, with R’HK $`=`$ $`4.81`$, indicating a rotation period of 13.6 d, and expected Doppler “jitter” of $``$7 m s$`^1`$. The other six G dwarfs have R’HK values similar or lower than the Sun, indicating rotation periods of 25 days or longer, and intrinsic Doppler “jitter” of 3 m s$`^1`$or less.

The H lines for the three K dwarfs reported in this paper are shown in Figure 2. All three of these stars are slow rotators. The rapid rotator HD 128311 (K0 V) is shown at the bottome for comparison. While even the slowly rotating K dwarfs show a slight line reversal in the core of the H line, they are clearly distinguished from rapid rotators like HD 128311.

The orbital parameters of the ten low mass companions are listed in Table 2, while the individual Keck Doppler velocity measurements are listed in Tables 3 through 12. The host stars are discussed below.

The Hipparcos Intermediate Astrometric Data (ESA 1997) have been analyzed with the orbits derived from the precision velocity data, using the technique outlined by Pourbaix & Arenou (2001) to constrain, or in three cases to solve for, the orbital inclination.

### 3.1 HD 4208

Based on Stromgren photometry, Eggen (1998) finds the metallicity of HD 4208 (G5 V) to be \[Fe/H\] = -0.21, in good agreeement with our estimate of -0.24 from spectral synthesis matched to our Keck HIRES spectra. The star is photometrically stable at the level of Hipparcos measurement error, $``$0.01 mag. Based on the $`BV`$ color, the Hipparcos derived absolute magnitude, and the metallicity, we estimate the mass of the primary to be 0.93 M. This star is slowly rotating and chromospherically inactive as indicated by the R’HK value.

Thirty–five Keck Doppler velocity observations have been made of HD 4208, spanning 4.9 years, as shown in Figure 3 and listed in Table 3. These observations cover two full orbital periods. The semiamplitude ($`K`$) of the Keplerian orbital fit is 18 m s$`^1`$, only the third extrasolar planet yet published with an amplitude less than 20 m s$`^1`$(Marcy et al. 2000; Fischer et al. 2002). The RMS of the velocity residuals to the Keplerian fit is 5.21 m s$`^1`$, slightly worse than the typical internal measurement error of 4.2 m s$`^1`$. The reduced $`\chi _\nu ^2`$ of the Keplerian fit is 1.33. Within measurement error the orbit is circular, thus joining 47 UMa (Butler & Marcy 1996; Fischer et al. 2002) and HD 27442 (Butler et al. 2001) as the only published planets in circular orbits beyond 0.2 AU.

With a semimajor axis $`a`$ = 1.67 AU, the maximum angular separation between planet and star for HD 4208 is 53 milli-arcsec. The orbit is similar to that of Mars in the Solar System. The companion is not detected in the Hipparcos astrometric data, thus constraining the orbital inclination to be larger than 1.4 deg.

Unlike most of the extrasolar planet candidates, HD 4208 is modestly metal poor relative to the Sun.

### 3.2 HD 114783

We estimate the metallicity of HD 114783 (K2 V) to be \[Fe/H\] = +0.33 based on uvby photometry. Strassmeier et al. (2000) report that the star is chromospherically inactive, in agreement with our Keck–derived value of log(R’HK) = -4.96. HD 114783 is photometrically stable at the level of Hipparcos measurement error, $``$0.01 mag.

A total of 37 Keck observations have been obtained between 1998 June and 2001 August, as shown in Figure 4 and Table 4. These observations cover slightly more than two full orbital periods of 501 days. The semiamplitude ($`K`$) of the Keplerian orbital fit is 27 m s$`^1`$, with an eccentricity of 0.10. The RMS to the Keplerian fit is 4.08 m s$`^1`$, with a reduced $`\chi _\nu ^2`$ of 1.48. There is 5$`\sigma `$ discrepant point near 2000.35. A periodogram of the Keplerian residuals reveal no other periodicities with a false alarm probability under 5%.

At a distance of 20.4 pc, HD 114783 is an attractive astrometric target. With a semimajor axis of 1.2 AU, the maximum separation between the planet and the host star is 65 milli-arcsec. The astrometric amplitude of the star due to the planet is 61 mas/$`\mathrm{sin}i`$, which should be easily within the capabilities of both ground and space–based interferometry within a few years. The orbital inclination is constrained to be greater than 1.26 degrees by the Hipparcos astrometric data.

### 3.3 HD 4203

HD 4203 (G5 V) is photometrically stable at the level of Hipparcos measurement error, and chromospherically quiet with a measured R’HK value of -5.13. At 77.8 pc, this is one of the most distant stars on this project. This star was added to the Keck project in 2000 July based on the suggestion of Laughlin (2000), who noted the star is metal rich (\[Fe/H\] = +0.22). ¿From interpolation of isochrones, Prieto & Lambert (1999) estimate the mass of the HD 4203 to be 0.98 M. This estimate does not take into account the extreme metallicity of the star. Based on its similarity to 51 Peg (Marcy et al. 1997), we estimate the mass of HD 4203 to be 1.06 M.

A total of 14 Keck observations have been obtained between 2000 July and 2001 October. These observations are listed in Table 5 and graphically displayed in Figure 5. The best–fit Keplerian orbit to these data has a period of 406 d, a semiamplitude (K) of 51 m s$`^1`$, and an eccentricity of 0.53. The minimum ($`M\mathrm{sin}i`$) mass of the companion is 1.6 MJUP. The RMS to the Keplerian fit is 3.97 m s$`^1`$, consistent with measurement uncertainty.

As only one orbit has been observed, the orbital period and amplitude remain somewhat uncertain. Only a few observations have been obtained near the velocity maximum. The next passage through maximum velocity will occur in Spring 2002. The orbital inclination is constrained to be greater than 0.44 deg by Hipparcos astrometry.

### 3.4 HD 68988

Based on Stromgren photometry, Laughlin (2000) finds the metallicity of HD 68988 (G2 V) to be \[Fe/H\] $`=`$ $`+0.36`$, consistent with the estimate of Feltzing & Gustafsson (1998), but somewhat higher than the +0.24 derived from matching spectral synthesis to our Keck HIRES spectra. HD 68988 is photometrically stable at the level of Hipparcos measurement error, $``$0.01 mag. We estimate the mass of the primary to be 1.2 M. This star is slowly rotating and chromospherically inactive as indicated by our measurement of the CaII line cores, log(R’HK) = -5.07. The age of the star is estimated to be approximately 6 Gyr based on the chromospheric diagnostic.

Thirteen Keck Doppler observations and six Lick observations have been made of HD 68988 spanning 1.4 years, as shown in Figure 6 and listed in Table 6. The orbital period is 6.276 d, the semiamplitude ($`K`$) is 187 m s$`^1`$, and the eccentricity is $`e`$ $`=`$ $`0.14`$. The RMS of the velocity measurements to the best–fit Keplerian is 4.36 m s$`^1`$, yielding a reduced $`\chi _\nu ^2`$ of 1.05. The minimum mass of this companion is 1.9 MJUP. In addition, there is a linear trend of -0.072 m s$`^1`$per day, indicating a second companion with an orbital period much greater than 4 years.

All of the published planets with orbital periods of less than 1 week111 HD 46375 (Marcy et al. 2000), HD 179949 (Tinney et al. 2001), HD 187123 (Butler et al. 1998), $`\tau `$ Boo (Butler et al. 1997), BD–103166 (Butler et al. 2000), HD 75289 (Udry et al. 2000, Butler et al. 2001), HD 209458 (Henry et al. 2000), 51 Peg (Mayor & Queloz 1995; Marcy et al. 1997), $`\upsilon `$ And b (Butler et al. 1997; Butler et al. 1999) have circular orbits. It is therefore surprising that the orbit of HD 68988 is markedly non-circular. This system is similar to the companion orbiting HD 217107 (Fischer et al. 1999) with its period of 7.12 d, $`M\mathrm{sin}i`$$`=`$ $`1.25`$ MJUP, and eccentricity of 0.134. HD 68988 and HD 217107 are both markedly metal rich.

The time scale for orbital circularization due to tidal dissipation within a planet is given by Marcy et al. (1997):

$$t_{circ}\frac{4Q}{63}\frac{m}{M}\frac{P_{orb}}{2\pi }\left(\frac{a}{R_p}\right)^5.$$
(1)
Assuming a radius 1.3 R$`_{\mathrm{JUP}}`$for the companion to HD 68988, this reduces to 5000Q/$`\mathrm{sin}i`$years, where Q is the tidal quality factor (Goldreich & Soter 1966). For Q values less than $``$ 106/$`\mathrm{sin}i`$, the orbit of HD 68988 should have circularized. The Q value for Jupiter is inferred to be, 105 $`<`$ Q $`<`$ 106 (Goldreich & Soter 1966; Yoder 1979; Yoder & Peale 1981). The Q value for solid planets is in the range of 100 to 1000 (Marcy et al. 1997). The planet orbiting HD 68988 is thus likely to be a gas giant. Three possibilities exist to explain the observed orbital eccentricity: the planet has a higher Q value than is inferred for Jupiter; it has only recently arrived in its current orbit; or it is being perturbed by another body.

### 3.5 HD 33636

Spectral synthesis matched to our Keck spectra yields \[Fe/H\] = -0.13 for HD 33636 (G0 V). The star is photometrically stable at the level of Hipparcos measurement error, and mildly active with log(R’HK) = -4.81. The Doppler velocity “jitter” associated with this level of activity for a G0 V star is $``$7 m s$`^1`$(Saar et al. 1998). The Hipparcos derived distance of HD 33636 is 28.7 pc.

A total of 21 observations of HD 33636 have been made at Keck between 1998 January and 2001 October. A further 11 observations were made at Lick with the 3-m Shane telescope and the 0.6-m CAT between 1998 Jan and 2001 Feb. These observations are listed in Table 7 and graphically displayed in Figure 7. The best–fit Keplerian orbit to the combined Keck and Lick data sets has a period of 1553 d, a semiamplitude (K) of 148 m s$`^1`$, and an eccentricity of 0.39. The minimum ($`M\mathrm{sin}i`$) mass of the companion is 7.7 MJUP. The RMS to the Keplerian fit is 8.69 m s$`^1`$, consistent with the expected Doppler “jitter” due to activity.

The Keck and Lick observations of HD 33636 cover nearly the full amplitude of the velocity variation, but slightly less than one orbital period, leading to a large uncertainty in the derived orbital period. Hipparcos astrometry is not able to place limits on the orbital inclination of this system because the orbital period is much longer than the duration of the Hipparcos mission.

### 3.6 HD 169822

Spectral synthesis matched to our Keck spectra yield \[Fe/H\] = -0.10 for HD 169822 (G5 V). The star is chromospherically quiet with R’HK $`=`$ $`4.97`$, and photometrically stable at the level of Hipparcos measurement error.

We have obtained a total of 22 Keck observations between 1999 July and 2001 October. The measured Doppler velocities are shown in Figure 8 and listed in Table 8. The best–fit Keplerian to these velocities yields a period of 293 d, an eccentricity of 0.48, and a semiamplitude ($`K`$) of 991 m s$`^1`$. The minimum ($`M\mathrm{sin}i`$) mass of the companion is 27.2 MJUP. The RMS to the Keplerian is 5.70 m s$`^1`$.

The Hipparcos catalog assigned a stochastic solution for HD 169822 because none of the investigated models gave a satisfactory solution at that time. The Hipparcos data has been re–investigated in light of the Doppler velocity signal, and now clearly yields a periodicity of 292.7 days, in excellent agreement with the Doppler velocity data. The new Hipparcos solution, taking into account the orbital motion, changes the measured parallax to this system from 37.04 mas to 31.2 mas, increasing the distance of this system from 27 to 32 pc. The Hipparcos–derived orbital inclination of 175 degrees yields a companion mass of 0.30 M. The companion is thus stellar, presumably with a spectral type around M3 V. With a semimajor axis of 0.84 AU, the maximum separation between the companion and the primary is 39 milli-arcsec. The M dwarf companion is expected to be $``$5.5 $`V`$ magnitudes fainter than the primary.

### 3.7 HD 184860

The metallicity of HD 184860 (K2 V), determined by spectral synthesis matched to our Keck spectra, is \[Fe/H\] $`=`$ $`0.13`$. The star is chromospherically quiet, with R’HK $`=`$ $`5.00`$. Hipparcos $`V`$ band photometric scatter is 0.016 mag.

A total of 21 Keck Doppler velocity of HD 184860 have been made between 1996 July and 2001 August. These observations are shown graphically in Figure 9 and in tabular form in Table 9. The best–fit Keplerian yields an orbital period of 693.0 d, a semiamplitude $`K`$ $`=`$ 1123 m s$`^1`$, eccentricity $`e`$ $`=`$ 0.67, and $`M\mathrm{sin}i`$$`=`$ $`32`$ MJUP. As no observations have been made near minimum velocity, the semiamplitude is poorly constrained. Minimum velocity will next occur in 2002 December.

With a semi–major axis of 1.44 AU, the companion to HD 184860 is separated from the primary by 47.5 milli-arsec. The minimum astrometric semiamplitude of the primary is $``$1.8 milli-arsec. Unfortunately, HD 184860 has a stellar companion separated by 5.1 arcsec. This stellar companion shows up in the Hipparcos Intermediate Astrometric Data, making it impossible to search for the astrometric signal that is caused by the brown dwarf candidate. It is thus not possible to derive an upper limit for the orbital inclination of the brown dwarf candidate. Moderate values of $`\mathrm{sin}i`$would render this an M dwarf.

### 3.8 HD 64468

Spectral synthesis matched to our Keck spectra yield \[Fe/H\] $`=`$ $`+0.00`$ for HD 64468 (K2 V). The star is chromospherically quiet with R’HK $`=`$ $`5.03`$, and photometrically stable within Hipparcos measurement error.

Thirteen observations of HD 64468 have been made at Keck between 1997 Jan and 2001 April. These observations are listed in Table 10 and graphically displayed in Figure 10. The best–fit Keplerian orbit to this data has a period of 161 d, a semiamplitude ($`K`$) of 5730 m s$`^1`$, and an eccentricity of 0.26. The minimum ($`M\mathrm{sin}i`$) mass of the companion is 0.13 M. As the companion is clearly stellar, we will drop this star from our observing program.

The orbital inclination of this system was derived from the Hipparcos astrometric data by freezing the orbital parameters from the Doppler velocities, yielding an inclination of 64 $`\pm `$13 degrees, and a true mass for the companion of 0.14 M. The companion is likely to be an M6 dwarf. With a semimajor axis of 0.56 AU, the maximum separation between the companion and the primary is 28 milli-arcsec. The M dwarf companion is expected to be $``$7.5 $`V`$ magnitudes fainter than the primary.

### 3.9 HD 35956A

Based on Stromgren photometry, Eggen (1998) reports the metallicity of HD 35956A (G0 V) to be \[Fe/H\] = -0.14, in reasonable agreeement with our photometric estimate of -0.21. The star is photometrically stable at the level of Hipparcos measurement error, $``$0.01 mag, and chromospherically quiet with R’HK = -4.92. The Hipparcos derived distance of this star is 28.9 pc. We estimate the mass of the primary to be 1.0 M.

We have obtained a total of 14 Keck observations between 1996 October and 2001 October. The measured Doppler velocities are shown in Figure 11 and listed in Table 11. The best–fit Keplerian orbit to these velocities yields a period of 1427 d, an eccentricity $`e`$ $`=`$ $`0.62`$, and a semiamplitude, $`K`$ = 3796 m s$`^1`$. The RMS to the Keplerian fit is 4.89 m s$`^1`$. The minimum ($`M\mathrm{sin}i`$) mass of the companion is 0.18 M. As this companion is stellar, we will drop this star from our target list.

Although the orbital period of this binary is about 50% longer than the duration of the Hipparcos mission, the astrometric signature of the companion was detected. Freezing the orbital elements from the precision velocities yields an orbital inclination of 78 $`\pm `$4 degrees. The true mass of the companion is thus only 2% larger than the minimum mass. The spectral type of the companion is expected to be roughly M6, with $`V`$ magnitude $``$16, separated from the $`V`$ $`=`$ $`6.7`$ primary by 0.1 arcsec.

### 3.10 HD 43587

HD 43587 is assigned a spectral type of F9 V by Simbad and G0.5 by Hipparcos. The $`BV`$ color of $`0.61`$ is consistent with G1 or G2 spectral type. Spectral synthesis matched to our Keck spectra yields Teff = 5795 K, \[Fe/H\] = -0.03, and $`V\mathrm{sin}i`$= 2.7 km s$`^1`$. The R’HK value measured from the H & K lines in the Keck spectra is -4.97. The star is photometrically stable at the level of Hipparcos measurement error. HD 43587 is thus a near solar twin.

Fourteen observations of HD 43587 have been made at Keck between 1996 Oct and 2001 April. These observations are listed in Table 12 and graphically displayed in Figure 12. The best–fit Keplerian to this data has a period of 33.7 yr, a semiamplitude (K) of 4323 m s$`^1`$, and an eccentricity of 0.80. The minimum ($`M\mathrm{sin}i`$) mass of the companion is 0.34 M. Given the long orbital period, we were fortunate to observe the extreme velocity drop of 1998. As the companion is clearly stellar, we will drop this star from our observing program. Hipparcos was unable to astrometrically detect the companion due to the short duration of the Hipparcos mission.

This star was observed 7 times during 5.3 yr around epoch 1987 by the CORAVEL Doppler velocity survey (Duquennoy & Mayor 1991), and found to be constant, with standard deviation of 0.36 km s$`^1`$. That constant velocity is consistent with our results. The present best–fit Keplerian predicts the last epoch of rapid velocity change would have been in 1964. During the epoch when CORAVEL measurements were made, the velocity of the star would have changed by only 0.2 km s$`^1`$per year. We expect that in retrospect, the CORAVEL velocities should reveal a trend just above the errors. The absolute velocity reported by Duquennoy & Mayor was 9.6 km s$`^1`$. We find that on 2 Dec 1997 (JD = 2450785) the velocity of HD43587 was 12.3$`\pm `$0.1 km s$`^1`$(Nidever et al. 2002). Thus, the measured change in the velocity of the star from 1987 to 2 Dec 1997 was +2.7 km s$`^1`$. In comparison, our orbital solution predicts a change of +2.4 km s$`^1`$during that 10–year period, in good agreement within errors.

With a semimajor axis a = 11.6 AU, the companion orbits about 0.6 arcseconds from the primary. For a minimum mass companion ($`\mathrm{sin}i`$= 1), the primary would be about 5.5 V magnitudes brighter than the companion. The value of $`\mathrm{sin}i`$is likely greater than 0.5, as smaller values would yield a companion with an expected difference in V magnitude of 3 or less, which would be detectable in our Keck spectra. This binary thus makes an interesting AO target.

## 4 Discussion

A total of 1,200 stars are currently being surveyed by the Keck, Lick, and Anglo–Australian precision velocity surveys. All three of these programs use the Iodine cell technique, and all three have demonstrated long term precision of 3 m s$`^1`$(Butler & Marcy 1997; Vogt et al. 2000; Butler et al. 2001). These three surveys have no selection bias against stars with brown dwarf companions.

A total of 48 substellar candidates ($`M\mathrm{sin}i`$$`<`$ 80 MJUP) have been uncovered from these surveys, of which 44 are planet candidates with plausible masses below the Deuterium burning limit ($`M\mathrm{sin}i`$$`<`$ $`12`$ MJUP). All of these companions have either been published or are currently submitted to referee journals, including 3 new companions from the AAT survey (Tinney et al. 2002; Jones et al. 2002).

Two of the four brown dwarf candidates from these surveys, HD 169822 and HD 164427 (Tinney et al. 2001), have been revealed by Hipparcos astrometry to be M dwarfs. Hipparcos is able to place upper limits on the masses of the remaining 49 companions (Pourbaix & Arenou 2001; Perryman et al. 1996).

The mass function of the 46 surviving substellar candidates from the Keck, Lick, and Anglo–Australian surveys is shown in Figure 13. All but 3 of these companions orbit within 3 AU. The mass function is flat and sparsely populated above 10 MJUP, and then begins rising abruptly below 10 MJUP.

Within 3 AU these surveys are complete for companions having $`M\mathrm{sin}i`$$`>`$ 10 MJUP. Incompleteness is greatest for the smallest mass bin, $`M\mathrm{sin}i`$$`<`$ 1 MJUP. Companions with $`M\mathrm{sin}i`$$`<`$ 1 MSAT are only detectable if they orbit within 1 AU. These selection effects strongly favor the detection of companions with $`M\mathrm{sin}i`$$`>`$ 10 MJUP at the expense of companions with $`M\mathrm{sin}i`$$`<`$ 1 MJUP. In spite of this, 15 companions having $`M\mathrm{sin}i`$$`<`$ 1 MJUP have been detected, while only two have been found with $`M\mathrm{sin}i`$$`>`$ 10 MJUP.

The abrupt and discontinuous jump in the mass function below 10 MJUP empirically sets the threshold between planets and brown dwarfs at $``$10 MJUP. Coincidentally the Deuterium-burning limit resides near this limit at 12 MJUP.

The two companions with the largest mass in Figure 13 are HD 168443c ($`M\mathrm{sin}i`$= 17 MJUP) and HD 184860 ($`M\mathrm{sin}i`$= 32 MJUP). As described in the previous section, Hipparcos can not put an upper limit on the mass of the companion to HD 184860 because of a nearby stellar companion. A moderate inclination of 23 degrees would be be sufficient to push this brown dwarf candidate into the M dwarf mass range. For the case of HD 168443c, Hipparcos astrometry is able to put an upper limit on the mass of the companion of 42 MJUP (Marcy et al. 2001a). Brown dwarf companions orbiting within 3 AU of main sequence stars are rare, with at most about 1 such companion per 500 stars. In contrast, planetary mass companions are common.

Relatively few planets have been found orbiting metal poor stars (Butler et al. 2000, Santos et al. 2001). Planets orbiting stars of less than solar metallicity (\[Fe/H\] $`<`$ $`+0.00`$) are indicated in Figure 13 by the cross hatched area. They constitute about 12% of the planets found in the Keck, Lick, and Anglo-Australian surveys (even though two–thirds of the field stars have metallicity less than solar). If planets are common around metal poor stars, these planets must either have smaller masses or larger orbital radii than the planets found to date.

We acknowledge support by NSF grant AST-9988358 (to SSV), NSF grant AST-9988087 and travel support from the Carnegie Institution of Washington (to RPB), NASA grant NAG5-8299 and NSF grant AST95-20443 (to GWM), and by Sun Microsystems. We thank the NASA and UC Telescope assignment committees for allocations of telescope time. This research has made use of the Simbad database, operated at CDS, Strasbourg, France.


text:
# What is the fate of the heavy electron at a quantum critical point?

## A Introduction.

Discoveries over the past decade have brought a new awareness of the importance of quantum critical points in condensed matter physics. A quantum critical point (QCP) is a zero-temperature instability between two phases of matter where quantum fluctuations develop long range correlations in both space and time. These special points exert wide-reaching influence on the finite temperature properties of a material. Systems close to quantum criticality develop a new excitation structure, they display novel thermodynamic, transport and magnetic behavior. They also have marked a predeliction towards the development of new kinds of order, such as anisotropic superconductivity. A dramatic example is provided by the cuprate superconductors. By doping with holes, these materials pass through one or more quantum phase transitions: from a Mott insulator to a

metal with a spin gap at low doping, and at higher doping a second QPT may occur when the spin gap closes  (Fig. 1). The singular interactions induced by quantum criticality may be the driving force for both the high temperature superconductivity and the anomalous metallic state above the spin gap temperature $`T^{}`$.

This paper discusses quantum criticality in heavy fermion materials. These materials offer a unique opportunity to study quantum criticality under controlled conditions. By the application of pressure or doping, heavy fermion materials can be tuned through a quantum critical point from a metallic antiferromagnet into a paramagnet (Fig. 2). Unlike the cuprate metals, here the the paramagnetic phase is a well characterized Fermi liquid, with heavy Landau quasiparticles, or “heavy electrons”. A central property of these quasiparticles, is the existence of a finite overlap “$`Z`$” between a single quasiparticle state, denoted by $`|\text{qp}^{}`$ and the state formed by adding a single electron to the ground-state, denoted by $`|e^{}=c_{𝐤\sigma }^{}|0`$. This quantity is closely related to the ratio $`m/m^{}`$ of the electron to quasiparticle mass,

$$Z=|e^{}|\text{qp}^{}|^2\frac{m}{m^{}}.$$
(1)
A wide body of evidence suggests that $`m^{}/m`$ diverges at a heavy fermion QCP, indicating that

$`Z0(PP_c).`$

The state which forms at the QCP is referred to as a “non-Fermi” or “singular Fermi liquid”.  By what mechanism does this break-down in the Landau quasiparticle occur?

## B Properties of the Heavy Fermion Quantum Critical Point

There is a growing list of heavy fermion systems that have been tuned to an antiferromagnetic QCP by the application of pressure or by doping (Table 1.). These materials display many common properties

Table. 1. Selected Heavy Fermion compounds with quantum critical points.    Compound $`x_c`$/$`P_c`$/$`H_c`$ $`\frac{C_v}{T}\mathrm{}`$? $`\rho T^a`$ Ref.  $`YbRh_2Si_{2x}Ge_x`$ $`x_c=0.1`$ $`Log\left(\frac{T_o}{T}\right)`$ (b) $`T`$   $`CeCu_{6x}Au_x`$ $`x_c=0.1`$ $`Log\left(\frac{T_o}{T}\right)`$ $`T+c`$   $`CeCu_{6x}Ag_x`$ $`x_c=0.09`$ $`Log\left(\frac{T_o}{T}\right)`$ $`T^{1.1}`$   $`CeNi_2Ge_2`$ $`P_c=0`$ $`Log\left(\frac{T_o}{T}\right)`$(a) $`T^{1.2}`$   $`U_2Pt_2In`$ $`P_c=0`$ $`Log\left(\frac{T_o}{T}\right)`$ $`T`$   $`U_2Pd_2In`$ $`P_c<0`$ ? $`T+c`$   $`CePd_2Si_2`$ $`P_c>0`$ ? $`T^{1.2}`$   $`CeIn_3`$ $`P_c>0`$ ? $`T^{1.5}`$   $`U_3Ni_3Sn_4`$ $`P_c>0`$ no ?   $`CeCu_2Si_2`$ $`P_c=0`$ no $`T^{1.5}`$   $`CeRu_2Si_2`$ $`H_c^{}=7.7T`$ $`Log\left(\frac{T_o}{T}\right)`$    $`UPt_3`$ $`H_c^{}=20T`$ $`Log\left(\frac{T_o}{T}\right)`$ $`T^{1.2}`$   $`Sr_3Ru_2O_7`$ $`H_c^{}=7.7T`$ $`Log\left(\frac{T_o}{T}\right)`$ $`T^{1.1}`$

(a) New data show a stronger divergence at lower temperatures, and $`\gamma AB\sqrt{T}`$ at intermediate temperatures.

(b) At low temperatures, $`\gamma `$ diverges more rapidly than $`Log\left(\frac{T}{T_o}\right)`$.

* Fermi liquid behavior in the paramagnet, as indicated by the emergence of a quadratic temperature dependence in the resistivity in the approach to the QPT $`\rho =\rho _o+AT^2`$  at ever lower temperatures.
* Divergent specific heat at the QCP, typically with a logarithmic temperature dependence,

$$\gamma (T)=\frac{C_v(T)}{T}=\gamma _0\mathrm{log}\left[\frac{T_o}{T}\right],$$
(2)
suggesting that the Fermi temperature vanishes and the quasiparticle effective masses diverge

$$T_F^{}0,\frac{m^{}}{m}\mathrm{}$$
(3)
at the QCP. Further support for this conclusion is provided by the observation that the quadratic coefficient $`A`$ of the resistivity grows, and may diverge in the approach to the quantum critical point.
* Quasi-linear resistivity

$$\rho T^{1+ϵ},$$
(4)
at the QCP with $`ϵ`$ in the range $`00.6`$. In critical $`YbRh_2Si_{2x}Ge_x`$, $`\rho T`$ over three decades.
* Non-Curie spin susceptibilities

$$\chi ^1(T)=\chi _0^1+cT^a$$
(5)
with $`a<1`$ observed in critical $`CeCu_{6x}Au_x`$ (x=0.1), $`YbRh_2Si_{2x}Ge_x`$ (x=0.1) and $`CeNi_2Ge_2`$.
* $`E/T`$ and $`H/T`$ Scaling. In critical $`CeCu_{6x}Au_x`$ and $`YbRh_2Si_{2x}Ge_x`$ the differential magnetic susceptibility $`dM/dH`$ exhibits $`H/T`$ scaling,

$$(dM/dH)^1=\chi _0^1+cT^ag[H/T],$$
(6)
where $`a0.75`$. Neutron measurements show $`E/T`$ scaling in the dynamical spin susceptiblity of critical $`CeCu_{6x}Au_x`$, throughout the Brillouin zone, parameterized in the form

$$\chi ^1(𝐪,\omega )=T^af(E/T)+\chi _0^1(𝐪)$$
(7)
$`F[x](1ix)^a`$. Scaling behavior with a single anomalous exponent in the momentum-independent component of the dynamical spin susceptibility suggests an emergence of local magnetic moments which are critically correlated in time at the quantum critical point.

Growing evidence suggests a strong connection between the heavy fermion QCP and “metamagnetism”. A number of strongly correlated electron systems, such as $`CeRu_2Si_2`$ and $`UPt_3`$, and $`Sr_3Ru_2O_7`$ exhibit a sudden rise in the magnetic polarization at a finite “critical” field. Properties characteristic of an heavy fermion QCP, (such as the logarithmic dependence of the specific heat and a quasi-linear resistivity) are seen to develop at such metamagnetic transitions. In all likelihood, this reflects the close vicinity to a ferromagnetic QCP (Fig 3), but why is the physics similar to that of a heavy fermion QCP?

## C Universality

Usually, the physics of a metal above its Fermi temperature depends on the detailed chemistry and band-structure of the material: it is non-universal. However, if the renormalized Fermi temperature $`T_F^{}(P)`$ can be tuned to become arbitrarily small compared with the characteristic scales of the material as one approaches a QCP, we expect that the “high energy” physics above the Fermi temperature $`T_F^{}`$ is itself, universal.

Quantum critical behavior implies a divergence of the long distance and long-time correlations in the material. Finite temperatures introduce the cutoff timescale

$$\tau _T=\frac{\mathrm{}}{k_BT}$$
(8)
beyond which coherent quantum processes are dephased by thermal fluctuations. Renormalization group principles imply that the quantum critical physics has an upper-critical dimension $`d_u`$. For $`d<d_u`$, the $`\tau _T`$ becomes the  correlation time $`\tau `$ of the system, so frequency dependent correlation functions and response functions take the form

$$F(\omega ,T)=\frac{1}{\omega ^\alpha }f(\omega \tau _T)=\frac{1}{\omega ^\alpha }f(\mathrm{}\omega /k_BT).$$
(9)
leading to $`E/T`$ scaling. By contrast, for $`d>d_u`$ the correlation time is sensitive to the details of the short-distance interactions between the critical modes, and in general $`\tau ^1T^{1+b}`$ Thus $`E/T`$ scaling with a non-trivial exponent strongly suggests that the underlying physics of the heavy fermion quantum critical point is governed by universal physics with $`d_u>3`$.

## D Failure of the Spin Density Wave picture

The commonly accepted picture of the heavy fermion QCP assumes the non-Fermi liquid behavior derives from Bragg diffraction of the electrons off a quantum-critical spin density wave (QSDW). The virtual emission of these soft fluctuations,

$$e^{}e^{}+\text{spin fluctuation}$$
(10)
generates a retarded interaction

$$V_{eff}(𝐪,\omega )=g^2\stackrel{\chi (𝐪,\omega )}{\stackrel{}{\left[\frac{\chi _0}{(𝐪𝐐)^2+\xi ^2\frac{i\omega }{\mathrm{\Gamma }_𝐐}}\right]}}$$
(11)
between the electrons, where $`\chi (𝐪,\omega )`$ is the dynamical spin susceptibility of the collective modes. The damping term $`i\omega /\mathrm{\Gamma }_𝐐`$ of the magnetic fluctuations is derived from the linear density of particle-hole states in the Fermi sea. $`\xi ^1`$ and $`\tau ^1=\mathrm{\Gamma }_𝐐\xi ^2`$ are the the inverse spin correlation length and correlation times respectively. In real space,

$$V_{eff}(r,\omega =0)\frac{e^{r/\xi }}{r}e^{i𝐐𝐫}$$
(12)
is a “modulated ” Yukawa potential whose range $`\xi (PP_c)^{\frac{1}{2}}\mathrm{}`$ at the QCP. Unlike a ferromagnetic QCP, the modulated potential only effects electron quasiparticles along “hot lines” on the Fermi surface, that are separated by the wave-vector $`𝐐`$ and satisfy $`ϵ_𝐤=ϵ_{𝐤+𝐐}`$. At a finite temperature, electrons within a momentum range $`\sqrt{T}`$ are affected by this critical scattering (Fig. 4.). This limits the ability of this singular potential to generate non Fermi liquid behavior.

There are then two major difficulties with the QSDW scenario for the heavy fermion QCP:

1. No breakdown of the Fermi liquid Away from the hot lines, the Fermi surface and Landau quasiparticles remain intact at the QCP. Thus the specific heat and typical quasiparticle mass do not diverge but exhibit a weaker singularity, $`C_V/T=\gamma _oA\sqrt{T}`$ in the QSDW picture.
2. No $`E/T`$ scaling  The quantum critical behavior predicted by this model has been extensively studied. In the interaction $`V_{eff}(𝐪,\omega )`$ the momentum dependence enters with twice the power of the frequency, so

$`\tau \xi ^z,(z=2).`$

In the renormalization group (RG) treatment time counts as $`z`$ space dimensions so the effective dimensionality is $`D_{eff}=d+z=d+2`$. The upper critical dimension is set by $`D_{eff}=4`$, or $`d_u=2`$, so 3D quantum spin fluctuations will not lead to $`E/T`$ scaling. In three dimensions, QSDW theory predicts that the scale entering into the energy dependent response functions should scale as $`T^{3/2}`$, with a non-universal prefactor.

### 1 Can a 2D spin fluid cure the difficulties?

One explanation of $`E/T`$ scaling and the logarithmically divergent specific heat is to suppose that the spin fluctuations form a quasi-two-dimensional spin fluid, lying at the critical dimension. Inelastic neutron scattering experiments on $`CeCu_{6x}Au_x`$, (x=0.1) support a kind of reduced dimensionality in which the critical scattering is concentrated along linear, rather than at point-like regions in reciprocal space. More recent data may support quasi-2D spin fluctuations at intermediate scales in $`CeGe_2Ni_2`$.

## E New Approaches

Unfortunately, quasi-two dimensionality still cannot explain the anomalous exponents in the $`E/T`$ scaling, and we are essentially forced to consider the possibility of a fundamentally new interacting fixed point. Two approaches have emerged:

### 1 Local Quantum Criticality

The momentum-independent scaling term in the inverse dynamic susceptibility (7) suggests that the critical behavior associated with the heavy fermion QCP contains some kind of local critical excitation. One possibility, is that this critical excitation is the spin itself, which would then presumably develop a slow power-law decay

$$S(\tau )S(\tau ^{})=\frac{1}{(\tau \tau ^{})^{2ϵ}},$$
(13)
where $`ϵ0`$ signals non- Fermi liquid behavior.

Si et al have developed this idea, proposing that the local  spin susceptibility $`\chi _{loc}=_\stackrel{}{q}\chi (\stackrel{}{q},\omega )|_{\omega =0}`$ diverges at a heavy fermion QCP, From (7),

$$\chi _{loc}(T)d^dq\frac{1}{(𝐪𝐐)^2+T^\alpha }T^{(d2)\alpha /2}$$
(14)
so a divergent local spin susceptibilty requires a spin fluid with $`d2`$. Si et al are thus motivated to propose that the non-trivial physics of the heavy fermion QCP is driven by the formation of a two-dimensional spin fluid. Si et al consider an impurity spin within an effective medium in which the local Weiss field $`H`$ has a critical power-spectrum (fig. 5.)

$$|H(\omega )|^2\chi _0^1(\omega )=\omega ^\gamma $$
(15)
where $`ϵ`$ is self-consistently evaluated using a dynamical mean-field theory, where $`q`$ dependence of self-energies is dropped. In principle, the method solves the dynamical spin susceptibility of the impurity $`\chi ^1(\omega )=\chi _o^1(\omega )+M(\omega )`$. This, in turn furnishes a “spin self-energy” $`M(\omega )`$ used to determine the spin susceptibility of the medium $`\chi ^1(\stackrel{}{q},\omega )=J(\stackrel{}{q})+M(\omega )`$.

Si et al find that a self-consistent solution is obtained for $`ϵ=1`$, if the spin-self energy contains a separate power-law dependence $`M(\omega )\omega ^\alpha `$ with an exponent $`\alpha =1/\rho \mathrm{\Gamma }`$ which is determined by the density of states $`\rho `$ and band-width $`\mathrm{\Lambda }`$ of the bond-strengths in the two-dimensional spin fluid. Although self-consistency requires a new power-law in the spin susceptibility, independent solutions of the impurity model have not yet shown that this feature is indeed generated by a critical Weiss field. This theory nevertheless raises many interesting questions:

1. Is the requirement of a two dimensional spin fluid consistent with experiment? For example- does the the cubic (and hence manifestly three dimensional) quantum critical material, $`CeIn_3`$ display a divergent specific heat?
2. If the spin-fluids are quasi-two dimensional, do we expect an ultimate cross-over to a three-dimensional QSDW scenario?
3. If $`\alpha `$ is non-universal, why are the critical exponents in $`CeCu_{6x}Au_x`$ and $`YbRh_2Si_{2x}Ge_x`$ so similar?
4. What stabilizes the local quantum criticality against intersite couplings?

### 2  Towards the critical Lagrangian.

If the heavy fermion QCP is a truly three-dimensional phenomenon, then a different approach is needed- we need to search for a new class of critical Lagrangian with $`d_u>3`$. On general grounds, the existence of a Fermi liquid in the paramagnetic phase tells us that it must find expression in terms of the quasiparticle fields $`\psi `$ in the Fermi liquid,

$$L=L_F[\psi ]+L_{FM}[\psi ,M]+L_M[M].$$
(16)
where $`L_F`$ describes the heavy Fermi liquid, far from the magnetic instability, $`L_M`$ describes the magnetic excitations that emerge above the energy scale $`T_F^{}(P)`$.

$`L_{FM}`$ describes the way that the quasiparticles couple to and decay into critical magnetic modes; it also determines the type of transformation which takes place in the Fermi surface which takes at the QCP. This last point follows because away from the QCP, magnetic fluctuations can be ignored in the ground-state, so that $`L_M0`$. In the paramagnetic phase, $`M=0`$ so $`L_{FM}0`$, but in the antiferromagnetic phase $`M0`$, i.e.

$`L_{\text{eff}}=\{\begin{array}{cc}L_F^{}[\psi ]& \text{paramagnet}\hfill \\ & \\ L_F^{}[\psi ]+L_{FM}[\psi ,M]& \text{a.f.m.}\hfill \end{array}`$

where the asterix denotes the finite renormalizations derived from zero-point fluctuations in the magnetization.

If the staggered magnetization is the fundamental critical field, then we are forced to couple the magnetic modes directly to the spin density of the Fermi liquid

$`L_{FM}^{(1)}`$ $`=`$ $`g{\displaystyle \underset{𝐤,𝐪}{}}\psi _{\stackrel{}{k}\stackrel{}{q}}^{}\stackrel{}{\sigma }\psi _\stackrel{}{k}\stackrel{}{M}_𝐪.`$ (17)
But once the staggered magnetization condenses, this leads directly back to a static spin density wave (Fig 4).

An alternative possibility is suggested by the observation that the magnetism develops spinorial character in the heavy fermi liquid. The Luttinger sum rule governing the Fermi surface volume $`V_{FS}`$ “ counts” both the electron density $`n_e`$ and the number of the number of local moments per unit cell $`n_{spins}`$ :

$$2\frac{𝒱_{FS}}{(2\pi )^3}=n_e+n_{spins}.$$
(18)
The appearance of the spin density in the Luttinger sum rule reflects the composite nature of the heavy quasiparticles, formed from bound-states between local moments and high energy electron states. Suppose the spinorial character of the magnetic degrees of freedom seen in the paramagnet also  manifests itself in the decay modes of the heavy quasiparticles. This would imply that at the QCP, the staggered magnetization factorizes into a spinorial degree of freedom $`\stackrel{}{M}(x)=z^{}(x)\stackrel{}{\sigma }z(x)`$, where $`z`$ is a two-component spin $`1/2`$ Bose field. “Spinorial magnetism” affords a direct coupling between the magnetic spinor $`z`$ and the heavy electron quasi-particles via an inner product, over the spin indices

$$L_{FM}^{(2)}=g\underset{𝐤,𝐪}{}[z_{𝐤𝐪\sigma }^{}\psi _{𝐤\sigma }\chi _𝐪^{}+\text{H.c}],$$
(19)
where conservation of exchange statistics obliges us to introduce of a spinless charge $`e`$ fermion $`\chi `$. This would imply that the composite heavy electron decays into a neutral “spinon”and a spinless charge e fermion $`e_\sigma ^{}s_\sigma +\chi ^{}`$. From this perspective, the heavy fermion QCP involves spin-charge separation, and the critical Lagrangian is a gauge theory.

To go beyond this general discussion we need to answer various questions:

* Can we connect the appearance of local criticality with an underlying gauge symmetry?
* What is the link between the heavy fermion QCP and metamagnetism?
* Is there a fundamental reason why the heavy fermion QCP exhibits $`H/T`$ scaling (e.g as opposed to $`H^\delta /T`$ scaling) ?
* What feature in the critical Lagrangian can push the upper-critical dimension above 3?

On the experimental front, Hall constant measurements may provide a good way to discern between the spin density wave and composite quasiparticle alternative. In the former, regions around the hot-line do not contribute to the Hall conductivity, and the change in the Hall constant is expected to evolve as the square of the staggered magnetization. By contrast, the composite fermion scenario leads to a much more rapid evolution: depending on whether the density of spinless fermions is finite at the QCP the Hall constant which either jumps or evolves linearly with the staggered magnetization.

$$\mathrm{\Delta }R_H\{\begin{array}{cc}M_𝐐^2,& \hfill \text{(vectorial )}\\ O(1),M_𝐐& \hfill (\text{spinorial})\end{array}$$
(20)
The only available Hall measurement at a QCP to date shows a change in sign takes place in the close vicinity of the QCP in critical $`CeCu_{1x}Au_x`$, it is not yet clear whether there is a discontinuity at the transition. This is clearly an area where more experimental input is highly desireable.

## F Summary

This paper has discussed the origin of the mass divergence at the heavy fermion quantum critical point, emphasizing that a quantum spin density wave picture can not explain the observed properties. The proposal of fundamentally new kinds of quantum critical point have been reviewed. This is clearly an area with a huge potential for progress both on the experimental, and theoretical front.