text:
# Discovery and Implications of a New Large-Scale Stellar Bar in NGC 5248

## 1 Introduction

Grand-design spiral galaxies generally host two long symmetric spiral arms which can be followed over a large azimuthal angle and dominate the optical disk, e.g., NGC 1300 (Sandage 1961), NGC 895, NGC 1566 (Elmegreen & Elmegreen 1995), M51, M100 (Elmegreen, Seiden, & Elmegreen 1989) and NGC 5248 (e.g., Patsis, Grosbol, & Hiotelis 1997). Within the framework of the spiral density wave (SDW) theory (e.g., Lin & Shu 1964) spiral arms in disk galaxies result from a traveling wave pattern that remains quasi-stationary in a frame of reference rotating around the galaxy center at a certain pattern speed. Spiral density waves can be excited by a companion or by a bar. The spiral waves can get swing amplified as they are sheared from leading to trailing (e.g., Goldreich & Lynden-Bell 1965) provided the disk is cold enough.

In this Letter, we focus on the driving mechanism for a SDW in the nearby ($`D`$=15 Mpc) well studied grand-design spiral NGC 5248 and discuss the implications of an extended stellar bar in this galaxy. Optical imaging (e.g., Fig. 1) of NGC 5248 reveals two bright, relatively symmetric spiral arms from about $`30\mathrm{}`$ (2.2 kpc) to $`80\mathrm{}`$ (6.0 kpc). NGC 5248 does not show evidence of any recent or ongoing interaction with its three faint irregular neighbors (UGC 8575, UGC 8614, and UGC 8629; Zaritsky et al. 1997) which have similar redshifts and are located more than $`30\mathrm{}`$ (135 kpc) away in projection. It is therefore likely that the grand-design spirals in NGC 5248 are driven by a massive large-scale bar, which is verified in the present work. The circumnuclear ring of Hii regions at a radius of $`5\mathrm{}`$ (375 pc) (e.g., Elmegreen et al. 1997) also suggests the presence of a bar since rings are commonly associated with the dynamical resonances of a bar (e.g., review by Buta & Combes 1996).

NGC 5248 is classified as SAB(rs)bc (de Vaucouleurs et al. 1991; hereafter RC3) and has been postulated to have a weak stellar bar with an ellipticity of 0.32 and a semi-major axis of $`22\mathrm{}`$ (1.6 kpc) at a position angle (PA) of $`110^{}`$ (e.g., Martin 1995). This oval feature is evident in optical images (e.g., Fig. 1). However, the extent and morphology of the grand-design spirals in NGC 5248 strongly suggests that the bar driving these spirals must be significantly longer and more massive than previously assumed. The two bright inner stellar spiral arms have prominent dust lanes on their inner (concave) sides (Fig. 1) out to a radius $`r70\mathrm{}`$ (5 kpc). Such a dust lane morphology is expected to exist only inside the corotation resonance (CR) of the arms, where young stars which form when gas is compressed by shocks, seen as dust lanes, can overtake the pattern From this argument alone, we expect that the CR of the spiral pattern in NGC 5248 is beyond $`70\mathrm{}`$. If the bar and spirals in this galaxy have the same pattern speed, the latter must be driven by a stellar bar whose semi-major axis is comparable in size to that of their common CR. Large-scale bars are expected to end at or near CR based on the studies of stellar orbits (Contopoulos & Papayannopoulos 1980), and on the analysis of shapes of offset dust lanes within bars (Athanassoula 1992). In order to test the existence of an extended bar in NGC 5248, we obtained a deep $`R`$-band image with a large field of view ($`\mathrm{\S }`$ 2). While this Letter focuses on the observational aspects of the large-scale bar in NGC 5248 and its consequences for galactic dynamics in nearby and high redshift galaxies, the subsequent paper (Jogee et al. 2002, hereafter Paper II) addresses gas dynamics in the inner few kpc of this galaxy, based on multi-wavelength observations and modeling.

## 2 Observations and Data Reduction

NGC 5248 was imaged through the Harris $`R`$-band filter for a total of 30 minutes using the Wide Field Camera on the 2.5 m Isaac Newton Telescope (INT) in La Palma in August 2001. The data frames had a plate scale of 0$`\stackrel{}{\mathrm{.}}`$33/pixel, a field of view (fov) of 11$`\stackrel{}{\mathrm{.}}`$3 $`\times `$ 22$`\stackrel{}{\mathrm{.}}`$5, and an average seeing of 1$`\stackrel{}{\mathrm{.}}`$8. Photometric calibration was performed by comparison with a shorter image of the galaxy in the same band, taken on the same night with the 1 m Jacobus Kapteyn Telescope. The $`R`$-band surface brightness (SB) goes down to 26.5 mag arcsec-2. Bias-subtracted and flat-fielded frames were obtained from the INT data reduction pipeline. Sky-fringing at the 1 % level is present, an effect known for the thinned EEV chip CCD4. We used the IRAF package for fixing bad pixels, cleaning cosmic rays, and combining frames into a final image. After masking out stars and faint background galaxies, isophotal analysis was performed with the “isophote” package in IRAF. We also refer to a $`K^{}`$-image obtained in May 2000 using the INGRID camera on the 4.2 m William Herschel Telescope. This image, and other data described in this paper, are fully detailed in Paper II.

## 3 The extended bar and faint outer disk

The $`R`$-band image (Fig. 1) reveals a hitherto unknown large-scale stellar bar-like feature embedded within a fainter, more circular, outer disk which is visible out to a radius of $`230\mathrm{}`$ (17.2 kpc). Below we describe the morphological, photometric, and dynamical properties of this bar. Earlier optical and NIR images lacked the sensitivity and/or field of view to clearly detect the end of the large-scale stellar bar and the faint outer disk. As a result, the bar was misinterpreted as an inclined outer disk. The new $`R`$-band image (Fig. 1) shows the actual outer disk and the full extent of the bar. In order to quantitatively describe the features in NGC 5248, we performed an isophotal analysis (Fig. 2) of the $`R`$-band and available $`K`$-band images. To minimize the effects of extinction, the $`K`$-band image of $`0\stackrel{}{\mathrm{.}}6`$ seeing is used out to a radius of $`50\mathrm{}`$. Beyond this point, its S/N ratio is low and we use the $`R`$-band image.

The bright optical stellar spiral arms extend along the leading ($`\mathrm{\S }`$ 4) edge of the bar and cross its major axis (PA $`135^{}`$) at around $`95\mathrm{}`$ (Fig. 1). This suggests that the bar has a semi-major axis of about $`95\mathrm{}`$ (7.1 kpc). Large-scale bars typically tend to end near the ultra harmonic resonance where chaos and orbital instability set in. The ratio of corotation to bar radii ($`R_{\mathrm{CR}}`$/$`R_{\mathrm{bar}}`$) is expected to lie in the range of $`1.2\pm 0.2`$ as found empirically by Athanassoula (1992). This corresponds to a CR of $`115\mathrm{}`$ (8.6 kpc) in NGC 5248 – a value supported by photometric evidence. As seen in Fig. 2, between $`60\mathrm{}`$ to $`115\mathrm{}`$, the projected ellipticity e rises to a maximum of 0.52 and the PA settles to a fairly constant value of $`135^{}`$, after some initial variation between $`120^{}`$$`135^{}`$. This small variation is due to the unusual prominence of the spiral arms within the bar. The peak in $`e`$ over a fairly constant PA, as observed in NGC 5248, is the behavior expected of bars (e.g., Knapen et al. 1995).

Between $`a`$ = $`115\mathrm{}`$$`150\mathrm{}`$, secondary faint spirals which stem off the bar at lower radii become particularly prominent. They dominate the fitted isophotes which show a varying e of 0.41–0.46 (Fig. 2). Such features are reminiscent of “plumes” which are seen in other barred systems such as NGC 1433 and may be secondary compression zones which arise near CR (Buta 1984). Between $`150\mathrm{}`$$`230\mathrm{}`$, the isophotes are dominated by the outer disk and very faint spirals which extend away from the bar and cover $`90^{}`$$`100^{}`$ in azimuth (Fig. 1). This leads to a low but varying $`e`$ (Fig. 2). The outer Lindlad resonance (OLR), where arms are expected to end, is probably around $`230\mathrm{}`$. The PA of the projected disk is relatively constant at $`105^{}`$. Assuming the outer disk is intrinsically circular, we derive a PA of $`105^{}\pm 2^{}`$ for the line of nodes and an inclination of $`40^{}\pm 4^{}`$. With these parameters, we performed a two-dimensional deprojection of the $`R`$-band and $`K`$-band images. The deprojected disk (Paper II) is circular with e below 0.02 (Fig. 2) while the deprojected bar is moderately strong with a peak e of 0.44. We note that the deprojection is somewhat uncertain because the intrinsic shape of the disk between CR and OLR is expected to be slightly elongated along the bar minor axis.

We now turn our attention to the inner regions ($`r50\mathrm{}`$) of NGC 5248. A very weak oval feature is present in the inner $`3\mathrm{}`$ radius of the $`K`$-band image (Fig. 3) where the deprojected e varies from 0.18 to 0.10 (Fig. 2). Two $`K`$-band spirals, extending from this oval, are lined with “hot spots” between $`3\mathrm{}`$$`9\mathrm{}`$. They cross the optically-visible starburst ring at $`5\mathrm{}`$ and between $`5\mathrm{}`$$`8\mathrm{}`$, they delineate the bright super star clusters seen in HST images (see Paper II). Near the bar major axis, the spirals fade out, and e falls from 0.20 to 0.11 between $`10\mathrm{}`$$`13\mathrm{}`$ (Fig. 2). Between $`17\mathrm{}`$$`26\mathrm{}`$, the stellar distribution looks smoother (Fig. 3) and the $`K`$-band spirals fade as they approach the bar major axis. Further out, between $`30\mathrm{}`$$`50\mathrm{}`$, the bright patchy outer $`K`$-band spirals (Fig. 3) dominate, leading to a higher $`e`$ (0.15–0.30) and a smoothly varying PA. The oval-shaped feature of radius $``$ $`22\mathrm{}`$ which is visible in the projected images (e.g., Fig. 1) has been interpreted as a primary bar of semi-major axis $`22\mathrm{}`$ (1.6 kpc) at a PA of $`110^{}`$ in earlier studies (e.g., Martin 1995). However, after deprojection, this feature has a peak $`e`$ of only 0.15 (Fig. 2). Therefore, the oval appearance of this feature in the sky plane seems to be largely due to projection effects. The dip in $`e`$ at around $`26\mathrm{}`$ is due to the K-band arms joining the oval at its minor axis. Uncertainties in the deprojection are unlikely to change this conclusion, namely, there appears to be no signs of a bar of size $`22\mathrm{}`$ in NGC 5248.

## 4 Discussion: Implications for barred spirals locally and at high redshifts

NGC 5248 shows a grand-design spiral morphology in stars, dust, and gas, as traced by $`B`$-band (Paper II), $`R`$-band (Fig. 1), $`K`$-band (Fig. 3), $`BI`$, H$`\alpha `$, and CO maps (Paper II). In $`B`$, $`R`$ and $`K`$ light, there are prominent stellar spiral arms between about $`30\mathrm{}`$$`80\mathrm{}`$ (2.2–6.0 kpc), on the leading edge of the large-scale stellar bar. The stellar arms have prominent dust spirals on their inner (concave) sides out to at least $`70\mathrm{}`$ (e.g., Fig. 1), as is expected inside the CR. The observed offset between the dust and stellar spirals is consistent with our understanding of gas flows in barred galaxies where offset dust lanes along the leading edges of the bar delineate the shocks which compress the gas, leading to the formation of massive young stars. Within this framework, the bar is rotating clockwise, and the patchy $`B`$, $`R`$, and $`K`$ stellar spiral arms immediately offset from the dust lanes are made of young stars. It is particularly noteworthy that the appearance of the $`K`$-band spirals in NGC 5248 is dominated by young stars. $`K`$-band light is generally assumed to trace old stellar populations, but it can be dominated by very young (8-10 Myr) supergiants in starburst regions (e.g., Knapen et al. 1995) Along the $`K`$-band arms in NGC 5248, there is good spatial coincidence between peaks in $`K`$-band emission and Hii regions as traced by H$`\alpha `$ images. The patchy morphology of the $`K`$-band spirals is in stark contrast to the smooth appearance of old stellar arms (e.g., Kennicutt & Edgar 1986). Taken together, these properties suggest that young stars dominate the $`K`$-band spirals.

While grand-design optical spiral arms in many strongly barred galaxies stem from the ends of a stellar bar, in NGC 5248, open spiral arms appear conspicuous within the moderately strong bar itself, delineating intense SF. The SF along the offset bar dust lanes in NGC 5248 is consistent with current understanding of gas flows and SF in barred galaxies. In strong bars which exhibit almost straight offset dust lanes delineating shocks, little SF is seen along the bar. Common examples are NGC 1300 and NGC 5383 (Tubbs 1982). Star formation is believed to be inhibited because of the strong shear in the postshock flow (e.g., Athanassoula 1992). In weaker bars, where the offset dust lanes are curved, the weaker shocks and shear can induce SF rather than inhibit it. For instance, the collapse of gas cores to form stars can be induced by weak shocks with speeds of order 20 to 30 km s-1 (e.g., Vanhala & Cameron 1998). It is thus common to see SF along weak bars, e.g., in M100 (Elmegreen et al. 1989), NGC 4254, and NGC 4303 (Koopmann 1997). NGC 5248 falls in this category as its large-scale bar is only moderately strong with curved offset dust lanes (Fig. 1), and a deprojected ellipticity of 0.44.

We note that there exists a class of weakly barred galaxies where bright optical spiral arms appear to emerge from a relatively short oval bar-like feature and extend to twice the radius of this bar (Elmegreen & Elmegreen 1995). We would have put NGC 5248 in this category if we had taken its primary stellar bar to be the oval feature of radius $`22\mathrm{}`$, as was long assumed to be the case. However, with our identification of the large-scale bar of radius $`95\mathrm{}`$, the bright optical spiral arms with concave dust lanes in NGC 5248 turn out to lie inside the bar. The short oval turns out to be primarily the result of projection effects ($`\mathrm{\S }`$ 3). It is therefore imperative to ask whether galaxies with bright extended optical spiral arms apparently driven by a short oval feature, are in fact hiding a more extended large-scale bar. NGC 1566 and ESO 111-110 (Buta 1995) are two such possible candidates.

While the majority of nearby spiral galaxies are barred (e.g., Grosbol 2002), a low fraction of bars ($`<`$ 10 %) has been reported in galaxies at redshifts $`z`$ 0.5–0.8 (e.g., Abraham et al. 1999). In this latter study, the authors apply the same bar identification method to nearby galaxies and to spirals at $`z`$ 0.5–0.8. Images which trace the rest-frame $`B`$-band light are used. The inclination $`i`$ is derived from an outer isophote which corresponds to an arbitrary fraction (1%) of the peak flux. After removing the projected ellipticity corresponding to $`i`$, a bar is identified based on the residual ellipticity in the inner bright regions of the spiral galaxy. The authors find a much lower bar fraction at intermediate redshifts than locally. It has been suggested that this result may be biased due to the small sample size or that bars are truly deficient at intermediate redshifts because of dynamically hotter disks or an enhanced efficiency in bar destruction (Abraham et al. 1999). Another possibility is that the method used is less efficient in identifying bars at intermediate redshifts than locally. This method relies on the $`B`$-band surface brightness ratio ($`R_1`$) between the central regions and the outer disk beyond the bar end. In a barred spiral with a bright central region and a faint outer disk beyond the bar end, $`R_1`$ is low, and this method can readily fail: the outer 1% isophote would lie inside the bar rather than in the faint, more circular outer disk and, consequently, the bar would be misinterpreted as an inclined disk. Bars may be missed more frequently at $`z`$ 0.5–0.8 than locally with this method because of two factors: evolution and imaging depth. If there is evolution in the surface-brightness profiles of spirals with redshift such that $`R_1`$ is on average lower at higher redshifts, then bars may be missed more frequently at $`z`$ 0.5–0.8 than locally. Detailed future studies with large samples are needed to quantify $`R_1`$, but our suggestion is not ruled out by existing data. Spiral galaxies at intermediate redshifts have bright blue central colors indicative of recent starbursts (Ellis, Abraham, & Dickinson 2001), and the outer disk/bar which have presumably formed recently may still be quite faint. Another factor is that images of spirals at $`z`$ 0.5–0.8 often do not reach the same sensitivity level as images of nearby galaxies.

We also discuss why a higher bar fraction at $`z`$ 0.5–0.8 is quite likely. It is well known that bars or $`m`$=2 modes can be easily excited in a dynamically cold disk either spontaneously or during tidal interactions and minor mergers. On the other hand, major mergers tend to produce highly disturbed systems which eventually evolve into elliptical-like systems (e.g., Barnes 1992) through violent relaxation. The fraction of spiral galaxies in pairs increases with redshift out to $`z`$ 0.75–1 (Le Fèvre et al. 2000), suggesting enhanced interaction between spiral galaxies. These interactions do not appear to involve major mergers since HST studies hint that sizes of spiral galaxies (e.g., Lilly et al. 1998) and the number counts of spirals with $`I<`$22 mag in HST fields (e.g., Glazebrook et al. 1995) do not show much evolution out to $`z`$ 1. Taken together, these results suggest that conditions out to $`z`$ 1 are favorable to the formation of bars in disk galaxies. Furthermore, if the recently discovered massive disk galaxy at $`z`$ =1.34 (Van Dokkum & Stanford 2001) is representative, it reveals that massive, organized, disks with ongoing SF are already in place by $`z`$ =1. The ring of SF at a $``$ 10 kpc radius may be associated with the dynamical resonances of a bar or alternatively could be an out-of-the-plane ring.

Support for this work was provided by a grant from the K. T. and E. L. Norris Foundation, and an NSF grant AST 99-81546 for S. J. and N. Z. S., and grants NAG 5-10823, HST GO-08123.01-97A, and WKU 516140-02-07 for I. S. We thank Thomas Augusteijn for INT service observations and Ron Buta for useful comments. Data were retrieved from the Isaac Newton Group (ING) archive. The William Herschel, Isaac Newton, and Jacobus Kapteyn Telescopes are operated on La Palma by the ING in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias.


format: text

text:
# Charge inversion in colloidal systems

## 1 Introduction

Charged colloids (macroions) might exhibit very surprising behaviors due to the Coulomb interaction and the two different length scales involved there: (i) microscopic (counterions) and (ii) mesoscopic (macroions). A very spectacular effect which has attracted great attention these last years is the phenomenon of overcharging . This situation occurs when the number of counterions in the vicinity of the macroion surface is so high that the macro-particle bare charge is overcompensated (charge inversion). This phenomenon has been observed experimentally by electrophoresis . On the theoretical side, a massive effort has been recently devoted to elucidate this striking phenomenon. Wigner crystal theories have been successfully applied to understand charge inversion  and quantitative agreement was found with simulations .

It turns out that all numerical and analytical methods, undertaken so far, neglect the discrete nature of the macroion charge distribution. But in realistic systems the macroion charge is not continuously distributed over the macro-particle, but rather carried by small microscopic ions. It is only very recently that a simulation study has been reported where the macroion charge discreteness has been considered .

In this paper, we concentrate on the effect of the macroion charge distribution discreteness. In particular, we study this effect on the overcharging properties. Besides we show that it is possible to have stronger overcharge in the discrete case than in the continuous case.

## 2 Computational details

### 2.1 Macroion charge discretization

The procedure is similar to the one used in a previous study . The macroion charge discretization is produced by using $`Z_m`$ monovalent microions of diameter $`\sigma `$ distributed randomly on the surface of the macroion (see Fig. 1). Then the macroion structural charge is $`Q=Z_me`$  where $`Z_m>0`$ and $`e`$ is the positive elementary charge. The counterions have a charge $`q=+Z_ce`$ where $`Z_c>0`$ stands for the counterion valence. The discrete colloidal charges (DCC) are fixed on the surface of the spherical macroion and are at a distance $`r_0`$ from the macroion center. Thus to produce a random discrete charge distribution on the surface we generated (uniformly) randomly the variables $`\mathrm{cos}\theta `$ and $`\phi `$ where $`\theta `$ and $`\phi `$ are the standard spherical angle coordinates. Excluded volume is taken into account by rejecting configurations leading to an overlap of microions. Phenomena such as surface chemical reactions , hydration, roughness  are not considered.

### 2.2 Molecular dynamics procedure

A MD simulation technique was used to compute the motion of the counterions coupled to a heat bath acting through a weak stochastic force W(t). The motion of counterion i obeys

$$m\frac{d^2𝐫_i}{dt^2}=_iU(𝐫_i)m\gamma \frac{d𝐫_i}{dt}+𝐖_i(t),$$
(1)
where m is the counterion mass, $`U`$ is the potential force having two contributions: (i) the Coulomb interaction and (ii) the excluded volume interaction, and $`\gamma `$ is the friction coefficient. Friction and stochastic force are linked by the dissipation-fluctuation theorem $`<𝐖_i(t)𝐖_j(t^{})>=6m\gamma k_BT\delta _{ij}\delta (tt^{^{}})`$. For the ground state simulations the stochastic force is set to zero.

Excluded volume interactions are taken into account with a pure repulsive Lennard-Jones potential given by

$$U_{LJ}(r)=4ϵ\left[\left(\frac{\sigma }{rr_0}\right)^{12}\left(\frac{\sigma }{rr_0}\right)^6\right]+ϵ,$$
(2)
if $`rr_0<r_{cut}`$, and $`U_{LJ}(r)=0`$ otherwise. For the microion-microion interaction (the microion being either a counterion and/or a DCC), whereas $`r_0=7\sigma `$ for the macroion-counterion interaction, and $`r_{cut}=2^{1/6}\sigma `$ is the cutoff radius. This leads to a macroion-counterion distance of closest approach $`a=8\sigma `$.

Energy and length units in our simulations are related to experimental units by taking $`ϵ=`$$`k_BT_0`$ (with $`T_0=298`$ K), and $`\sigma =3.57`$ Å respectively.

The pair electrostatic interaction between any pair ij, where i and j denote either a DCC, a counterion or the central charge (for the case of the continuous surface charge distribution), reads

$$U_{coul}(r)=k_BT_0l_B\frac{Z_iZ_j}{r},$$
(3)
where $`l_B=e^2/4\pi ϵ_0ϵ_rk_BT_0`$ is the Bjerrum length (here $`l_B=10\sigma `$) describing the electrostatic strength.

The macroion and the counterions are confined in a spherical impenetrable cell. The macroion is held fixed and is located at the center of the cell. To avoid image charge complications, the permittivity $`ϵ_r`$ is supposed to be identical within the whole cell (including the macroion) as well as outside the cell.

## 3 Results

The simplest way to quantify overcharging in salt-free environment is to compute the total electrostatic energy of the system in the ground state (i. e., $`T=0`$K) as a function of the number of overcharging counterions $`n`$. Although this approach is very idealized, it allows an insight into the mechanisms governing overcharging.

In what follows we are going to compare systematically (for a macroion charge $`Z_m=60`$) the overcharging occurring in the continuous case (CC), as is conventionally done, and in the discrete case (with DCC ions). Fig. 1 shows a typical starting equilibrium configuration with a discrete macroion charge distribution where the system is globally neutral (before adding excess overcharging counterions).

Fig. 2(a) shows overcharging energy profiles for $`Z_m=60`$ and $`monovalent`$ counterions ($`Z_c=1`$). It is observed that overcharging is stronger in the continuous case than in the discrete case as one could expect. More precisely the gain in energy as well as the maximal number of stabilizing overcharging counterions (corresponding to a minimum in the $`En`$ curve) are higher when the surface macroion charge is continuous. Indeed, having in mind that overcharging is enhanced by counterion ordering, it is clear that in the discrete case where strong ionic pairing occurs (between DCC site and counterion \- see Fig. 1) the counterion ordering is lowered and thus leads to a weaker charge inversion. However since overcharging occurs in the discrete case, the excess counterions still experience an attractive interaction with their neighboring dipoles (ionic pairs). This ion-dipole interaction is the driving force for the overcharging in the discrete case with monovalent counterions. Nevertheless it has been shown elsewhere  that by decreasing the mean separation between DCC sites (i. e., increasing macroion charge density) the difference in overcharging between discrete and continuous macroion charge distribution decreases.

When we consider $`multivalent`$ counterions the situation is qualitatively different. Indeed, one observes that for $`Z_c=`$ 2 and 3 overcharging is even stronger in the discrete case \[see Figs. 2 (b-c)\]. This is a rather counter-intuitive phenomenon, since ionic pairing is still effective. The crucial difference between monovalent and multivalent counterion systems is that overcharging occurs with ionic pairing for multivalent counterions (since free monovalent DCC sites are unbound), whereas for monovalent counterions only ion-dipole interaction can take place (since in the neutral state $`each`$ DCC site is already bound with a monovalent counterion). Obviously for these systems under consideration depicted in Figs. 2 (b-c), the ionic pairing energy is a fundamental quantity that controls the overcharging strength. More precisely, the stronger the ionic bonding the more efficiently it overcomes the $`excess`$ overcharge ”self-energy”, and therefore the higher the overcharging for multivalent counterion systems.

Beside of the ionic bonding strength there is another important parameter controlling overcharging: the surface macroion charge density. When the macroion charge density is high enough (not reported here) then the $`E(n)`$ curves become quasi-identical for discrete and continuous charge distributions. In terms of length scales, the relevant parameter controlling charge reversal with discrete distribution is the ration $`\rho `$ between the mean distance between DCC and the distance separation of the ionic pair. For low $`\rho `$ the effects of discretization are very important whereas for $`\rho 1`$ one recovers the continuous limit.

## 4 Conclusions

We carried out MD simulations to elucidate the effect of discrete macroion charge distribution on the overcharging. We showed that overcharging is still effective in the discrete case. For multivalent counterions and a sufficiently low macroion charge density, overcharging can even be stronger in the discrete case than in the continuous case. Future works will address aqueous solutions at finite temperature as well as the presence of salt-ions.


text:
# Structure of Vortices in Two-component Bose-Einstein Condensates

## Abstract

We develop a three-dimensional analysis of the phase separation of two-species Bose-Einstein condensates in the presence of vorticity within the Thomas-Fermi approximation. We find different segregation features according to whether the more repulsive component is in a vortex or in a vortex-free state. An application of this study is aimed at describing systems formed by two almost immiscible species of rubidium-87 that are commonly used in Bose-Einstein condensation experiments. In particular, in this work we calculate the density profiles of condensates for the same conditions as the states prepared in the experiments performed at JILA \[ Matthews, et al., Phys. Rev. Lett. 83, 2498 (1999)\].

Vortices in Bose-Einstein condensates (BECs) have thoroughly been studied from a theoretical point of view since the first experiments on condensates were carried out, due to their direct connection to superfluidity. For a review of these issues see for example , and for the most recent works references in . But until the experiment of Matthews et al.  was performed, no evidence about their existence had ever been reported. This experiment was based on a method proposed by Williams and Holland  to create vortices in two-component BECs. They trapped atoms of 87Rb in two different hyperfine states $`|F=1,m_f=1`$ and $`|F=2,m_f=1`$, henceforth denoted as $`|1`$ and $`|2`$, respectively, and created vortices either in the $`|1`$ or $`|2`$ component, while the other species remained without rotation. Because of the relation between intra- and inter-particle scattering lengths, these states have the property of being almost immiscible. We shall see that this fact has important consequences on the structure of the vortices.

In connection with this experiment interesting theoretical work has recently been developed  in order to describe different properties of many types of two-component condensate systems. García-Ripoll and Pérez García, by means of a robust formalism, analyzed the stability and dynamics of a variety of two-species condensates, taking into account the effect of the interchange of species in the system, and also the variation of the number of particles. In particular, in Ref.  they set the relative number of particles equal in both species and varied the total number of particles, and in Ref.  they pursued the research by also varying the relative composition using a 2D model. However, in these works  little is said about the structure of the stationary states and no comments are made about segregation of species. In a recent work, Chui et al.  computed the energy of different configurations within the Thomas-Fermi Approximation (TFA). They varied the number of particles and interchanged the species in order to derive conclusions about the relative stability. But, although they used a 3D model and considered phase separation in the mixture, a hypothesis was made regarding the form of the interface between species that we find is not consistent with our results and apparently also with new available experimental data . We shall comment on this later.

In 1996, in a pioneering work, Ho and Shenoy  constructed an algorithm to determine the density profiles of binary mixtures of alkali atoms within the TFA. We develop here a different procedure which allows one to obtain conclusions about the general features of segregation through simple calculations.

The aim of the present article is to analyze and perform a classification of the structure of condensates near phase separation, when two almost immiscible species are involved. Our starting point for calculating the stationary states are the Gross-Pitaevskii equations in the TFA :

$`\mu _1\mathrm{\Psi }_1`$ $`=`$ $`(V_1+N_1G_{1,1}|\mathrm{\Psi }_1|^2+N_2G_{1,2}|\mathrm{\Psi }_2|^2)\mathrm{\Psi }_1,`$ (1)
$`\mu _2\mathrm{\Psi }_2`$ $`=`$ $`(V_2+N_1G_{1,2}|\mathrm{\Psi }_1|^2+N_2G_{2,2}|\mathrm{\Psi }_2|^2)\mathrm{\Psi }_2,`$ (2)
where $`N_i`$ denotes the number of atoms of species $`i`$, and $`G_{k,l}=u_{k,l}U`$, with $`u_{k,l}`$ the relative interaction strengths between species $`k`$ and $`l`$. We set the most repulsive component $`|1`$ in the $`\mathrm{\Psi }_1`$ state, and fix $`u_{1,1}=1`$ $`(>u_{2,2})`$, so $`U=\frac{4\pi \mathrm{}^2a}{M}`$, $`a`$ being the scattering length of the $`|1`$-species and $`M`$ the atom mass.

For simplicity, we consider the system in a spherically symmetric trap and make a change of variables according to $`\sqrt{\frac{M}{2}}w_0\stackrel{}{r}\stackrel{}{r}`$, with $`w_0`$ the trap angular frequency. Thus the effective potentials $`V_i`$ written in cylindrical variables read

$$V_i=r^2+z^2+\frac{\kappa _i}{r^2}.$$
(3)
The first two terms correspond to the trapping potential, while the third includes the centrifugal term , $`\kappa _i=\frac{m_i^2\mathrm{}^2w_0^2}{4}`$, with $`m_i`$ the number of quanta of circulation associated with the wave function for species $`i`$.

Let us first analyze what happens when $`Du_{1,1}u_{2,2}u_{1,2}^2=0`$. Physically the system starts undergoing phase separation. Mathematically, for both $`|\mathrm{\Psi }_i|^2>0`$, Eqs. (2) become linearly dependent. So to have a compatible set of equations the following condition must hold:

$$\frac{\mu _1}{G_{1,1}}\frac{V_1}{G_{1,1}}=\frac{\mu _2}{G_{1,2}}\frac{V_2}{G_{1,2}}.$$
(4)
If we define $`\beta _2=u_{1,2}/u_{1,1}`$, the above equation can be rewritten as

$$\frac{\mu _2\beta _2\mu _1}{1\beta _2}r^2z^2+\frac{\beta _2\kappa _1\kappa _2}{(1\beta _2)r^2}=0.$$
(5)
This equation defines different families of surfaces that we shall call $`𝒮`$$`{}_{}{}^{i}{}_{0}{}^{}`$, where $`i`$ labels the family. As was mentioned above, this is the only region of space where both $`\mathrm{\Psi }_i`$ may be non-zero, so these surfaces necessarily determine the interfaces between species. Such families are characterized by the values of

$$A=\frac{\mu _2\beta _2\mu _1}{1\beta _2},B=\frac{\beta _2\kappa _1\kappa _2}{(1\beta _2)},$$
(6)
(note that $`D=0`$ implies $`\beta _2<1`$).

Depending on the values of the above parameters, three topologically different kinds of surfaces exist, namely, spheres, cylinder-like surfaces, and toroids. In Fig. 1 we show a partition of the parameter space $`(A,B)`$ indicating the domains corresponding to each family. In order to classify these surfaces let us study the $`z(r)`$ curves defined by Eq. (5). If $`z(r)`$ intersects the $`z`$ axis, the associated surfaces $`𝒮`$$`{}_{}{}^{a}{}_{0}{}^{}`$ are spheres, and this occurs for $`B=0`$ and $`A>0`$. If $`z(r)`$ has only one root the corresponding surfaces $`𝒮`$$`{}_{}{}^{b}{}_{0}{}^{}`$ are cylinder-like, and this takes place when $`B>0`$. Finally, if $`z(r)`$ has two roots the related surfaces $`𝒮`$$`{}_{}{}^{c}{}_{0}{}^{}`$ are toroids, and it is easy to verify that the conditions to be fulfilled are $`A^2+4B>0`$, $`A>0`$, and $`B<0`$. The remaining region of the parameter space has no associated surfaces.

To perform a displacement of a point within the parameter space two mechanisms are possible. On the one hand, the $`A`$ value may be changed by varying the number of particles in each species, which will cause a change in the chemical potentials. A variation of this type does not change the family, but the shape of the surface is modified. For example, for $`B>0`$ one can get a straighter “cylinder” by moving from positive to negative $`A`$ values, and, in terms of the number of particles, this corresponds to a decrease in the $`N_2/N_1`$ ratio. On the other hand, $`B`$ can only take discrete values and it is modified by changing the vorticity.

It is worthwhile mentioning that the only possible realization for $`B=0`$, (and hence for $`𝒮`$$`{}_{}{}^{a}{}_{0}{}^{}`$), within reasonable values of $`m_i`$, is a vorticity free condensate $`m_1=m_2=0`$. In such a case, regarding the location of each species, it is easy to conclude that the energy minimum is obtained when the less repulsive component lies inside the sphere of radius $`\sqrt{A}`$, while the other species forms a spherical shell around it. In the presence of vorticity ($`B0`$), the location of each species can be predicted from energetic considerations. The state with greater vorticity should be located in the region of space that excludes $`r=0`$, because of the centrifugal term.

From this simple analysis we may state that close to phase separation the presence of vorticity dramatically changes the segregation structure.

Now, let us consider the general case $`D0`$. The solution of Eqs. (2 ) can be easily obtained and has the following different expressions depending on whether there exists any overlap between the wave functions of both species. a) In the region where only one wave function is non-vanishing ( $`|\mathrm{\Psi }_i|^20`$ and $`|\mathrm{\Psi }_k|^2=0`$, for $`ik`$) Eqs. (2) are decoupled and the solution is

$$|\mathrm{\Psi }_i|^2=\left[\mu _ir^2z^2\kappa _i/r^2\right]/(G_{i,i}N_i).$$
(7)
b) In the region where both wave functions are non-vanishing, $`|\mathrm{\Psi }_i|^2>0`$, one obtains

$`|\mathrm{\Psi }_1|^2`$ $`=`$ $`\left[{\displaystyle \frac{\mu _1\beta _1\mu _2}{(1\beta _1)}}r^2z^2+{\displaystyle \frac{\kappa _2\beta _1\kappa _1}{r^2(1\beta _1)}}\right]A_1u_{2,2},`$ (8)
$`|\mathrm{\Psi }_2|^2`$ $`=`$ $`\left[{\displaystyle \frac{\mu _2\beta _2\mu _1}{(1\beta _2)}}r^2z^2+{\displaystyle \frac{\kappa _1\beta _2\kappa _2}{r^2(1\beta _2)}}\right]A_2u_{1,1},`$ (9)
where $`\beta _1=\frac{u_{1,2}}{u_{2,2}}`$ and $`A_i=\frac{(1\beta _i)}{UN_iD}`$.

We call $`𝒮`$1 ( $`𝒮`$2) the surface defined by equating the expression inside the square bracket in $`|\mathrm{\Psi }_1|^2`$ ($`|\mathrm{\Psi }_2|^2`$) to zero. It is interesting to notice that the expressions mentioned above are of the same type as the one that defines $`𝒮`$0 (cf. Eq. (5)), and so a similar classification can be applied to these surfaces. The surfaces $`𝒮`$1 and $`𝒮`$2 define the boundary of the coexisting region of the $`|1`$ and $`|2`$ species. Depending on the value of the interparticle interaction strength, the regions where the overlapping of the waves function occurs exhibit different features. For $`D>0`$ and $`u_{1,2}`$ small enough, the coexisting region can be very large, as may be observed in the figures shown in Ref.. When increasing $`u_{1,2}`$ this region becomes narrower. At $`D=0`$, it is easy to verify that $`𝒮`$1 $`=`$ $`𝒮`$2 $`=`$ $`𝒮`$0, and hence the coexisting region is reduced to a single surface. When $`u_{1,2}`$ increases making $`D<0`$, this surface splits into two distinct surfaces, forming thus a new narrow 3D overlapping region, due to mutual repulsion of the species.

To have the problem completely solved one has to calculate the chemical potentials. In order to obtain $`\mu _1`$ and $`\mu _2`$ the normalization condition for the wave functions given by Eqs. (7) and (9) must be used. This leads to a system of two coupled integral equations that must be solved numerically. We have accomplished this in order to describe the states obtained in the JILA experiment , choosing parameters equal to those of the experiment. Namely, the trapping potential has an angular frequency $`w_0=2\pi \nu _{\text{trap}}`$, being $`\nu _{\text{trap}}=7.8`$ Hz. The interaction strengths relative to the $`|1`$ component of the 87Rb species are $`u_{1,1}=1,u_{1,2}=.97`$ and $`u_{2,2}=.94`$. And the number of particles of each species is $`N_1=N_2=4`$$`\times `$$`10^5`$. Note that for these values of the interaction strengths, $`D=9`$$`\times `$$`10^4`$ is close to the phase separation value, so the region between interfaces is very narrow, and hence $`𝒮`$1 $``$ $`𝒮`$2 $``$ $`𝒮`$0. Thus hereafter, in order to identify the type of segregation of species encountered, we will make use of the classification we performed for $`𝒮`$0.

For completeness, we also include calculations for the two species without vorticity, although the experimental realization consisted only of one vortex and one non-vortex state in each component.

In Fig. 2 we show the density contours of each wave function in the (r,z) plane for the following three configurations. (i) A vortex-free system $`(m_1,m_2)=(0,0)`$ yields chemical potentials in units of $`\mathrm{}w_0`$ as $`\mu _1=24.50`$ and $`\mu _2=24.06`$. It may be seen that in this case the segregation evolves as $``$ $`𝒮`$$`{}_{}{}^{a}{}_{0}{}^{}`$. The more repulsive component is on a spherical shell around the other, as expected. (ii) For $`(m_1,m_2)=(1,0)`$ we obtained $`\mu _1=24.54`$ and $`\mu _2=24.04`$. The vortex segregates in a ring surrounding the non vortex state. It is easy to verify that the species become separated by an interface of the type $`𝒮`$$`{}_{}{}^{b}{}_{0}{}^{}`$. The core size is about $`9\mu `$m and the whole vortex has a radius of $`26\mu `$m. This reproduces very well the experimental data. (iii) When the vortex is in the less repulsive component $`(m_1,m_2)=(0,1)`$, we obtained $`\mu _1=24.50`$ and $`\mu _2=24.15`$. In this case the vortex is segregated inside the non-vortex state. The corresponding surface that describes such segregation is the toroid $`𝒮`$$`{}_{}{}^{c}{}_{0}{}^{}`$. The inner and outer radii of the torus are $`5\mu `$m and $`18\mu `$m, respectively. However, this does not correspond to what was obtained at the first stage of the experiment, which looks quite similar to the pattern described in (ii). Mainly, in the experiment, the radius of the vortex was approximately $`26\mu `$m, while according to our calculations it should be $`18\mu `$m. Taking into account that after vortex creation, the vortex immediately shrank, a possible and likely explanation for the discrepancy is that the system was not created near the equilibrium configuration but moved towards it. The case with a vorticity of $`(m_1,m_2)=(1,0)`$ is created much closer to equilibrium . It is important to notice that in every case, even when the vortex configuration is out of equilibrium, the whole condensate is almost spherical with the same radius value of $`26\mu `$m.

In Fig. 3 we display the density profiles as a function of the radial coordinate at $`z=0`$, for the same three configurations described in Fig. 2. It may be seen in all of the cases that the shape of the total density is like a parabola, a characteristic that also holds for any other direction considered. Regarding the densities of each one of the components separately, they abruptly go to zero in the narrow region between species. It is worthwhile to mention that this behavior would not be so sharp for real systems and some interpenetrating region would exist .

With respect to the hypothesis made by Chui et al. , they assumed for condensates including vorticity in a symmetric trap that one component is confined inside a ball and the other one in a spherical shell around it; and that the species with vorticity have a small healing length sized vortex core. In the context of our article this means that they worked as if segregation of species, including vortices, always evolves as $``$ $`𝒮`$$`{}_{}{}^{a}{}_{0}{}^{}`$, which is not consistent with our results. In a very recent experimental work  Anderson et al. have included a figure (Fig. 3) in which a two-component vortex is viewed along its axis direction and along an orthogonal one. Views of a vorticity free condensate are also included. The in-trap pictures of Ref.  cannot rule out a healing length core because such a feature would be too small ( $`0.7\mu `$m ) to resolve. So in this experiment Chui’s condensates should look like a vortex free condensate. As it may be seen in the above mentioned figure the images of a system with and without vorticity are quite different. In particular, we want to note that in the vortex case there is a visible evidence of a cylinder type segregation. However to be conclusive about these statements more experimental information is required.

In Refs.  the authors employed, for describing the stationary states, wave functions expanded on a finite subset of the harmonic oscillator basis. This means that they did not work with exact solutions; unfortunately in these articles there is no precise information about the structure of the states. However, looking at both the domains of the wave function (Fig. 2) and the form of the density itself (Fig. 3), it seems that they should have taken into account a large number of elements of this basis to describe similar configurations.

It is well known that the TFA turns out to be very accurate when systems with a large number of particles are considered, so the main features of our analysis could be regarded as a guide to check the accuracy of different approximations.

Finally, we want to mention that very interesting experimental  and theoretical  work dealing with off-centered vortices in two-species condensates, is forthcoming.

In conclusion, we think that our analysis of states within the TFA could be helpful both from the experimental and theoretical viewpoints. In the first case, to prepare a configuration as close as possible to a stationary state, and in the second, to make a suitable choice for the set of wave functions in which such kinds of states are to be expanded.

We gratefully acknowledge useful conversations with B. P. Anderson.


text:
# Towards A Consistent Modeling Of Protein Thermodynamic And Kinetic Cooperativity: How Applicable Is The Transition State Picture To Folding and Unfolding?

## I INTRODUCTION

Our physical knowledge of protein folding is measured by the extent to which current ideas of elemental polypeptide interactions are capable of reproducing experimental data. Tremendous experimental progress has been made in the recent past.1-6 During the same time, investigations of simplified self-contained polymer models have contributed much physical insight.7-16 These models are successful in physically rationalizing many general features of proteins in terms of polymer properties, building a foundation for future advances. To move forward in theoretical development, it is necessary to recognize what common protein properties have not been predicted by models to date and target our efforts towards rectifying such deficiencies. A prime example is the current lack of chain models that can quantitatively reflect the extreme kinetic and thermodynamic cooperativity of small single-domain proteins.5,17 This highlights our insufficient understanding of protein energetics even at a “big-picture” level, suggesting that as heteropolymers natural proteins may be quite special in some respects.

We have recently investigated the severe constraints imposed on protein polymer models by the experimental observations of calorimetric and other hallmarks of thermodynamic two-state cooperativity.18-20 These cooperativity requirements appear to be more stringent than most other generic protein properties studied so far. A variety of flexible heteropolymer models with additive residue-based contact energies are able to explain significant aspects of the folding process.7-16 But such additive models — at least for the several examples evaluated to date — are found to be insufficient to satisfy the thermodynamic cooperativity requirements, even though deviations from proteinlike thermodynamics can be lessened in some instances by enhancing interaction heterogeneity through using larger numbers of letters and repulsive energies in model alphabets.19 As far as thermodynamic cooperativity is concerned, in scenarios tested thus far, we find that proteinlike thermodynamics can arise from nonadditive cooperative contributions that originate from an interplay between local conformational preferences and (mostly nonlocal) interactions that favor formation of protein cores.18,20

In the folding literature, an intimate correspondence between protein thermodynamics and kinetics has figured prominently in theoretical,7,21-23 modeling,24-27 and interpretative2,4 discourses. Therefore, we ask: Given that a heteropolymer model has already been constrained to satisfy a set of proteinlike thermodynamic properties, to what extent proteinlike kinetic behavior would follow automatically? For instance, would such a model be sufficient for two-state kinetics as observed for many small single-domain proteins5? More generally, what improvement in proteinlike kinetics would such a model enjoy over other models that are now known to be thermodynamically less cooperative?

In analyses of folding/unfolding kinetics experiments, free energy profiles are used extensively to provide useful insight2,4,28-33 and as a picturesque device to summarize data. However, other than the folding free energy and rate measurements themselves, independent experimental techniques to accurately define and determine such profiles are currently lacking. Moreover, many of these profiles were proposed without considering explicit chain representations. Therefore, the applicability and generality of their implied physical pictures remain to be ascertained. Coarse-grained protein chain models are well suited to shed light on this fundamental issue because they allow for broad conformational sampling. Free energy profiles in coarse-grained models can be obtained directly from chain population distributions, without regard to (and therefore independent of) kinetic rates. It follows that a rigorous evaluation of the applicability of transition state theory to protein folding can be conducted by comparing the transition-state-predicted rates and the actual kinetic rates in these models. We study one such model below.

## II A MODEL FOR THERMODYNAMIC COOPERATIVITY

The present analysis is based on a thermodynamically cooperative 55mer lattice protein model that folds to a ground state with a four-helix core (Fig. 1). The intrachain interaction scheme includes additive 5-letter contact energies,20,34 repulsive interactions disfavoring left-handed helices and sharp turns at the end of a helix, as well as cooperative “native hydrogen bond burial” terms20 (c.f. refs. 35, 36). The total energy $`E`$ is defined by Eq. (1) in ref. 20. Although “native-centric” interactions were introduced to enhance thermodynamic cooperativity in the present model, unlike the usual Gō construction, they are not necessary for recognizing the ground state. This is because the general, non-native-centric terms in the model (all terms in $`E`$ except the “native hydrogen bond burial” terms) are sufficient to provide global favorability to the proteinlike four-helix ground state. We note that several other studies36-39 have also emphasized cooperative interactions in protein folding; and nonadditive aspects of hydrophobic effects are being explored.40-45 As we have emphasized,20 although the present model is useful for exploring the issues at hand, it should be regarded as tentative, partly because it does not provide an explicit account of other possible physical origins of protein thermodynamic cooperativity such as sidechain packing.19,46,47 Furthermore, in view of the current lack of definitive understanding of hydrogen bonding energetics (see discussion and references in ref. 20), the cooperative “native hydrogen bond burial” energy in the present model should be broadly interpreted as representative of a general favorable coupling between local conformational preference and formation of proximate tertiary contacts, the physical mechanisms of which remain to be further elucidated.

In addition to the formulation in ref. 20, here we introduce a parameter $`ϵ`$ to model protein behaviors at different intrachain interaction strengths, such that the effective energy of a conformation with energy $`E`$ is equal to $`ϵE`$, hence its Boltzmann weight equals $`\mathrm{exp}(ϵE/k_BT)`$, where $`k_BT`$ is Boltzmann’s constant times absolute temperature. It follows that the partition function $`Q=_Eg(E)\mathrm{exp}(ϵE/k_BT)`$, where $`g(E)`$ is the number of conformations with energy $`E`$, and $`g(E)`$ is estimated by a parameter-space Monte Carlo histogram technique.20 The formulation in ref. 20 corresponds to $`ϵ=1`$. Because of the peculiar and significant temperature dependence of the solvent-mediated interactions in real proteins, varying $`ϵ/k_BT`$ serves better as a model for how effective intrachain interactions are modulated at constant temperature by denaturant concentration48 or denaturant activity37 than for how Boltzmann weights changes with temperature.11,49-51 Here, as a first approximation, denaturant effects are simply taken to be uniform over different interaction types. Models with different denaturant effects on different interaction types remain to be explored.

As we have discussed from a general polymer perspective,19 matching theoretical considerations with the experimental practice of calorimetric baseline subtractions necessitate a multiple-conformation native state that entails considerable fluctuations beyond small-amplitude vibrations. Here we further investigate the implications of native-state conformational diversity. To that end, we study different definition of “native” and “denatured” states by assigning different values for a “transition” energy $`E_t`$ demarcating the native and denatured ensembles, such that $`P_\mathrm{N}`$ $`=_{EE_t}g(E)\mathrm{exp}(ϵE/k_BT)/Q`$ is the fractional native population and $`P_\mathrm{D}=1P_\mathrm{N}`$ is the fractional denatured population (Fig. 1). For each of these definitions to be tested, a range of energies is spanned by the native state, except for the special case when $`E_t`$ is chosen to be equal to the ground-state energy. Remarkably, despite the differences in the definitions of “native” and “denatured” states, the thermodynamic ($`\mathrm{\Delta }G=0`$) transition midpoints of the different $`E_t`$’s in Fig. 1 are very similar, all at $`ϵ/k_BT1.56`$.

Stabilities of different denatured ensembles relative to that of the ground state are shown in Fig. 1 (dashed curves). These quantities correspond roughly to native state hydrogen exchange (HX) free energies,52-54 (see also ref. 55), for it is reasonable to expect that certain amides become exposed and exchangeable when the effective energy of a conformation is above a certain threshold. Our results share the same general trend as that observed in these experiments,52-54 suggesting that some of the fluctuations observed by HX may be considered as part of a multiple-conformation native state.18,19,56 Figure 1 indicates that linear extrapolation of $`\mathrm{\Delta }G`$ from the transition region to the strongly native regime (more negative $`ϵ/k_BT`$) is only valid for the free energy difference between the set of “fully unfolded” open conformations and the ground state (top dashed curve).

## III FREE ENERGY PROFILES AND CHEVRON PLOTS

Consistent with the model’s thermodynamically two-state character,20 its energy distribution is bimodal under denaturing conditions and moderately native conditions (Fig. 2), although the distribution becomes one-sided or “downhill”7,57,58 under strongly native conditions (e.g., when $`ϵ/k_BT=2.0`$). We emphasize that here the determination of $`\mathrm{ln}P(E)`$ is entirely independent of any kinetic consideration. Therefore, the present $`\mathrm{ln}P(E)`$ function reflects the actual thermodynamics of the model. As such, its physical origin is fundamentally different from free energy profiles that have been empirically constructed or postulated to fit rate data. We can therefore use these $`\mathrm{ln}P(E)`$’s to assess the transition state picture, with $`E`$ as the reaction coordinate. Different reaction coordinates have been used in other investigations.59,60

We employ standard Metropolis Monte Carlo dynamics19,22 to explore physically plausible kinetic scenarios, using the number of attempted moves as the model time. This approach has been proven useful7-11,13,61 despite its obvious limitations.11,16,50 The present set of elementary chain moves consists of end flips, corner flips, crankshafts,59-61 rigid rotations,11 and local moves that transform two turns of a right-handed helix among its three possible orientations while holding its two end monomers fixed. The relative frequencies of attempting these moves are 2.3%, 27%, 60.6%, 10%, and 0.1%, respectively. Some chain moves can lead to large changes in energy, hence movements along the model free energy profile need not be continuous (inset of Fig. 2). Therefore, it is more justified to regard the transition state as a region rather than a single highest point along this particular reaction coordinate.22,26,60 The group of conformations represented by the shaded area in Fig. 2 clearly serves the role of a transition state because all conformational interconvertions between the native and unfolded sides of the population distribution must pass through one or more conformations in the shaded area.

Figure 3 reports simulated mean first passage times59,60 for a range of intraprotein interaction strength on both sides of the transition midpoint, in a format identical to typical experimental chevron plots.29,31,32,49,51,62 We explore a variety of kinetic folding and unfolding criteria by monitoring the time it takes for the chain to first cross several different “finish lines.” This results in an appreciable variation in apparent rates under strongly native conditions (Fig. 3, more negative $`ϵ/k_BT`$). Similar effects may be operative when multiple experimental probes are used to monitor kinetics.63-65

The trajectory in Fig. 4 (upper panel) depicts the model’s heuristically “two-state” behavior at the transition midpoint. For the two chain properties shown, native and denatured parts of the trajectory can be easily discerned, with very little time spent in between; strongly suggesting that the kinetics is first order. Fluctuations in $`E`$ is considerable within the native (low $`E`$) part of the trajectory, underscoring the utility and necessity of a multiple-conformation native state (see Fig. 1 and below). Another facet of the native-denatured interconversion is provided by the $`R_g`$ trace. Consistent with a recent kinetic $`R_g`$ measurement on a small protein,66 it shows that the chain undergoes sharp kinetic transitions between a native state that has minimal fluctuations in $`R_g`$ and a denatured state that spans a wide range of $`R_g`$’s.

A more quantitative test introduced by Gutin et al.67 is performed in the lower panel of Fig. 4. It indicates that folding kinetics is essentially first order (i.e., single-exponential) for an extended range of model intraprotein interaction strength, covering moderately native conditions ($`ϵ/k_BT1.80`$) through conditions that are less favorable to folding (less negative $`ϵ/k_BT`$), although deviations from single-exponential behavior occur under strongly native conditions in the model ($`ϵ/k_BT<1.85`$). Using the same technique, unfolding kinetics (Fig. 3) is found to be essentially single-exponential for the entire range of unfolding $`ϵ/k_BT`$ investigated (detailed results not shown). We have confirmed these conclusions by analyzing first passage time (FPT) distributions as in ref. 68 at several $`ϵ/k_BT`$’s, paying special attention to folding kinetics under moderately native conditions that are not far from the onset of drastic chevron rollover and non-single-exponential behavior (Fig. 3). For example, we have obtained the logarithmic FPT distribution at $`ϵ/k_BT=1.72`$ by binning 1,080 simulated trajectories into time slots of $`10^6`$, and found that 98% of these trajectories can be fitted by a single exponential with a correlation coefficient $`r=0.95`$. If one assumes that the unit model time needed for each elementary chain move corresponds roughly to a real time scale of $`10^{11}`$$`10^9`$ sec (ref. 69), the fastest model folding rate in Fig. 3 is in the order of $`10^2`$$`10^4`$ sec-1.

The contrast between the present model and its corresponding Gō model is intriguing. We have shown that the Gō model in Fig. 4 is thermodynamically significantly less cooperative,20 yet it folds faster than our model. This scenario of a negative correlation between folding speed and thermodynamic cooperativity may bear on the issue of folding rate overestimation in folding theories that use Gō-like potentials.23 It also raises a more basic question as to whether and when the Gō prescription is sufficiently adequate for capturing minimal frustration21 mechanisms in real proteins.

The present model is proteinlike in that it predicts a mild chevron rollover concomitant with single-exponential folding kinetics, consistent with experiments.29,31,62 But the drastic chevron rollover (with an appreciable decreasing folding rate with increasing native conditions, and non-single-exponential folding kinetics) predicted for strongly native conditions in the present model ($`ϵ/k_BT<1.85`$ in Fig. 3) has not been documented for real proteins. This suggests that such conditions, which coincide with downhill folding7,57,58 in the present model (see above), may not be realizable. If so, this is not surprising. Native stability can be arbitrary high in the model ($`ϵ/k_BT`$ can be arbitrarily large), but for real proteins native stability is limited by the actual chemistry at zero denaturant. It follows that the experimental zero-denaturant situation for most proteins most likely corresponds to $`ϵ/k_BT>1.80`$ in Fig. 3. It would be interesting to explore whether special experimental situations corresponding to the very strongly folding conditions in the model can be found for some proteins such that similar drastic chevron rollovers can be observed.

## IV ASSESSING THE TRANSITION STATE PICTURE

Consistent with the thermodynamics of our model, kinetic rate of folding to the ground state and of unfolding to an open state meet at the approximate transition midpoint determined thermodynamically in Fig. 1 (c.f. the lower “V” in Fig. 3). Interestingly, a similar consistency is seen by matching rates of crossing the peak of the free energy profile in the inset of Fig. 2 in the folding and unfolding directions (upper “V” in Fig. 3). Near the transition midpoint, folding rates defined by crossing several different finish lines at low but non-ground-state energies are very close to one another, and are only slightly faster than the rate of folding to the ground state (Fig. 3). This implies that, under midpoint to moderately native conditions, kinetics is rapid once the folding chain has cleared the shaded transition state region in Fig. 2 and proceeds on to the native side. But the folding rates for different finish lines are very different under strongly native conditions ($`ϵ/k_BT2.0`$), indicative of glassy dynamics (Figs. 3 and 4).

Despite the essentially single-exponential and heuristically “two-state” kinetics discussed above, the folding/unfolding kinetics of the present model differs from the strictly two-state variety observed for an increasing number of small single-domain proteins4,5,24,25,70 such as a 64-residue form of chymotrypsin inhibitor 2. Figure 3 shows that folding rates under moderately to strongly native conditions are slower than that required for such strictly thermodynamically and kinetically two-state proteins (inclined dashed-dotted line). In fact, Fig. 3 is reminiscent of experimental chevron plots with rollovers.31,32,71 Examples include wildtype barnase29 and ribonuclease A.62 Hence we believe the present lattice construct may serve as a tool for better understanding the folding kinetics of these proteins.

How much kinetic information can be inferred from free energy profiles such as those in Fig. 2? The conventional transition state picture of protein folding2,28,30 stipulates that

$$\mathrm{rate}=F\mathrm{exp}\left(\frac{\mathrm{\Delta }G^{}}{k_BT}\right),$$
$`(1)`$
where $`\mathrm{\Delta }G^{}`$ is the activation free energy for the kinetic process in question. We call $`\mathrm{\Delta }G^{}/k_BT`$ the transition-state exponent. $`F`$ is the pre-exponential front factor11 or prefactor,23 which depends on solvent viscosity (not considered here) but is often taken to be insensitive to intraprotein interaction strength.30 Figure 5 examines whether this picture applies to the present model. It does so by investigating the dependence of $`F`$ on $`ϵ/k_BT`$. For the sake of generality, several physically motivated $`P(E)`$-based definitions of “transition state,” “folded state” and “unfolded state” are evaluated. In the tests conducted here, “transition state” is defined by either the shaded area in Fig. 2 ($`38E30`$) or $`E=34`$ (peak of barrier); “folded state” is defined by the ground state only ($`E=52`$) or by $`E34`$ (left of the barrier); and “unfolded state” is defined by the approximate position of the denatured free energy minimum ($`E=4`$), by the bulk of the open conformations ($`E>10`$), or by $`E34`$ (right of the barrier). Transition-state exponents for folding and unfolding are then computed, respectively, by the \[transition/unfolded\] and \[transition/folded\] population ratios at the given $`ϵ/k_BT`$ (Fig. 5).

The scaling of $`F`$ with respect to $`ϵ/k_BT`$ is found to be not sensitive to these variations in definition (Fig. 5, upper panel). Our results show that the simple transition state picture does not apply to folding in this model. For the quasi-linear part of the chevron plot in Fig. 3, the relationship between $`\mathrm{\Delta }G^{}/k_BT`$ and $`\mathrm{ln}(\mathrm{MFPT})`$ is approximately linear (Fig. 5), with slope $`1.5`$ (filled symbols) or $`1.9`$ (open symbols). This implies that $`F`$ has approximately the same functional form as the activation factor in this regime, but with an exponent of opposite sign, viz., $`F\mathrm{exp}(\mu \mathrm{\Delta }G^{}/k_BT)`$, where $`\mu `$ $`0.33`$ for folding to the ground state (filled symbols). On the other hand, unfolding appears to be well described by the simple transition state picture. The corresponding slope for unfolding $`1.0`$, hence $`F`$ constant. Folding and unfolding front factors are approximately equal near the thermodynamic midpoint (Fig. 5, lower panel), reflecting the fact that in that region folding and unfolding rates are essentially equal (Fig. 3). The role of front factors in understanding folding rates has recently been emphasized by Portman et al.23 The pattern in Fig. 5 is similar in some respects to the results of a recent model study by Nymeyer et al.22 Using a different reaction coordinate for 2-, 3-letter and Gō 27mer models, these authors found approximate linear relations with non-unity slopes between kinetically simulated rates and rate quantities deduced from free energy profiles, although they did not consider a broad range of overall interaction strengths as that in chevron plots.

Figure 5 shows that $`F`$ for folding decreases with increasingly native conditions. $`F`$ may be identified as an effective diffusion coefficient. It corresponds to an internal friction term arising from the impediment to motion imposed by the chain segments on one another.2,7,23,60,71-74 Plaxco and Baker72 have experimentally investigated internal friction in protein folding, and concluded insightfully that internal friction effects are limited for strictly two-state proteins. But the functional form they considered is different from the novel approximate exponential form $`F\mathrm{exp}(\mu \mathrm{\Delta }G^{}/k_BT)`$ discovered here for the quasi-linear part of chevron plots and the part with a mild rollover (corresponding to a limited range of $`ϵ/k_BT`$) for proteins that are not kinetically two-state. In the present model, when a wider range of $`ϵ/k_BT`$ is considered (Fig. 5, lower panel), the folding $`\mathrm{ln}F`$ reveals its nonlinear character, a feature anticipated by energy landscape theory21,23 and consistent with a pioneering simulation study of Socci et al.60

## V THE CALORIMETRIC CONNECTION: WHAT IS THE NATIVE STATE?

A thermodynamic consideration of the model’s free energy profiles (Fig. 2) and the above kinetic analysis suggest that a natural way to define the “native” and “denatured” states is to have their demarcation line at $`E_t=34`$. Figure 6 investigates the calorimetric implications of different $`E_t`$ choices. As a first test of principles, in this section we take the effective intraprotein energies as temperature independent. Our conclusions are not expected to be significantly altered by the incorporation of proteinlike temperature-dependent effective interactions.18

We have argued that empirical calorimetric baseline subtractions correspond essentially to an operational definition of the native and denatured ensembles.19 The demarcation energy (or enthalpy) between the native and denatured states may be ascertained by matching empirical baselines19 to the nonlinear “formal two-state” baselines of Zhou et al.75 Figure 6 shows that the baselines from empirically extrapolating the native and denatured tails of the heat capacity curve of the present model20 essentially coincide with the formal two-state baselines for $`E_t=34`$, implying that by adopting such empirical baselines $`E_t=34`$ is effectively adopted. In that case, the two chain conformations shown in Fig. 6 would belong to the native state and therefore may be regarded as sitting “below the calorimetric baseline.”19 These chains may model the sparsely populated conformations revealed by native-state HX.52-54 A multiple-conformation native state view is supported by a recent experimental observation of pretransitional conformational changes in ribonuclease A.76

## VI DISCUSSION

The logic of the present analysis is premised on a comparison between simulated folding/unfolding rates and transition-state predictions based on independently obtained free energy profiles. The conventional transition state picture of folding posits a weak or nonexistent dependence of the front factor on a protein’s intrachain interaction strength. In the present model, which exhibits a chevron rollover, we find that the conventional picture holds approximately for unfolding but not for folding. In particular, for the quasi-linear part of the folding arm of the chevron plot, the folding front factor adopts a re-scaled form of the exponential factor, harboring an exponent opposite in sign to that of the activation term. These findings are consistent with internal friction and diffusive folding dynamics ideas from energy landscape theory. They suggest that simple transition state theory with a constant front factor may not be generally applicable in the presence of a chevron rollover, even if the kinetics is apparently first order. In this regard, future single-molecule measurements77 of folding time distributions may provide important insight into the physics underlying approximate single-exponential folding kinetics, since these measurements may detect small deviations from first-order kinetics11 (e.g., a possible small non-exponential tail in the distribution) that would otherwise be difficult to ascertain from traditional measurements of ensemble-averaged folding rates.

Chevron rollovers have been rationalized by dead-time discrete intermediates29 and by movements of the transition-state peak on broad activation barriers.31 We have not been able to detect these proposed features in our model free energy profiles. Instead, the present results offer an alternate rationalization in terms of diffusive dynamics and an interaction-dependent folding front factor. It follows that, in general, analyses that focus exclusively on free energy profiles may be incomplete. Inasmuch as chevron rollovers are a manifestation of an interaction-dependent front factor, as suggested here, experimental observations of significant mutational effects on rollover behavior29 imply that mutations can have a significant effect not only on the free energy profile itself, but also on front factors not afforded by such profiles.

The present account of salient features of chevron rollover and native-state HX in terms of an essentially continuous energy distribution (Fig. 3) is similar in spirit to the recent idea that these features may originate from a “burst phase continuum.”32,33 However, the burst phase continuum view is based on postulated free energy profiles, not free energy profiles derived from models with explicit chain representations. Further effort will be required to elucidate the relationship between the burst phase continuum and the present chain-based perspectives, as there are apparent differences between the two. For instance, the present study suggests that some of the states detected by native-state HX are on the native side of the conformational distribution (Figs. 3 and 6) rather than on the denatured side as envisioned by the burst phase continuum scenario.

In summary, we emphasize that while the current study proposes a new physical rationalization for chevron rollover, it does not by itself rule out other mechanisms. Chevron rollovers in real proteins may arise from a combination of effects. Obviously, the generality of the present interaction-dependent front factor scenario should be further tested using model proteins with non-helical native topologies as well as using geometrically more realistic off-lattice continuum models.26,78

As for the relationship between generic features of folding/unfolding kinetics and thermodynamics of small globular proteins, the qualitative agreement between Fig. 3 and typical chevron rollover plots for real proteins supports the idea that proteinlike thermodynamics necessarily lead to proteinlike folding/unfolding kinetics. A case in point is the folding kinetics of a set of 20-letter model sequences reported by Gutin et al.67 Our test calculations show that random sequences of this particular 20-letter alphabet with additive contact energies are not calorimetrically two-state (data not shown). Although much useful insight has been gained from them10,14,68 (see also ref. 16), recent calculations19 indicate that even some designed sequences in this 20-letter model are thermodynamically less cooperative than the present model.20 Apparently, the folding kinetics of these 20-letter model sequences are less proteinlike as a result, in that their folding rates decrease when native stability is increased from the transition midpoint (Figs. 2, 5, 7 and 8 in ref. 67). This is because the maximum folding rates and the onset of drastic rollover in these 20-letter models occur around the thermodynamic transition midpoint, rather than under strongly native conditions as in the cooperative model studied here (Fig. 3 of the present work). Thus, the chevron trend predicted by these 20-letter models under transition midpoint through moderately folding conditions is opposite to that observed experimentally,49,51 because experiments almost invariably show an increasing folding rate when native stability is increased from the transition midpoint.

The present model’s kinetics is proteinlike but not two-state. In this respect, it is reassuring that the exercise here fares no worse than Nature’s. This is because a protein’s calorimetric two-state cooperativity,79 such as that of hen lysozyme,63,71 barnase,29 and ribonuclease A,62 is no guarantee for two-state kinetics.4,5 However, the present exercise does suggest that additional or alternate interaction mechanisms have to be discovered to account for the strictly two-state behavior of many small single-domain proteins. In that regard, it would be interesting to investigate the connection between the strictly two-state proteins’ apparently nonglassy kinetics80 and the possibility that their front factors might be minimally sensitive to the variation in intrachain interaction strength.

Acknowledgments

We thank Alan Davidson, Julie Forman-Kay, Yuji Goto, Bob Matthews, and Tetsuya Yomo for helpful discussions, and David Baker, Martin Gruebele, Walid Houry, Kevin Plaxco, and Peter Wolynes for their critical reading of an earlier draft of this paper and their very useful suggestions. This work was supported by Medical Research Council of Canada grant no. MT-15323 and a Premier’s Research Excellence Award (Ontario). H. S. C. is a Canada Research Chair in Biochemistry.


text:
# A Near-Infrared Wide-Field Proper Motion Search for Brown Dwarfs

## 1 Introduction

Although the sub-stellar initial mass function (IMF) has been studied in a range of environments such as star-forming clusters (Luhman et al. 1998; Luhman 2000; Najita et al. 2000), young open clusters (Bouvier et al. 1998; Barrado y Navascues et al. 2001) and the field (Reid et al. 1999), the IMF of low-mass companions is not well understood, especially at “wide” ($`>`$ 100 AU) separations. Radial velocity searches around solar-type main sequence stars (e.g., Mayor & Queloz 1995; Marcy & Butler 1996) have produced few confirmed brown dwarfs at separations $`<`$ 3 AU. Fewer than 0.5% of their sample have brown dwarf companions at those separations. A coronagraphic search for companions in the range 40-100 AU (Oppenheimer et al. 2001) produced only one brown dwarf, GJ 229B (Nakajima et al. 1995), well below the 17-30% multiplicity observed for all stars (Reid & Gizis 1997). Other types of surveys, such as high spatial resolution space-based observations (Lowrance et al. 1999; Lowrance et al. 2000) and ground-based adaptive optics (Els et al. 2001), have also resulted in discoveries of low-mass stellar and sub-stellar companions. However, the frequency of stellar and sub-stellar companions at close separations remains distinctly different, resulting in the idea that there is a “brown dwarf desert”.

To date there has been only one systematic search for brown dwarf companions at wide separations and with a volume-limited sample (Simons et al. 1996; hereafter, SHK). This was mainly a color-based search around M dwarfs within 8 pc of the Sun and did not turn up any new brown dwarfs, although, given the surprisingly blue colors of GJ 229B, cool brown dwarfs with intermediate $`J`$-$`K`$ colors may have been overlooked in the survey.

Proper motion searches for companions have been used for many years to identify low-mass objects (e.g., van Biesbroeck 1961) and offer a less biased way of finding low-mass companions than color-based surveys. Therefore, we have conducted the planned second epoch survey of the SHK sample, in order to identify low-mass companions to M dwarfs at wide separations out to over 1000 AU. The choice of M dwarf primaries is significant: Reid & Gizis (1997) and Reid et al. (1999) show that the distribution of mass ratios for a sample of 80% M dwarfs has a peak at $`q=`$ 0.95, where $`q`$ is the ratio of the secondary mass to the primary mass. They conclude that their sample shows a distinct bias towards approximately equal-mass systems and that the mass function for stellar companions is different from the IMF of field stars. If these conclusions extend to brown dwarf masses, M primaries may harbor more sub-stellar companions than other stellar types. On the other hand, Reipurth & Clarke (2001) suggest that brown dwarfs have been ejected by dynamical interactions during the star formation process and cannot accrete enough mass to become stars. In this case M dwarf primaries may not be accompanied by such companions except in a multiple M dwarf systems with a correspondingly large gravitational potential.

Thus, our proper motion search around one spectral class of primaries fills a unique niche in the search for low-mass stellar and brown dwarf companions. We describe the data acquisition and reduction in $`\mathrm{\S }`$ 2 and discuss the results of the survey in $`\mathrm{\S }`$ 3.

## 2 Observations and Data Reduction

### 2.1 Sample Selection

Our sample is identical to Henry’s (1991) list of M dwarfs within 8 pc of the Sun. The M dwarfs were chosen initially from the Second Catalog of Nearby Stars (Gliese 1969) and its updates (Gliese & Jahreiss 1979) along with other additions from more recent literature (e.g., LHS 292; see SHK for more details). The sample consists of 75 M dwarf primaries with MV$``$ 8.0 mag, trigonometric parallaxes $``$ 0$`\stackrel{}{\mathrm{.}}`$125 and declinations north of $`30\mathrm{°}`$. SHK observed 66 of these systems, discarding three due to confusion toward the galactic plane (GJ 701, GJ 729, and GJ 752), leaving a total of 63 systems.

Of the original 75 systems, we observed 74 targets. Table 1 contains a list of the targets, their proper motions, the dates of observation and those of SHK, and other relevant parameters. The median distance of the M dwarf targets is 5.8 pc, corresponding to a median search radius of 1480 AU in the present survey. Over the complete distance range, the search radius varies from 800-2100 AU. Figure 1 shows the distribution of total proper motions of our objects between the first and second epoch observations. Due to discoveries of new objects and to the measurement of more accurate parallaxes, the SHK list is no longer a complete volume-limited sample. Table 1 has five objects from the original sample whose updated parallaxes move them beyond the 8 pc limit (GJ 185, GJ 623, GJ 686, GJ 1230, and GJ 884) and separately lists four objects whose redetermined parallaxes or recent discoveries (e.g., G 180-060; Ducourant et al. 1998) place them within the survey criteria.

### 2.2 Imaging

First epoch images (SHK) were obtained at the University of Hawaii’s 24-inch telescope between 1991 August and 1992 August with a facility 256 x 256 NICMOS camera in both the $`J`$ and $`K^{}`$ bands with a scale of 2$`\stackrel{}{\mathrm{.}}`$0 pixel-1. Exposure times were typically 1 hour, and images were processed using conventional techniques. A custom program searched for point sources above a 3 $`\sigma `$ detection level and performed photometry on all sources using a 10$`\mathrm{}`$ aperture.

Between 1998 April and 2000 December, $`J`$-band images of the same fields were taken with the PISCES camera (McCarthy et al. 2001) at the Bok 2.3 m telescope on Kitt Peak. PISCES has an 8$`\stackrel{}{\mathrm{.}}`$5 diameter field-of-view and a 0$`\stackrel{}{\mathrm{.}}`$5 pixel-1 plate scale at this telescope. Nine 30 second exposures were obtained, centered on a program M dwarf, with a 10$`\mathrm{}`$ dither between each exposure. All images were corrected for quadrant cross-talk effects known to be present in HAWAII arrays (McCarthy et al. 2001). The images were then dark-subtracted, flat-fielded, masked for hot pixels, corrected for geometric distortion, and combined with standard IRAF tasks. The flat-field was produced through a median combination of the dithered science frames. The flat-field is predictably poor in the region of the bright M star, which was allowed to saturate the detector. However, the flat-fielding does not affect the astrometry of the field, changing the position of the M dwarf by less than half a pixel in multiple test cases. Also, because accurate photometry of the objects has already been carried out by SHK, the resulting flat-field is adequate for this survey.

Figure 2 shows a sample fully reduced field (GJ 752) along with the identical field observed by SHK. The field-of-view sizes are almost perfectly matched, except that PISCES has a circular inscribed field. The SHK images show an internal reflection due to the optics in the camera to the lower right of the M dwarf primary that is not in the second epoch set. The SHK survey has a limiting magnitude of mJ$``$ 16.5 and is sensitive to companions down to 40 MJupiter assuming an age of 5 Gyr at the mean distance of the survey, 5.8 pc, based on the models of Burrows et al. (1997). We use 5 Gyr, following the findings of Henry (1991) on this sample; however, because age dating M dwarfs is difficult, 5 Gyr may not be an accurate average age. If the M dwarfs are instead 1 Gyr old, the survey is sensitive to 16 MJupiter. The sensitivity of the PISCES images matches or exceeds that of the SHK images in all cases with mJ$``$ 17.0.

Each second epoch image is compared to its matching first epoch image using an IRAF script originally designed to identify supernovae in nearby galaxies (Van Dyk et al. 2000). Using input coordinates of identical objects in the two frames, the script matches the pixel scales of the two cameras, accounts for any differences in geometric distortion, matches the point spread functions for the two images, and subtracts them. Typically, 10-15 background stars from each M dwarf field were used as input references for the program. The total number of background stars ranges from 50 to over 1000 sources for the crowded fields near the galactic plane. The mean total proper motion of the M dwarfs between epochs is $`8\stackrel{}{\mathrm{.}}8`$ ($``$ 18 pixels on the PISCES camera). Moving objects are revealed by adjacent positive and negative images. Companion objects would have proper motion vectors identical to the M dwarf primaries which have been accurately measured by Hipparcos (Perryman et al. 1997). Objects with such large proper motions can easily be detected by visually examining the subtracted images.

Figure 3 shows the detection of the known low-mass companion van Biesbroeck 10 (VB 10, GJ 752B). Other known wide companions were also detected as indicated with an asterisk in Table 1. These results demonstrate the reliability of the subtraction method. GJ 570A, known to have a T dwarf companion at a separation of 258$`\stackrel{}{\mathrm{.}}`$3 (Burgasser et al. 2000), was not included in the sample because the primary is a K4 dwarf.

## 3 Results

No new low-mass stellar or brown dwarf companions were detected in this 8 pc sample of M dwarfs. The same conclusion was reached by SHK from a $`J`$-$`K^{}`$ color search. However, as many as nine wide (3600 $``$$`\mathrm{\Delta }`$$``$ 120 AU) companions have recently been detected around nearby (9.6-39 pc) stars using the Two Micron All-Sky Survey (2MASS) database (Kirkpatrick et al. 2000, 2001; Burgasser et al. 2000; Wilson et al. 2001). Seven of these new objects are common proper motion companions. In general, it is difficult to compare the 2MASS results with the present survey because the 2MASS primary stars have uncertain ages and generally higher luminosities. Based on initial 2MASS results, Gizis et al. (2001) estimate that $``$ 1% of primaries with masses 0.6-1.5 M (MV$`<`$ 9.5) have wide ($``$ 1000 AU) L dwarf companions and that the frequency of all wide brown dwarf companions is 5-13 times greater. Extending this analysis to our M dwarf sample of 63 objects with masses between 0.08-0.6 M, we would expect to detect between 3 (5% of our sample) and 9 (13%) brown dwarf companions. The apparent difference between our results and Gizis et al. (2001) might be resolved if the frequency of brown dwarf companions depends strongly on primary mass and orbital separation. This possibility could be tested either by systematic data mining of the 2MASS survey or by extending the present PISCES survey to other spectral types.

J. L. H. thanks Chien Peng, Rose Finn, and Roelof de Jong for valuable discussions concerning near-infrared data reduction and Michael Meyer for comments on early drafts of this paper. The HAWAII detector and electronic controller for PISCES were purchased by the Air Force Office of Scientific Research through grant F49620-96-1-0285. The PISCES camera was also briefly supported by the National Science Foundation through grant NSF 96-23788. D. A. S. was supported by the Gemini Observatory, which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the Particle Physics and Astronomy Research Council (United Kingdom), the National Research Council (Canada), CONICYT (Chile), the Australian Research Council (Australia), CNPq (Brazil) and CONICET (Argentina).