text:
# Implications of the discrepancy between proton form factor measurements

## I INTRODUCTION

The proton electromagnetic form factors, $`G_E`$ and $`G_M`$, parametrize deviations from a point particle in elastic electron-proton scattering, and are related to the charge and magnetization distribution of the proton. The form factors depend only on $`Q^2`$, the square of the four-momentum transfer, and until recently it was believed that the electric and magnetic form factors showed approximate scaling, i.e. nearly identical $`Q^2`$ dependence Walker et al. (1994). More recent Jefferson Lab measurements Jones et al. (2000); Gayou et al. (2001, 2002) utilized the polarization transfer technique to measure the ratio $`G_E/G_M`$ and found that $`G_E`$ decreases more rapidly than $`G_M`$ at large $`Q^2`$. The polarization transfer measurements are more precise at high $`Q^2`$, and significantly less sensitive to systematic uncertainties than the Rosenbluth separation measurements. However, the two techniques disagree significantly even in the region where both yield precise results.

At the present time, it is not known why the techniques give different results. The systematic uncertainties of the polarization transfer measurements, primarily spin transport and backgrounds, have been carefully studied Punjabi et al. . A detailed global analysis of the cross section measurements Arrington (2003) does not show any inconsistencies in the cross section datasets, or yield any likely candidate to explain the discrepancy. To resolve the discrepancy, a systematic error in the cross section would have to have a significant dependence on the virtual photon polarization $`ϵ`$, $`ϵ^1=1+2(1+Q^2/4M_p^2)\mathrm{tan}^2(\theta _e/2)`$, where $`M_p`$ is the proton mass and $`\theta _e`$ is the electron scattering angle. Such a systematic error would have to yield a (5–7)% $`ϵ`$ dependence in the cross section, roughly linear in $`ϵ`$, in order to resolve the discrepancy.

There appear to be two possibilities: either a fundamental flaw in the Rosenbluth or polarization transfer formalism, or an error in either the cross section or polarization transfer measurements. Recent works have suggested that additional radiative correction terms, related to two-photon exchange corrections, may lead to an error in determining the form factors from the measured cross sections Blunden et al. (2003); Guichon and Vanderhaeghen (2003); Rekalo and Tomasi-Gustafsson (2003). If the two-photon exchange mechanism, or some other correction that is neglected in the cross section extraction, is the source of the discrepancy, then the form factors extracted from a Rosenbluth separation of cross section data will not represent the underlying structure of the proton, but they will parametrize the elastic electron-proton cross section in the usual one-photon approximation. Conversely, the true form factors will not yield the correct cross sections, and will thus give incorrect results if used as a parameterization of the elastic cross section in data analysis.

If the two-photon exchange term explains the discrepancy, then the polarization transfer result will relate to the true form factors, assuming that the two-photon exchange has a much smaller effect on the polarization transfer than on the Rosenbluth extractions. However, the existing polarization transfer experiments Jones et al. (2000); Gayou et al. (2001, 2002) have extracted the ratio $`G_E/G_M`$, rather than the individual form factors. To extract the form factors, these data must be combined with cross section measurements to determine the absolute magnitudes of $`G_E`$ and $`G_M`$. If the two-photon exchange correction modifies the cross sections from those calculated from the underlying form factors, then it is not possible to consistently combine the two kinds of measurements without some assumption about the two-photon exchange correction.

In this paper, we present two extractions of the proton form factors. From a global analysis of cross section measurements, we extract the “Rosenbluth form factors”. From a combined analysis of cross section and polarization transfer data, with a “minimal” assumption about the nature of the two-photon exchange corrections, we extract the “Polarization form factors”. If two photon corrections are the source of the discrepancy, then the Rosenbluth form factors will parametrize the elastic cross section, and are therefore useful as input to analysis or simulations that require the electron-proton cross section. The Polarization form factors will provide the true form factors, which relate to the underlying structure of the proton. These form factors are often described as the Fourier transformations of the charge and magnetization distributions of the proton in the Breit frame, although relativistic effects and the fact that each value of $`Q^2`$ corresponds to a different Breit frame lead to substantial theoretical difficulties in extracting charge and magnetization distributions Kelly (2002).

## II GLOBAL ROSENBLUTH ANALYSIS

The “Rosenbluth form factors” are determined from a global fit to elastic electron-proton cross section measurements. The details of the fitting procedure are described in Ref. Arrington (2003). For the present analysis we include more recent Jefferson Lab measurements of elastic scattering Christy et al. ; Dutta et al. (2003); Niculescu (1999), as well additional datasets to constrain the low $`Q^2`$ behavior Borkowski et al. (1974); Murphy et al. (1974); Simon et al. (1980, 1981) to the datasets used in Ref. Arrington (2003); not . In addition, we include all of the high $`Q^2`$ data, up to 30 GeV2, while the previous analysis was limited to 8 GeV2. The older data have updated radiative corrections, and the small-angle data from Walker et al. Walker et al. (1994) are excluded, as described in Arrington (2003). The form factors are fit to the following form:

$$G_E(Q^2),G_M(Q^2)/\mu _p=[1+p_2Q^2+p_4Q^4+\mathrm{}+p_{2N}Q^{2N}]^1,$$
(1)
where $`\mu _p`$ is the magnetic dipole moment of the proton and $`Q^2`$ values are in GeV2. Reasonable fits are achieved for $`N3`$. Note that this is a different functional form than used in previous fits Bosted (1994); Arrington (2003); Brash et al. (2002), which used polynomials in $`q=\sqrt{Q^2}`$. The polynomial in $`q`$ is a very general form, with adequate flexibility to reproduce the data, but does not have the proper behavior as $`Q^20`$.

The fit is quite insensitive to the order of the polynomial above $`N=6`$, except for $`G_E`$ at large $`Q^2`$. For $`Q^2`$ above 6 GeV2, fits with nearly identical $`\chi ^2`$ values can have $`G_E/G_M`$ either rise or fall dramatically with $`Q^2`$. This is a result of the reduced sensitivity to $`G_E`$ and the limited $`ϵ`$ coverage for $`Q^2`$ values above 6 GeV2. To avoid unreasonable behavior in the region where $`G_E`$ is unconstrained by data, we keep the ratio $`G_E/G_M`$ fixed for all $`Q^2`$ values above 6 GeV2. This leads to a fit for $`G_E`$ which is continuous, but not smooth, at $`Q^2=6`$ GeV2. Because $`G_E`$ has relatively little contribution to the total cross section at these momentum transfers, the cross section extracted is still quite smooth, and the value of $`G_E`$ at large $`Q^2`$ values has little effect on the cross section, as long as the fit is constrained to avoid $`|\mu _pG_E||G_M|`$.

The normalization factor for each data set is allowed to vary along with the parameters of the fitting functions for $`G_E`$ and $`G_M`$. The total $`\chi ^2`$ from the cross section measurements and normalization factors is:

$$\chi _\sigma ^2=\underset{i=1}{\overset{N_\sigma }{}}\frac{(\sigma _i\sigma _{\mathrm{fit}})^2}{(d\sigma _i)^2}+\underset{j=1}{\overset{N_{\mathrm{exp}}}{}}\frac{(\eta _j1)^2}{(d\eta _j)^2},$$
(2)
where $`\sigma _i`$ and $`d\sigma _i`$ are the cross section and error (excluding normalization uncertainties) for each of the $`N_\sigma `$ data points, $`\eta _j`$ is the fitted normalization factor for the $`j`$th dataset, and $`d\eta _j`$ is the normalization uncertainty for that dataset. We fit to 470 data points ($`N_\sigma `$=443, $`N_{exp}`$=27) with 39 parameters (six parameters each for the electric and magnetic form factors, and 27 normalization parameters).

The result of the global fit to the cross section data is shown in Fig. 1. The fit yields a total $`\chi ^2`$ of 326.7 for 431 degrees of freedom, yielding a reduced $`\chi ^2`$, $`\chi _\nu ^2=\chi ^2/N_{dof}`$, of 0.758. This yields an unreasonably high confidence level, indicating that the quoted uncertainties of the measurements are too large. As was observed in the previous fit Arrington (2003), the majority of the datasets, 20 out of 27, have values of $`\chi _\nu ^2<1`$, indicating that most of the experiments were overly conservative in estimating their uncertainties. Table 1 lists the parameters for the Rosenbluth form factors. The fit includes cross sections for $`Q^2`$ values from 0.005 to 30 GeV2, and should be valid over this range, though the separation of $`G_E`$ and $`G_M`$ is only well constrained by the data for $`Q^2<6`$ GeV2.

The normalization factors were generally smaller than the quoted scale uncertainties of the experiments ($`\chi ^2`$=18.0 for 27 normalization factors). The average normalization factor is 0.65%, and the RMS normalization factor is 2.7%. The normalization factors are very close to those obtained in the previous global fit Arrington (2003). The average normalization factor differs by approximately 0.5%, and the individual normalization factors differ by less then 1% for 18 of the 20 experiments. Because the previous fit excluded data below $`Q^2`$=0.6 GeV2, the agreement indicates that it is the self-consistency, rather than the form factor constraint at $`Q^2=0`$, that dominates the determination of the normalization factors.

We can test the self-consistency of the individual data sets by comparing the global fit to the results of single-experiment extractions of $`G_E`$ and $`G_M`$. By comparing only the single-experiment extractions, we avoid the potentially large and correlated uncertainties that arise from the relative normalization of different datasets. Comparing the ratio $`G_E/G_M`$ from the fit to the individual experiments, taken from Refs. Christy et al. ; Dutta et al. (2003) and the reanalysis of older experiments presented in Ref. Arrington (2003) yields $`\chi ^2=45.3`$ for 50 data points ($`\chi ^2=17.8`$ for the 20 points above $`Q^2=1.5`$ GeV2).

We can estimate the uncertainties in the form factors by performing direct Rosenbluth separations in several $`Q^2`$ bins using the full dataset, with normalization factors determined from the global fit. For each $`Q^2`$ bin, the data are scaled to the average $`Q^2`$ value of the data points in that bin, using the global fit as the scaling function. $`Q^2`$ bins were chosen so that there are at least three data points in the bin, the $`ϵ`$ range covered is at least 0.3, and the correction for scaling each point to the average $`Q^2`$ value was $`<`$10% (typically $`<2`$%). The scaling was also done using the fits of Refs. Bosted (1994) and Brash et al. (2002), shown in Fig. 1. Varying the scaling procedure changed the ratios by $`1`$%, except for the very highest (lowest) $`Q^2`$ points, where the change in $`G_E`$ ($`G_M`$) was as much as 3%, but was still much smaller than the uncertainty in the extracted form factor.

Figure 2 shows the fits to $`G_E`$ and $`G_M`$, along with the direct Rosenbluth separation points, using the normalization factors from the fit. Except for the very low $`Q^2`$ values, typical uncertainties on $`G_M`$ are $``$1%, increasing to $``$2% for $`Q^2`$=10 GeV2 (8% for $`Q^2`$=30 GeV2). At low $`Q^2`$, the experimental uncertainties become quite large, but the constraint on the behavior as $`Q^20`$ yields a much smaller uncertainty on the fit. For $`G_E`$, the uncertainties are (1–2)% at low $`Q^2`$, but are (5–10)% for intermediate $`Q^2`$ values (2–4 GeV2), and grow rapidly as $`Q^2`$ increases. Note that the uncertainties in $`G_E`$ and $`G_M`$ are highly anti-correlated, due to the way the form factors are separated from the cross section measurements. This can be seen in the anti-correlation of the deviation of the points from the fits in Fig. 2. Thus, the uncertainty on the cross sections extracted from this parametrization is not just the sum of the uncertainties in the contributions from $`G_E`$ and $`G_M`$. Up to $`Q^2`$$``$4 GeV2, there is a large body of cross section measurements with point-to-point uncertainties of $``$1%. Because the normalization factors are determined in the fit, and the residual uncertainty in the normalization is small, the absolute cross sections should be known to better than 2%. Above $`Q^2`$=4 GeV2, the number of data points decreases, and the uncertainties in the cross sections grow, reaching 10% at $`Q^2`$=25 GeV2.

Even with the uncertainty related to the discrepancy between Rosenbluth and polarization transfer, this fit yields a precise parametrization of the elastic cross section in the one-photon exchange formalism. While these may not be the underlying form factors of the proton (e.g., if there are missing radiative correction terms), this is still the appropriate parametrization to use as input to a calculation or analysis that requires the elastic cross section. Using the form factors derived from the polarization transfer technique will not yield the correct cross section, even in a combined analysis of Rosenbluth and polarization transfer such as performed in Refs. Brash et al. (2002); Arrington (2003). More importantly, an inconsistent combination of cross section and polarization transfer results can magnify the error. Combining a parametrization of $`G_M`$ from a Rosenbluth analysis with the form factor ratios measured in polarization transfer decreases $`G_E`$, and thus decreases the total cross section, relative to the best fit to the cross section data, without allowing a corresponding increase in $`G_M`$. This leads to form factors which give cross sections that are (4–10)% below the measured cross sections at large $`ϵ`$ over a large $`Q^2`$ range ($`0.1<Q^2<15`$ GeV2).

## III COMBINED POLARIZATION AND ROSENBLUTH ANALYSIS

While the Rosenbluth form factors yield the best parametrization for the cross section in the usual one-photon exchange picture, the ratio does not agree with the ratio extracted from the polarization transfer technique. For the larger $`Q^2`$ values, the polarization transfer technique is less sensitive to knowledge of the kinematics, radiative correction, and other systematic uncertainties that are important in the Rosenbluth separation.

If this discrepancy is related to a problem in the cross section data, then the polarization transfer will yield the true ratio of the form factors, but has to be combined with cross section data to obtain both $`G_E`$ and $`G_M`$. We present in this section a combined analysis of the polarization transfer and cross section data, which will yield the “Polarization form factors”.

In order to obtain a consistent extraction of the form factors, we must make an assumption about the nature of the discrepancy. We assume that the difference comes from a common systematic error in the cross section measurements. Analyses of this discrepancy Arrington (2003); Guichon and Vanderhaeghen (2003) indicate that there must be an $`ϵ`$-dependent correction of (5–7)%, roughly linear in $`ϵ`$, for $`1<Q^2<6`$ GeV2.

In the combined analysis, we apply a linear, $`ϵ`$-dependent correction of 6%, to all datasets. This is the ‘minimal’ assumption necessary to make the two techniques consistent, to the extent that a correction that was not linear in $`ϵ`$, or which modified only some of the datasets, would have to be larger. A correction that is nearly linear in $`ϵ`$ and fairly $`Q^2`$-independent is consistent with the form for the two-photon exchange term in the analysis of Ref. Blunden et al. (2003), although the size of the correction in Ref. Blunden et al. (2003) is only $``$2%, less than half the size necessary to explain the discrepancy.

We repeat the fit from Sec. II, but with cross sections modified by the linear $`ϵ`$ dependence, and with the polarization transfer data included in the fit, as described in Arrington (2003). The correction to the cross section could either lower the cross section at large $`ϵ`$ values, or increase it at small $`ϵ`$ values:

$`\sigma _{c1}=\sigma _0(10.06ϵ)`$ (3)
$`\sigma _{c2}=\sigma _0(10.06(ϵ1))=0.94\sigma _{c1}`$ (4)
The first correction is consistent with the form from Ref. Guichon and Vanderhaeghen (2003), while the second is consistent with the behavior of Blunden et al. (2003). The second form was chosen for the main fit because the correction is small at large $`ϵ`$ (small $`\theta _e`$), where comparisons of positron to electron scattering from SLAC Mar et al. (1968) set fairly tight limits on the size of two-photon exchange.

The Polarization form factors, from the combined fit to the cross section and the 26 polarization transfer data points from Refs. Jones et al. (2000); Gayou et al. (2001, 2002) is shown in Fig. 3. The fit yields a total $`\chi ^2`$ of 391.6 for 457 degrees of freedom, $`\chi _\nu ^2=0.857`$, including the additional $`\chi ^2`$ contribution for the polarization transfer data (Eq. (8) of Ref. Arrington (2003)). Table 2 gives the fit parameters for the Polarization form factors.

The fit was also performed with the correction of Eq. 3. This leads to an overall rescaling of all of the cross sections by 6%, relative to the correction of Eq. 4. However, this does not yield a simple rescaling of the form factors, because each dataset has a normalization factor that is determined in the fit, and because the form factors are constrained to reproduce the charge and magnetic moment at $`Q^2=0`$. While a two-photon correction of this size for large $`ϵ`$ values would appear to be ruled out by the SLAC positron-proton measurements Mar et al. (1968), an $`ϵ`$-dependent systematic other than two-photon exchange could also resolve the discrepancy. However, this fit yields a much worse $`\chi ^2`$ value: 575.1 for 457 degrees of freedom, and so we choose to apply Eq. 4 for our combined fit.

Note that the result of the combined fit (Table 2) will not reproduce the measured elastic cross section in the one-photon exchange formalism; it will reproduce the modified cross sections of Eq. 4. Therefore, the Polarization form factors should not be used to model elastic electron-proton cross section measurements. However, if the ‘minimal’ assumption that a correction consistent with the form of Eq. 4 explains the discrepancy, this should yield a consistent extraction of the underlying form factors of the proton.

## IV CONCLUSIONS

Form factors extracted using the Rosenbluth technique provide a parametrization of the deviation of the elastic electron-proton cross section from the point-scattering cross section. If the cross section has additional corrections, such as two-photon exchange terms, that are not being taken into account, then the Rosenbluth extraction does not yield the true proton form factors that relate to the structure of the proton. In this case, $`G_E`$ must be extracted from the polarization transfer measurements, which yield $`G_E/G_M`$, and the cross section data must be utilized to determine $`G_M`$. While we cannot know how to properly combine the polarization transfer and cross section data until we understand the cause of the discrepancy, the uncertainties in $`G_M`$ that arise from this problem are much smaller than those in $`G_E`$. The same holds true if there is some other correction or combination of corrections to the cross section other than the two-photon exchange (e.g. Coulomb corrections Higinbotham (2003)). It is of course possible that the discrepancy is due to a problem with the polarization transfer data or technique rather than the cross section data. If so, then the Rosenbluth form factors represent both the correct cross section and the correct nucleon structure. However, there do not appear to be any obvious candidates for problems in the technique, and the experiment should be less prone to systematic uncertainties than the Rosenbluth extractions.

We have presented two extractions of the proton electromagnetic form factors. The Rosenbluth form factors come from a global Rosenbluth extraction of the form factors from electron-proton elastic scattering measurements. The Polarization form factors come from a combined fit to the cross section and polarization transfer data, under the assumption that the discrepancy between the techniques is caused by a linear, $`ϵ`$-dependent correction to the cross sections. The Rosenbluth form factors give a global parametrization of the elastic electron-proton scattering cross section in the one-photon exchange approximation. Even if there is a correction to the cross sections, neglected in the one-photon exchange formalism, this parametrization will yield the correct cross sections in the one-photon approach. Under the above assumption of an unknown correction to the cross sections, the Polarization form factors yield the underlying form factors, but will not reproduce cross sections, and will therefore yield incorrect results if used as input for an analysis that requires the elastic cross section, as was observed in an analysis of quasielastic scattering from nuclei Dutta et al. (2003).

Additional data will help shed light on the origin of the discrepancy. An improved “Super-Rosenbluth” separation measurement Arrington et al.  completed at Jefferson Lab in 2002 will yield a precise extraction of $`G_E/G_M`$, and will determine if the discrepancy can be explained by experimental problems in the Rosenbluth extractions. A new polarization transfer experiment Perdrisat et al. , approved to run at Jefferson Lab, will provide an independent confirmation of the existing polarization transfer results, as well as extending the measurements to higher $`Q^2`$ values. If sufficiently improved calculations or direct measurements of the two-photon exchange corrections become available, we should be able to determine if they are responsible for the discrepancy, and if so, remove the current uncertainty in combining cross section and polarization transfer measurements.

If the discrepancy is explained by two-photon corrections or some other effect on the cross sections, and we have reliable calculations for these effects, then the cross section data can be combined with the polarization transfer data to extract the form factors without ambiguity. These form factors will represent the underlying structure of the proton and provide a useful parametrization of the elastic electron-proton cross section, as long as the effect is properly accounted for. Until the discrepancy is well understood, however, both sets of form factors are necessary, and it is important to use form factors that are (1) extracted consistently from the cross section and/or polarization transfer data, and (2) appropriate for the problem being addressed.

###### Acknowledgements.

This work was supported by the U. S. Department of Energy, Nuclear Physics Division, under contract W-31-109-ENG-38.


format: text

text:
# 𝑑-Mott phases in one and two dimensions

## Abstract

We use exact diagonalization to determine the spectrum of reduced Hamiltonians based on renormalization group flows to strong coupling. For the half-filled two-leg Hubbard ladder we reproduce the known insulating $`d`$-Mott groundstate with spin and charge gaps. For the saddle point regions of the two-dimensional Hubbard model near half-filling we find a crossover to a similar strong coupling state, which truncates the Fermi surface near the saddle points. At lower scales $`d`$-wave superconductivity appears on the remaining Fermi surface.

In recent years the application of renormalization group (RG) methods to interacting fermions in low dimensions has made great progress. In one dimension (1D) where the Fermi surface (FS) is a discrete set of points, the method has a long history going back to the pioneering work of Menyhárd and Sólyom solyom . In two dimensions (2D) the FS is a continuous curve and the discrete set of coupling constants is replaced by a 4-point vertex which is a function of three independent momenta. Zanchi and Schulz zanchi  developed a numerical method to evaluate the functional renormalization group flow of the scattering vertex based on a break-up of the FS into $`N`$ patches. The resulting coupled differential equations can be integrated numerically to obtain the RG flow. This method has been applied to the Hubbard model on a square lattice with nearest-neighbor (n.n.) and next-nearest neighbor (n.n.n.) hopping using a Wilsonian scheme which successively integrates out high energy modes down to an infrared cutoff zanchi ; halboth ; honerkamp . The RG method has the advantage that it allows an unbiased treatment of competing instabilities including even Stoner ferromagnetism, when the recent reformulationtflow , with the temperature as the flow parameter, is used.

A serious limitation to RG schemes occurs when the scattering vertices flow to strong coupling as the RG scale is lowered. Although the one-loop flow may still be a good approximation even when the interaction strength has become comparable to the kinetic energy, self energy corrections honerkamp  typically become important shortly after that and higher order terms or non-perturbative effects might not be negligible any longer. Thus these RG flows cannot be extrapolated to scales below the divergence. Generally those are the most interesting cases as divergent RG flows signal a breakdown of the perturbative normal state. In the special case of 1D, bosonization allows a complete analysis of the strong coupling regime. In higher dimensions, in cases where a single instability dominates the RG flow, a mean field approximation can be introduced to describe the strong coupling phase. This can be corrected for fluctuations around the mean field to obtain a satisfactory description of the ordered phase. Yet in some especially interesting cases, there can be several mutually reinforcing instabilities so that it is not clear how to formulate a mean field approximation.

In this Letter we propose a new and more general approach to determine the strong coupling phase that emerges on scales below divergent RG flows. The starting point is to recognize that a mean field approximation based on factorization of the interaction term in the Hamiltonian, can also be expressed as an exact solution of a reduced Hamiltonian e.g. BCS theory was formulated as an exact solution for a reduced Hamiltonian in which only terms describing the scattering of Cooper pairs with zero total momentum were kept BCS . Note these are the most divergent terms in the RG flow of the 4-point vertex. Thus we can regard the mean field approximation as contained in the solution of a reduced Hamiltonian whose form is derived from the diverging RG flow. In those cases mentioned above the reduced Hamiltonian is more complex and cannot be simply solved analytically. However it is amenable to numerical diagonalization and this is the approach we take here. In particular in the study of a n.n.n. Hubbard model on a square lattice as a prototypical model for underdoped cuprates, Furukawa et al.furukawa , and Honerkamp et al.honerkamp , found a mutual reinforcement of antiferromagnetic (AF) and $`d`$-wave pairing tendencies when the Fermi energy was close to the van Hove singularities at the saddle points (SP) at $`(\pm \pi ,0),(0,\pm \pi )`$. They pointed out that the same mechanism is at work in the RG flows in the half-filled 2-leg Hubbard ladder (2LHL). At the heart of the mutual reinforcement is the geometrical fact that in these cases a significant fraction of particle pairs that experience strong Umklapp scattering in the $`(\pi ,\pi )`$ channel, e.g. two particles close to the SP $`(\pi ,0)`$, have small total momentum and thus couple into the Cooper channel as well. Since both the AF $`(\pi ,\pi )`$-scattering and the large momentum transfer pair scattering in the $`d_{x^2y^2}`$-channel are repulsive, a mutual reinforcement occurs. In the half-filled 2LHL, this works down to arbitrarily low scales. In 2D, the overlap between the channels decreases for low scales such that a finite interaction strength is needed for this mechanism. The strong coupling behavior of the ladder has been fully analyzed using Density Matrix Renormalization Group (DMRG)noack  and by bosonization by Lin et al.fisher  who named the groundstate a $`d`$-Mott phase. It shows the unusual behavior of enhanced correlations simultaneously in the AF, $`d`$-wave pairing and $`d`$-density wave (or orbital AF) channels and is intimately related to resonating-valence bond (dRVB) states at large $`U`$. Note all correlations remain strictly short range due to the presence of both spin and charge gaps. Clearly this is a case of a nonperturbative groundstate which cannot be even qualitatively described by a mean field approximation based on long range order.

We analyze first the case of the 2LHL which is a good test case since the groundstate is highly nontrivial and turn then to the n.n.n. Hubbard model on a square lattice. The limitation of Lanczos algorithms to a relatively small number of $`\stackrel{}{k}`$-points means that in 2D a representation of the whole FS is impossible. However in the RG flows of the 2D n.n.n. Hubbard model, the dominant divergent scattering processes appear in the vicinity of the SPs. This allows us to restrict our reduced model, derived from the RG flow, to processes connecting a discrete set of $`\stackrel{}{k}`$ points located near the two SPs. If this discrete set of $`\stackrel{}{k}`$-points is restricted to just the two SPs, then the RG flow reduces to a small number of coupled differential equations. Two different fixed points were found by Dzyaloshinskiidzyaloshinskii  and by Lederer et al.lederer . The latter found a divergent flow to strong coupling when the next to leading order particle-hole processes are included. The 2D RG flows have a similar form for the dominant processes connecting the SP regions. Our numerical results for this case reveal both spin and charge gaps at the saddle points and show enhanced correlations in both AF and $`d`$-wave pairing channels. They confirm the similarity to the $`d`$-Mott phase of the 2LHL, in line with the earlier interpretations.

Let us briefly describe the numerical approach. An extended discussion and results for the 2LHL and other 1D systems will be published elsewhere amlOnMethod . Our scheme is based on the numerical analysis of the reduced Hamiltonian $`H_\mathrm{\Lambda }`$ obtained from the RG flow close to the critical scale $`\mathrm{\Lambda }_c`$. $`H_\mathrm{\Lambda }`$, acting on states in a $`\mathrm{\Lambda }`$-shell around the Fermi points or surface, is discretized using a finite number of $`\stackrel{}{k}`$-points (illustrated in Fig. 1 for two Fermi points). Generically it has the form

$`H_\mathrm{\Lambda }`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{k},\sigma }{}}\epsilon (\stackrel{}{k})c_{\stackrel{}{k},\sigma }^{}c_{\stackrel{}{k},\sigma }`$
$`+{\displaystyle \frac{\lambda }{2\mathrm{\Omega }}}{\displaystyle \underset{\genfrac{}{}{0pt}{}{\stackrel{}{k}_1,\stackrel{}{k}_2,\stackrel{}{k}_3}{\sigma ,\sigma ^{}}}{}}V_\mathrm{\Lambda }(\stackrel{}{k}_1,\stackrel{}{k}_2,\stackrel{}{k}_3)c_{\stackrel{}{k}_3,\sigma }^{}c_{\stackrel{}{k}_4,\sigma ^{}}^{}c_{\stackrel{}{k}_2,\sigma ^{}}c_{\stackrel{}{k}_1,\sigma }`$
where $`\epsilon (\stackrel{}{k})`$ denotes the kinetic energy, $`\lambda `$ is a global coupling constant, $`\mathrm{\Omega }`$ the total volume, $`V_\mathrm{\Lambda }(\stackrel{}{k}_1,\stackrel{}{k}_2,\stackrel{}{k}_3)`$ the discretized coupling function and $`\stackrel{}{k}_4=\stackrel{}{k}_1+\stackrel{}{k}_2\stackrel{}{k}_3`$ (modulo umklapp) due to momentum conservation. The momentum dependence of the interaction function $`V_\mathrm{\Lambda }(\stackrel{}{k}_1,\stackrel{}{k}_2,\stackrel{}{k}_3)`$ is determined by the RG couplings at scale $`\mathrm{\Lambda }`$:

$$V_\mathrm{\Lambda }(\stackrel{}{k}_1,\stackrel{}{k}_2,\stackrel{}{k}_3)=\stackrel{~}{g}[\mathrm{\Lambda }](\text{Patch}(\stackrel{}{k}_1),\text{Patch}(\stackrel{}{k}_2),\text{Patch}(\stackrel{}{k}_3)).$$

Here $`\stackrel{~}{g}[\mathrm{\Lambda }]`$ denotes the ratios of the diverging couplings and $`\text{Patch}(\stackrel{}{k})`$ is the patch index of $`\stackrel{}{k}`$-point $`\stackrel{}{k}`$. Note that in the reduced BCS Hamiltonian, only interactions $`V_\mathrm{\Lambda }(\stackrel{}{k}_1,\stackrel{}{k}_2,\stackrel{}{k}_3)`$ with $`\stackrel{}{k}+\stackrel{}{k}_2=0`$ are kept. In contrast to this, here in order to describe the short-range $`d`$-Mott state, a whole set of scattering processes with nonzero $`\stackrel{}{k}+\stackrel{}{k}_20`$ must be kept in $`H_\mathrm{\Lambda }`$. An analogous width for $`V_\mathrm{\Lambda }(\stackrel{}{k}_1,\stackrel{}{k}_2,\stackrel{}{k}_3)`$ for momentum transfers $`(\pi ,\pi )`$ is necessary. To incorporate this essential feature we assume that all scattering amplitudes with the same patch indices take the same value. The global coupling constant $`\lambda `$ is in principle determined by the the initial conditions and the final scale $`\mathrm{\Lambda }`$. Here however we use it as a parameter in order to investigate the interplay between kinetic energy and scattering processes and to analyze the effects of the finite number of $`\stackrel{}{k}`$-points.

In the next step the discretized Hamiltonian is diagonalized using a Lanczos algorithm. Due to momentum conservation (modulo umklapp scattering) a number $`N_k20`$ of orbitals can be treated. In this way energy gaps, static and dynamical correlation functions are easily accessible, similar to standard real space diagonalization algorithms. As a future development the Lanczos scheme might be replaced by a DMRG algorithm, thus possibly allowing more $`\stackrel{}{k}`$-points for a higher resolution on the FS.

We illustrate the method with an application to the two leg Hubbard ladder at half filling. We choose $`t_{}<2t_{}`$ in order to have 4 Fermi points. The RG equations have been derived and analyzed in fisher . In particular for general repulsive initial interactions the couplings flow towards a fixed ray, i.e. a reduced Hamiltonian with fixed ratios among the diverging scattering vertices.

We characterize the groundstate of the half-filled 2LHL first by calculating different gaps: the spin gap: $`\mathrm{\Delta }_{\text{Spin}}=E_0(N_e,S^z=1)E_0(N_e,S^z=0)`$, the single particle gap: $`\mathrm{\Delta }_{\text{1P}}=[E_0(N_e1,1/2)+E_0(N_e+1,1/2)]/2E_0(N_e,0)`$ and the two particle gap: $`\mathrm{\Delta }_{\text{2P}}=[E_0(N_e2,0)+E_0(N_e+2,0)]/2E_0(N_e,0)`$. Results are shown in the left panel of Fig. 2 for different system sizes and choices of $`\lambda `$. One nicely sees that all gaps remain finite, even for $`N_k\mathrm{}`$. Interestingly the spin gap and the two particle gap are exactly degenerate, while the single particle gap is slightly larger. This can be understood as an approximate realization of the $`SO(8)`$ symmetry of the fixed ray, discovered in Ref. fisher . In our discretized scheme only a $`SO(5)`$ symmetry seems to be manifest.

Now we consider the correlation functions. We investigate particle-hole and particle-particle response with both $`s`$-wave and $`d`$-wave formfactors. The $`\stackrel{}{Q}=(\pi ,\pi )`$ charge density wave structure factors are defined as follows: $`𝒮_{CDW}=\frac{1}{N_k}\left|_{\stackrel{}{k},\sigma }f(\stackrel{}{k})c_{\stackrel{}{k}+\stackrel{}{Q},\sigma }^{}c_{\stackrel{}{k},\sigma }|\mathrm{\Psi }\right|^2,`$ where $`f(\stackrel{}{k})=1`$ for the s-wave channel and $`f(\stackrel{}{k})=\mathrm{exp}(ik_y)`$ for the $`d`$-wave channel. The structure factors for the spin density wave and the pairing correlations are defined similarly. The results for the case $`\lambda /\mathrm{\Lambda }=2`$ are shown in the right panel of Fig. 2, normalized to the respective expectation value in the noninteracting $`N_\stackrel{}{k}=8`$ groundstate. The response for the standard $`(\pi ,\pi )`$ spin density wave, the $`(\pi ,\pi )`$ $`d`$-charge density wave (orbital AF) and the $`d`$wave pairing correlations are equally enhanced. However the finite size scaling does not support long range order, in complete agreement with the results obtained by bosonization fisher .

The $`N`$-patch RG analysis for the n.n.n. 2D Hubbard model revealed a saddle point regime honerkamp  where the FS is close to the SPs but sufficiently curved such that the Brillouin zone diagonals are not nested. Then the leading flow is given by the scattering processes between the two SP regions. The rest of the FS, which we will call arcs, contributes less to the flow and will be considered later. First let us focus exclusively on the SP regions. Then the flow between the SPs is the one described by Lederer et al.lederer  and Furukawa et al.furukawa . There are four different types of processes, $`g_1,\mathrm{}g_4`$, depicted in Fig. 3 (a). In principle the scattering vertex $`V_\mathrm{\Lambda }(\stackrel{}{k}_1,\stackrel{}{k}_2,\stackrel{}{k}_3)`$ depends on the precise location of the wave vectors $`\stackrel{}{k}_i`$ in the SP regions. As a simplification we assume that they have the same value. To a certain extent this can be justified by the N-patch RG. Note however that in the one-loop flow the couplings change smoothly from arc to SP regions.

We implement the two patch model using the two discretizations shown in Fig. 3 (b). From the RG approachfurukawa  the ratios of the couplings close to the critical scale are known: $`g_1:g_2:g_3:g_4=0:1:2.2:1`$. Using this we calculate energy gaps and correlation functions in the same way as for the 2LHL case discussed above. The gap structure and structure factors are shown in Fig. 4. The sSDW and dSC channels are enhanced in agreement with the divergences found in the 1-loop RG flow furukawa  and the dCDW is also enhanced. The results show a striking analogy to the $`d`$-Mott phase of the half-filled 2LHL. This clearly supports the idea that the initial divergence leads to the truncation of the FS with spin and charge gaps opening near the SP. We examine next the behavior of the remaining arcs of the non-nested FS as the RG scale is lowered further. These arcs couple to the SP regions in the Cooper channel.

Although there is a charge gap for particle pairs $`s_{\stackrel{}{p},s}`$, $`s_{\stackrel{}{p},s^{}}`$ in the SP regions, quasiparticle pairs with zero total momentum $`a_{\stackrel{}{k},s}`$, $`a_{\stackrel{}{k},s^{}}`$ from the arcs can be scattered there as virtual excitations. The coupling can in principle be derived from the $`N`$-patch RG, its action is ($`k=(i\omega ,\stackrel{}{k})`$)

$$S_{AS}=\frac{V_{AS}T^3}{2\mathrm{\Omega }}\underset{k,p}{}g_\stackrel{}{k}\overline{a}_{k,s}\overline{a}_{k,s^{}}g_\stackrel{}{p}s_{p,s^{}}s_{p,s}.$$

Here we focus on the $`d_{x^2y^2}`$ channel, $`g_\stackrel{}{k}=\mathrm{cos}k_x\mathrm{cos}k_y`$. Formally integrating out the SP regions gives rise to an induced pairing interaction on the arcs

$$S_{AA}^{\mathrm{ind}}=\frac{V_{AS}^2\chi _S^dT^3}{2\mathrm{\Omega }}\underset{\genfrac{}{}{0pt}{}{k,k^{}}{s,s^{}}}{}g_\stackrel{}{k}\overline{a}_{k,s}\overline{a}_{k,s^{}}g_\stackrel{}{k}^{}a_{k^{},s^{}}a_{k^{},s}.$$

Here, $`\chi _S^d`$ is the static $`d`$-wave pairing susceptibility per volume of the SP regions. $`\chi _S^d`$ remains finite at all temperatures $`T`$, as the pair correlations in the $`d`$-Mott state decay exponentially. However the induced attraction enhances the weak $`d`$-wave attraction within the arcs, and eventually, at low enough temperatures the arcs will undergo a superconducting transition. This mechanism resembles that proposed by Geshkenbein et al. geshkenbein  with preformed pairs at the saddle points. As the SP regions are insulating, the superfluid weight will be that of the arcs alone. With a very similar reasoning for the $`(\pi ,\pi )`$-processes in the spin channel we can understand how the system develops AF order when the arcs become more nested at sufficiently small doping.

We may further ask what criterion determines whether the superfluid density is reduced in the groundstate. The answer lies in the strength of the mutual reinforcement of Cooper and $`(\pi ,\pi )`$-channels at the SPs. If the $`(\pi ,\pi )`$-’almost-Cooper’ processes enhance the $`d`$-wave pair scattering at the SPs more than the coupling in the Cooper channel to the arcs, we can expect that the system favors the reduced superfluid density. In this situation the energy gain through the mutual reinforcement - manifest in spin and charge gap - is larger than the loss of superconducting condensation energy by the partial gap to Cooper pairs in the SP regions. If the mutual reinforcement is not strong enough compared to other processes, a $`d`$-Mott or equivalently dRVB condensate will not appear as the initial instability. Instead a long range order such as ordinary $`d`$-wave superconductivity on the whole FS at larger doping, when the FS moves away from flat SP regions, or AF order at small doping - sacrificing the partial spin gap - when the nesting over the the whole FS becomes too strong.

Summarizing we have determined a reduced fermionic Hamiltonian that contains the essential ingredients of a $`d`$-Mott state. For the half-filled two-leg Hubbard ladder we have shown that exact diagonalization of this Hamiltonian reproduces the known spin and charge gaps. For the 2D n.n.n. Hubbard model we have described how such a Hamiltonian captures the leading scattering processes between the saddle points, and hence how a $`d`$-Mott condensate can truncate the Fermi surface near the saddle points. Finally we have outlined how superconducting and antiferromagnetic states can arise as additional instabilities on the remaining parts of the Fermi surface, so that the $`d`$-RVB phase can be viewed as an unstable fixed point as proposed by Anderson anderson . In 2D our picture still relies on several approximations and assumptions, but we hope that the present work may provide some guidance for further research.

###### Acknowledgements.

We thank D. Poilblanc and M. Salmhofer for stimulating discussions. Financial support through the Swiss National Science Foundation (A.L.) and the German research foundation (DFG, C.H.) is acknowledged.


text:
# A search for disks around massive young stellar objects

## 1. Introduction

The existence of circumstellar disks around low-mass protostars is now well-established (Mundy, Looney & Welch 2000). However, the situation for their massive counterparts remains unclear. Indirect observational evidence for circumstellar disks around massive young stellar objects (YSOs) comes from the presence of bipolar outflows, polarization patterns and the presence of ionized winds extending in an equatorial direction (Hoare 2002). In one case, a model of the 7-mm emission shows that the data are consistent with a circumstellar disk with a diameter of 130 AU (Shepherd, Claussen & Kurtz 2001). Molecular emission line signatures of rotation have mostly been confined to larger and therefore non-circumstellar scales (e.g. Zhang, Hunter, & Sridharan 1998).

The question of whether or not circumstellar disks are present around massive stars is key in current theoretical debate. One important question is whether high-mass star formation is simply a scaled-up version of low-mass star formation where accretion onto the central star occurs via a disk. Direct observation of disks around massive YSOs is key to testing the theories for how massive stars form. It is therefore important to determine the circumstellar mass distribution around luminous YSOs.

## 2. Observations

We have used the Berkeley-Illinois-Maryland Association (BIMA: see Welch et al. 1996) millimetre array to observe a sample of 8 massive young stellar objects: GL 490, Cep A, NGC7538-IRS1, NGC7538-IRS9, G35.2$``$0.7N, GGD27-IRS, W75N and GL 2591. The observations were made in continuum bands at 2.7 mm (107 GHz) for all sources (except GGD27-IRS which was observed at 95 GHz) and 1.4 mm (216 GHz) for GL 490, Cep A and NGC7538-IRS1 only. Data were conventionally calibrated using a nearby phase calibrator, while fast-switching was employed in the more extended configurations (see Looney, Mundy & Welch 2000). In three sources (Cep A, NGC7538-IRS1 and G35.2$``$0.7N), the 107-GHz methanol maser was sufficiently strong to apply a self-calibration solution (the 95-GHz methanol maser was used for GGD27-IRS). The configurations used were A and B (at 2.7 mm) and B and C (at 1.4 mm), yielding an angular resolution varying from 0.3 arcsec (for Cep A, GL 490 and NGC7538-IRS1) to approximately 1 arcsec at 2.7 mm, which translates to a radius of 105 to 350 AU at the distance to our nearest source (Cep A). At 1.4 mm the angular resolution was $``$1.5–2.0 arcsec. The typical sensitivity was 3–8 mJy bm-1 at 2.7 mm, rising to 20–30 mJy bm-1 at 1.4 mm.

## 3. Results

### 3.1. GL 490

Figure 1 shows the 2.7-mm and 1.4-mm continuum emission towards GL 490. The 2.7-mm emission arises from an elongated structure with dimensions 0.7$`\times `$0.5 arcsec2, corresponding to 700$`\times `$500 AU2 at a distance of 1 kpc. The emission peaks on the position of the VLA source with an extension 0.3-arcsec to the south-east. It is not clear whether this represents part of a circumstellar disk around GL 490 or if it marks the location of a secondary companion. However, we note that the high-resolution $`H`$-band speckle-interferometry image of Hoare, Glindemann & Richichi (1996) shows no evidence for multiplicity reinforcing the disk interpretation. The position angle of the 2.7-mm emission is in good agreement with the larger structures observed by Mundy & Adelman (1988) and more recently by Schreyer et al. (2002).

An image made using only the longest baselines yields a weak point source centred on the VLA position. The flux of this source is in good agreement with the predicted flux of the wind contribution at 2.7 mm indicating that the flux is dominated by the ionized wind on scales of 0.1 arcsec (100 AU) or less.

On the other hand, the 1.4-mm emission is dominated by dust, although our resolution is not high enough to resolve any structure. The image shown in Fig. 1b is consistent with a point source within the 2.4$`\times `$1.1 arcsec2 beam.

### 3.2. Cep A

Cep A2 is the well-known radio-jet source in the Cep A region. The 2.7-mm continuum shown in Figure 2a is slightly elongated along the direction of the jet with FWHM dimensions of 0.45$`\times `$0.28 arcsec2 corresponding to 315$`\times `$200 AU2 at a distance of 700 pc. The emission at a resolution of 0.3 arcsec remains dominated by the ionized jet although images made using the more compact configurations reveal the presence of dust emission from the envelope. At 1.4 mm (Fig. 2b) the emission is dominated by Cep A2 although there appears to be some extended emission which encompasses the methanol masers to the west of Cep A2 with a further extension to the SE. This extended emission may represent the remnant flattened core from which Cep A2 formed. A Gaussian fit to the dust emission yields a slightly extended source with deconvolved dimensions of 0.8$`\times `$0.5 arcsec2 (560$`\times `$350 AU2). Weak dust emission is also detected from two of the Cep A3 group of sources (3b and 3c: Garay et al. 1996).

While the water masers mostly cluster around Cep A2 itself, the methanol masers show no clear association with any particular source. Our BIMA observations clearly resolve the two 107-GHz maser spots of Mehringer et al. (1997) separate from Cep A2. The stronger maser coincides with one of the 12.2 GHz methanol masers, both spatially and spectrally.

### 3.3. Other sources

NGC7538-IRS1 is barely resolved at 2.7 mm: the ‘south spherical’ source of Gaume al. (1995) is discernible. The presence of this source complicates the spectrum making it difficult to disentangle the contributions from the dust and the ionized gas.

The strong 95-GHz methanol maser is detected towards GGD 27-IRS and resolved to be a single maser spot 9.2 arcsec north-north-east of the jet source. This position places the maser along the axis of the jet and it may arise in one of the jet shocks. The flux of GGD 27 is somewhat uncertain but is probably $``$50 mJy.

Not all of our mapped sources had strong compact emission at 2.7 mm. NGC7538-IRS9 and G35.2$``$0.7N were only weakly detected at the highest resolutions, and GL 2591 was not detected at all to a 3-$`\sigma `$ sensitivity limit of 20 mJy in a 1-arcsec beam at 2.7 mm.

## 4. Analysis: dust emission around massive YSOs

Figure 3 plots the spectral energy distribution from radio through to millimetre wavelengths for four of our sources. Radio measurements were taken from the literature while the 2.7 and 1.4-mm points (filled squares) are from our data. The radio spectral index has been used to predict the contribution from the ionized wind at millimetre wavelengths and is shown by the straight lines in Fig. 3. In each case it is clear that dust emission has been detected even at 2.7 mm.

Accurate mass estimates require detailed modelling, but here we present simple estimates of the total mass assuming that the dust emission is optically thin and is at a uniform temperature of 75 K (see, e.g., Mundy & Adelmann 1988). We take the dust emissivity to be the same as that used by Looney et al. (2000) for low-mass sources. Note that the mass estimates roughly scale linearly (but inversely) with temperature. In all cases the mass estimate is of order 1–4 M within a region $``$1 arcsec in diameter (see Table 1).

While it would be a nice result to be able to claim that we have detected disks around these massive YSOs, in practice our resolution and sensitivity are not high enough to determine whether the emission comes from flattened disk-like structures or from a larger-scale envelope. At 2.7 mm images of both GL 490 and Cep A2 made with just the very longest baselines of the A configuration show unresolved sources ($`<`$0.2 arcsec FWHM) which are consistent with emission from the ionized stellar wind. For GL 490 the wind contribution of $``$25 mJy allows us to place an upper limit on the mass within a radius of 0.15′′ (or 150 AU) of $``$0.5 M (assuming the same dust properties as Table 1). In order to separate the dust emission on scales smaller than 0.2 arcsec, high-resolution observations at 1.4-mm (or shorter wavelengths) are required. Even then, discovery of a disk can not be claimed until a kinematic signature is observed. We are currently examining our line data but so far have not detected a clear rotation signature.

## 5. Conclusions: the search for disks goes on

We have presented results from a search for dusty disks around a sample of massive YSOs. We have found that while dust is clearly present around these sources, the dust emission at 2.7 mm is too weak to separate from the free-free emission from the ionized wind. Our observations at 1.4-mm do not have sufficient resolution to probe the geometry, but do show that several solar masses of dusty material exist on scales of a 1000 to 2000 AU from the central source. The exception is GL 490 which we clearly resolve into an elongated structure with dimensions 0.7$`\times `$0.5 arcsec2, or 700$`\times `$500 AU2. Although it is not possible for us to determine unambiguously whether this structure represents a circumstellar disk or if we are beginning to resolve a binary system, we note that the near-infrared speckle imaging of Hoare et al. (1996) revealed no evidence for multiplicity towards GL 490.

Future observations at 2.7 mm will need significantly higher sensitivity and resolution. A more promising approach is to employ observations at 1.4 mm or shorter wavelengths where the wind contribution will be negligible. In this respect the completion of CARMA (http://www.mmarray.org/) at its high-altitude site will mark the beginning of a new era in the search for disks around massive YSOs.

### Acknowledgments

BIMA is funded by NSF grant AST-9981289.

## References

Garay, G.M., Ramírez, S., Rodríguez, L.F., Curiel, S., Torrelles, J.M. 1996, ApJ, 459, 193

Gaume, R.A., Goss, W.M., Dickel, H.R., Wilson, T.L., Johnston, K.J. 1995, ApJ, 438, 776

Hoare, M.G. 2002, in ASP Conf. Ser. Vol. 267, The Earliest Stages of Massive Star Birth, ed P.A. Crowther (San Francisco: ASP), 137

Hoare, M.G., Glindemann, A., Richichi, A. 1996, in The Role of Dust in the Formation of Stars, ESO Proc, ed H.U. Käufl & R. Siebenmorgen (Berlin: Springer-Verlag), 35

Looney, L.W., Mundy, L.G., Welch, W.J. 2000, ApJ, 529, 477

Mehringer, D.M., Zhou, S., Dickel, H.R. 1997, ApJ, 475, L57

Mundy, L.G., Adelman, G. 1988, ApJ, 329, 907

Mundy, L.G., Looney, L.W., Welch, W.J. 2000, in Protostars and Planets IV, ed V.G. Mannings, A.P. Boss & S.S. Russell (Tucson: University of Arizona Press), 355

Schreyer, K., Henning, T., van der Tak, F.F.S., Boonman, A.M.S., van Dishoeck, E.F. 2002, A&A, 394, 561

Shepherd, D.S., Claussen, M.J., Kurtz, S.E. 2001, Science, 292, 1513

Welch, W.J., et al. 1996, PASP, 108, 93

Zhang, Q., Hunter, T.R., Sridharan, T.K. 1998, ApJ, 505, L151


text:
# The Dual Dust Chemistry - Binarity Connection

## 1. Prologue

WC Wolf-Rayet central stars of planetary nebulae (\[WC\] central stars of PNe) are H-deficient stars that exhibit strong ionic emission lines of helium, carbon and oxygen from their dense stellar winds. Amongst the coolest central stars in this group are CPD–56o8032 (the nucleus of the PN He 3-1333) and the central star of He 2-113 (both classified as \[WC10\], Crowther et al. 1998). Cohen et al. (1986) found the mid-infrared KAO spectra of both these objects to show very strong unidentified infrared bands (UIBs – usually attributed to polycyclic aromatic hydrocarbons, PAHs). Indeed, both the nebular C/O ratios (De Marco, Barlow & Storey 1997) and the ratio of UIB luminosity to total IR dust luminosity (Cohen et al. 1989) for these two objects are the largest known. It was therefore a major surprise when mid- and far-IR Infrared Space Observatory spectra of these two objects showed the presence of many emission features longwards of 20 $`\mu `$m that could be attributed to crystalline silicate and water ice particles (Barlow 1997, Waters et al. 1998a, Cohen et al. 1999), indicating a dual dust chemistry, i.e. the simultaneous presence of both C-rich dust and O-rich dust. The dual dust chemistry phenomenon in PNe appears to show a strong correlation with the presence of a late WC (\[WCL\]) nucleus – four out of six \[WC8-11\] nuclei studied by Cohen et al. (2002) showed similar dual dust chemistries. In the context of a single star scenario, this would point to a recent transition (within the last $``$1000 yr) between the O- and the C-rich surface chemistries. However, the probability of finding a post-AGB object that had recently changed from an O-rich to a C-rich surface chemistry due to a third dredge-up event should be very low indeed. An alternative scenario envisages these systems as binaries (Waters et al. 1998a, Cohen et al. 1999, 2002), in which the O-rich silicates are trapped in a disk as a result of a past mass transfer event, with the C-rich particles being more widely distributed in the nebula as a result of recent ejections of C-rich material by the nucleus.

## 2. A Dusty Disk around CPD–56o8032

CPD–56o8032 exhibits a variable optical lightcurve (Pollacco et al. 1992, Jones et al. 1999), with minima which are $``$1 mag deep and have a 5-year pseudo-period (Cohen et al. 2002). This was interpreted by Cohen et al. as precession of a dusty disk. We present an updated light-curve for CPD–56o8032 in Fig. 1.

De Marco et al. (2002) presented HST/STIS spatially-resolved spectroscopy of CPD–56o8032, obtained during its third visual light minimum, which revealed the stellar near UV and blue optical continuum to be split into two spectra separated by about 0.10 arcsec. Once other possibilities were excluded (e.g., instrumental problems or the presence of two stars), the simpler explanation is that CPD–56o8032’s light is seen reflected by the upper and lower rim of a dusty disk or torus seen close to edge-on. If so, CPD–56o8032 could be similar to the Red Rectangle (HD44179; Osterbart et al. 1997), whose binary central star light is seen indirectly, reflected by the rims of an edge-on disk. In this binary scenario, CPD–56o8032’s light-curve could also be interpreted as CPD–56o8032’s orbital motion taking it in and out of alignment with a denser part of the disk, such that its brightness is modulated in step with the binary period. An alternative scenario that is more in line with the asymmetric lightcurve declines, conceives a dust clump in orbit with the disk passing in front of CPD–56o8032 every 5 years. The companion cannot be bright enough to contaminate substantially its spectrum, which is reasonably fit with a single Wolf-Rayet model atmosphere (De Marco & Crowther 1998). CPD–56o8032’s putative companion could plausibly be a low mass main sequence star, such as a K dwarf.

Interestingly, the Red Rectangle, also a binary, is the only system to have a double dust chemistry and a central star which is not hydrogen-deficient. From this we could suggest that the Red Rectangle’s A supergiant central star (the other is thought to be an unseen white dwarf \[Men’shchikov et al. 2002\]) might be on its way to becoming hydrogen deficient.

The conjecture that \[WC\] stars are binary systems is in harmony with the fact that binarity is generally known to promote mass-loss. For instance, all the massive Wolf-Rayet stars in the metal-poor Small Magellanic Cloud are binaries. This has been explained with the fact that the low metallicity of the SMC leaves its stars with atmospheres with relatively low opacities, too low to develop the dense Wolf-Rayet winds. Hence the only massive Wolf-Rayet stars possible in the SMC are those where a companion has facilitated mass-loss.

A direct detection of binarity via radial velocity monitoring (see Section 3), is unlikely in the case of \[WC\] central stars, because of their intrinsically variable (e.g. Balick et al. 1996), broad emission lines.

## 3. Preliminary Results of a Central Star Radial Velocity Survey

There are increasing indications that binary-star processes are intimately related to the ejection of many, or possibly even most, PNe. The evidence includes: the fact that $``$10% of PN nuclei are found to be very close binaries (periods of hours to a few days; Bond & Livio 1990, Bond 2000) through photometric monitoring; population-synthesis studies suggesting that these may be just the short-period tail of a much larger binary population extending up to orbital periods of several months (Yungelson et al. 1993, Han et al. 1995); and the prevalence of highly non-spherical morphologies among PNe.

The photometric search technique does not work for binaries with periods of more than a few days, since it relies on proximity effects. We have therefore carried out radial-velocity monitoring of a sample of PN nuclei, in order to search for the anticipated population of binaries with periods up to a few months. If they do exist, there will be new implications for the evolution of binary populations, the origin of compact binaries (CVs, SN Ia progenitors), and even the question of whether single stars can produce visible ionized PNe at all.

This program was started in May 2002 at the 3.5-m WIYN telescope. Despite being plagued by very bad weather and some instrument problems, we report that 8 of the 18 central stars for which we have more than five data points have clear radial velocity variability. The independent, but similar program of Pollacco et al. (these proceedings) returned a similar radial velocity variable central star fraction. However, it will take additional observations for us to have robust statistics, and to determine orbital periods needed to show conclusively that the velocity variations are due to binary motion and not, for instance, to wind variability.

### Acknowledgments.

OD gratefully acknowledges Janet Jeppson Asimov for financial support. OD and MC acknowledge support from NASA grant HST-GO-08711.05-A. We are grateful to the WIYN telescope team.

## References

Balick, B., Rodgers, B., Hajian, A., Terzian, Y., Bianchi, L. 1996, AJ111, 843 Barlow, M. J. 1997, ApSS, 255, 315

Bond, H. E. 2000, in Asymmetrical planetary nebulae II, eds. J. H. Kastner, N. Soker & S. Rappaport (ASP Conference Series 199), 115

Bond, H. E., & Livio, M. 1990, ApJ, 355, 568

Cohen, M., Allamandola, L., Tielens, A. G. G. M., Bregman, J., Simpson, J. P., Witteborn, F. C., Woodenm D. & Rank, D. 1986, ApJ, 302, 736

Cohen, M., Tielens, A. G. G. M., Bregman, J. D., Allamandola, L. J., Wooden, D., & Jourdain de Muizon, M. 1989, ApJ, 341, 246

Cohen, M., Barlow, M. J., Sylvester, R. J., Liu, X.-W., Cox, P., Lim, T., Schmitt, B. & Speck, A. K. 1999, ApJ, 513, L135

Cohen, M., Barlow, M. J., Liu, X.-W., Jones, A. F. 2002, MNRAS, 332, 879

Crowther, P. A., De Marco, O., & Barlow, M. J. 1998, MNRAS, 296, 367

De Marco, O., Barlow, M. J. & Storey, P. J. 1997, MNRAS, 292, 86

De Marco, O., & Crowther, P. A. 1998, MNRAS, 296, 419

Han, Z., Podsiadlowski, P., & Eggleton, P. 1995, MNRAS, 272, 800

Jones, A., Lawson, W., De Marco, O., Kilkenny, D., van Wyk, F., & Roberts, G. 1999, The Observatory, 119 76

Men’shchikov, A. B., Schertl, D.; Tuthill, P. G., Weigelt, G., Yungelson, L. R. 2002, A&A, 393, 867

Osterbart, R., Langer, N., & Weigelt, G. 1997, å, 325, 609

Pollacco, D. L., Kilkenny, D., Marang, F., van Wyk, F., & Roberts, G. 1992, MNRAS256. 669

Waters, L. B. F. M., et al. 1998, å331, L61

Yungelson, L. et al. 1993, ApJ, 418, 794


text:
# The Thermal Stability of Mass-Loaded Flows

## 1 Introduction

The thermal instability of radiative media was examined by Field (F1965 (1965)), who considered both the presence and absence of thermal conduction, and derived the growth rates of isobaric and isentropic perturbations. The first numerical calculations of catastrophic cooling in shock heated gas were performed by Falle (F1975 (1975), F1981 (1981)). Radiative shocks were shown to exhibit a global overstability by Langer et al. (L1981 (1981)), and have since been extensively examined (e.g., Chevalier & Imamura CI1982 (1982); Imamura, Wolff & Duisen IWD1984 (1984); Gaetz, Edgar & Chevalier GEC1988 (1988); Blondin & Cioffi BC1989 (1989); Strickland & Blondin SB1995 (1995); Walder & Folini WF1996 (1996)).

The interaction of radiative flows with cold embedded clouds is known to significantly modify these flows (see, e.g., Pittard, Dyson & Hartquist PDH2001 (2001), and references therein). However, there has yet to be an investigation into how mass-loading may affect the stability properties of such flows. This is the aim of this paper.

In Sec. 2 we perform a linear stability analysis of a static medium, in thermal equilibrium, undergoing conductively-driven mass-loading. In Sec. 3 we present numerical models of mass-loaded radiative shocks to examine the suppression of thermal instability in them. We finish in Sec. 4 with a discussion on the conditions necessary for mass-loading to suppress the thermal instability in a typical planetary nebula, and provide numerical estimates for the Helix nebula.

## 2 Instability in a uniform medium, including the effects of conductively-driven mass-loading

The dynamics are governed by the standard hydrodynamic equations for plane parallel flow,

$$\frac{\rho }{t}+\frac{(\rho v)}{z}=q_0\left(\frac{T}{T_0}\right)^\lambda q_0,$$
(1)
$$\frac{(\rho v)}{t}+\frac{(\rho v^2)}{z}+\frac{P}{z}=0,$$
(2)
$`{\displaystyle \frac{1}{2}}{\displaystyle \frac{(\rho v^2)}{t}}+{\displaystyle \frac{1}{(\gamma 1)}}{\displaystyle \frac{P}{t}}`$ (3)
$`+{\displaystyle \frac{}{z}}\left({\displaystyle \frac{1}{2}}\rho v^3+{\displaystyle \frac{\gamma }{(\gamma 1)}}vP\right)=\rho ,`$
together with an equation of state,

$$P=\frac{R}{\mu }\rho T.$$
(4)
We have assumed that the clouds are all at their equilibrium radius so that if $`T<T_0`$ they grow and that if $`T>T_0`$ they evaporate (McKee & Cowie MC1977 (1977)). $`q_0`$ has the dimensions of a mass-loading rate per unit volume.

The equilibrium state is characterized by $`\rho =\rho _0`$, $`T=T_0`$, $`v=0`$, and $`(\rho _0,T_0)=0`$. Assuming perturbations of the form

$$a(z,t)=a_1\mathrm{exp}(nt+ikz),$$
(5)
we find the linearized equations for the perturbations to be

$$n\rho _1+\rho _0ikv_1=q_0\lambda T_1/T_0,$$
(6)
$$n\rho _0v_1+ikP_1=0,$$
(7)
$$\frac{n}{\gamma 1}+\frac{\gamma }{\gamma 1}ikv_1P_0=\rho _0_TT_1\rho _0_\rho \rho _1,$$
(8)
and

$$P_1\frac{R}{\mu }\rho _1T_0\frac{R}{\mu }\rho _0T_1=0,$$
(9)
where $`_\rho (/\rho )_T`$ and $`_T(/T)_\rho `$ are evaluated for the equilibrium state. We are left with 4 variables ($`\rho _1,P_1,T_1`$, and $`v_1`$) in 4 equations. The resulting dispersion relation is

$`n^3+n^2c\left(k_T+{\displaystyle \frac{k_m}{\gamma }}\right)+nc^2\left(k^2+{\displaystyle \frac{k_\rho k_m}{\gamma }}\right)`$ (10)
$`+{\displaystyle \frac{c^3k^2}{\gamma }}(k_Tk_\rho +k_m)=0.`$
The adiabatic speed of sound is $`c=(\gamma P_0/\rho _0)^{1/2}`$, and we have introduced the wavenumbers

$$k_\rho =\frac{\mu (\gamma 1)\rho _0_\rho }{RcT_0},$$
(11)
$$k_T=\frac{\mu (\gamma 1)_T}{Rc},$$
(12)
and

$$k_m=\frac{\gamma \lambda q_0}{\rho _0c}.$$
(13)
By introducing the non-dimensional variables

$`y={\displaystyle \frac{n}{kc}},\sigma _\rho ={\displaystyle \frac{k_\rho }{k}},\sigma _T={\displaystyle \frac{k_T}{k}},`$
$`\sigma _m={\displaystyle \frac{k_m}{k}},`$ $`\sigma _T^{}=\sigma _T+\sigma _m,`$ (14)
we can write the dispersion equation in the form

$$y^3+y^2\sigma _T^{}+y(1+\sigma _\rho \sigma _m)+(\sigma _T+\sigma _m\gamma \sigma _\rho )/\gamma =0.$$
(15)
The coefficient in the y term can be removed by the introduction of the variable

$$y^{}=\frac{y}{(1+\sigma _\rho \sigma _m)^{1/2}}.$$
(16)
The dispersion relation then becomes

$$y^3+\frac{\sigma _T^{}}{(1+\sigma _\rho \sigma _m)^{1/2}}y^2+y^{}+\frac{\sigma _T+\sigma _m\gamma \sigma _\rho }{\gamma (1+\sigma _\rho \sigma _m)^{3/2}}=0.$$
(17)
This equation is now in the same form as Eq. 18 in Field (F1965 (1965)).

The growth of the isobaric instability (which Field refers to as a condensation mode) requires

$$\sigma _T\sigma _\rho <\gamma \sigma _m.$$
(18)
Since $`\sigma _m`$ is positive by definition, mass-loading always acts to reduce the growth rate of this instability mode (as Field found for conduction). However, unlike the corresponding equation including conduction, Eq. 18 is independent of $`k`$, so mass-loading stabilizes all wavelengths equally effectively against isobaric perturbations. Rewriting Eq. 18, we find that the instability is suppressed if

$$\frac{q_0}{\rho _0}>\frac{(\gamma 1)\mu }{\gamma \lambda R}\left(\frac{\rho _0_\rho }{T_0}_T\right).$$
(19)
This is roughly equivalent to requiring that the mass-loading timescale be less than the cooling timescale of the medium.

The growth of the isentropic instability requires

$$\sigma _T+\frac{\sigma _\rho }{\gamma 1}<\frac{\gamma }{\gamma 1}\sigma _m\sigma _\rho (\sigma _T+\sigma _m).$$
(20)
Mass-loading is again always stabilizing, but is not equally effective at all wavelengths in suppressing isentropic perturbations. Rewriting Eq. 20 we find that the critical wavenumber above which perturbations are stabilized is

$$k_\mathrm{c}=\sqrt{\frac{k_mk_\rho (k_T+k_\rho )}{(\gamma 1)k_T+k_\rho }}.$$
(21)
and that the isentropic instability is suppressed if

$$q_0>\frac{k^2c^4}{\gamma ^2(\gamma 1)\lambda _\rho }\frac{(\gamma 1)_T+\rho _0_\rho /T_0}{_T+\rho _0_\rho /T_0}.$$
(22)
Over the temperature range $`5\times 10^5<T<5\times 10^7\mathrm{K}`$, a good approximation is $`=\rho \mathrm{\Lambda }T^{1/2}`$ (e.g., Kahn K1976 (1976)). This leads to greatly simplified versions of Eqs. 19 and 22, which respectively become

$$\frac{q_0}{\rho _0}>\frac{3}{2}\frac{(\gamma 1)\mu }{\gamma \lambda R}\rho _0\mathrm{\Lambda }T_0^{3/2},$$
(23)
and

$$q_0>\frac{2k^2c^4}{\gamma ^2(\gamma 1)\lambda \mathrm{\Lambda }T_0^{1/2}}\left[1\frac{1}{2}(\gamma 1)\right],$$
(24)
and to the relationship $`k_\rho /k_T=\sigma _\rho /\sigma _T2`$.

## 3 Hydrodynamical calculations

The nature of the overstability of radiative shocks is known to depend on the temperature dependence of the local cooling rate. For a power-law dependence (e.g., $`\rho T^\alpha `$), the system is overstable for values of $`\alpha `$ below some critical value, $`\alpha _{\mathrm{cr}}`$. Systems with $`\alpha >\alpha _{\mathrm{cr}}`$ are stable. Previous numerical work has shown that $`\alpha _{\mathrm{cr}}0.4`$ (e.g., Imamura et al. IWD1984 (1984); Strickland & Blondin SB1995 (1995)), in good agreement with the linear stability analysis of Chevalier & Imamura (CI1982 (1982)).

To investigate the effect of mass-loading on the stability of a radiative shock, we have computed 1D numerical calculations for a Mach 10 shock with the condition that $`\alpha =0.0`$. Our computations were initialized in a similar fashion to that presented by Strickland & Blondin (SB1995 (1995)) for the case of an isolated planar shock, and were performed using the same hydrodynamical code, VH-1 (see Blondin et al. BKFT1990 (1990)).

Briefly, we assumed a steady-state radiative shock in the centre of the grid, with pre-shock flow from the left grid boundary, and cold, dense, post-shock material exiting the right grid boundary. The most common downstream condition for similar calculations in the literature is that of a reflecting wall, but a continuous outflow has several benefits. For instance, feedback between the cold dense post-shock layer and the cooling gas is properly treated, and the problem more closely resembles the common occurence of an isolated shock in the interstellar medium. We also repeat the cautionary note of Strickland & Blondin (SB1995 (1995)) that the boundary conditions play an important role in determining the stability of the shock. Therefore, it is desirable to place the shock well away from grid boundaries. We refer the reader to Strickland & Blondin (SB1995 (1995)) for a fuller description of the code and initial conditions. In our runs we assumed $`\gamma =5/3`$.

In the left panel of Fig. 1 we show the development of the overstability for a Mach 10 flow with $`\alpha =0.0`$ and $`q_0=0.0`$ (i.e. no mass-loading). The pre-shock density and velocity are set to unity. At the beginning of the simulation, corresponding to the top of the plot, the flow is initialized as closely as possible to the steady state solution for a radiative shock with no mass-loading. The overstability is excited from weak perturbations produced by the numerical mapping of the steady state solution onto the computational grid. Linear growth was observed for about the first 40% of the run, after which the amplitude of the oscillations saturate and the system enters the nonlinear regime. The behaviour shown in the left panel of Fig. 1 is in excellent agreement with Fig. 3 of Strickland & Blondin (SB1995 (1995)) where the overstability for a Mach 40 flow is displayed.

In the right panel of Fig. 1 we show the evolution of a mass-loaded radiative shock with $`q_0=150.0`$ and $`\lambda =5/2`$. The process of mass-loading increases the density and thermal pressure, and decreases the temperature and velocity of the flow. The pressure increase causes the shock position to initially expand away from the dense cooling layer. However, the enhanced density leads to more rapid cooling, and thermal pressure support is soon lost. This causes the shock to reverse its direction of motion and to fall towards the cold dense layer. The shocked gas is then repressurized and this cycle is repeated. With each cycle, the amplitude of the oscillations decreases as the mass-loading damps this instability.

With $`\alpha =0.0`$, $`=1`$, $`\mu /\mathrm{R}=1`$ and the post-shock values $`\rho =4`$, $`T=0.193`$, Eq. 19 shows that we require $`q_0>13`$ for mass-loading to begin to damp the overstability. In practice, we find that we need values appreciably larger than this for strong damping, although the oscillation amplitude does decrease slightly for values of $`q_0`$ near 13. There are a couple of reasons why we need higher values of $`q_0`$ than suggested by Eq. 19. The principal reason is due to the fact that adding mass to the shocked gas decreases the cooling timescale while simultaneously increasing the mass-loading timescale. Hence, while the initial addition of mass to the immediate post-shock flow may be at a rate such that the mass-loading timescale is less than the cooling timescale, the difference between the two timescales will diminish with time and may even be reversed.

Secondly, the large dynamic range in temperature which exists for high Mach number shocks (e.g., for a Mach number of 40, the immediate post-shock temperature is almost 3 orders of magnitude greater than the ambient temperature) leads to a smaller value of $`T_{\mathrm{avg}}/T_{\mathrm{ps}}`$ than is the case for lower Mach number shocks. Here $`T_{\mathrm{avg}}`$ is the average temperature over the cooling length of the shock, and $`T_{\mathrm{ps}}`$ is the immediate post-shock temperature. We find that as we reduce the shock Mach number, the value of $`q_0`$ needed for significant damping more closely matches that from Eq. 19.

While we have not been concerned with the ionization state of the gas in this model, it is worth noting that the radiative cooling rate can be greatly enhanced when neutral gas is introduced into hot plasma, and that the ionization energy can also be significant under such conditions (Slavin et al. SSB1993 (1993)).

## 4 Discussion

As a relevant example, the above results can be applied to planetary nebulae (PNe). Such objects often display clumps which appear to mass-load the shocked wind of the central star. The clumps were presumably formed by instabilities in the atmospheres of the red supergiant stars, prior to the evolution of the central star into a hot object with a fast wind. In the following we derive relationships for the number density of clouds and their total mass, with the condition that mass injection from the clouds is able to suppress the thermal instability.

Since the central stars of PNe appear to fall into 3 groups (e.g., Kudritzki et al. K1997 (1997)), with winds spanning a large range of mass-loss rate ($`\dot{M}10^910^6\mathrm{M}_{}\mathrm{yr}^1`$), we choose not to focus on an individual nebula. Instead, we note that the typical thermal pressure in PNe, wind-blown-bubbles, and starburst superwinds is comparable to the pressure in clumps embedded within them, $`P/k10^7\mathrm{cm}^3\mathrm{K}`$. In PNe, this must balance the ram pressure from the wind of the central star. The post-shock number density is then $`n=4\times 10^7v_7^2k/\mu (P/k)_7=55v_7^2(P/k)_7\mathrm{cm}^3`$, where $`v_7`$ is the wind speed in units of $`1000\mathrm{km}\mathrm{s}^1`$ and $`(P/k)_7`$ is $`P/k`$ in units of $`10^7\mathrm{cm}^3\mathrm{K}`$. In the last step, and throughout the remainder of these calculations, we set $`\mu =10^{24}\mathrm{g}`$. The post-shock temperature, $`T=1.36\times 10^5v_7^2\mathrm{K}`$.

To suppress the thermal instability by mass-loading, Eq. 19 gives $`q_0>10^{32}v_7^7(P/k)_7^2\mathrm{g}\mathrm{cm}^3\mathrm{s}^1`$, where we have adopted $`\mathrm{\Lambda }=1.33\times 10^{19}/\mu ^2`$ and $`\alpha =\frac{1}{2}`$ (cf. Kahn K1976 (1976)), and substituted for $`T`$ and $`n`$. The mass evaporation rate from a single clump is $`\dot{m}=2.75\times 10^4T^{5/2}R_{\mathrm{pc}}(30/\mathrm{ln}\mathrm{\Lambda })\mathrm{g}\mathrm{s}^1`$ (Cowie & McKee CM1977 (1977)), where $`T`$ is the temperature of the ambient surroundings, and $`R_{\mathrm{pc}}`$ is the radius of the clump in parsecs. Since $`\mathrm{ln}\mathrm{\Lambda }30`$, $`\dot{m}1.9\times 10^{17}v_7^5R_{\mathrm{pc}}\mathrm{g}\mathrm{s}^1`$, and the required number density of clumps is $`n_\mathrm{c}=q_0/\dot{m}>1.6\times 10^6v_7^{12}R_{\mathrm{pc}}^1(P/k)_7^2\mathrm{pc}^3`$.

Observations of the Helix nebula (NGC 7293) indicate that the density within the clumps, $`N_\mathrm{c}10^6\mathrm{cm}^3`$, while they are unresolved on spatial scales of $`10^3\mathrm{pc}`$ (O’Dell & Handron DH1996 (1996)). Adopting $`R_{\mathrm{pc}}=10^3`$, $`v_7=2`$, and $`(P/k)_7=1`$ yields a number density of clumps of $`n_\mathrm{c}>4\times 10^5\mathrm{pc}^3`$. Since the radius of the Helix nebula is $`0.2\mathrm{pc}`$, we then expect $`1600`$ clumps, with a total mass of $`0.1\mathrm{M}_{}`$. While these estimates are in good agreement with observational inferences from the Helix nebula, small changes in $`v_7`$ can greatly affect these values. For the evaporation to be smooth on a spatial scale of $`1/k_\mathrm{m}`$, we require that the clumps are smaller than $`L=2\pi /k_\mathrm{m}`$, i.e. $`L<1.5\times 10^2v_7^6(P/k)_7^1\mathrm{pc}`$. With the above parameters, $`L<1\mathrm{pc}`$, and is consistent with the derived number density of clumps.

###### Acknowledgements.

We would like to thank the referee, John Raymond, for timely and constructive comments. JMP would also like to thank PPARC for the funding of a PDRA position. This research has made use of NASA’s Astrophysics Data System Abstract Service.


text:
# Carbon Nanotube Based Bearing for Rotational Motions

Also at \]SPEC, CEA Saclay, F-91191 Gif-sur-Yvette, France.

## Abstract

We report the fabrication of a nanoelectromechanical system consisting of a plate rotating around a multiwalled nanotube bearing. The motion is possible thanks to the low intershell friction. Indeed, the nanotube has been engineered so that the sliding happens between different shells. The plate rotation is activated electrostatically with stator electrodes. The static friction force is estimated at $`210^{15}`$ N/Å2.

Microfabricated motors have been one of the most studied microelectromechanical systems in the 80s and 90s Trimmer ; Fan ; Tai ; Mehregany . Many efforts have been invested in improving performance characteristics such as the rotational speed. It has been shown that the characteristics are mostly limited by the fabrication process through the low reproduction accuracy of the design geometry Mehregany . For example, microfabrication cannot realize regular enough surfaces at the contact between the rotor and the bearing. These irregularities are important sources of energy loss and reduce considerably the rotational speed.

Multiwalled carbon nanotubes (MWNTs) consist of several nested cylindrical shells. Electron transmission microscopy has shown the high degree of perfection at the atomic level of their structure Dresselhaus . Their structural perfection and geometry make MWNTs attractive as elements enabling a rotational motion for nanoelectromechanical systems (NEMS) Roukes ; Craighead . Fig. 1a shows the operating principle of such a MWNT element where one or more inner shells can slide with respect to some outer shells.

The utilization of MWNTs as bearings for rotational motion in NEMS has been further motivated by recent experiments that show that the intershell friction can be low in MWNTs Cummings ; Yu . In Ref.Cummings  the static friction force has been experimentally estimated to be lower than $`610^{15}`$ N/Å2, which gives a very low friction torque due to the thin diameters of the shells. Easy sliding motion between shells has also been predicted theoretically Charlier ; Crespi ; Saito . Interestingly, corrugation against sliding of incommensurate and free of disorder shells is predicted to be extremely small Crespi .

We report the realization of a NEMS where a metal plate rotates around a MWNT axle. The suspended MWNT is clamped between two anchor pads. The MWNT has been engineered with the electrical-breakdown technique Collins1 ; Collins2 ; Bourlon , as shown in Fig. 1(a), to allow access to the inner shells on which is fixed the metal plate. The metal plate is shown to rotate due to surface tension forces of a liquid drop that evaporates on the substrate. The rotational motion can also be activated electrostatically with two stator electrodes. In this way, the plate does not achieve complete rotations but can be positioned in directions lying between the two stator electrodes.

Our rotational NEMSs are fabricated through a four-step process. The MWNTs are synthetised by arc-discharge evaporation and carefully purified Bonard . From a dispersion in dichloroethane the nanotubes are dispersed onto a 1000 nm oxidized Si wafer. The nanotubes are imaged with atomic force microscopy (AFM). Thick MWNTs of diameter between 15 and 25 nm are selected because they are mechanically more robust. Their positions are recorded with respect to metallic marks. The nanotubes are then connected with two electrodes using electron beam lithography (Fig. 1(b)). These electrodes will be used both as conducting electrodes and as anchor pads that hold the nanotube bearing. The electrodes, which are 3 $`\mu `$m wide and consist of 10 nm Cr and 60 nm Au, have to be wide and thick enough such that they stay horizontal and suspended when the SiO2 is etched away at the end of the process.

Fig. 1(c) shows that several outer shells are removed in between the electrodes in order to have access to the inner ones. These outer shells are removed with the electrical-breakdown technique Collins1 ; Collins2 ; Bourlon  that we review briefly in this paragraph. The bias voltage applied on the tube is continuously increased while the current is recorded. When a sharp step in the current is detected the bias is quickly reduced to 0. This sharp step is usually a current reduction of $``$20 $`\mu `$A which corresponds to the oxidation of the outermost shell. Importantly, the shell is then removed along with most of the tube portion that is not covered by the electrodes. This leaves plenty of room for access to the inner shells.

The electrical-breakdown technique is repeated until several shells are removed to obtain a sliding motion with a lower friction. The sliding will occur through a self-selection process between the shells with the most perfect surfaces that offer the least resistance to motion. However, we cannot repeat this process until we are left with one or a very few shells as the structure is not robust enough to sustain the drying step at the end of the process. We have found that the samples with inner rotating shells of diameter around 15 nm work well.

Fig. 1(d) shows the structures that are fabricated in a second lithography step. A 500 nm long plate is attached above the rotating inner shells. The plate is asymmetrically positioned with respect to the tube so that the longer section can be electrostatically attracted to one of the two stator electrodes that are fabricated during the same fabrication step. These 200 nm wide electrodes are designed to be narrow so that they will be deposited on the substrate in the etching step making the plate rotation between the two electrodes easier.

The last step of the process consists of etching 500 nm of the SiO2 with BHF Nygard . The wafer is rinsed in DI and in ethanol. It is then dried under nitrogen flow. The ethanol is warmed to just below the boiling point in order to reduce the surface tension. This reduces the risk of the structures breaking down during drying.

Figure 2 shows two examples of obtained structures at the end of the fabrication process. Interestingly, the plate is not horizontal but has already rotated. The motion has been induced by surface-tension forces of ethanol when drying the sample. Such forces are well known in NEMS fabrication to be important and to deform suspended structures. In contrast, it has been shown that the plate stays horizontal for samples where outer shells of the MWNT are not removed Williams . This indicates that the sliding between the moving and the fixed elements of the NEMS occurs between MWNT shells. Note that the rotation has been observed on most but not all of the samples. In the cases where the plate has not moved, some lithography-resist is observed between the plate and the rest of the structure that may block the motion.

We now show that the plate can be electrostatically driven by stator electrodes. Fig. 3 shows top down images of the device made with a scanning electron microscope (SEM) that has been modified in order to electrically access the different electrodes. In Fig. 3(a) the plate is oriented toward the stator electrode on the right side. A bias voltage $`V`$ is applied between the plate and the stator electrode on the left side in order to rotate the plate in this direction. The motion does not happen at once but occurs in the following manner. $`V`$ is increased continuously from zero while the other stator electrode and the anchor pads are grounded. The plate stays immobile until 40 V when the plate rotates suddenly to an almost vertical position (Fig. 3(b)). Another angular displacement is observed at 49 V (Fig. 3(c)). The bias has to be increased to 59 V to reach the final position, where the plate faces the biased electrode (Fig. 3(d)).

Importantly, the plate stays in this direction when the bias is turned off. This shows that the rotation is not enabled through the torsion of the tube which would return the plate to the initial position and which occurs for non engineered MWNTs Williams . We note that the plate stays fixed in this direction even if 100 V is applied on the right side electrode. This electrode lies further away from the tube and it is not surprising that the associated electrostatic torque is not enough large to initiate the rotation. Optimization of the design geometry should solve this.

We estimate now the static friction force. Using a finite element method program and the geometry of the plate and the stator electrode together with $`V=40`$ V, we have calculated the electrostatic energy as a function of the plate position. We deduce that the electrostatic couple $`10^{16}`$ Nm. This suggests that the static friction force is $`210^{15}`$ N/Å2 considering that sliding happens on a cylindrical surface that is 1 $`\mu `$m long and 18 nm in diameter from AFM. This is a low friction value compared with previously reported values, which shows that the electrical breakdown technique is suitable for the fabrication of such NEMSs.

Here we discuss briefly the rotational motion of the plate though a systematic investigation is left open for further work. The plate may block at some specific angles due to variations of the friction between the shells. Friction variation can result from the displacement of some disorder centers situated along the tube or some dangling bonds at the extremities left over from the shell oxidation. Another explanation might be the variation of the capacitance between the plate and the bias electrode. The blockage might then occur when the electrostatic torque becomes lower than the friction torque due to corrugation in the intershell interaction.

After this work was completed, related results were reported in Ref. Fennimore . A third stator electrode was used in Ref. Fennimore  that has allowed for the complete rotation of the plate around the axle. In contrast to our devices, the electrical-breakdown technique has not been used to engineer the MWNT. The plate has been attached directly on the pristine MWNT. The rotational motion is then obtained by removing one or more outer shells in the region between the plate and the anchors. The removal method is to mechanically fatigue and to shear the outer shells by successive application of large stator voltages.

We wish to acknowledge D. Esteve, P. Joyez, H. Camon and N. Fabre for discussions, B. Placais, J.M. Berroir and C. Delalande for support and P. Morfin for technical assitance. LPMC is CNRS-UMR8551 associated to universities Paris 6 and 7. The research has been supported by the DGA and ACN programs.