text:
# Direct observation of the phonon energy in a Bose-Einstein condensate by tomographic imaging

\[

## Abstract

The momentum and energy of phonons in a Bose-Einstein condensate are measured directly from a time-of-flight image by computerized tomography. We find that the same atoms that carry the momentum of the excitation also carry the excitation energy. The measured energy is in agreement with the Bogoliubov spectrum. Hydrodynamic simulations are performed which confirm our observation.

\]

In a Bose-Einstein condensate (BEC), the long wavelength excitations are phonons, characterized by a linear dispersion relation, as well as a wave function consisting of atoms with positive and negative momenta. Both the dispersion relation and the wave function were measured indirectly, by Bragg spectroscopy  .

A two-photon Bragg process coherently imparts an energy $`\mathrm{}\omega `$ and momentum $`\mathrm{}k`$, determined by the frequency difference and the angle between the Bragg beams, to the condensate . In Bragg spectroscopy, the two-photon energy $`\mathrm{}\omega `$ is varied, and the response of the BEC to the two-photon Bragg transition is observed   . Thus, the directly measured energy is the photon energy. Calculating the excitation energy requires knowledge of the process by which the photons impart momentum and energy to the BEC. Analogously, excitations in superfluid $`^\text{4}`$He were observed indirectly, by measuring the momenta and the energies of scattered neutrons .

The phonon energy is much larger than that of a free particle with the same momentum, and consists mostly of interaction energy between the phonon and the condensate. When the trapping potential is turned off, all the interaction energy is transformed into kinetic energy during a short acceleration period .

In this letter we present a direct measurement of this energy from time-of-flight (TOF) images of the released atoms. Previously the ground state interaction energy of the condensate was also measured from TOF absorption images . Since absorption images provide only the density of the cloud integrated along the axis of the absorption beam, whereas energy measurements require knowledge of the full density distribution, fits to model functions of the density were needed . Here we use computerized tomography of the TOF absorption images, in which the cylindrical symmetry of the cloud is used to reconstruct the full density profile from a single absorption image.

Several measurements are performed using this technique. First, the ground state interaction energy of a released, unexcited BEC is measured and compared to the value extracted from a fit to the radial size of the absorption image. Second, the density distribution of a Bragg-excited condensate is reconstructed. A clear separation between a released-phonon cloud and a condensate cloud, which is not clear in the original absorption image, is visible. The energy and momentum added to the cloud by the Bragg pulse are measured from the reconstructed image. In this way the excitation energy is found without sweeping through frequencies, and by looking directly at images of the atoms. Finally, the energy carried by the released phonons is measured by observing the distribution of the total energy between the two clouds. Thus, we extract the excitation energy from a single absorption image. Hydrodynamic simulations of the phonon release process are performed and found to be consistent with our observations.

In a previous experiment the interaction energy between two identical condensates, with a large relative velocity, was measured by their interference pattern . For this case a one-dimensional treatment was sufficient.

The relation between $`E_k`$, the energy, and $`\mathrm{}k`$, the momentum, of an excitation in a homogenous condensate has the Bogoliubov form $`E_k=\sqrt{E_r\left(E_r+2gn\right)}`$ , where $`n`$ is the condensate density, $`g`$ is the coupling constant, given by $`g=4\pi \mathrm{}^2a/m`$, $`m`$ is the mass of the atoms, $`a`$ is the scattering length, and $`E_r=\left(\mathrm{}k\right)^2/2m`$ is the corresponding free-particle energy. According to the Feynman relation, the probability of creating an excitation with momentum $`\mathrm{}k`$ is proportional to $`E_r/E_k`$ . Excitations with $`k<\xi ^1`$ , where $`\xi =\sqrt{\mathrm{}^2/2mgn}`$ is the healing length, are called phonons whose dispersion relation is approximated by the linear relation $`E_k\sqrt{gn/m}k`$ .

For an inhomogeneous condensate, we use the local density approximation (LDA) according to which the system behaves locally as a part of a homogenous condensate. The excitation energy is then defined by the average of the local excitation energies over the entire condensate ,

$$E_k=\frac{1}{NS_k}ϵ_k\left(r\right)\frac{E_r}{ϵ_k\left(r\right)}n\left(r\right)𝑑r$$
(1)
where $`ϵ_k\left(r\right)`$ is the excitation energy corresponding to $`n\left(r\right)`$, $`S_k={\displaystyle \frac{1}{N}}{\displaystyle \frac{E_r}{ϵ_k\left(r\right)}}n\left(r\right)𝑑r`$ is the structure factor of the condensate, and $`N`$ is the number of atoms.

In our experiment a nearly pure BEC of approximately 10$`^\text{5}`$ $`^{\text{87}}`$Rb atoms in the $`|F,m_f=|2,2`$ ground state, is formed in a QUIC type magnetic trap . The thermal component of the cloud is measured to be 5% or less. The trap is elongated and cylindrically symmetric, with radial and axial trapping frequencies of $`220`$ Hz and $`25`$ Hz, respectively. The radial and axial Thomas-Fermi (TF) radii of the condensate are thus $`R_{TF}=3`$ $`\mu m`$ and $`Z_{TF}=27`$ $`\mu m`$ . $`\xi ^1=4.3`$ $`\mu m^1`$ averaged over the entire condensate.

The two Bragg beams, produced by the same laser, are detuned 6.5 GHz from the 5S1/2, F=2 $``$ 5P3/2, F’=3 transition. The beams pass through two acousto-optic modulators that control the frequency difference between them. Bragg pulses of duration varying between 1 and 3 msec are applied to the condensate. The angle between the two beams is varied to produce excitations with varying momenta along the $`z`$ axis (to ensure cylindrical symmetry). The beams intensities are adjusted such that the number of atoms in the excitation is always less than 25% of the total number of atoms in the condensate. The frequency difference between the two beams is varied to produce excitations of different strengths. In the range of $`k`$’s measured, the dominant broadening of the Bragg resonance is due to the inhomogeneous density in the condensate , and is always narrower than the finite-time broadening of our pulses. Since the pulse frequency spectrum is broad compared to the excitation line width, excitations are formed at every point in the condensate . Thus, the excitation energy is given, according to Eq.(1), by the average energy added per excitation, i.e. the energy added to the cloud by the pulse, divided by the number of excitations.

After the Bragg pulse, the magnetic trapping potential is rapidly shut off, and repulsion between atoms causes a short acceleration period, during which all the interaction energy is transformed into kinetic energy. After 38 msec of free expansion the cloud is imaged by an on-resonance absorption beam, perpendicular to $`z`$. Figure 1a shows an absorption image for excitations with $`k\xi =0.96`$. Two overlapping clouds are seen. The left cloud corresponds to the condensate, and the right cloud to the released phonons. The right cloud is clearly larger in the radial direction than the left cloud, reflecting the interaction energy of the phonons.

In the TOF image, the atom’s position is proportional to their momentum. Therefore, if the density distribution $`n(r,z)`$ after free expansion is known, we can directly calculate the energy and momentum of the atoms. The absorption picture only provides $`n(r,z)`$ integrated along the absorption beam axis $`f(y,z)=n(r,z)𝑑x`$. However, since the cloud is cylindrically symmetric we reconstruct $`n(r,z)`$ from a single absorption image using computerized tomography. According to the Fourier slice theorem, $`n(r,z)`$ is reconstructed by evaluating

$$n(r,z)=\frac{1}{\left(2\pi \right)^2}F(\kappa _y,z)J_0\left(\kappa _yr\right)𝑑\kappa _y$$
(2)
where $`F(\kappa _y,z)`$ is the 1D Fourier transform of$`f(y,z)`$ along the y direction, and $`J_0`$ is the zeroth order Bessel function. The center of coordinates $`(r=0,z=0)`$ is determined by a Gaussian fit to the BEC cloud. Figure 1b shows the reconstructed $`n(r,z)`$ of Fig. 1a. The released-phonon cloud is seen to have a shell-like shape. A clear separation between the two clouds, which is not visible in Fig. 1a, is now seen. This separation is caused by the released repulsion energy between the two clouds.

From the reconstructed $`n(r,z)`$, the energy of the atoms is calculated as,

$$E=\frac{\pi m}{\tau ^2}n(r,z)\left[r^2+\left(zz_{in}\left(z\right)^2\right)\right]r𝑑r𝑑z$$
(3)
where $`\tau =38`$ msec is the time of flight, and $`z_{in}(z)`$ is the in-trap z position of an atom with position z after TOF ($`r_{in}`$ is approximately $`0`$). For the BEC cloud we take $`z_{in}\left(z\right)=\left(z/Z_{final}\right)Z_{TF}`$, where $`Z_{final}`$ is the parabolic radius in the $`z`$ direction of an integrated parabola, fitted to the absorption image . For the released-phonon cloud, we have no theory for $`z_{in}(z)`$. Thus, we take $`z_{in}(z)=0`$.

First, the energy of an unperturbed BEC is measured. To reduce background we measure the BEC energy within an ellipsoid with radii obtained from a parabolic fit to the absorption image. The average energy per atom is measured in 12 images of unperturbed BEC’s to be $`E_{BEC}/h`$ $`=640`$ Hz, which is approximately equal to the value of $`E_{BEC}/h=590`$ Hz, obtained by the more conventional technique of fitting an integrated parabola to the same absorption images.

The energy of a Bragg-excited BEC is then measured in the combined volume of two ellipsoids, centered around the BEC cloud, and the released-phonon cloud (see Fig. 1b). The size of the second ellipsoid is enlarged to include all the atoms in the released-phonon cloud. A reference image with no excitations is subtracted in the area of the released-phonon cloud to exclude the energy carried by the thermal atoms. The measured momentum of the cloud, divided by $`\mathrm{}k`$, yields the number of excitations. During the pulse the average energy per atom in the condensate cloud $`E_{BEC}`$ decreases by $``$12%. We take the time-averaged value $`E_{BEC}/h=600`$ Hz to be the relevant value.

The filled circles in Fig. 2 are the measured energy of the atoms in both clouds as a function of the measured number of excitations for $`k\xi =0.96`$. The slope of these points is $`E_k/h=1675\pm 290`$ Hz and corresponds to the energy of a single excitation. It is in agreement with $`E_k/h=1780`$ Hz, calculated from the Bogoliubov dispersion relation in the LDA and also with the measured value of  $`E_k/h=1750\pm 20`$ Hz. Figure 3 shows the excitation energies measured this way at four different momenta. The measured energies agree well with the Bogoliubov excitation spectrum for our measured $`E_{BEC}`$ in the LDA with no fit parameters, indicated by a solid line.

In the reconstructed $`n(r,z)`$ of Fig. 1b, the condensate cloud and the released-phonon cloud are distinct and separate. Figure 4a and 4b show the linear density profile in each of the clouds, along the z and the r directions respectively, for $`k\xi =0.96`$. The released-phonon cloud is significantly larger than the condensate cloud in the r direction.

The condensate cloud seems not to be distorted by the interaction with the phonons. This may imply that each phonon is repelled by a fraction of the condensate that is much heavier than the phonon. Since the condensate cloud is much heavier than the phonon cloud, this would suggest that most of the excitation energy is carried by the released-phonon cloud. In order to check this hypothesis, we measure the energy of each cloud separately. The open circles in Fig. 2 are the measured energy in the released-phonon cloud only, as a function of the measured number of excitations. From a linear fit, the slope of these points is $`2390\pm 90`$ Hz. When the average energy of an atom in the condensate cloud is subtracted from this value , we get an energy of $`1840\pm 100`$ Hz per excitation in the released-phonon cloud, which is consistent with the Bogoliubov value $`E_k/h=1780`$ Hz. This indicates that all of the excitation energy is indeed carried by the atoms in the released-phonon cloud.

The filled circles in Fig. 5 are the excitation energies, measured from the slope of the energy of the released-phonon cloud versus the number of excitations. Except for the $`k\xi =0.32`$ point, the measured points agree with the Bogoliubov spectrum in the LDA, indicated by a solid line.

Given that all the excitation energy is carried by the released-phonon cloud, this energy can be extracted from a single image. This technique requires the collection of far less data than does the indirect measurement, where an entire spectrum is needed. The open circles in Fig. 5 are the excitation energies extracted from the image with the largest number of excitations, for each $`k`$. These points agree with the filled circles measured from the slope.

During the release process some of the interaction energy is transformed into motion in the $`z`$ direction, causing the condensate to recoil. The zero of momentum is therefore not at the center of the condensate cloud, introducing a systematic error in our procedure. For the higher $`k`$ values the interaction energy is transformed mainly into motion in the radial direction so this systematic error is negligible. For the smaller $`k`$ values a larger fraction of the energy is transformed into motion in the $`z`$ direction. This may explain the low energies measured at the $`k\xi =0.32`$ point.

We perform hydrodynamic simulations of the release of the excited cloud . Excitations at different initial positions in the condensate are chosen with a probability $`{\displaystyle \frac{E_r}{ϵ_k\left(r\right)}}n\left(r\right)`$ . The condensate is then released and the excitations propagate with an initial momentum $`\mathrm{}k`$. The force at each time step is equal to the gradient of the interaction energy in the condensate at the current position of the excitation, multiplied by a factor, accounting for the excitation’s initial interaction energy. Since the excitations are a small perturbation, the condensate is assumed to expand as would an unperturbed condensate . The calculated average energy of the excitations minus the average energy of the BEC, after 38 msec of free expansion, is shown by the dotted line in Fig. 5. The simulations are seen to be consistent with the observation that the released-phonon cloud carries all of the excitation energy.

In conclusion, the interaction energy of phonons in a Bose-Einstein condensate is directly measured by computerized tomography of TOF images. The excitation energy is extracted from a single absorption image. This measurement of the excitation energy relies solely on images of the atoms, and does not require any information about the photon energy. Furthermore, no assumptions are made with respect to the shape of the cloud, except for cylindrical symmetry. The measured energy is found to agree with the predicted Bogoliubov spectrum in the LDA. After the cloud is released, the excitation energy is seen to be carried by the same atoms that carry the excitation momentum. Hydrodynamic simulations of the release process are performed, and are found to be consistent with our experimental findings.

Computerized tomography of cold atomic clouds could also be used for the observation of non-trivial density profiles, such as the $`s`$-wave scattering shell of two colliding clouds, vortices, solitons, or the fringe pattern of two interfering condensates. In cases of cylindrical symmetry, computerized tomography requires only a single absorption image.


format: text

text:
# Scaling Analysis of Surfactant Templated Polyacrylamide Gel Surfaces

## I Introduction

Polyacrylamide (PAAm) gels are used extensively for separations of biological macro-molecules. Because of their intrinsically broad distribution of pore sizes they are very useful for a wide range of separations for protein isolation and sequencing of single-stranded DNA chiari95 . In particular, protein separations achieved by electrophoresis on gels templated with surfactant micelles have recently been reported rill96 , and the possibility of improvements of the separation properties, when pores of a particular size are templated within the gels, have been proposed rill98 . Extensive work has also been done on the possibility of using such templated gel matrices in gel permeation chromatography (GPC) patterson00 .

Since its introduction in 1986 binnig86 , atomic force microscopy (AFM) has been widely used for recording images of many surfaces with resolution at small length scales. Recent developments in AFM technology have enabled its use, not only in high-resolution profiling of surface morphology and nanostructure, but also in determining local mechanical properties and in surface compositional mapping of heterogeneous samples. These applications have been extensively discussed by Magonov et al. magonov97 ; magonov98 . Despite such developments, there have been only a few reports in the literature on the direct observation of polymer gel surfaces by AFM suzuki96 ; suzuki97 ; maaloum98 . Here, an AFM study of the surface morphologies of templated PAAm hydrogels and their scaling properties are presented.

We use micelles of tetradecyltrimethyl ammonium bromide (TTAB) surfactant to template pores in PAAm. Aqueous solutions of surfactant, acrylamide, and N’,N’-methylene bisacrylamide (cross-linker) are mixed at various concentrations that allow the surfactant molecules to self-assemble into micelles. The solutions are then polymerized, after which the surfactants are removed by soaking in distilled water to produce templated hydrogels. A schematic representation of the templating process is shown in Fig. 1. In previous applications of TTAB templated PAAm gels for GPC patterson00 , mean pore sizes, exclusion limits, and mechanical integrity were determined. It was shown that templated gels yield improved resolution for the separation of globular proteins from 10,000 to 150,000 Da molecular weight. Electrophoretic separations achieved on templated PAAm were also shown to be consistent with the pore-size distributions observed by GPC analysis. Templated PAAm gel matrices have the potential to improve separation efficiency over a wide range of globular protein sizes.

In this paper we present a quantitative study of the structural changes due to the templating process, using scaling analysis of AFM images of the surfaces of templated and untemplated gels. A wide variety of image-analysis tools are available with the Digital Instruments AFM software crmdi3000 . However, these are not well suited for scaling analysis of soft surfaces. Hence, we investigated the surface structure over a wide range of TTAB concentrations using our own, custom-designed surface-analysis software. This scaling analysis appears to be a robust tool to analyze and characterize surfaces imaged by AFM.

The rest of this paper is organized as follows. In Sec. II we describe the experimental procedure, apparatus, and results. In Sec. III we discuss the scaling-analysis methods. This section is subdivided into two-dimensional (Sec. III.1) and one-dimensional (Sec. III.2) methods. In Sec. IV we discuss the results of the scaling analysis of the AFM images, and in Sec. V we present our conclusions.

## II Experimental

PAAm gels were formed at room temperature by free-radical polymerization with chemical initiators from a mixture of monomer, cross-linker, buffer, and surfactant. The gels were cast between two glass slides separated by $`1.5`$ mm spacers. The glass slides were ordinary precleaned microscope slides (Fisher Scientific). They were cleaned with detergent and rinsed in distilled water and air dried before use. AFM images (not shown here) of these slides revealed an root-mean-square (rms) roughness of less than a nanometer over several square micrometers. Acrylamide (Electrophoresis Grade, Fisher Scientific) solutions were polymerized in the presence of a TBE buffer (45 mM tris-borate, 1 mM EDTA) at pH $`8.3`$. The chemical initiators were ammonium persulfate (Ammonium Peroxydisulfate, Certified A.C.S., Fisher Scientific) and N, N, N’, N’-tetramethylethylenediamine (TEMED, Fisher Scientific). The concentration of acrylamide was kept constant at $`40\%`$ by weight, and that of the cross-linker, N’,N’-methylene bisacrylamide (Bis-acrylamide, Electrophoresis Grade, Fisher Scientific) was held at $`7\%`$ of the acrylamide concentration in all samples. Various amounts of the surfactant, TTAB (Sigma-Aldrich), were added to form the templated gels. The amount of surfactant added was so that the final surfactant concentration ranged from $`0\%`$ to $`40\%`$ by weight in steps of $`10\%`$. The presence of surfactant micelles did not prevent the formation of PAAm gels. The rates of polymerization were, however, altered by the presence of surfactants. Higher concentration of surfactants in the pre-gel solution reduced the rate of acrylamide polymerization. All pre-gel solutions were clear, homogeneous solutions with viscosities dependent on the amount of surfactant. The resulting templated gels were all optically clear before surfactant removal.

After polymerization, the gels were peeled off the glass plates and soaked in distilled water to remove the surfactants by diffusion. More than $`98\%`$ of the surfactants diffused out of the gel upon soaking for about two days. Several methods, including Raman spectroscopy have been used to quantify the removal of the surfactant templates by free diffusion into water patterson00 . During the surfactant-removal process, gels with surfactant concentrations of more than $`20\%`$ transformed into white homogeneous materials. This suggests a structural change in the templated gels, involving a length scale comparable to the wavelength of visible light. This further indicates that the structural change is associated with kinetic restructuring of the gel during and after the removal of micelles.

The gel surfaces were imaged using AFM (Digital Instruments, D3000) with a silicon nitride probe (Digital Instruments) of spring constant 0.32 N/m. Scans were performed in contact mode under water (HPLC grade). The set-point voltage was kept at the lowest possible value to minimize the risk of damage to the gel. The inherent soft nature of the gels under investigation can cause difficulties in imaging and affect the reproducibility of the images. In order to overcome this difficulty, we devised a sample holder which helped provide the required stability and hence improved the reproducibility of the images. This sample holder consists of a stainless-steel block with a well in the center, into which the sample is placed. The sample is then covered using a thin steel cover slip with a hole in the center for imaging. The cover slip is held in place with mounting screws (see Fig. 2).

After the surfactants were removed by soaking, three different locations on each gel surface were selected for imaging, and a set of centered, zoomed images were obtained for each location. The image sets were created by scanning the surface from small scan sizes up to larger sizes, in sequence. The images were all $`L\times L`$ squares with sides $`L`$=1.25, 2.5, 5.0, 10.0, and 20.0 $`\mu `$m, respectively. This process yielded a set of successively zoomed-out images of the same location.The lateral drift of the sample was very small. Figure 3 shows one such zoomed set, and the excellent reproducibility of the images justifies the use of the cover slip holder. To further test the reproducibility, drift, and surface distortion, two images were obtained for each scan size, first in the forward scanning direction, and next in the reverse scanning direction, before changing to the next scan size (see Fig. 4). Some distortion along the direction of the tip motion is noticeable. All images have $`512\times 512`$ pixel resolution, regardless of scan size, and all imaging was performed at room temperature.

This telescoping scanning procedure produced a set of ten images for each location and a total of thirty images for each gel. With this procedure, five different gels were examined with surfactant concentrations of 0, 10, 20, 30, and 40$`\%`$, yielding a total set of 150 images (75 forward scans and 75 reverse scans). Many of the images were too noisy for analysis, and all 150 images were examined for quality in a blind procedure in which two of the authors, who had not themselves acquired the data, were separately shown all 150 images. They indicated whether each image should be retained for analysis. If an image was indicated as too noisy by either individual, that image was excluded from analysis. This procedure yielded a total of 106 images for analysis, as detailed in Table 1.

The gel surfaces were also imaged using tapping mode under water. A single set of ten images as described above was obtained for gels with surfactant concentrations of $`0\%`$ and $`40\%`$. The surface morphologies were statistically indistinguishable from those obtained using contact mode. Since tapping mode is rather time consuming and did not provide any apparent advantage over contact mode, all the data presented here were collected in contact mode.

Figure 5 shows examples of AFM images ($`L=10.0\mu \mathrm{m}`$) of gel surfaces with each of the five different template concentrations. The gels exhibit marked statistical differences, most notably that the height range of the surface depends on both the surfactant concentration and the length scale of the scan. This is clearly seen in the three-dimensional visualizations of the gel surface, shown in Fig. 6. Typical scan lines for gel surfaces at various surfactant concentrations are shown in Fig. 7. On the 10 $`\mu `$m length scale, the untemplated gels have height ranges on the order of a few nanometers, while the 40$`\%`$ templated gels have height ranges on the order of a few hundred nanometers. AFM images and Raman spectroscopy of glass plates before and after polymerization and gel removal indicated that there was no polymer left behind on the glass. This suggests that the gel surfaces imaged by AFM are relaxed, free surfaces, rather than fracture surfaces.

The surface analysis routines available with the D3000 only quantify changes in the overall surface width (often called “the surface roughness”). In contrast, our statistical analysis shows self-affine scaling behavior at small length scales, which crosses over to a length-scale independent behavior at large length scales, as discussed in Sec IV. As we are interested in the behavior of the surface on a wide range of length scales, no image filtering or processing was performed, with the exception of flattening (described in Sec. III.1). It is our opinion that scaling analysis of filtered images could be highly suspect since the essence of most image processing techniques is to alter the data differently on different length scales. All analysis was performed with our own code, written in C, as discussed in Sec. III.

## III Analysis Methods

We here present a well-defined procedure to characterize surfaces, based on the statistical analysis methods described in Refs. barabasi95 ; yang93 . One common measure of surface morphology is the rms height of the surface, which is sometimes referred to as the “rms roughness” (referred to here as the rms surface width),

$$w=\sqrt{\frac{_{i=1}^{M^2}h_i^2}{M^2}\left[\frac{_{i=1}^{M^2}h_i}{M^2}\right]^2},$$
(1)
where $`h_i`$ is the $`i`$th value of the height in an image consisting of $`M\times M`$ pixels. (For all of our images, $`M=512`$.) However, such a drastic simplification of the surface structure to only one single number ignores nearly all of the morphological information contained in the image. Reporting the surface width without specifying further details can be meaningless. The value of the width is highly dependent, not only on the image processing methods used fang97 , but also on the length scale of observation iwasaki93 . We therefore present a comprehensive scaling analysis of the gel surfaces which goes beyond the simple single $`w`$ description.

The surface width defined by Eq. (1) for a $`L\times L(\mu \mathrm{m}^2)`$ image of $`M\times M`$ pixels varies with $`L`$. For a variety of different surfaces the width is found to exhibit different length dependences with two limits:

$$w\{\begin{array}{cc}L^\alpha \hfill & Ll_\times \hfill \\ w_{\mathrm{sat}}\hfill & Ll_\times \hfill \end{array},$$
(2)
where the exponent $`\alpha `$ is the Hurst exponent (in this study equal to the roughness exponent meakin98 ), $`w_{\mathrm{sat}}`$ is a constant saturated surface width, and $`l_\times `$ is a cross-over length scale between the two behaviors. In general, $`\alpha `$, $`l_\times `$, and $`w_{\mathrm{sat}}`$ depend on the details of the surface, such as the surfactant concentration in the present study.

The power law in Eq. (2) is a consequence of self-affine scaling behavior, in which the surface is statistically invariant under the transformation barabasi95

$$\begin{array}{ccc}xbx,& yby,& hb^\alpha h,\end{array}$$
(3)
where $`x`$ and $`y`$ are distances in the two image-scan directions, $`h`$ is the height, and $`b`$ is a dimensionless re-scaling factor. Self-affine behavior is different from self-similar behavior, which involves invariance under the isotropic scale-change transformation barabasi95  $`xbx`$, $`yby`$, and $`hbh`$. This distinction between self-affine and self-similar structures is often overlooked, and the term “fractal” is used somewhat freely for both in the literature.

If the surface were ideally self-affine, then the power-law behavior in Eq. (2) would continue on all length scales, and there would be neither a saturation width, $`w_{\mathrm{sat}}`$, nor any identifiable length scale such as $`l_\times `$ or a typical feature size. For a surface whose heights are described by an ideal random walk, $`\alpha =1/2`$ and $`l_\times \mathrm{}`$. For most surfaces, measured values of $`\alpha `$ are between $`1/2`$ and $`1`$ maloy92 ; bouchaud93 ; delaplace99 , even for such different surfaces as brittle rocks ($`\alpha 0.78`$brown85  and CuCl islands on CaF2 ($`\alpha 0.84`$tong94 . As discussed in Sec. IV, our estimates of $`\alpha `$ for the different gel surfaces all fall within the range $`[0.5,1.0]`$, and $`w_{\mathrm{sat}}`$ and $`l_\times `$ vary strongly with the surfactant concentration.

In the following subsections, we discuss several different methods used to examine the scaling behavior of the gel surfaces. Both two-dimensional analysis methods that analyze the image as a whole, and one-dimensional methods that analyze the image one scan line at a time were examined.

### III.1 Two-Dimensional Methods

The AFM records the vertical surface height, $`h(x,y)`$, as a function of the fast ($`x`$) and slow ($`y`$) scanning directions. As already discussed, $`L`$ dependent surface widths can be obtained from $`L\times L(\mu \mathrm{m}^2)`$ images. After imaging at a sufficiently large number of scan sizes (different $`L`$’s), $`\alpha `$, $`w_{\mathrm{sat}}`$, and $`l_\times `$ can be obtained. However, the number of images required for such an analysis is large. By noting that each image contains information about length scales from $`L/M`$ (where $`M`$ is the number of pixels) up to $`L`$, one can perform analyses over a wide range of length scales from a single AFM image. Several variations of two-dimensional analysis methods were studied. We discuss a spectral method based on the Fourier transform of the image yang93 , a box-counting method from subsections of images, and widths calculated from the same subsections of the image barabasi95 . All of the images were flattened using the AFM software crmdi3000  and not processed any further. The flattening procedure fits a straight line to each scan line by the method of least squares, which is then subtracted from the data. The fitting and subtracting procedure was applied, independently, to each scan line in the image. After flattening, the average height and slope in the $`x`$-direction (the fast scan direction) are zero.

For the spectral methods, we computed the structure factor, $`S(\stackrel{}{k})|_\stackrel{}{k}(h(x,y))|^2`$, which is the absolute square magnitude of the complex Fourier transform of the surface height. Here, $`\stackrel{}{k}`$ is the wave vector in the $`x,y`$ plane. For large $`k`$, the structure factor should scale as $`S(k)k^{(2+2\alpha )}`$ yang93 , a generalized Porod’s Law porod83 . The presence of anisotropy in the AFM images produces ridges along the $`k_x`$ and $`k_y`$ directions as shown in Fig. 8. This made circular averaging unreliable, and hence the results from spectral methods were suspicious. The presence of the ridges seems to be a general feature of AFM imaging, as has been observed for hard materials as well fang97 . Removal of these ridges (which have a non-zero width) by a technique such as that suggested in Ref. fang97  would require a priori knowledge of their functional form.

For both the box-counting and width methods, the $`L\times L(\mu \mathrm{m}^2)`$ image of $`M\times M`$ height values (pixels) was recursively subdivided into non-overlapping images of $`l\times l(\mu \mathrm{m}^2)`$ with $`(M/n)\times (M/n)`$ pixels, where $`l=L/n`$ $`\mu `$m and $`nM/2`$ is an integer denoting the level of subdivision. For each value of $`l`$, $`(L/l)^2`$ different sub-images were obtained. Non overlapping sub-images were used to ensure statistical independence.

A box-counting method was tried following Ref. barabasi95 . The minimum number of cubic boxes of side $`l`$ needed to span the height range of a sub-image of $`l\times l\mu \mathrm{m}^2`$ was calculated. However, the height range was nearly always less than $`l`$, yielding the total number of boxes to be one, regardless of the length scale. This would seem to imply $`\alpha =1`$. Thus the box-counting method is entirely insensitive to the scaling behavior of the gel surfaces.

The rms width for the $`j`$-th sub-image of size $`l\times l`$,

$$w_{lj}=\sqrt{\frac{_{i=1}^{(M/n)^2}h_i^2}{(M/n)^2}\left[\frac{_{i=1}^{(M/n)^2}h_i}{(M/n)^2}\right]^2},$$
(4)
and its average, $`w_l`$, over all sub-images of size $`l`$ were calculated. Typical results are shown in Fig. 9. For a discussion of the effects of changing the order of the averaging and the square root, see Ref. BUEN02 .

Among the several two-dimensional analysis methods used to study the scaling properties of the gel surfaces, the calculation of $`w_l`$ was found to provide the most reliable results. The other methods have inherent disadvantages as discussed in the preceding paragraphs. As we shall explore in the next subsection, the scanning process introduces noise between the scan lines, and for a more accurate analysis, one-dimensional methods should be used.

### III.2 One-Dimensional Methods

If there is noise in the measurement process we would expect to measure $`h_{\mathrm{measured}}(x,y)=h(x,y)+h_{\mathrm{noise}}`$, where the noise adds an effective offset. If we further assume that $`h_{\mathrm{noise}}`$ is a random walk in time, noise correlations in time would appear as correlations in space. Furthermore, if the noise varies sufficiently fast compared to the scanning velocity, $`h_{\mathrm{noise}}`$ at the beginning of each scan line will be completely uncorrelated. Thus for a perfectly flat surface with $`h(x,y)=0`$, each scan line would appear to be a random walk with $`\alpha =1/2`$, with little or no correlation between scan lines. Motivated by these considerations and the extremely small height ranges in the gels with 0$`\%`$ and $`10\%`$ surfactant, we constructed a one-dimensional analysis procedure which analyzes the images one scan line at a time, where a scan line is a line in the fast ($`x`$) scan direction.

One-dimensional surface widths $`w_l`$ were calculated in a manner equivalent to that defined in Eq. (4) with $`j`$ now representing a sub-interval in the fast scan direction instead of a sub-image and $`(M/n)^2`$ replaced by $`(M/n)`$. The $`w_{lj}`$, for all $`j`$ and all scan lines, were then averaged to get the $`w_l`$ for the image. The surface width calculations at the smallest image sub-interval depends strongly on the definition of $`w_l`$ (see Eq. (4)). When applied to a straight line ($`\alpha =1`$), the standard one-dimensional estimate of $`w_l`$ used here overestimates $`\alpha `$ by as much as $`20\%`$ for sub-intervals that contain less than four data points. The results of the one-dimensional surface width analysis are presented in Table 2 and in Sec. IV.

While surface-width calculations contain information about all length scales $`L`$, calculation of correlation functions yield information at a specific value of relative separation. We define the height-height correlation function of a scan line at fixed $`y`$ as

$$C(r,y)=h(x,y)h(x+r,y)_xh(x,y)_x^2,$$
(5)
where $`r`$ is a relative displacement in the $`x`$ direction, and $`_x`$ denotes an average over all values of $`x`$ in the scan line. Note that after flattening, the second term is zero, and that $`C(0,y)`$ is the rms width of the scan line, $`w(y)`$. Although the flattening was necessary to compensate for non-horizontal mounting of the gel and/or gradients in the gel thickness, it introduces correlations at the length scale $`L`$. Consequently, correlation functions for flattened images are meaningful only for lengths, $`lL/2`$. For scaling analysis, it is more convenient to deal with an alternative form of the correlation function, which is sometimes referred to as the variance of increments feder88  or the increment correlation function,

$$H(r,y)=|h(x+r,y)h(x,y)|^q_x.$$
(6)
For the special case of $`q=2`$,

$$H(r,y)=[h(x+r,y)h(x,y)]^2_x=2[w(y)^2C(r,y)].$$
(7)
The correlation functions in the slow ($`y`$) scan direction, $`C(x,r)`$ and $`H(x,r)`$, can be analogously defined.

To improve the statistical accuracy of the increment correlation function measurements, multiple scan lines must be used. Thus, we define the first and second moments as

$$\begin{array}{c}H_{\mathrm{fast}}(r)=H(r,y)=_{j=1}^p\left[_{i=1}^MH(r,y=iL/M)\right]/pM,\hfill \\ H(r,y)^2=_{j=1}^p\left[_{i=1}^MH(r,y=iL/M)^2\right]/pM,\hfill \end{array}$$
(8)
where $`_{j=1}^p`$ runs over the $`p`$ different scan lines. The standard error of $`H(r,y)`$ is then $`\sigma _{\mathrm{fast}}=\sigma (r,y)=\sqrt{\frac{H(r,y)^2H(r,y)^2}{pM1}}`$. Again, the moments and error in the slow scan direction can be analogously defined.

The increment correlation function for $`q=2`$ is expected to scale as

$$H_{\mathrm{fast}}(r)\{\begin{array}{cc}r^{2\alpha }\hfill & rl_\times \hfill \\ 2w_{\mathrm{sat}}^2\hfill & rl_\times \hfill \end{array}.$$
(9)
In Fig. 10 we see power-law scaling for small $`r`$ and an apparent cross-over to saturated values at large $`r`$. When the noise-induced offsets between scan lines are small, $`H_{\mathrm{slow}}(r)`$ should have the same scaling as $`H_{\mathrm{fast}}(r)`$. However, if the signal-to-noise ratio is small, we expect $`H_{\mathrm{slow}}(r)`$ to measure only the noise with $`\alpha =1/2`$. Figure 11 shows plots of $`H_{\mathrm{slow}}(r)`$ and $`H_{\mathrm{fast}}(r)`$ for $`0\%`$ and $`40\%`$ templated gel surfaces. For the untemplated gel surfaces the fluctuations due to noise in the slow scan direction are comparable to the signal. Thus $`H_{\mathrm{slow}}(r)`$ does not follow $`H_{\mathrm{fast}}(r)`$. For the $`40\%`$ templated gel surfaces, the random offsets between scan lines are relatively smaller, but not negligible. Hence $`H_{\mathrm{slow}}(r)`$ is much closer to $`H_{\mathrm{fast}}(r)`$.

By evaluating the increment correlation function and plotting $`H_{\mathrm{fast}}(r)`$ on a log-log scale, the roughness exponent $`\alpha `$, the limiting surface width $`w_{\mathrm{sat}}`$, and the lateral cross-over length $`l_\times `$ can be determined from the slope in the small-$`r`$ linear regime (see Fig. 10), the limiting value of $`H_{\mathrm{fast}}(r)`$ for large $`r`$, and the inflection point of the logarithmic derivative of $`H_{\mathrm{fast}}(r)`$ (see Fig. 12), respectively.

## IV Results and Discussion

As shown in the previous sections, among the two-dimensional analysis methods, only the $`l`$-dependent surface width $`w_l`$ gives meaningful results for this system. The $`l`$-dependence of this quantity was presented in Fig. 9 and is essentially similar to that of the one-dimensional surface-width (not shown) and the one-dimensional increment correlation function shown in Fig. 10. All three quantities show a self-affine scaling region on short length scales, and a saturated plateau on large length scales. These two scaling regimes are separated by the crossover length $`l_\times `$. Both the crossover length and the saturation value increase with increasing surfactant concentration.

The values of the roughness exponent $`\alpha `$ extracted from the self-affine regime for the surface width measurements (both two-dimensional and one-dimensional), along with those calculated from the one-dimensional increment correlation function with $`q=2`$, are listed in Table 2. Using $`q2`$ in Eq. (6) yielded the same values of $`\alpha `$ as obtained with $`q=2`$. This indicates that the gel surfaces are perfectly self-affine on these short length scales, rather than multi-affine meakin98 ; BUEN02 ; BARA91 . A plot of the estimated values of $`\alpha `$ vs the surfactant concentration is shown in Fig. 13. The estimates obtained from the two-dimensional surface width measurements and from the one-dimensional increment correlation function agree well for $`20\%`$, $`30\%`$, and $`40\%`$ templated surfaces. We believe that the lack of agreement for $`0\%`$ and $`10\%`$ templated surfaces is associated with the low signal-to-noise ratio in the images for these gels. The vertical resolution limit of the AFM ($`0.1\mathrm{nm}`$ crmdi3000 ) contributes to the uncertainty since the height variations are only about 1 nm for the lowest surfactant concentrations.

The values for $`\alpha `$, determined from the one-dimensional surface-width measurements, are consistently higher than those obtained from the other two measurements. As discussed in Sec. III.2, the roughness exponents obtained using the one-dimensional surface width for intervals with only a few data points are consistently too high, and this may be at least part of the explanation for the discrepancy. The values of $`\alpha `$ for $`20\%`$, $`30\%`$, and $`40\%`$ templated surfaces obtained from two-dimensional surface-width and one-dimensional correlation-function measurements are clustered between 0.80 and 0.89. However, measurements of surface-roughness exponents from relatively small data sets, such as the present ones, tend to have systematic errors, as well as random ones, and the calibration curves for synthetic data sampled over 512 points, given in Ref. CAST01 , indicate that the values obtained here may be underestimated by up to $`10\%`$. Thus, based on our results it is not unreasonable to conclude that all of these surfaces have the same roughness exponent in the range 0.85 to 1.0, irrespective of the surfactant concentration. A roughness exponent in this range is consistent with previous results for N-isopropylacrylamide gel surfaces by Suzuki et al. suzuki97 . However, as discussed in the next paragraph, the magnitude of the height variations depends strongly on the surfactant concentration, as do the crossover length $`l_\times `$ and the saturation surface width $`w_{\mathrm{sat}}`$.

The limiting value of the surface width on large length scales, $`w_{sat}`$, is found to increase by almost two orders of magnitude as the surfactant concentration is increased from 0% to 40%. This behavior is evident from Figs. 7, 9(a), 10(a), and 14. The measurements of $`w_{sat}`$ are listed in Table 2 and displayed vs surfactant concentration in Fig. 14. Similarly, but less spectacularly, $`l_\times `$ increases by about a factor of two, as shown in Table 2 and Fig. 15.

As the surfactants are removed, the transparency of some of the gels change. The gels that contained no or little surfactant stay clear. On the other hand, gels that had higher concentrations of surfactant (more than $`20\%`$) become opaque and white as the surfactant diffuses out. Using the theory of light scattering by refractive index fluctuations, an expression for the spectral density of light scattered can be obtained. Reasonable functional forms for the correlation function can then be used to obtain the turbidity or inverse attenuation length. A detailed derivation is provided in the Appendix. The turbidity for a wavelength $`\lambda =550`$ nm in the visible region is shown as a function of the correlation length $`\xi `$ in Fig. 16. The two curves presented correspond to an isotropic exponential and Gaussian correlation function, respectively. If we associate the correlation length $`\xi `$ with the saturation value of the rms surface width, $`w_{\mathrm{sat}}`$, we notice that the turbidity increases by about six orders of magnitude as $`w_{\mathrm{sat}}`$ changes from 1 nm for the $`0\%`$ templated gel to 80 nm for the $`40\%`$ templated gel. These changes in the turbidity are consistent with a change from a transparent to an opaque material. Our justification for the association between $`\xi `$ and $`w_{\mathrm{sat}}`$ is discussed below.

The experimental results presented here indicate that while the templated gel surfaces are self-affine on small to moderate length scales, with $`\alpha `$ in the range 0.85 to 1, both the crossover length and the saturation value of the surface width increase significantly with the bulk concentration of surfactant. From Figs. 14 and 15 it is clear that most of these increases occur only after the surfactant concentration exceeds about 20%. This indicates that the dramatic changes seen in the surface structure are reflections of a sharp transition in the bulk structure, most likely a percolation transition at which the surfactant structure changes from a spatial distribution of individual micelles to a percolating network. These different structures apparently do not seriously affect the mechanical properties of the gel as long as the surfactants are in place. However, once they are soaked out, the resulting percolating void network is likely to weaken the gel and lead to large relaxations as the elastic polymer network readjusts its configuration. While the percolation threshold is a non-universal quantity and thus depends on the particular system, for site percolation on standard three-dimensional lattices it lies roughly in the range from 20 to 30 percent by volume STAU92 . This is consistent with the range in which the surface structure and the opacity of the templated gels undergo their strongest changes. Previous work on the phase diagrams of alkyltrimethylammonium surfactants jonsson88  suggests a transition from micelle to hexagonal phase at around $`35\mathrm{wt}\%`$ TTAB at room temperature.

Further support for the percolation hypothesis is provided by recent computer simulations of a two-dimensional frustrated spring network with a free surface and a variable density of holes BUEN01 ; BUEN02 , which were inspired by the experimental results presented here. In those simulations the structure of the free surface remains relatively constant, with a self-affine structure on short length scales and a finite saturation width on large length scales, as long as the volume fraction of holes is far below the percolation limit. As the hole concentration approaches the percolation threshold, both the the saturation width and the crossover length increase in a manner qualitatively similar to our experimental observations. The saturation width, in particular, appears to follow the same power-law divergence as the bulk correlation length over a wide range of hole concentrations approaching the percolation threshold. These observations form our justification for the association of $`w_{\mathrm{sat}}`$ with the bulk correlation length in the estimate of the change in the turbidity, discussed above.

## V Conclusion

We have obtained AFM images of relaxed, free surfaces of PAAm hydrogels templated with TTAB surfactant and performed a numerical scaling analysis of the observed surface structure. In doing so, we have developed methods that can be used not only for our data, but also to perform scaling analysis of any surface imaged by an AFM. We showed that the calculation of the two-dimensional surface width and one-dimensional correlation functions gave consistent, useful results. Two-dimensional box-counting and spectral methods were unsatisfactory, while results from one-dimensional surface-width measurements provided over-estimates at small length scales. The one-dimensional increment correlation function in the fast scan direction, $`H_{\mathrm{fast}}(r)`$, is simple to calculate and gives a reliable characterization of rough surfaces. By calculating this correlation function we have determined that the surfaces of the templated gels are self-affine over a significant range of small and intermediate length scales, with a roughness (Hurst) exponent in the range 0.85 to 1. This self-affine regime is terminated by a cross-over length, above which the surface width is scale independent. The increases in the cross-over length and the saturation width of the surface with increasing surfactant concentration coincides with the observation of a large increase in the scattering of visible light. We propose that the observed effects are caused by a percolation transition in the bulk surfactant system, from isolated micelles to a percolating network. Further investigations of the bulk structure of the gels at different surfactant concentrations are planned.

## Acknowledgments

We thank R. L. Rill, G. M. Buendía, P. Meakin, and S. Hong for useful discussions, and E. Lochner for his assistance with the AFM and the development of the new sample holder. G. Brown, K. Park, and S. Wang are acknowledged for their comments on the manuscript. We thank B. van de Burgt for the Raman spectroscopy results and analysis, and A. Beheshti is acknowledged for his help in testing the validity of our data-analysis code. Supported in part by NSF grants No. BES-951381 and DMR-9981815, and by Florida State University through the Center for Materials Research and Technology and the School of Computational Science and Information Technology.

## Appendix A

The theory of light scattering by refractive-index fluctuations leads, in the Born approximation, to the following expression for the spectral density of light scattered from a volume $`V`$, a distance $`R`$ away in the far-field zone berne00 :

$$I(\stackrel{}{k})=I_0\frac{V}{R^2}\frac{(2\pi )^4\mathrm{sin}^2\psi }{2\lambda ^4}𝒮^{(\delta n)}(\stackrel{}{k}),$$
(10)
where

$$𝒮^{(\delta n)}(\stackrel{}{k})=_Vd^3rG(\stackrel{}{r})\mathrm{exp}[i\stackrel{}{k}\stackrel{}{r}].$$
(11)
Here $`I_0`$ is the intensity of the incident light, $`\lambda `$ is the wavelength of the scattered light, $`\psi `$ is the polarization angle, $`\stackrel{}{k}`$ is the scattering wave vector, and $`G(\stackrel{}{r})=\delta n(\stackrel{}{r}_0)\delta n(\stackrel{}{r}_0+\stackrel{}{r})/(\delta n)^2`$ is the normalized correlation function for the refractive-index fluctuations $`\delta n=nn`$.

We can calculate the total scattering intensity at a particular wavelength $`\lambda `$ by using the relation $`|\stackrel{}{k}|=\frac{4\pi }{\lambda }\mathrm{sin}(\frac{\theta }{2})`$, where $`\theta `$ is the scattering angle, and integrating $`I(\stackrel{}{k})`$ over all angles, i.e.,

$$I(\lambda )=I_0\frac{V}{R^2}\mathrm{\Lambda }^1(\lambda ,\xi )$$
(12)
where the turbidity or inverse attenuation length is

$$\mathrm{\Lambda }^1(\lambda ,\xi )=𝑑\theta 𝑑\varphi \mathrm{sin}\theta \frac{1}{2}(1+\mathrm{cos}^2\theta )𝒮^{(\delta n)}(\stackrel{}{k}),$$
(13)
and $`\xi `$ is the bulk correlation length. The term in the integral of the form $`(1+\mathrm{cos}^2\theta )/2`$ corresponds to unpolarized incident light johnson94 . Assuming an isotropic exponential correlation function, $`G(r)=\mathrm{exp}[r/\xi ]`$, we find

$$\mathrm{\Lambda }^1(\lambda ,\xi )=\frac{4\pi ^3(\delta n)^2}{\lambda }\frac{a^2+2}{a}\left[\frac{a^2+2}{a^2+1}\frac{2\mathrm{ln}(a^2+1)}{a^2}\right],$$
(14)
where $`a=4\pi \xi /\lambda `$. Alternatively, using a Gaussian correlation function, $`G(r)=\mathrm{exp}[(r/\xi )^2]`$, we find

$$\mathrm{\Lambda }^1(\lambda ,\xi )=\frac{2\pi ^{7/2}(\delta n)^2}{\lambda a^3}\left\{(648a^2+a^4)(64+8a^2+a^4)\mathrm{exp}[a^2/4]\right\}.$$
(15)
A reasonable value for the mean refractive-index fluctuation between the PAAm and water regions in a gel is $`(\delta n)^2=0.01`$. The resulting turbidity, as given by Eqs. (14) and (15), is shown in Fig. 16.


text:
# Resolved nuclear CO(1-0) emission in APM08279+5255: Gravitational lensing by a naked cusp?

## 1 Introduction

Identified serendipitously in a search for high latitude carbon stars, the $`z=3.9`$ broad absorption line quasar APM08279+5255 is coincident with an IRAS source with a flux of 0.95Jy at 100$`\mu `$m (Irwin et al. 1998). Observations with SCUBA reveal that APM08279+5255 possesses a significant submillimetre flux of 75mJy at 850$`\mu `$m (Lewis et al. 1998), implying a bolometric luminosity of $`5\times 10^{15}L_{}`$. Imaging reveals that APM08279+5255 is not point-like, but rather is extended over a fraction of an arcsecond with a structure indicative of gravitational lensing (Irwin et al. 1998). The composite nature of APM08279+5255 was confirmed in adaptive optics (AO) images obtained by Ledoux et al. (1998), with the system appearing as a pair of point-like images separated by 0.4 arcsec. Observations with NICMOS on the Hubble Space Telescope (Ibata et al. 1999) and AO images taken with the Keck telescope (Egami et al. 2000) also uncovered a fainter third component between the brighter two. Gravitational lens models derived from these observations suggest that the quasar continuum source has been magnified by $`90`$.

Using IRAM, Downes et al. (1999) detected emission in CO(4-3) and CO(9-8), revealing the presence of warm circumnuclear gas in APM08279+5255. Papadopoulos et al. (2001) were able to search for CO(1-0) and CO(2-1) in this system using the Very Large Array (VLA). Both were clearly detected associated with the quasar nucleus, as well as a more extended component located several arcsecs from the quasar images. Using locally established values of the CO-to-H2 ratio, this lone cloud represents $`10^{11}M_{}`$ of cold and/or subthermally excited gas.

In this paper, we present an analysis of nuclear CO(1-0) emission in APM08279+5255 using VLA 111The VLA is operated by the National Radio Astronomy Observatory, which is a facility of the National Science Foundation, operated under cooperative agreement by Associated Universities, Inc. at high spatial resolution (0.3 arcsec). The CO appears as a partial ring of $``$0.6 arcsec diameter. These data suggest a total revision in the gravitational lens model for this source, with the new model involving a ‘naked cusp’, which naturally accounts for the observed odd-number of images. They also imply that the nuclear CO must be spatially extended on a scale of at least 400 pc, making this the second source in which gravitational lensing can be used as a ‘telescope’ to explore sub-kpc scale structure of molecular gas in the AGN host galaxy.

## 2 Observations and Results

Observations of the CO(1-0) emission from APM08279+5255 were made in March and April 2001 using the VLA in the B (10 km) configuration. A total of 20 hours were spent on the source. Due to limitations with the VLA correlator, we chose to observe in continuum mode using two 50 MHz bandwidth IFs, each with two polarizations. One IF was centered on the redshifted CO(1-0) line, corresponding to a frequency of 23.465GHz, while the second IF was tuned away from the emission line to a frequency of 23.365GHz. This observing set-up maximizes sensitivity since the effective bandwidth (45 MHz, 575 km/s) closely matches the obseved CO linewidth (Downes et al. 1999), but sacrifices velocity information (Carilli, Menten, & Yun 1999). The source 3C286 was used for absolute gain calibration. The rms noise on the final image is 35$`\mu `$Jy beam-1 with an effective spatial resolution of FWHM=$`0.39^{\prime \prime }\times 0.28^{\prime \prime }`$ with a major axis position angle of $`70^o`$. Fast switching phase calibration was employed with a calibration duty cycle of 130 seconds (Carilli & Holdaway 1999). The phase stability for all the observations was excellent, such that the array was easily phase coherent over the calibration cycle time. This was demonstrated by imaging the phase calibrator (0824+558) with a similar calibration cycle time, and by imaging the radio continuum emission from APM08279+5255 itself.

The images of the line and continuum emission from APM08279+5255 are shown in Figure 1. The CO line image was generated by subtracting from the visibility data a Clean-Component model of the continuum emission made from the off-line data. The peak surface brightness for the continuum emission is $`0.24\pm 0.035`$ mJy beam-1, with a total flux density $`0.41\pm 0.07`$ mJy. The corresponding numbers for the line emission are $`0.183\pm 0.035`$ mJy beam-1 and $`0.39\pm 0.09`$ mJy. The latter corresponds to a velocity-integrated line flux of $`0.22\pm 0.05`$ Jy km/sec, consistent with $`0.15\pm 0.045`$ Jy km/sec deduced for the inner $`1`$arcsec by Papadopoulos et al. (2001)

The continuum emission is extended north-south, with a position angle and spatial extent as expected based on the three component optical gravitational lens. An attempt was made to decompose the continuum map into three point-like images, centred upon the positions derived from the analysis of HST images (Ibata et al. 1999). Given the image resolution, it was not possible to cleanly separate the north-most images (A+C) which have a separation of $``$0.15 arcsecs. The relative flux of (A+C)/B$``$1.5, quite similar to the optical ratio presented by Ibata et al. (2001). The line emission, on the other hand, shows distinct curvature away from the continuum position angle, and at the 3$`\sigma `$ surface brightness level it appears as an almost complete ring with a diameter of about 0.6 arcsec. While the map of the line emission is of low signal to noise, the ring-like structure is apparent in each of several days worth of observations. Also, its size is in agreement with that of the CO(4-3) and CO(9-8) emission (Downes et al. 1999) and the CO(2-1) emission (Papadopoulos et al. 2001). Hence, these observations reveal the gross features of the CO(1-0) emission, although more observations are required to uncover the finer details.

The observed CO(1-0) brightness temperature of the ring (averaged over $``$575 km/sec) is 1.4$`\pm `$0.8 K, which corresponds to an emitted brightness temperature of $`T_b`$7 K at z$``$3.9. A lower limit for the magnification factor of the CO(1-0) emission can be derived assuming the warm gas emitting the high-J CO lines (Downes et al. 1999) to be also emitting the J=1-0. This gas phase is optically thick with $`T_{kin}200K`$, in agreement with the inferred dust temperatures (Lewis et al. 1998), which then yields a velocity-averaged filling factor of $`f7/200=0.035`$. Since differential lensing will “boost” a compact warm region at the expense of a more extended and possibly sub-thermally excited gas phase that emits the CO J=1-0, the true value of f will be larger.

## 3 Gravitational Lensing

The lensing galaxy has yet to be identified in APM08279+5255. Hence, only the quasar positions and magnitudes are available to constrain any lensing model. Other than APM08279+5255, all other lens systems possess an even-number of images. Theoretically, any non-singular mass distribution should produce an odd number of lensed images (Burke 1981), and the ubiquity of even numbers of images has been used to limit the core radius in lensing systems, as small cores result in the demagnification of one of the images. Due to the brightness of the central source, however, models for APM08279+5255 have required the opposite, a very circular model with a large core / shallow cusp (Lewis et al. 1999; Egami et al. 2000; Munoz et al. 2001). We further explore the lens model of APM08279+5255 in light of the observations presented here.

The concordance between the optical images of APM08279+5255 and the radio continuum at 3.6cm (Ibata et al. 1999) and 23GHz demonstrates they come from coincident regions with a similar scale size. Associated with the active core of APM08279+5255, these regions are smaller than the gravitational lensing caustics. The quite different morphology displayed in the CO image, however, indicates that this emission arises in a larger region, distinct from the continuum radiation. With its extended nature, the CO emitting region can lie under different parts of the caustic network. The CO source, however, can not be significantly more extended than the caustic network, as the resulting image would not show the ring structure seen in Figure 1. Therefore, for a given model, the size of the ring image provides a probe of the scale size of the emitting region.

To form ring-like images, the source must be extended and cover a substantial fraction of the inner caustic structure of the lens; the resulting image is appears as a ring, following the outer critical line in the image plane (see Kochanek, Keeton & McLeod 2001). Examining the critical line structure in all previously published models for APM08279+5255 (Ibata et al. 1999; Egami et al. 2000; Munoz et al 2001), the resulting structure for the Einstein ring of a source centred upon the quasar nucleus should be quite circular, with a radius of $`0.2`$ arcsecs, passing through the two brighter quasar images; such images for extended source can be seen in the models of Egami et al. (2000). This is quite different to the CO structure displayed in Figure 1, with the CO emission clearly showing a roughly east-west extension, with the hole in the ring occurring $`0.5`$ arcsecs from image A. Using the model derived from the HST data (Ibata et al. 1999), we explored the image configurations a range of sizes and shapes, centred upon the quasar source, for the CO emitting region. None reproduced the observed image structure and we reject this previous model.

Naked cusps occur in highly elliptical systems, such as flattened disks, when the inner diamond caustic extends outside the elliptical caustic. A source inside this extension produces three roughly colinear images of similar brightness (Maller, Flores & Primack 1997; Bartelmann & Loeb 1998; Keeton & Kochanek 1998; Moller & Blain 1998; Blain, Moller & Maller 1999; Bartelmann 2000). While gravitational lensing by spiral galaxies has been observed (e.g. B 1600+434; Koopmans, Bruyn & Jackson 1998) no observed lensed quasar system has so far been associated with a naked cusp.

To examine the possibility that the observed image configuration is due to gravitational lensing by a naked cusp, we constructed a simple mass model. Following Batelmann and Loeb (1998), this consists of a truncated flattened disk in a spherical halo. APM08279+5255 has brightened considerably since its discovery (Lewis, Robb & Ibata 1999; E. Ofek priv. comm., see http://wise-obs.tau.ac.il/$``$eran/LM), potentially due to the effects of gravitational microlensing. Comparing the images obtained in 1998 (Ibata et al. 1999; Egami et al. 2000) and those obtained in 1999 (Munoz et al. 2001), the relative image brightness have changed appreciably, with image B changing from being 78% of image A to only 50% in the latter epoch; such behaviour is extremely suggestive of the action of gravitational microlensing, although longer term monitoring is required to confirm this. Hence, the relative image brightnesses cannot be used to constrain any mass model. Given the sparsity of constraints, we choose to find a model that can recover the image positions, while providing a reasonable description of the observed CO emission. As can be imagined, the parameter space is large, so a range of models that reproduced the quasar positions were chosen and then modeling of the CO emission was undertaken by-eye. With this, therefore, we do not claim that the model presented here is unique, only consistent with the general form of data.

In our chosen model, the disk is highly inclined, presenting a projected axis ratio of 0.25, and the disk truncated at 8$`h^1`$kpc 222$`\mathrm{\Omega }=1`$ and $`\mathrm{\Lambda }=0`$ is assumed throughout, and possess a core radius of 0.065$`h^1`$kpc. The rotation velocity of the disk is 200km$`s^1`$. Figure 2 presents the source and image plane for this model, with the elliptical and diamond caustic, and corresponding critical lines apparent. In this model, the quasar images are not substantial magnified, with a total magnification of $`7`$, with the intrinsic source of APM08279+5255 being correspondingly luminous, $`L_{bol}7\times 10^{14}L_{}`$. While extreme, this value is not necessarily outrageous as the unlensed quasar HD1946+7658 possess an intrinsic luminosity of $`4\times 10^{14}L_{}`$ (Hagan et al. 1992), and APM08279+5255 may be a member of this very luminous class of quasars. It must be conceded, however, that the non-uniqueness of the lens model translates into uncertainty in the model magnification and a true determination of the intrinsic properties of APM08279+5255 require models derived from better observational constraints.

The CO source is taken to be have an circular surface brightness distribution centred upon the quasar position. One important aspect of the results presented herein is that APM08279+5255 becomes the second system for which the gravitational lens can be used to study structure on sub-kpc scales in the molecular gas associated with the AGN host galaxy, the first system being the Clover Leaf quasar, H1413+117 (Yun et al. 1998; Kneib, Alloin, & Pello, R. 1998). For APM08279+5255, the observations and lens model require the CO to be distributed on a scale covering a substantial fraction of the caustics in the image plane, but not too large to lose the ring structure. The lower limit to the CO source size based on the modeling is $`400h^1`$ pc, while a rough upper limit is $``$1 kpc. With this model, the CO(1-0) has been magnified by a factor of $`2.53`$. Like the Clover Leaf, we find that the spatial extent and mass of the molecular gas in APM08279+5255 are comparable to those seen in nearby nuclear starburst galaxies (Sanders and Mirabel 1996; Downes and Solomon 1998).

Downes et al. (1999) determined a CO source size of $`(80135)h^1`$pc for the estimated magnification factors of $`207`$. This size is much smaller than the one calculated above, while their magnification factors are larger. Their analysis is based upon modeling of CO emission in the gravitationally lensed ultraluminous galaxy IRAS F10214+4724 (Downes, Solomon & Radford 1995), whose image is clearly an extended arc-like feature which possesses an essentially linear magnification. Such a simple model is probably a poor representation of the lensing in APM08279+5255. Additionally, Downes et al. (1999) assumed that the velocity filling factor is unity, substantially larger than the value derived in Section 2; as the intrinsic source radius in their model scales inversely with this value and the magnification factor is proportional to it. For f$``$0.35, a value well within our estimated range, the Downes et al. (1999) model yields an upper limit for the intrinsic source size and a lower limit for the magnification factor that our similar to ours.

Accounting for the influence of gravitational lensing, the velocity-integrated CO(1-0) flux density implies that the nuclear content of molecular gas in APM08279+5255 is $`10^{10}h^2M_{}`$, assuming a CO-to-H2 conversion factor of $`1(\mathrm{M}_{}\mathrm{km}\mathrm{sec}^1\mathrm{pc}^2)^1`$ which is typical for starburst/kinematically violent, UV-intense environments of gas-rich, IR-ultraluminous systems (Downes & Solomon 1998). Assuming the CO is in a rotating disk, the dynamical mass can be calculated from the radius of $`500`$ pc set by the lens modeling, and using a rotational velocity of 350 km sec-1 set by the observed of line velocity HWHM = 250 km sec-1 and assuming a disk inclination angle of 45o (Downes et al. 1999). The implied dynamical mass is $`1.5\times 10^{10}`$ M within $``$500 pc of the nucleus, consistent with the value derived from the CO flux. Hence it appears likely that the molecular gas mass makes a significant, and perhaps dominant, contribution to the total mass within a few hundred parsecs of the nucleus in APM08279+5255, unless the nuclear CO disk is close to face-on. A similar conclusion has been reached for most nearby nuclear starburst galaxies (Downes and Solomon 1998).

## 4 Conclusions

This paper has presented resolved images of nuclear CO(1-0) emission in the gravitationally lensed BAL quasar APM08279+5255. While the continuum emission is found to be well aligned with the optical quasar images, the CO(1-0) is more extended, with a broken ring-like appearance. Such a structure is consistent with the action of gravitational lensing, with the continuum emission occurring on the scale of the quasar core, while the CO(1-0) arises from a larger region and is differentially magnified. The three-image nature of APM08279+5255 has posed a problem for lens modeling, as an extremely large, flat core is required to produce the central image. Such three image configurations are a nature consequence of gravitational lensing by a flattened potential which can produce naked cusps. Modeling of the CO(1-0) emission supports this hypothesis, although a deficit in constraints implies that the model is not unique. An immediate prediction of this model is that the lensing galaxy, whose position could be revealed by observing below the Lyman limit for this system $`(\stackrel{<}{}4400\mathrm{\AA })`$, hence removing the glare from the quasars, should be offset $`0.5`$arcsec from the quasar image, rather than lying behind the quasar images.

Currently, our CO images of APM08279+5255 are of limited signal-to-noise. However, with further integration a detailed map of the CO image can be made. As this region will be free from the effects of microlensing, and as its extended nature provides many more constraints (Kochanek, Keeton & McLoed 2001), such imaging has the potential to provide a more accurate model of the lensing in APM08279+5255 than from the quasar images.


text:
# The brown dwarf desert as a consequence of orbital migration

## 1 Introduction

Recent surveys have demonstrated a high abundance of brown dwarfs in young clusters (Martín et al. 2000), star forming regions (Béjar et al. 2001), and the field (Kirkpatrick et al. 1999, 2000). Brown dwarfs are also reasonably common in wide binaries (Gizis et al. 2001). The sole well established exception is in close (semi-major axis $`a\stackrel{<}{}4\mathrm{AU}`$) binaries. The same radial velocity surveys that have been so successful in finding massive extrasolar planets show that brown dwarfs are rarely close binary companions to Solar-type stars (Marcy & Butler 2000; Halbwachs et al. 2000). This brown dwarf desert supports the conventional belief that massive planets and brown dwarfs form in distinctly different ways, and also suggests that some aspect of the formation or early evolution of brown dwarfs differs from that of stars.

One possibility for explaining the desert is to postulate that the formation mechanism for brown dwarfs is entirely different from that of either stars or giant planets. Reipurth & Clarke (2001), for example, suggest that brown dwarfs are objects whose growth towards stellar masses was curtailed by ejection from small multiple systems (see also Reipurth, Clarke & Delgado-Donate 2001). In their model, the desert arises because substellar objects that were not ejected continued to accrete, eventually reaching stellar masses.

In this paper, we investigate a less radical (and more limited) possibility. We show that the absence of brown dwarfs in close binaries with Solar-type stars could be due to orbital migration within an evolving protoplanetary disc. This is the same process that is invoked to explain the presence of massive planets at small orbital radii (Lin, Bodenheimer & Richardson 1996), and which may also play a role in the orbital evolution of cataclysmic variables (Taam & Spruit 2001) and supermassive black hole binaries (Gould & Rix 2000). The critical assumption is that brown dwarfs form contemporaneously with the star, and thus become embedded within a young and relatively massive protoplanetary disc. Under these conditions, we show that migration efficiently clears out the desert by forcing the brown dwarfs into mergers with the star. This is in contrast to the evolution of massive planets, which have a chance of escaping the same fate by forming later, when the disc is close to being dispersed and cannot drive migration through to merger (Trilling, Lunine & Benz 2001; Armitage et al. 2001).

## 2 Brown dwarf migration

We assume that close brown dwarfs companions form at essentially the same time as the star, either via fragmentation during cloud collapse, or via disc instabilities during the earliest protostellar stages (e.g. Bonnell & Bate 1994). The brown dwarf companions open up a gap in the protoplanetary disc, and their subsequent fate then depends upon the evolution of the coupled disc-brown dwarf system. If the disc mass is large enough, relative to the brown dwarf mass, then the brown dwarf migrates in radius as the disc evolves as if it were a fluid element in the gas. This can be inward or outward, depending upon the sign of the radial velocity in the disc. If the disc mass is too low, on the other hand, the brown dwarf will act as a dam in the disc, holding up accretion and suffering minimal orbital migration (Syer & Clarke 1995). We show later that the first regime, where migration is rapid, is probably appropriate for brown dwarfs around Solar mass stars, while the no-migration outcome is likely if the star has very low mass.

### 2.1 Numerical methods

The interaction of a low mass secondary with particulate or gaseous discs has been extensively studied (Goldreich & Tremaine 1980; Lin & Papaloizou 1986). Based on this work, we model the interaction using a one dimensional treatment, in which we solve for the coupled evolution of the disk surface density $`\mathrm{\Sigma }(r,t)`$ and brown dwarf semi-major axis $`a`$. The governing equation, including mass loss in a disc wind at a rate $`\dot{\mathrm{\Sigma }}_w(r)`$, is (Lin & Papaloizou 1986; Trilling et al. 1998),

$$\frac{\mathrm{\Sigma }}{t}=\frac{1}{r}\frac{}{r}\left[3r^{1/2}\frac{}{r}\left(\nu \mathrm{\Sigma }r^{1/2}\right)\frac{2\mathrm{\Lambda }\mathrm{\Sigma }r^{3/2}}{(GM_{})^{1/2}}\right]\dot{\mathrm{\Sigma }}_w.$$
(1)
where $`M_{}`$ is the stellar mass. The first term on the right-hand side describes the diffusive evolution of a disk in which the angular momentum transport can be parameterized via a kinematic viscosity $`\nu `$. The second term describes how the disk responds to the torque from the secondary, which is approximated as a fixed function of radius $`\mathrm{\Lambda }(r,a)`$. We adopt,

$`\mathrm{\Lambda }`$ $`=`$ $`{\displaystyle \frac{q^2GM_{}}{2r}}\left({\displaystyle \frac{r}{\mathrm{\Delta }_p}}\right)^4r<a`$
$`\mathrm{\Lambda }`$ $`=`$ $`{\displaystyle \frac{q^2GM_{}}{2r}}\left({\displaystyle \frac{a}{\mathrm{\Delta }_p}}\right)^4r>a`$ (2)
where $`q=M_{\mathrm{BD}}/M_{}`$ is the mass ratio of the binary, and $`\mathrm{\Delta }_p`$ is given in terms of the disc scale height $`h`$ by,

$$\mathrm{\Delta }_p=\mathrm{max}(h,|ra|).$$
(3)
The rate of migration of the brown dwarf is,

$$\frac{\mathrm{d}a}{\mathrm{d}t}=\left(\frac{a}{GM_{}}\right)^{1/2}\left(\frac{4\pi }{M_{\mathrm{BD}}}\right)_{r_{\mathrm{in}}}^{r_{\mathrm{out}}}r\mathrm{\Lambda }\mathrm{\Sigma }dr.$$
(4)
Tests show that the detailed form of the torque function is unimportant for determining the migration rate once a gap has been opened. We have used the same form adopted by Armitage et al. (2001) to study massive planet migration.

The use of equation (1) implicitly assumes that (a) the brown dwarf orbit remains circular, and (b) that there is no ongoing accretion across the gap. Current theoretical understanding does not permit a definitive answer to either question. Existing simulations suggest that the assumption of circular orbits ought to be reasonable for massive planets, and perhaps for low mass brown dwarfs, but that for more massive brown dwarfs ($`M_{\mathrm{BD}}\stackrel{>}{}20M_{\mathrm{Jupiter}}`$) or stars the interaction with the disc may drive eccentricity growth (Artymowicz et al. 1991; Papaloizou, Nelson & Masset 2001). Ongoing accretion is likely to be substantially reduced at the relatively large mass ratios appropriate to brown dwarf secondaries (Bryden et al. 1999), but the exact level is rather uncertain.

### 2.2 Disc model

For the disc, we use a variant of the Clarke, Gendrin & Sotomayor (2001) model, which combines viscous evolution with mass loss at large radii. The model is motivated by observations of disc photoevaporation in Orion (Johnstone, Hollenbach & Bally 1998), though the same process may also occur at a lower level for relatively isolated stars (Shu, Johnstone & Hollenbach 1993). Including mass loss from the disc enables the model to reproduce the observed rapid transition between accreting Classical T Tauri stars and non-accreting Weak-Lined systems (Clarke, Gendrin & Sotomayor 2001). We take,

$$\nu =1.75\times 10^{13}\left(\frac{r}{1\mathrm{AU}}\right)^{3/2}\mathrm{cm}^2\mathrm{s}^1$$
(5)
consistent with a $`\mathrm{\Sigma }r^{3/2}`$ surface density profile. We start the runs with a steady-state surface density profile (accretion rate constant with radius), and assume $`h=0.05r`$ at all radii.

Appropriate values for the mass loss are poorly known for the situation, relevant here, where the star is not part of a rich cluster such as Orion. We assume that mass is lost from the disc outside a critical radius $`r_{\mathrm{crit}}=5\mathrm{AU}`$, with a radial scaling $`\dot{\mathrm{\Sigma }}_wr^1`$. For most of our runs, we normalise $`\dot{\mathrm{\Sigma }}_w`$ such that the total mass loss rate integrated to 25 AU is $`5\times 10^9M_{}\mathrm{yr}^1`$. To gauge the sensitivity of the results to the mass loss prescription, we also run models with the mass loss rate reduced to $`5\times 10^{10}M_{}\mathrm{yr}^1`$.

Figure 1 shows the evolution of the disc mass for initial masses of $`0.1M_{}`$ and $`10^2M_{}`$. We concentrate on the higher initial disc mass case, which is likely to be representative of the disc around a Solar mass star. A disc mass of the order of $`0.1M_{}`$ is comparable to the most massive discs inferred from mm wavelength observations (Osterloh & Beckwith 1995), and can also be justified theoretically – it corresponds to the marginally gravitationally unstable state expected as an endpoint of the disc formation process (e.g. Lin & Pringle 1990). For this model, a break in the evolution occurs at around 3 Myr, when the disc wind becomes dominant. The subsequent evolution is consistent both with the inferred time-dependence of the accretion rate for T Tauri stars (Hartmann et al. 1998), and with the evolution of the disc fraction with cluster age (Haisch, Lada & Lada 2001).

Operationally, we solve equation (1) using an explicit method on a non-uniform grid. The runs described below all used 300 mesh points, with an inner radius of 0.075 AU and an outer radius of 33.3 AU. At $`R_{\mathrm{in}}`$ a zero-torque boundary condition ($`\mathrm{\Sigma }=0`$) was used, while at $`R_{\mathrm{out}}`$ we set the radial velocity $`v_r=0`$.

In our current implementation, brown dwarfs reaching the inner boundary at $`R_{\mathrm{in}}`$ are declared as mergers. We do not include any ‘stopping mechanism’, sometimes invoked to slow or halt migration of extrasolar planets at very small radii. This is because the evidence for such mechanisms is weak. Although there are plausible reasons why migration might stall at small radii – for example due to an inner magnetospheric cavity in the disc (Königl 1991) – the observed radial distribution of extrasolar planets can be matched satisfactorily without including any stopping mechanism (Trilling et al. 2001; Armitage et al. 2001).

### 2.3 Migration-driven mergers

Figure 2 shows the migration history of brown dwarfs embedded in protoplanetary discs with initial masses of $`0.1M_{}`$ and $`10^2M_{}`$. The brown dwarf mass was $`0.04M_{}`$, i.e. substantially less than the disc for the more massive initial disc model. Disc evolution in the initially more massive disc drives the brown dwarf into merger with the star for all initial orbital radii $`a\stackrel{<}{}10\mathrm{AU}`$. Migration is especially rapid for brown dwarfs with initial orbits inside $`2\mathrm{AU}`$. These objects are swept inside $`0.1\mathrm{AU}`$ within less than a Myr after the start of the calculation. Reducing the rate of mass loss from the outer disc by an order of magnitude tends to increase the extent of migration, as it means that there is a larger reservoir of gas exterior to the brown dwarf orbit at late times. For orbits with initial semi-major axis $`a\stackrel{<}{}10\mathrm{AU}`$, however, the change is relatively small.

As expected, this rapid migration does not occur if the disc mass is much smaller than the brown dwarf mass. For a disc with an initial mass of only $`10^2M_{}`$, the angular momentum reservoir of the brown dwarf is large enough to prevent substantial orbital evolution. Only a modest migration occurs before the disc wind disperses the disc and freezes the brown dwarf into its final orbit. The low mass runs provide a qualitative idea of the evolution around very low mass stars, where the initial disc mass is likely to be smaller than the mass of any brown dwarfs present. If we assume that the initial disc mass is proportional to the stellar mass, then the results suggest that substantial migration would still occur for moderately less massive stars than the Sun, with masses $`0.5M_{}`$. For stars with masses of $`0.10.2M_{}`$, however, very little migration is expected.

Migration will efficiently empty the desert of brown dwarfs with initial semi-major axis $`a\stackrel{<}{}8\mathrm{AU}`$, provided that the initial disc mass is of the order of $`0.1M_{}`$. The only brown dwarfs that will remain at radii observable with radial velocity surveys are those that have migrated inwards from greater distances, before becoming stranded at small separations by the dispersal of the disc. Continuity demands that such objects must exist, but their numbers are generically likely to be small. In most models the disc gas is outflowing at large radii (Basu 1998; Hartmann et al. 1998) and tending to drive outward, rather than inward, migration. Moreover, as is evident from the divergence of the tracks in Fig. 2, migration leads a dilution in the number of brown dwarfs per unit logarithm in radius. Surveys have discovered some brown dwarfs at small radii – for example the outer companion (at $`3\mathrm{AU}`$) of the HD 168443 system (Marcy et al. 2001) – so the observations require only a substantial rather than a complete depletion of the population. Quantitatively, Gizis et al. (2001) estimate that if brown dwarfs followed the same distribution of separations as stars, at least ten times the observed number would be expected with separations less than 3 AU. There is considerable uncertainty in this estimate, but it nonetheless suggests that a theory for the desert needs to produce roughly an order of magnitude depletion in the population of short period brown dwarfs.

To estimate the range of initial orbital radii that will lead to surviving short-period brown dwarfs, we have run a series of models in which brown dwarfs were placed at larger initial radii in a $`0.1M_{}`$ disc. The results of these calculations are also shown in Fig. 2. For our disc model, binaries with initial orbital radii in the range $`10\mathrm{AU}\stackrel{<}{}a\stackrel{<}{}14\mathrm{AU}`$ are progenitors of observable systems in which the final orbital radius is $`a<4\mathrm{AU}`$. Smaller initial separations lead to merger, as discussed already, while larger separations fail to migrate inwards into the observable window prior to disc dispersal.

To use these results to estimate the net depletion of the desert, we need to know the initial distribution of brown dwarf orbital radii. If we assume that the formation of binaries that include a brown dwarf is scale free (i.e. equal numbers per logarithmic interval of separation), then migration leads to a net depletion by an order of magnitude in the $`0.1\mathrm{AU}<a<4\mathrm{AU}`$ range of orbital radii. Alternatively, we could assume that the initial radii of the binaries mirrored that observed for solar-type main sequence binaries generally. Duquennoy & Mayor (1991) parameterize the distribution with period $`P`$ (in days) as,

$$f(\mathrm{log}P)e^{\left[\frac{(\mathrm{log}P\overline{\mathrm{log}P})^2}{2\sigma _{\mathrm{log}P}^2}\right]},$$
(6)
where $`\overline{\mathrm{log}P}=4.8`$ and $`\sigma _{\mathrm{log}P}=2.3`$. This distribution is an increasing function over the range of initial radii of interest here, and thus we estimate a rather smaller depletion of brown dwarfs using this as the initial distribution. Between $`0.1\mathrm{AU}<a<4\mathrm{AU}`$, we obtain a depletion by a factor of 6.

### 2.4 Rotation rates of merger remnants

The best observational test of this model would be the detection of a significantly larger brown dwarf frequency among close pre-main-sequence binaries, in which accretion is still ongoing. There may also be evidence for mergers within the distribution of stellar rotation rates. Unlike gas, which can be accreted via a magnetosphere with essentially zero net change in stellar specific angular momentum (Königl 1991; Armitage & Clarke 1996), a merging brown dwarf arrives with the Keplerian angular momentum near the stellar surface. The ratio of the brown dwarf to stellar angular momentum is,

$$\frac{L_{\mathrm{BD}}}{L_{}}=\frac{1}{k^2}\left(\frac{M_{\mathrm{BD}}}{M_{}}\right)\left(\frac{\mathrm{\Omega }_K}{\mathrm{\Omega }_{}}\right)$$
(7)
where $`\mathrm{\Omega }_{}`$ is the stellar angular velocity, $`\mathrm{\Omega }_K`$ the Keplerian angular velocity at the stellar surface, and $`k^20.2`$ is the radius of gyration for the star (assumed fully convective). For an initial rotation period $`P7\mathrm{dy}`$, appropriate to many Classical T Tauri stars (Bouvier et al. 1993; but see also Stassun et al. 1999), $`L_{\mathrm{BD}}/L_{}1`$, so the merger leads to substantial spin-up. The results outlined above suggest that of the order of 10% of stars might have suffered such mergers with brown dwarfs. If the spin-down time-scale subsequent to merger is relatively long (of the order of a Myr), we would then expect a small fraction (a few percent) of Classical T Tauri stars to be rapid rotators as a consequence of mergers.

## 3 Summary

We have shown that inward orbital migration of brown dwarfs, within an evolving protoplanetary disc, can account for a low frequency of brown dwarfs as close binary companions to stars with masses $`M_{}M_{}`$. Migration depletes the initial frequency of brown dwarf companions at all radii where there is significant viscous evolution of the protoplanetary disc. This region could extend out to several tens of AU. Brown dwarfs with initially larger orbital radii $`a10^2\mathrm{AU}`$ would be unaffected. In our specific model, migration leads to the destruction of all brown dwarfs with initial orbital radii $`a\stackrel{<}{}10\mathrm{AU}`$, via mergers with the star which cause significant stellar spin-up. Orbits with $`0.1\mathrm{AU}<a<4\mathrm{AU}`$ are partially replenished by brown dwarfs migrating inwards from still greater radii, but the net brown dwarf frequency at these radii is still reduced by a factor of 5-10.

The main prediction of the model is that brown dwarfs in close orbits ought to be up to an order of magnitude more common amongst the youngest pre-main-sequence stars (less than a Myr), as compared to the main sequence. For migration to occur, we also require that the initial disc mass be at least comparable to the mass of the brown dwarf. No brown dwarf desert is thus expected around the lowest mass hydrogen-burning stars ($`0.10.2M_{}`$), or indeed around other brown dwarfs, whose discs would have been too feeble to drive significant brown dwarf migration.