text:
# 1 Introduction

## 1 Introduction

It is commonly observed that under the influence of a flow of traffic, dirt roads develop regular corrugations of longitudinal “pitch,” or wavelength, between 0.5 and 1 m, and amplitude up to 50 mm . At first glance such an instability of the road surface might seem counterintuitive, as one might guess that a flow of traffic would exert downward forces on the surface which tend to smooth and compact the road bed, thereby suppressing pattern development rather than promoting it. Yet irregular, rough roads do not in fact heal themselves, but instead become progressively rougher and more corrugated under a flow of cars. Similar phenomena also occur, though sometimes on different length and time scales, on paved roads, on railroad rails , and on the rollers used to calendar paper . The purpose of this paper is to investigate this phenomenon using a simple, tractable physical model.

An early attempt at explaining road corrugation is due to Relton , who proposed that the underlying instability mechanism is a “relaxation oscillation,” caused essentially by stick-slip dynamics. According to this view, a moving wheel pushes grains ahead of it. The grains pile up in front of the wheel and form a heap. When the heap grows large enough, the wheel sticks momentarily, and then slips, running over the heap and leaving it behind as a ridge. For a given uniform speed, this stick-slip process will be fairly periodic, and so will generate equidistant ridges.

A different picture was provided by Mather , who argued that the origin of the road instability is the bouncing motion of the wheel, caused by random irregularities on the ground. When the bouncing occurs, the car is projected upward along a certain angle and is airborne for a brief time. When it then strikes the ground, the car creates a crater and the motion then repeats itself. According to this picture, it is not the piling up of grains ahead of the vehicle that is responsible for the instability, but the impact stress of the vehicle on the ground. The wavelength of the resulting corrugations will be determined by the competition between the typical distance the car flies over the ground and the size of the crater generated by the impact stress, which should in turn depend on the hardness of the ground and the relaxation time of the ejected grains. Mather’s picture is similar to other surface instabilities involving granular materials, in particular the ripple patterns in wind blown sand , where ejected grains are carried away by the wind and land in a place far from the ejection point. In such a nonlocal model, what sets the wavelength of the ripple is the ratio of the flux of the grains to the appropriately scaled saltation length, i.e., the distance that an ejected grain is carried by the wind.

The model we present in this paper builds on Mather’s picture, but we ignore any nonlocal transport of grains along the road, and we assume that the cars, and more specifically their wheels, generally do not lose contact with the road. Instead, we model the cars simply as masses attached to damped springs. We assume that the downward contact forces exerted on the road by the wheels causes a permanent downward deflection of the road surface, by either erosion or plastic deformation. In either case we model the effect of the contact force on the road as a simple proportionality between the deformation at any point on the road and the force exerted on that point, with a phenomenological “softness” parameter as the proportionality constant. We also include a diffusive relaxation that tends to smooth the road surface, which may come about, for instance, as a result of rain. We find that the diffusive relaxation is a stabilizing effect, as one would expect, but that its presence or absence generally has no qualitative effect on the corrugation phenomena. What is important, however, is the phenomenon of hardening. We expect that the softness parameter, and also the diffusion coefficient if it is present, will decrease as the road compacts, so that the road becomes less susceptible to the passage of more cars once the first several have compacted it and perhaps produced corrugations. This turns out to be a stabilizing effect, narrowing the parameter range in which corrugations occur.

Since we are modeling the cars, which in reality are complicated mechanical systems, by simple damped harmonic oscillators, the question quickly arises as to what the appropriate natural frequency might be. Given the observed pitch of road corrugations of 0.5 to 1 m, and the presumed speeds in the range of, say, 10 to 25 m/s of the cars that produce them, we expect that the important oscillation modes must have natural frequencies of about 0.02 to 0.1 s. This is one or two orders of magnitude shorter than the frequency of the car body bouncing on its suspension, so that is unlikely to be the relevant oscillation; indeed drivers normally adjust to the roughness of the road they are on and slow down to avoid such bouncing. This suggests that the relevant mode may be that of the wheel attached to the suspension. Alternatively, the important oscillation may be an elastic deformation of the wheel itself . This could account for an observed difference in the pitch of corrugations produced by hard and soft tires .

## 2 Model

Imagine a flux of cars traveling with some average horizontal velocity $`v_x`$ along a road surface whose height above some arbitrary zero level is given as a function of time $`t`$ and position $`x`$ along the road by $`H(x,t)`$; see Fig. 1. The cars are supported on springs, such that the natural angular frequency of vertical oscillation of the cars is $`\omega _0`$. Further, we assume that the springs are damped with damping constant $`b`$. Let $`Z(x,t)`$ be the height of a moving car, relative to a zero level chosen such that $`H(x,t)Z(x,t)`$ is the amount by which the springs are compressed. In a frame of reference moving horizontally with the car, the vertical component of the car’s equation of motion comes simply from Newton’s second law,

$`M{\displaystyle \frac{d^2}{dt^2}}Z(x,t)+b{\displaystyle \frac{d}{dt}}[Z(x,t)H(x,t)]=`$
$`M\omega _0^2[H(x,t)Z(x,t)].`$ (1)
The time derivatives here are total derivatives; that is, the relevant value of $`x`$ is time-dependent, since this equation follows the vertical motion of a single car. We may convert this into an equation for the height of the car at a fixed $`x`$ by replacing the total time derivative $`d/dt`$ with $`/t+v_x(/x)`$. Thus in a reference frame that is fixed to the road bed, the vertical equation of motion is

$`M\left({\displaystyle \frac{}{t}}+v_x{\displaystyle \frac{}{x}}\right)^2Z+b\left({\displaystyle \frac{}{t}}+v_x{\displaystyle \frac{}{x}}\right)(ZH)`$
$`+M\omega _0^2(ZH)=0.`$ (2)
Note that we have neglected any possible variation in the horizontal velocity $`v_x`$ of the cars.

We must now write an evolution equation for the height $`H(x,t)`$ of the road surface. We assume first that the road surface sinks at a rate which is proportional to the downward force on the road. This may be due to compaction of the road bed under the surface or to ejection of loose grains at the surface as a car passes; in either case we expect that the proportionality constant between the downward force on the road and its rate of sinking will decrease as the road sinks. That is, the road should harden as more and more cars pass over it. In addition, we assume that there is a “diffusive” relaxation process which tends to even out any roughness in the road surface. This can result from the action of wind or rain, and may also contain a contribution from the passage of the cars, as the loosely connected grains at the surface are fluidized by the flow of cars. Again, we expect the effective diffusion coefficient to decrease as cars pass and the road hardens. With these assumptions, the equation of motion for the road height reads

$$\frac{H}{t}=D(H)\frac{^2H}{x^2}a(H)\left[Mg+M\omega _0^2(HZ)\right],$$
(3)
where $`Mg`$ is the weight of a car. We are neglecting the geometrical distinction between $`H/t`$ and the rate of motion of the road surface normal to itself, and also the fact that the compressive force is not really vertical when $`H`$ is not constant. These are satisfactory approximations provided the vertical scale of the surface corrugations is much smaller than the horizontal scale. However, these effects should be included in any model of road corrugation which is nonlinear in the amplitude of the corrugation.

The proportionality factor $`a(H)`$ above represents the softness of the road, that is, its responsiveness to compressing forces. It should include the flux of cars as a multiplicative factor, as the rate of road sinking should also depend on the density of the traffic it sustains. The compression term should be replaced by zero whenever the quantity in brackets is negative, since that represents the situation in which the cars are airborne, so that the compressive force on the road would really be zero rather than negative. Requiring that the road harden as it compacts means that $`a`$ (and probably also $`D`$) should decrease as $`H`$ decreases, so we want $`a(H)`$ to be positive and an increasing function of $`H`$. A physically acceptable form of $`a(H)`$ is shown in Fig. 2.

Before proceeding, we first nondimensionalize the equations of motion for $`Z`$ and $`H`$. We choose the time scale to be $`1/\omega _0`$, the horizontal length scale to be $`v_x/\omega _0`$, and the vertical length scale to be $`g/\omega _0^2`$. The equation of motion for the cars then becomes

$$\left(\frac{}{t}+\frac{}{x}\right)^2Z+2\mathrm{\Gamma }\left(\frac{}{t}+\frac{}{x}\right)(ZH)+ZH=0,$$
(4)
and the equation for the road surface is

$$\frac{H}{t}=\mathrm{\Delta }(H)\frac{^2H}{x^2}\alpha (H)(1+HZ),$$
(5)
where the new dimensionless parameters are given by

$`\mathrm{\Gamma }=b/2M\omega _0,\mathrm{\Delta }(H)=(\omega _0/v_x^2)D(H),`$
$`\alpha (H)=M\omega _0a(H).`$ (6)
The parameter $`\mathrm{\Gamma }`$ is a property of the cars alone; the cars’ springs are underdamped for $`\mathrm{\Gamma }<1`$. Since the instability is controlled by the interplay of essentially two mechanisms, ejection of grains from the surface or compaction of the road by the flux of cars and subsequent relaxation of the road surface due to diffusion, $`\mathrm{\Gamma }`$ and $`\mathrm{\Delta }`$ will control the dynamics of the surface. As we will see, the hardening of the road also plays an important role in the development of the instability.

## 3 Linear Stability Analysis

We may solve the system of partial differential equations (4) and (5) numerically, and we will discuss this below, but useful information regarding the behavior of the system may also be extracted from an approximate linear stability analysis. The first step in the analysis is to recognize that as time increases, the spatially averaged values of $`Z(x,t)`$ and $`H(x,t)`$ will tend to decrease, reflecting the gradual settling of the road bed. We denote this spatially averaged, i.e. spatially independent, quiescent solution by $`H_0(t)`$ and $`Z_0(t)`$. Substituting this solution into the equations of motion, we find that it satisfies

$$\ddot{Z}_0+2\mathrm{\Gamma }(\dot{Z}_0\dot{H}_0)+Z_0H_0=0.$$
(7)
and

$$\dot{H}_0=\alpha (H_0)(1+H_0Z_0),$$
(8)
The full solution may now be written as

$$H(x,t)=H_0(t)+h(x,t),Z(x,t)=Z_0(t)+z(x,t),$$
(9)
in which the functions $`h(x,t)`$ and $`z(x,t)`$ will carry information about pattern development. Substituting these forms of $`H(x,t)`$ and $`Z(x,t)`$ into (4) and (5), using (7) and (8) for the time evolution of $`Z_0`$ and $`H_0`$, and expanding to first order in $`h`$ and $`z`$, we obtain the following linearized equations of motion:

$$\left(\frac{}{t}+\frac{}{x}\right)^2z+2\mathrm{\Gamma }\left(\frac{}{t}+\frac{}{x}\right)(zh)+zh=0,$$
(10)
$$\frac{h}{t}=\mathrm{\Delta }\frac{^2h}{x^2}\alpha (hz)\beta h,$$
(11)
where $`\alpha `$ and $`\mathrm{\Delta }`$ are evaluated at $`H=H_0`$, and $`\beta `$ is given by

$$\beta =(1+H_0Z_0)\alpha ^{}(H_0)=\frac{d\text{ln}\alpha (H_0)}{dt}.$$
(12)
The last equality here follows from (8). Note that $`\beta `$ should be positive, since we expect that $`\alpha `$ will increase with $`H`$ and decrease with time.

If $`\alpha `$ and $`\mathrm{\Delta }`$ are nontrivial functions of $`H`$, then it is not simple to solve the linearized equations (10) and (11), because the coefficients in the latter equation are functions of time. However, we may perform an approximate stability analysis by regarding $`\alpha `$, $`\beta `$, and $`\mathrm{\Delta }`$ as constants, at least for short time intervals, and calculating the linear growth rates that obtain when those parameters have their current values. Thus we will determine time dependent growth rates and phase velocities, provided we can ascertain the time dependent forms of $`\alpha `$ and $`\mathrm{\Delta }`$.

We proceed by assuming both $`h`$ and $`z`$ to be proportional to $`\mathrm{exp}(ikx+\sigma t)`$, where the linear growth rate $`\sigma `$ is in general complex. As usual, $`Re[\sigma ]`$ is the exponential growth or decay rate of the amplitude of a perturbation with wave number $`k`$, and $`Im[\sigma ]`$ is related to the phase velocity $`c`$ of that mode through $`c=Im[\sigma ]/k`$. This substitution yields

$$[(\sigma +ik)^2+2\mathrm{\Gamma }(\sigma +ik)+1]z=[1+2\mathrm{\Gamma }(\sigma +ik)]h$$
(13)
and

$$\sigma h=\mathrm{\Delta }k^2h\alpha (hz)\beta h.$$
(14)
Eliminating $`h`$ between these equations gives the stability relation

$$[(\sigma +ik)^2+2\mathrm{\Gamma }(\sigma +ik)+1](\beta +\sigma +\mathrm{\Delta }k^2)+\alpha (\sigma +ik)^2=0,$$
(15)
from which we can find the growth rates and phase velocities of spatially sinusoidal perturbations in terms of their wave numbers.

It is possible to determine the stability boundary for the model, at least parametrically, from (15). That is, we could determine the locus in parameter space on which the real part of the linear growth rate $`\sigma `$, as a function of $`k`$, has a global maximum at height $`Re[\sigma ]=0`$. However, it is much more instructive to simplify the already approximate problem further by taking $`\alpha `$ and $`\mathrm{\Delta }`$ to be small, as we expect to be the case once the road has had a chance to harden sufficiently. To set the stage for this calculation, imagine setting $`\alpha =0`$ in (15). The cubic equation for $`\sigma `$ would then factor. Two of the solutions would always have negative real parts, since $`\sigma +ik`$ for these two solutions would be a root of a quadratic with positive coefficients; these represent decaying perturbation modes. The third would be $`\sigma =\beta \mathrm{\Delta }k^2`$, so this mode too would always be stable unless $`\beta `$ is also small, on the same order as $`\alpha `$ or smaller. Note that small $`\alpha `$ does not necessarily imply that $`\beta `$ must be small, because $`\beta `$ is the logarithmic derivative of $`\alpha `$.

Since we are interested in finding parameter ranges for which corrugations grow, we now focus on the case where $`\alpha `$, $`\beta `$, and $`\mathrm{\Delta }`$ are all small. As we just saw, only one of the three solutions for $`\sigma `$ can possibly be positive; to leading order in the small parameters it is given by

$$\sigma =\beta \mathrm{\Delta }k^2+\alpha k^2\frac{(1k^2)2i\mathrm{\Gamma }k}{(1k^2)^2+4\mathrm{\Gamma }^2k^2}.$$
(16)
We may locate an approximate stability boundary by seeking combinations of parameters for which the real part of this $`\sigma `$ has a maximum, as a function of wave number $`k`$, at $`Re[\sigma ]=0`$. Thus we solve the equations $`Re[\sigma ]=0`$ and $`d(Re[\sigma ])/dk=0`$ to obtain the critical values of $`\beta /\alpha `$ and $`\mathrm{\Delta }/\alpha `$ parametrically in terms of $`k`$. The results are shown for several values of the damping parameter $`\mathrm{\Gamma }`$ in Fig. 3. As is clear from (16), both $`\beta `$ and $`\mathrm{\Delta }`$ represent stabilizing effects, so for a given $`\mathrm{\Gamma }`$, the flat road surface is stable above and to the right of the stability boundary curve for that $`\mathrm{\Gamma }`$. The stability boundaries move down and left as $`\mathrm{\Gamma }`$ is increased, so damping of the car’s springs also tends to stabilize the flat road surface. One can show explicitly that the stability boundary reaches $`\mathrm{\Delta }=0`$ at $`\beta /\alpha =[4\mathrm{\Gamma }(1+\mathrm{\Gamma })]^1`$, and it reaches $`\beta =0`$ at $`\mathrm{\Delta }/\alpha =[4\mathrm{\Gamma }(1\mathrm{\Gamma })]^1`$ for $`\mathrm{\Gamma }<1/2`$ or $`\mathrm{\Delta }/\alpha =1`$ for $`\mathrm{\Gamma }>1/2`$. Moreover, for small $`\mathrm{\Gamma }`$ – that is, when the oscillations of the car springs are very slightly damped – the stability boundary is given approximately by the line $`\mathrm{\Delta }+\beta =\alpha /4\mathrm{\Gamma }`$.

To find the wave numbers of the most important perturbations, we first note from (16) that the parameter $`\beta `$ only enters the expression for the growth rate additively. As a result, the wave number $`k_{max}`$ at which the real part of $`\sigma `$ has its maximum is independent of $`\beta `$; it is given explicitly by

$$\frac{\mathrm{\Delta }}{\alpha }=\frac{(1k_{max}^2)^24\mathrm{\Gamma }^2k_{max}^4}{[(1k_{max}^2)^2+4\mathrm{\Gamma }^2k_{max}^2]^2}.$$
(17)
Fig. 4 shows $`k_{max}`$ as a function of $`\mathrm{\Delta }/\alpha `$ for several values of the damping constant $`\mathrm{\Gamma }`$. (Below and to the right of the dotted curve, the maximum growth rate is negative for any positive value of $`\beta `$, so the flat road surface can only be unstable when $`\mathrm{\Delta }/alpha`$ is to the left of the dotted curve.) We see that the wave number of the most rapidly growing mode decreases with increasing damping $`\mathrm{\Gamma }`$; it also decreases with increasing $`\mathrm{\Delta }/\alpha `$, although the dependence on $`\mathrm{\Delta }/\alpha `$ is weak when $`\mathrm{\Gamma }`$ is small. In fact, for small $`\mathrm{\Gamma }`$ we always have $`k_{max}=1\mathrm{\Gamma }+O(\mathrm{\Gamma }^2)`$. In no case do we find $`k_{max}`$ to be above 1. Recall that $`k`$ is the wave number of the perturbation in units of $`\omega _0/v_x`$. Thus the fact that $`k_{max}`$ is always less than 1 means that the wavelength of the most rapidly growing perturbation is always greater than $`(2\pi /\omega _0)v_x`$, the distance a car travels in one period of its natural oscillation.

From the imaginary part of $`\sigma `$ we obtain the drift velocity $`c`$ of a perturbation of wave number $`k`$ (in units of $`v_x`$),

$$c=Im[\sigma ]/k=\frac{2\alpha \mathrm{\Gamma }k^2}{(1k^2)^2+4\mathrm{\Gamma }^2k^2}.$$
(18)
Note the factor $`\alpha `$: since $`\alpha `$ is small, the drift velocity is small compared to $`v_x`$, the speed of the cars. Also, since $`\alpha `$ contains a factor of the number of cars passing per unit time, the drift is actually a fixed distance per passing car, rather than per unit time. The drift velocity is positive, so that corrugations are pushed in the direction that traffic is moving as they develop. It reaches its maximum value of $`\alpha /2\mathrm{\Gamma }`$ at $`k=1`$.

## 4 Numerical Results

Since the linear stability analysis of the preceding section is only approximate, we have found it useful to check its results by carrying out numerical simulations of the full model given by Eqns. (4) and (5). We choose a wave number $`k`$, choose a sinusoidal form for the initial $`H(x)`$ and $`Z(x)`$ with small amplitudes, and then integrate the equations on an interval of length $`2\pi /k`$ with periodic boundary conditions. We the monitor the amplitudes and phase velocities of $`H`$ and $`Z`$ as time passes. The logarithmic derivative of the amplitudes with respect to time gives the instantaneous exponential growth rate $`\sigma (t)`$.

For the non-constant $`\alpha `$ cases we use the simple ansatz

$$\alpha (H)=\alpha _0\mathrm{exp}(ϵH),ϵ>0,$$
(19)
which has the qualitative form shown in Fig. 2 and is mathematically tractable. The parameter $`ϵ`$ controls the rate at which the soil hardens. We may solve for the time evolution of the quiescent road level $`H_0`$ provided $`\alpha `$ is small, for then $`Z_0`$ varies on a much shorter time scale than $`H_0`$, and so $`Z_0`$ should always remain near $`H_0`$. Then (8) reduces to

$$\dot{H}_0=\alpha (H_0)=\alpha _0\mathrm{exp}(ϵH_0).$$
(20)
This can be integrated immediately to yield $`\mathrm{exp}(ϵH_0)=ϵ\alpha _0(t+t_0)`$, where $`t_0`$ is a constant of integration. From this we get

$$\alpha =\frac{1}{ϵ(t+t_0)},$$
(21)
and from the definition (12) of $`\beta `$ we then get

$$\beta =\frac{1}{t+t_0}=ϵ\alpha ,$$
(22)
so that the combination $`\beta /\alpha `$ that appears in the approximate linear stability analysis is just the constant $`ϵ`$.

Figs. 5a and 5b show typical numerical and theoretical determinations of the growth rate $`Re[\sigma (t)]`$ and phase velocity $`c(t)=Im[\sigma (t)]/k`$ for $`\alpha `$ given by (19) above and $`\mathrm{\Delta }`$ constant, for various values of the hardening parameter $`ϵ=\beta /\alpha `$. As expected, for a given value of $`k`$ the agreement is quite strong; this continues to hold over a range of $`\alpha _0`$ and $`ϵ`$.

It is clear from the the linear stability analysis that if the diffusion parameter $`\mathrm{\Delta }`$ remains finite as $`\alpha `$ decreases to zero, or more generally if $`\alpha /\mathrm{\Delta }`$ goes to zero as the road compacts, then we will eventually reach a situation in which diffusion dominates the dynamics. From that point on, any corrugation in the road will decay. On the other hand, we certainly expect that as the road compacts and its surface hardens, the hardening should inhibit the lateral transfer of material which we have modeled as diffusion, as well as further compaction of the road. We see, then, that in order to get nontrivial corrugations on the road surface, we must have $`\mathrm{\Delta }`$ decrease at least as rapidly as $`\alpha `$ as the road compacts.

To investigate the possibility of generating corrugation patterns, we examine the simplest nontrivial case, namely where $`\mathrm{\Delta }`$ and $`\alpha `$ have the same $`H`$-dependence, $`\mathrm{\Delta }(H)=\mathrm{\Delta }_0\mathrm{exp}(ϵH)`$ and $`\alpha (H)=\alpha _0\mathrm{exp}(ϵH)`$. As above, this ansatz leads to $`\beta /\alpha `$ being constant – specifically, equal to $`ϵ`$ – and all three parameters $`\alpha `$, $`\beta `$, and $`\mathrm{\Delta }`$ decreasing as $`t^1`$ at long times $`t`$. From (16) we then see that the linear growth rate $`\sigma `$ is also proportional to $`t^1`$ for large $`t`$. Figures 6a and 6b show this slow decrease of the real and imaginary parts of $`\sigma (t)`$ for a few single mode solutions, arrived at by numerical solution and stability analysis.

Since $`\sigma `$ is the logarithmic time derivative of the amplitude of a perturbation, we see that the amplitude itself grows or decays algebraically in the long-time limit:

$$\frac{d\mathrm{ln}A}{dt}=\sigma =\frac{\sigma _0}{t+t_0}A(t+t_0)^{\sigma _0}.$$
(23)
Thus the growth rate of the corrugation amplitude becomes small at long times in this model, but the corrugation does not reach a true steady state.

It is clearly of interest to determine whether restabilized steady states exist. However, such states would depend on a host of nonlinear effects which are omitted in the model embodied in Eqns. (4) and (5). While we do not believe our model as presented is capable of predicting steady states, we have observed a remarkable feature in some of our numerical calculations. We have investigated cases in which $`\alpha `$ and $`\mathrm{\Delta }`$ are constants. These cases neglect the hardening effect which, as we saw above, tends to stabilize the flat road surface; the parameter $`\beta `$ vanishes. The linear stability analysis is then exact, and we expect to find purely exponential growth or decay of $`h`$ and $`z`$. Our numerical calculations show that modes that are stable according to the linear stability analysis do indeed decay in amplitude. Modes which are linearly unstable grow until their amplitudes are so large that we have $`1+H(x,t)Z(x,t)<0`$ for part of the cycle. When this happens, (5) suggests – unphysically – that the compression of the road is negative. What is happening here is that the cars are losing contact with the road. To avoid having the road surface spontaneously rise when an airborne car passes over it, we set $`\alpha =0`$ whenever we are in this situation.

Such a detachment in fact does occur in real situations, and has been termed “bouncing” . In a laboratory experiment involving two rotating disks that are in contact with each other under static compression, the two disks lose contact and bounce against each other when the amplitude of the corrugation along the perimeter of the disks exceeds about 1/3 the static compression. In the case of a car moving on a corrugated surface, such bouncing will kick the car into the air, but the car will quickly land in a different place, a mechanism suggested by Mather . Such a bouncing motion should involve a local flux of the cars and a saltation function that relates the landing position to the start-off position, as in the case of wind-blown sand . We note that such a nonlocal behavior would be difficult to account for in the model completely, but our mechanism of setting $`\alpha `$ locally to zero captures some of its flavor. Our expectations are that the modes that grow large enough to be subject to the local removal of the the force term in (5) will obviously have their growth rates reduced from those predicted by the linear stability analysis in the absence of such a force removal, and that this reduction in growth rate may be sufficient to generate steady states with finite amplitudes. To test this intuition, we numerically solve the system using single cycles of sinusoidal modes $`k`$, in periodic systems whose lengths are a single wavelength. The initial disturbances are chosen with amplitudes large enough to cause $`1+H(x,t)Z(x,t)`$ to be negative on occasion.

Provided a given mode is unstable according to the linear stability analysis, our numerical analysis indicates that it evolves toward a quasi-steady state in which the amplitude and phase velocity $`c`$ of the road corrugations remain fixed, even while the average height of the road bed continues to sink. Fig. 7 shows the steady state amplitudes of $`h`$ as a function of $`k`$ for $`\mathrm{\Gamma }`$, $`\mathrm{\Delta }`$, and$`\alpha `$ as given in the figure caption. For fixed $`k`$, the steady state amplitudes increase with increasing $`\alpha `$. Thus we see for this choice of parameters that softer roads (i.e., larger $`\alpha `$) support larger disturbances in the steady state than do harder ones. This result is intuitively satisfying, as we expect the ease of deformation of the material to correlate with the amplitude of the disturbance it supports. The phase velocities of the quasi-steady states exhibit particularly interesting behavior. Fig. 8a shows the observed quasi-steady state phase velocities, plotted as symbols, and an “envelope” in the $`kc`$ plane. If there were no step function in the term in $`\alpha `$, the phase velocities would lie in the region bounded by the envelope. The top curve would be the velocity curve for the largest $`\alpha `$ used (0.0055) and the bottom curve would be the velocity curve for the smallest $`\alpha `$ used (0.0040). The step function in the compactivity has forced the phase velocities to “collapse” so that they seem to lie along or near a single curve in the $`kc`$ plane. Examination of the quasi-steady state velocities on a finer scale, as shown in Fig. 8b, indicates a persistence of the kind of variation of phase velocity with softness we have typically seen (that is, the softer the road, or equivalently the larger $`\alpha `$ is, the greater the phase velocity, all other things being equal), though this variation exists on a much finer scale in the steady state case than in the case in which the step function is not invoked.

Even though we have identified interesting steady state behavior in the solution of the differential equations in this limiting case of $`\alpha `$ being constant, we must reiterate that setting $`\alpha `$ to a constant is certainly unphysical, and that many nonlinear effects have been omitted from our model equations. Moreover, our model does not account carefully for what happens when contact with the road and cars is lost; it merely turns off the coupling between cars and road in the equation of motion for $`H(x,t)`$, and furthermore assumes the evolution of $`Z(x,t)`$ remains well described by Eq. (4). Thus we cannot expect that our nonlinear calculations will speak to the behavior of real roads. Nonetheless, we find it remarkable that simply setting $`\alpha =0`$ when $`1+H(x,t)Z(x,t)`$ is negative seems to contain the growth of the unstable modes.

## 5 Conclusions and Prospects

In this work, we have explored the origin of the corrugation instability of dirt roads subjected to a constant flux of cars. We have presented a simple phenomenological model for the evolution of the road surface and carried out a linear stability analysis to uncover the gross features of the instability. ¿From our approximate stability analysis and its numerical verification, we find that the dynamical processes of diffusion and compression in a road bed are sufficient to generate instabilities that can give rise to pattern formation. Small-amplitude corrugations on the road surface grow when the diffusion parameter $`\mathrm{\Delta }/\alpha `$ and the hardening parameter $`\beta /\alpha `$ lie below the stability boundary (for the appropriate damping coefficient $`\mathrm{\Gamma }`$ for the cars) shown in Figure 3. From the definitions (6) of $`\alpha `$ and $`\mathrm{\Delta }`$ in terms of the original, dimensional softness and diffusion coefficients $`a(H)`$ and $`D(H)`$, we see that the dimensionless combination $`\mathrm{\Delta }/\alpha `$ is inversely proportional to the square of $`v_x`$, the horizontal speed of the cars across the road. Thus we find that there is a critical car speed above which the flat road surface is unstable. Also, since $`a(H)`$ is proportional to the flux of cars, there is a critical flux for any given speed, above which the flat road is unstable. Both diffusive relaxation and hardening of the surface are stabilizing effects. The wave number of the most rapidly growing mode is given implicitly in dimensionless terms by (17), which is plotted in Fig. 4.

Our model is clearly schematic; a quantitative theory of road corrugations would need to incorporate a number of effects we left out of our calculations. We have not attempted to account, for instance, for the geometric distinction between $`h/t`$ and the normal velocity of the road surface as it compacts, nor for the fact that the force exerted by a moving car on a corrugated road surface is not purely vertical. For these reasons alone we have not carried out any nonlinear analysis on our model – the model is explicitly not valid beyond linear order in the amplitude of the incipient corrugation. We have also ignored all details of the physical processes by which the road compacts, so it is not clear how well our phenomenological picture of compaction being proportional to applied vertical force or impulse can describe the dynamics of a real road. Our model for the cars is the simplest possible, neglecting the multiplicity of oscillation modes of real cars and the Hertzian nature of the contact force between the tires and the road. In particular, the fact that real vehicles have two or more pair of tires separated by fixed distances may introduce a new, relevant length scale into the problem. Finally, a quantitative theory would have to incorporate distributions of vehicle sizes, weights, oscillation frequencies, speeds, and the like.

Despite the fact that most nonlinear effects are omitted from our model, we find particularly interesting the numerical result that simply setting the contact force to zero when the equations would make it negative can lead to an apparently steady state. At least this results in a drastic slowing of the dynamics once a certain road configuration has been established. Continuing work will focus on determining the selection mechanism of such “frozen” states, and their response to local perturbations.


format: text

text:
# Surface soft phonon and the √3×√3↔3×3 phase transition in Sn/Ge(111) and Sn/Si(111).

\[

## Abstract

Density Functional Theory (DFT) calculations show that the reversible Sn/Ge(111) $`\sqrt{3}\times \sqrt{3}3\times 3`$ phase transition can be described in terms of a surface soft phonon. The isovalent Sn/Si(111) case does not display this transition since the $`\sqrt{3}\times \sqrt{3}`$ phase is the stable structure at low temperature, although it presents a partial softening of the $`3\times 3`$ surface phonon. The rather flat energy surfaces for the atomic motion associated with this phonon mode in both cases explain the experimental similarities found at room temperature between these systems. The driving force underlying the $`\sqrt{3}\times \sqrt{3}3\times 3`$ phase transition is shown to be associated with the electronic energy gain due to the Sn dangling bond rehybridization.

\]

A soft phonon is a vibrational mode of a crystalline material whose frequency decreases as T falls, eventually reaching zero. At this point the crystal is unstable in relation to the corresponding atomic displacements and undergoes a transition to a lower symmetry phase. Typical examples are the high–low transition in quartz, the ferroelectric transition in $`BaTiO_3`$ and the ferroelastic transition in $`SrTiO_3`$ . On semiconductor surfaces, different mechanisms, like Peierls transitions or Charge Density Waves, have been proposed as the origin of two dimensional phase transitions. However, no soft phonon transition has been unambiguously identified in these systems so far.

In this letter we consider the Sn/Ge(111) $`\sqrt{3}\times \sqrt{3}3\times 3`$ phase transition, which is a prototypical example of a T-induced reversible change of symmetry at a semiconductor surface. The driving force underlying this phase transition has been under intense debate since its discovery on the Pb/Ge(111) system . Different physical mechanisms have been proposed to play an important role in this transition: surface Fermi wavevector nesting or electron correlations, leading to the formation of a Charge Density Wave  at low T; dynamical fluctuations of an underlying $`3\times 3`$ structure, giving rise to the observation, on average, of a $`\sqrt{3}\times \sqrt{3}`$ symmetry at room T ; and the interaction between the $`3\times 3`$ periodicity and Ge–substitutional defects that act as nucleation centers. Recently, several groups  have also considered the isovalent Sn/Si(111) system and have found strong similarities and some striking differences with the Ge case. The electronic structure is very similar in both cases, with two major components in the Sn $`4d`$ core level spectra and two surface bands close to the Fermi level , associated with two different Sn–dangling–bonds. At variance with the Sn/Ge(111), where a transition to a $`3\times 3`$ pattern is clearly observed at low T, the Sn/Si(111) surface shows a $`\sqrt{3}\times \sqrt{3}`$ pattern both in LEED and STM, for temperatures as low as 70 K , even in the presence of a significant density of Si–substitutional defects . Moreover, the two surface bands mentioned above are not as clearly resolved at low T as they are in the Sn/Ge(111) case. By analyzing theoretically the Sn/Ge and Sn/Si systems, we conclude that the $`\sqrt{3}\times \sqrt{3}3\times 3`$ transition is due to a surface soft phonon. The frequency of this mode in Sn/Ge(111) goes to zero for a $`\overline{k}`$–vector corresponding to a reciprocal lattice vector of the $`3\times 3`$ periodicity. For Sn/Si(111) we find that the $`\sqrt{3}\times \sqrt{3}`$ is the stable structure at low T, and therefore this system should not display the $`\sqrt{3}\times \sqrt{3}3\times 3`$ transition. The comparison of the Sn/Ge and Sn/Si systems with their corresponding H-covered cases shows that the driving force for this transition is the electronic energy gain associated with the surface band splitting induced by the $`3\times 3`$ distortion.

The atomic displacements associated with the $`\sqrt{3}\times \sqrt{3}3\times 3`$ transition in Sn/Ge(111) provide the clue for the soft mode responsible for this transition. The signature of the $`3\times 3`$ structure is the upward displacement of one of the three Sn-adatoms, and the corresponding downward movement of the other two . Figure 1 and Table I show how these upward and downward movements are interconnected: as one of the Sn atoms ($`Sn_1`$) moves upwards, the three Ge nearest neighbours ($`A_1`$) follow this motion by moving towards it (with both upward and in–plane displacements) . These Ge displacements force corresponding in-plane movements of the Ge (labelled D2) towards the position of $`Sn_1`$. Since the Ge–D2 atoms are also bonded to first–layer Ge–atoms ($`B_1`$,$`C_1`$) linked with the two other Sn atoms ($`Sn_2`$ and $`Sn_3`$), these are forced to move downwards.

Does the Sn/Si(111) surface present a similar $`3\times 3`$ structure at low T?. We have studied this possibility by means of DFT calculations . In this analysis, we start from a $`\sqrt{3}\times \sqrt{3}`$ structure (in a $`3\times 3`$ unit cell), where the 3 Sn atoms are equivalent, we select one of them ($`Sn_1`$) and force it to move in the direction perpendicular to the surface. For each of these displacements, the other Sn atoms ($`Sn_2`$ and $`Sn_3`$) and all the semiconductor atoms are allowed to relax (up to the fifth layer) to their zero force positions under the constraint of the $`Sn_1`$ displacement. These calculations were performed for both the Sn/Si and Sn/Ge systems (see Figure 2). The energy vs. $`Sn_1`$ displacement is very flat for Sn/Ge(111); the ground state corresponds to a $`3\times 3`$ geometry with a $`Sn_1`$ displacement of $`0.18`$ Å, as found in previous calculations  and in agreement with the experimental evidence . The Sn/Si(111) surface, however, does not present any stable $`3\times 3`$ distortion, i.e., the $`\sqrt{3}\times \sqrt{3}`$ is the stable structure at low T.

The inset of Fig. 2 shows another interesting finding: $`Sn_2`$ and $`Sn_3`$ move together in the opposite direction to the constrained motion of the $`Sn_1`$ atom, for both the Sn/Ge(111) and Sn/Si(111) surfaces. These results suggest a correlated up/down motion of the Sn atoms  that, at a sufficiently low temperature, is frozen into the $`3\times 3`$ structure in the case of Ge (but not in the case of Si). As reflected in Table I, the displacements associated with this motion are essentially localized in the tetrahedra formed by a Sn atom and the three Ge of the first layer bonded to it. In order to analyze the corresponding normal modes associated with the up/down motion of the tetrahedra, we have considered a force constant model which includes effective interactions between first nearest neighbour tetrahedra:

$$F_i=\alpha z_i+\beta \underset{j,nn}{}(z_jz_i)$$
(1)
where $`F_i`$, the force acting on each tetrahedron, is proportional to its center–of–mass displacement (with respect to the $`\sqrt{3}\times \sqrt{3}`$ structure), $`z_i`$, and the differences in displacement between nearest neighbours, $`z_jz_i`$.

The force constants $`\alpha `$ and $`\beta `$ can be determined from the calculations presented above. These calculations provide the accurate information required for this purpose, since they are focused precisely on the deformations we are analyzing here. Imposing that $`F_2`$ and $`F_3`$ are zero in the above equations yields a relation between the force constants and the displacements $`z_1`$ and $`z_2(=z_3)`$: $`3\beta /\alpha =z_2/(z_2z_1)`$. For small displacements, ($`z_2`$) is around 0.5$`z_1`$ for Ge, while it is only 0.35 $`z_1`$ in Si. Therefore, we obtain $`3\beta /\alpha =1/3`$ for Ge and $`3\beta /\alpha =0.23`$ for Si.

The phonon dispersion relation of this mode is easily found, using eqn. 1, to be:

$$M\omega ^2=(\alpha +6\beta )\beta \underset{j}{}cos(\overline{k}\overline{R}_j)$$
(2)
where $`\overline{k}`$ is the momentum parallel to the surface, $`\overline{R}_j`$ the coordinates of the six first nearest neighbours, and $`M`$ the tetrahedron mass. In this approach, the atoms outside the tetrahedra are treated adiabatically, i.e. we neglect their mass and assume them to follow “instantaneously” the displacement imposed by the motion of the tetrahedra. This means that the forces between, say, the Ge-D2 atoms (see fig. 1) and the atoms in the tetrahedra are automatically included in the force constants of eqn. 1 and in the phonon dispersion relation of eqn. 2. We have estimated the error introduced by this adiabatic approximation in the phonon frequencies by comparing the kinetic energies of the atoms whose inertia is neglected in our approach with that of the atoms in the tetrahedra. This yields errors of 7% and 16 % for the Si and Ge cases respectively. These values set the order of magnitude of the accuracy that our model presents for the phonon frequencies.

Figures 3a and 3b show the phonon dispersion curves for Ge and Si, compared with the projected phonon bulk band structure. Notice that we find a zero-frequency mode at the $`\overline{K}^{^{}}`$ point for Ge. The corresponding $`\overline{k}`$–vector defines the new 3$`\times `$3 periodicity associated with this soft mode. For Si, the phonon dispersion relation shows only a minimum at the $`\overline{K}^{^{}}`$ point, with a phonon energy of 5.5 meV. These results are in agreement with our previous comments on the stability of the Sn/Si(111)-$`\sqrt{3}\times \sqrt{3}`$ structure, as shown in Figure 2.

The phonon branches shown in Fig. 3 provide an explanation for the main similarities and differences between the Si and Ge cases. At low T, while for Ge we can expect a phase transition to the 3$`\times `$3 structure, there is no stable 3$`\times `$3 phase for Si. On the other hand, at high T (say, for $`k_BT>>6`$ meV, room T is a high T limit) we can expect Ge and Si to be alike. This can be understood by considering the energy curves in Fig. 2: at room T the Sn atoms vibrate with a large amplitude in both surfaces (with displacements between -0.2 and 0.3 Å for Ge and -0.2 and 0.2 Å for Si). In this vibration, Sn atoms are correlated with their Sn–nearest–neighbours. This is true even in the Si case, because although the 3$`\times `$3 structure is not an stable minimum, due to the fact that the vibration associated with the 3$`\times `$3 structure is a minimum of the phonon dispersion, the atomic motion can be expected to be dominated by this mode. For the Sn/Ge(111)–$`\sqrt{3}\times \sqrt{3}`$ phase, anharmonic effects renormalize the frequency of the soft-phonon mode to a non-zero value, as shown schematically in Figure 3.

The atoms involved in this vibration spend most of the time around the turning points of their classical trajectories, where the atom velocity goes to zero. At these turning points, the electron occupation of the Sn dangling bonds is different: atoms in “up” positions have a fully occupied dangling bond, while those at “down” positions present only a partial occupation of their dangling bonds. This explains the double peak observed in the XPS Sn $`4d`$ core level spectra of both systems, as originally proposed in the dynamical fluctuations model to explain the Sn/Ge(111) $`3\times 3\sqrt{3}\times \sqrt{3}`$ transition. In this model, Sn atoms display correlated up/down vibrations which present, at RT, large amplitudes, but keep memory of the underlying 3$`\times `$3 phase. The soft phonon shown in fig. 3 is responsible for this atomic motion and provides the physical mechanism for the dynamical fluctuation model.

We have also addresed the problem of which driving force is softening the surface mode associated with the 3$`\times `$3 structure. As Si and Ge are isovalent atoms and the electronic structure of both interfaces is very similar, one might expect the difference in the relative stability of the two structures to be related to the different elastic properties of the two materials. In order to single out the elastic contribution from the one associated with the Sn–dangling–bond rehybridization  and corresponding surface–band splitting, we have studied the surface phonons of these surfaces with hydrogen atoms saturating the Sn–dangling–bonds. As shown in Figs. 3a and 3b, the dispersion curves for the H–saturated surfaces are rather different from the Sn–surfaces; in particular they display no “softening” of the mode at $`\overline{K}^{^{}}`$. This different behaviour is related to the Sn–dangling–bonds rehybridization.

In order to quantify these effects, we have considered the relaxed Sn–surfaces discussed above (see Fig. 2) for two different $`Sn_1`$ displacements, namely, 0.0 Å(the $`\sqrt{3}\times \sqrt{3}`$ geometry) and 0.18 Å(the “3$`\times `$3” geometry). The elastic contribution is calculated, in each case, as the difference in energy between the two geometries with H–atoms saturating the Sn–dangling–bonds. Substracting this elastic contribution from the total energy difference for the Sn–surfaces, the energy contribution due to the Sn–dangling–bonds rehybridization is obtained. These calculations provide the following results: $`E_{elastic}(Ge)=156`$ meV and $`E_{electronic}(Ge)=158`$ meV; $`E_{elastic}(Si)=131`$ meV and $`E_{electronic}(Si)=112`$ meV . Therefore, the main difference between the Sn/Ge and Sn/Si cases is the different hybridization between Sn and the surface atoms. A comparison of the surface bands for both cases confirms this result: the Ge case presents in the 3$`\times `$3 geometry a larger band splitting, 152 meV as compared with 93 meV for Si. This effect can be related to the different sizes of Si and Ge . When Sn is closer in size to the semiconductor atom the surface rehybridization is stronger and the surface gains more energy by means of the surface band splitting caused by the 3$`\times `$3 distortion.

In conclusion, we have shown via state-of-the-art DFT calculations that the reversible $`\sqrt{3}\times \sqrt{3}3\times 3`$ phase transition in Sn/Ge(111) is associated with a surface soft mode. The $`3\times 3`$ structure is not stable in the Sn/Si(111) case, although it is a minimum in the phonon dispersion, and consequently this surface should not display the $`\sqrt{3}\times \sqrt{3}3\times 3`$ transition. Comparing the Sn/Ge and Sn/Si cases, and considering the corresponding H-covered surfaces, we have found that the driving force underlying the $`\sqrt{3}\times \sqrt{3}3\times 3`$ phase transition lies in the electronic energy gain due to the Sn dangling bond rehybridization.

This work has been partly funded by the spanish CICYT under contract No.PB-97-0028. We acknowledge helpful discussions with V.R. Velasco and J.M. Pérez-Mato.


text:
# Not enough stellar mass Machos in the Galactic halo

## 1 Research context

The search for gravitational microlensing in our Galaxy has been going on for a decade, following the proposal to use this effect as a probe of the dark matter content of the Galactic halo . The first microlensing candidates were reported in 1993, towards the lmc  and the Galactic Centre  by the eros, macho and ogle collaborations.

Because they observed no microlensing candidate with a duration shorter than 10 days, the eros1 and macho groups were able to exclude the possibility that more than 10% of the Galactic dark matter resides in planet-sized objects .

However a few events were detected with longer timescales. In their two-year analysis , the macho group estimated, from 6-8 candidate events towards the lmc, an optical depth of order half that required to account for the dynamical mass of a “standard” spherical dark halo111 $`4\times 10^{11}\mathrm{M}_{}`$ within 50 kpc; the typical Einstein radius crossing time of the events, $`t_E`$, implied an average mass of about 0.5 M for halo lenses . Based on two candidates, eros1 set an upper limit on the halo mass fraction in objects of similar masses , that is below that required to explain the rotation curve of our Galaxy222 Assuming the original two eros1 candidates are microlensing events, they would correspond to an optical depth six times lower than that expected from a halo fully comprised of machos..

The second phase of the eros programme was started in 1996, with a ten-fold increase over eros1 in the number of monitored stars in the Magellanic Clouds. The analysis of the first two years of data towards the Small Magellanic Cloud (smc) allowed the observation of one microlensing event  also detected by . This single event, out of 5.3 million monitored stars, allowed eros2 to further constrain the halo composition, excluding in particular that more than 50 % of the standard dark halo is made up of $`0.010.5\mathrm{M}_{}`$ objects . In contrast, an optical detection of a halo white dwarf population was reported , compatible with a galactic halo full of white dwarfs.

Very recently, the macho group presented an analysis of 5.7 year light curves of 10.7 million stars in the lmc  with an improved determination of their detection efficiency and a better rejection of background supernova explosions behind the lmc. They now favour a galactic halo macho component of 20% in the form of 0.4 M objects. Within a few days, the detection of a halo white dwarf population at the level of a 10% component was also reported by . Simultaneously, the eros2 group presented its results from a two-year survey of 17.5 million stars in the lmc . One eros1 microlensing candidate, eros1-lmc-2, was seen to vary again, 8 years after its first brightening, and was thus eliminated from the list of microlensing candidates. Two new candidates were identified (eros2-lmc-3 and 4). Because this is much lower than expected if machos are a substantial component of the galactic halo, and because these two new candidates do not show excellent agreement with simple microlensing light curves, eros chose to combine these results with those from previous eros analyses, and to quote an upper limit on the fraction of the galactic halo in the form of machos.

In this article, we describe an update on the eros2 lmc data, an analysis of the three-year light curves from 25.5 million stars. While the sensitivity is improved, the main conclusions are unchanged compared to . One of the two-year candidates was seen to vary in the third season and was thus rejected. Three new candidates have been detected. We combine these eros2 lmc results with those of previous independent eros analyses, and derive the strongest limit obtained thus far on the amount of stellar mass objects in the Galactic halo.

## 2 Experimental setup and LMC observations

The telescope, camera, telescope operation and data reduction are as described in . Since August 1996, we have been monitoring 66 one square-degree fields in the lmc, simultaneously in two wide passbands. Of these, data prior to May 1999 from 39 square-degrees spread over 64 fields have been analysed. In this period, two thirds of the fields were imaged about 210 times in average; the remaining third were imaged only about 110 times. The exposure times range from 3 min in the lmc center to 12 min on the periphery; the average sampling is once every 4 days (resp. 8 days).

## 3 LMC data analysis

The analysis of the lmc data set was done using a program independent from that used in the smc study, with largely different selection criteria. The aim is to cross-validate both programs (as was already done in the analysis of eros1 Schmidt photographic plates ) and avoid losing rare microlensing events333 We have checked that the present program finds the same smc candidate as reported in .. The analysis is very similar to that reported in . We only give here a list of the various steps, as well as a short description of the differences with respect to our two-year analysis. A detailed description of the analysis will be provided in .

We first select the 6% “most variable” light curves, a sample much larger than the number of detectable variable stars. This subset of our data is “enriched” in genuine variable stars444 We monitor our selection efficiency with Monte-Carlo simulated variable star and microlensing light curves., but also and mainly in photometrically biased light curves, i.e. those of stars especially sensitive to the observing conditions, such as stars very close to nebulosities or to bright stars. Working from this “enriched” subset, we apply a first set of cuts to select, in each colour separately, the light curves that exhibit significant variations. We first identify the baseline flux in the light curve - basically the most probable flux. We then search for runs along the light curve, i.e. groups of consecutive measurements that are all on the same side of the baseline flux. We select light curves that either have an abnormally low number of runs over the whole light curve, or show one long run (at least 5 valid measurements) that is very unlikely to be a statistical fluctuation. We then discard light curves with a low signal-to-noise ratio by requiring that the group of 5 most luminous consecutive measurements be significantly further from the baseline than the average spread of the measurements. We also check that the measurements inside the most significant run show a smooth time variation.

The second set of cuts compares the measurements with the best fit point-lens point-source constant speed microlensing light curve (hereafter “simple microlensing”). They allow us to reject variable stars whose light curves differ too much from simple microlensing, and are sufficiently loose not to reject light curves affected by blending, parallax or the finite size of the source, and most cases of multiple lenses or sources. We also require that the fitted time of maximum magnification lie within the observing period or very close to it, and that the fitted timescale is shorter than 300 days. The latter cut is equivalent to requiring that the baseline flux of the star is observed for at least a few months; this is necessary in any analysis using this baseline flux. At this stage of the analysis, all cuts have been applied independently in the two passbands.

After this second set of cuts, stars selected separately in the two passbands represent about 0.01% of the initial sample; almost all of them are found in two thinly populated zones of the colour-magnitude diagram. The third set of cuts deals with this physical background. The first zone contains stars brighter and much redder than those of the red clump; variable stars in this zone are rejected if they vary by less than a factor two or have a very poor fit to simple microlensing. The second zone is the top of the main sequence, where the selected stars, known as blue bumpers , display variations that are almost always smaller than 60% of the base flux or at least 20% lower in the visible passband than in the red one. These cannot correspond to simple microlensing, which is achromatic; neither can they correspond to microlensing plus blending with another unmagnified star, as it would imply blending by even bluer stars, which is very unlikely. We thus reject all candidates from the second zone exhibiting one of these two features (see Fig. 1).

Compared to the analysis in , two new cuts are introduced to reject other types of variable stars that were not present in the two-year analysis. The first one is aimed at stars which have a roughly constant luminosity for some time, then vary typically over one or two months to reach a new constant level. We cannot yet conclude whether these are physical variable stars or some kind of instrumental problem. The second cut is aimed at novæ and supernovæ. It rejects light curves which have a rise time significantly smaller than the decline time; it is not applied to events with a timescale longer than 60 days, in order not to reject microlensing phenomena with parallax effects, that also show an asymmetry.

The final cuts are simply tighter cuts on the fit quality, applied to both colours (whereas similar previous cuts were applied independently in each passband), and a requirement that the observed magnification be at least 1.40 .

The tuning of each cut and the calculation of the microlensing detection efficiency are done with simulated simple microlensing light curves, as described in . For the efficiency calculation, microlensing parameters are drawn as follows : time of maximum magnification $`t_0`$ uniformly within the observing period $`\pm 150`$ days, impact parameter normalised to the Einstein radius $`u_0[0,2]`$ uniformly, and timescale $`t_E[1,400]`$ days uniformly in $`\mathrm{ln}(t_E)`$. All cuts on the data were also applied to the simulated light curves.

Only four candidates remain after all cuts. Of the two candidates presented in , eros2-lmc-3 is still a member of this list, while eros2-lmc-4 was seen to vary at least twice in the third season and was thus rejected. There are three new candidates, numbered 5 to 7. Their light curves are shown in Figs. 2 and 3; microlensing fit parameters are given in Table 1. Although the candidates pass all cuts, agreement with simple microlensing is not excellent. In particular, eros2-lmc-5 is dubious : it has a bad fit to simple microlensing and is located in an atypical region of the colour-magnitude diagram. The geometric mean of the candidates timescales is about 32 days, including that of the eros1 candidate lmc-1.

The microlensing detection efficiency of this analysis, normalised to events with an impact parameter $`u_0<1`$ and to an observing period $`T_{\mathrm{obs}}`$ of three years, is summarised in Table 2. The main source of systematic error is the uncertainty in the influence of blending. Blending lowers the observed magnifications and timescales. While this decreases the efficiency for a given star, the effective number of monitored stars is increased so that there is partial compensation. This effect was studied with synthetic images using measured magnitude distributions . Our final efficiency is within 10% of the “naive” sampling efficiency. Compared to the efficiency in , the present one is improved for the longest and shortest durations, but slightly lower for average durations around 50 days. This is largely explained by the fact that we have included in the present analysis stars in external lmc fields that were sampled less frequently.

## 4 Limits on Galactic halo MACHOs

eros has observed six microlensing candidates towards the Magellanic Clouds, one from eros1 and four from eros2 towards the lmc, and one towards the smc. As discussed in , and further in , we consider that the long duration of the smc candidate together with the absence of any detectable parallax, in our data as well as in that of the macho group , indicates that it is most likely due to a lens in the smc555 Alternatively, it can be argued that, if due to a galactic halo lens, this event corresponds to a heavy lens . . For that reason, the limit derived below uses the five lmc candidates. (The limit corresponding to all six candidates would be about a factor 1.13 times the limit shown, for masses larger than 0.01 solar mass.)

The limits on the contribution of dark compact objects to the Galactic halo are obtained by comparing the number and durations of microlensing candidates with those expected from Galactic halo models. We use here the so-called “standard” halo model described in  as model 1, but have checked that we obtain similar results for other reasonable halo models. The model predictions are computed for each eros data set in turn, taking into account the corresponding detection efficiencies ( and Table 2 above), and the four predictions are then summed. In this model, all dark objects have the same mass $`M`$; we have computed the model predictions for many trial masses $`M`$ in turn, in the range \[$`10^8\mathrm{M}_{}`$, $`10^2\mathrm{M}_{}`$\].

The method used to compute the limit is as in . We consider two ranges of timescale $`t_E`$, within or outside the interval I $`=[7.5;190]`$ days. This interval was chosen as follows. We first determine the average mass corresponding to the mean duration of the five lmc candidates, at about $`0.2\mathrm{M}_{}`$. We then compute the expected distribution of microlensing timescales for this average mass and check that the observed spread in timescales for the candidates is compatible with the width of this distribution. This means that our candidates are compatible with the hypothesis that their spread in mass contributes very little to the width of the timescale distribution. The interval I is then chosen as a symmetrical interval in $`\mathrm{ln}(t_E)`$ that contains 99% of the timescale distribution for halo machos of $`0.2\mathrm{M}_{}`$. Of course, all five lmc candidates have timescales well within the interval I.

We can then compute, for each mass $`M`$ and any halo fraction $`f`$, the combined Poisson probability for obtaining, in the four different eros data sets taken as a whole, zero candidate outside I and five or less within I. For any value of $`M`$, the limit $`f_{\mathrm{max}}`$ is the value of $`f`$ for which this probability is 5%. Whereas the actual limit depends somewhat on the precise choice of I, the difference (smaller than 5%) is noticeable only for masses around 0.01 and $`10\mathrm{M}_{}`$. Furthermore, we consider our choice for I to be a conservative one.

Figure 4 shows the 95% C.L. exclusion limit derived from this analysis on the halo mass fraction, $`f`$, for any given dark object mass, $`M`$. The solid line corresponds to the five lmc candidates; it is the main result of this article. This limit rules out a standard spherical halo model fully comprised of objects with any mass function inside the range $`[10^710]M_{}`$. In the region of stellar mass objects, where this result improves most on previous ones, the new lmc data contribute about 73% to our total sensitivity (the smc and eros1 lmc data contribute 10% and 17% respectively). The total sensitivity for $`t_E=50`$ days, that is proportional to the sum of $`N_{}T_{\mathrm{obs}}ϵ(t_E=50\mathrm{d})`$ over the four eros data sets, is about 3.2 times larger than that of  and two thirds that of .

## 5 Discussion

After nine years of monitoring the Magellanic Clouds, eros has a meager crop of five microlensing candidates towards the lmc and one towards the smc, whereas about 30 events are expected towards the lmc for a spherical halo fully comprised of $`0.4\mathrm{M}_{}`$ objects. Moreover, some of the candidates cannot be considered excellent. These candidates were obtained from four different data sets analysed by independent, cross-validated programs. So, the small number of observed events is unlikely to be due to bad (and overestimated) detection efficiencies.

This allows us to put strong constraints on the fraction of the halo made of objects in the range \[$`10^7\mathrm{M}_{}`$, $`10\mathrm{M}_{}`$\], excluding in particular at the 95% C.L. that more than 40% of the standard halo be made of objects with up to $`1\mathrm{M}_{}`$. The preferred value quoted in , $`f=0.2`$ and $`0.4\mathrm{M}_{}`$, is consistent with our limit as can be seen in Fig. 4. (The upper part - about 25% - of the domain allowed by  is excluded by the limit we report here.)

There are several differences which should be kept in mind while comparing the two experiments. First, eros uses less crowded fields than macho with the result that blending is relatively unimportant for eros. (Were eros results to be corrected for blending, the detection efficiency would increase slightly and the reported limit would be stronger.) Second, eros covers a larger solid angle (64 deg2 in the lmc and 10 deg2 in the smc) than macho, which monitors primarily 15 deg2 in the central part of the lmc. The eros rate should thus be less contaminated by self-lensing, i.e. microlensing of lmc stars by dimmer lmc objects, which should be more common in the central regions. The importance of self-lensing was first stressed in . Third, the macho data have a more frequent time sampling.

The results from eros and macho are apparently consistent, but the way they are interpreted is different. macho reports a signal and considers the contamination of its sample as low or null. eros2 quotes an upper limit and does not claim its sample to be background-free. The position of the lenses along the line of sight, halo or Magellanic Clouds, is also an issue. macho has compared the spatial distribution of its candidates across the face of the lmc and observes a better agreement with the halo hypothesis than with a specific model of the lmc. On the other hand, because the eros stars are spread over a wider field, the fact that the eros sample corresponds to a lower central value of the event rate (about twice lower than that of macho) is compatible with an interpretation where a noteable fraction of the events are due to self-lensing. The small number of eros candidates precludes at present any definitive conclusion on that topic.

It seems likely that the single most important input to the question of the position of the lenses will come from the comparison of the microlens candidates samples towards the smc with those towards the lmc. Because the two lines of sight are rather close (about 20 degrees apart), the timescale distributions of microlensing candidates towards the two Clouds should be nearly identical if lenses belong to the galactic halo. Also, the event rates should be comparable, although the ratio is more halo model dependent. At present, eros has analysed two seasons of smc data  and macho has not yet presented its detection efficiency towards the smc. From the published eros efficiencies, and assuming that the macho efficiencies towards the smc are similar to those towards the lmc, it can be expected that the completed experiments will have gathered between five and ten microlenses towards the smc. This should allow a significant comparison of the timescales (see also the discussion in ).

Finally, let us mention that, given the scarcity of our candidates and the possibility that some observed microlenses actually lie in the Magellanic Clouds, eros is not willing at present to quote a non zero lower limit on the fraction of the Galactic halo comprised of dark compact objects with masses up to a few solar masses.

Acknowledgements We are grateful to D. Lacroix and the staff at the Observatoire de Haute Provence and to A. Baranne for their help with the MARLY telescope. The support by the technical staff at ESO, La Silla, is essential to our project. We thank J.F. Lecointe for assistance with the online computing. We thank all non-eros members who participated in data taking.


text:
# Comment on “Efficient, multiple-range random walk algorithm to calculate the density of states”

\[

. \]

In their recent Letter Wang and Landau have proposed a very effective MC method of producing evenly distributed histograms for the degeneracy $`g(E_i)`$ of the energy levels of different discrete spin models. The method suffers, however, from a drawback. In the course of their calculations Wang and Landau have to reduce their multiplication factor $`f`$ to a value $`f_{final}`$ which differs from 1 only by $`10^8`$. At this stage the algorithm comes very close to entropic sampling. Indeed when $`g(E_i)`$ has the value $`10^8`$ then the modified density of states (DOS) $`g(E_i)g(E_i)f_{final}`$ in the new algorithm and the modified DOS: $`g(E_i)g(E_i)+1`$ in entropic sampling are identical. Thus the new method has the definite advantage of rapidly building a rough estimator for the DOS but in the end the work of refining it to useful accuracy is as time consuming as entropic sampling.

We have ourselves used a similar method as Ref.  to calculate histograms of the degeneracy of $`g(E_i)`$ of three dimensional Ising models and of the degeneracy $`g(E_0,M_i)`$ of 2d-Ising models for a fixed value $`E_0`$ of the energy as a function of the magnetization $`M_i`$. We have also produced histograms for the three dimensional ANNNI-model. All these data are unpublished.

In our approach the acceptance probability of a spin flip is also governed by the DOS which has been accumulated so far. In this respect our approach is identical to that of Wang and Landau, but it differs from Ref.  in two important aspects: (i) the histogram $`g(E_i)`$ of the DOS is updated (multiplied by $`f`$) only every tenth step, (ii) the transition observables $`T(E_i,\mathrm{\Delta })`$ are also evaluated by recording all attempts to move from an energy $`E_i`$ to $`E_i+\mathrm{\Delta }`$ with $`\mathrm{\Delta }=0,\pm 4,\pm 8`$ for the 2d-Ising model. On taking the logarithm one obtains the entropy $`S(E_i)=\mathrm{ln}g(E_i)`$ from the DOS and the inverse temperature $`\beta =\mathrm{ln}\left[T(E_i+\mathrm{\Delta },\mathrm{\Delta })/T(E_i,\mathrm{\Delta })\right]/\mathrm{\Delta }`$ from the transition observables. It turns out that the latter data suffer much less from statistical errors than the former. In contrast to Ref.  we do not vary $`f`$ – it is kept constant during the simulation. In the examples which follow we have added $`\eta =0.004`$ viz. $`\eta ^{}=0.002`$ to the entropy $`S(E_i)`$ after every 10 attempted spin flips. This leads to multiplication factors $`f=e^\eta `$ of 1.004 viz. 1.002.

It turns out that the error one makes is roughly proportional to $`\eta `$, i.e., $`\mathrm{ln}g_\eta (E_i)S(E_i)+\eta f(E_i)+K_\eta `$ where $`\eta f(E_i)`$ is this error and $`K_\eta `$ is an unimportant constant. A very good approximation for $`S(E_i)`$ can therefore be determined from two MC runs with different values of $`\eta `$ and $`\eta ^{}`$: $`S(E_i)\left[\eta \mathrm{ln}g_\eta ^{}(E_i)\eta ^{}\mathrm{ln}g_\eta (E_i)\right]/\left(\eta \eta ^{}\right)+K`$. Fig. 1 has been produced for the $`32\times 32`$ 2d-Ising model. It shows the second derivative $`d^2s(\epsilon )/d\epsilon ^2`$ of our simulated data. Here $`s=S/N`$ and $`\epsilon =E/N`$ are the entropy and the energy per spin. These data scatter around the exact data of Beale. The difference to these data is shown in the inset. Gaussian smoothing over 20 channels reduces the maximum deviation to 0.005.

As the entropy is determined up to an unimportant constant, the derivatives of the entropy are more assertive. Therefore one should not compare the entropy itself, but its derivatives with the exact results. We have chosen the second derivative for this comparison because it enters into the microcanonically defined specific heat $`c(\epsilon )=s^2(\epsilon )/s^{\prime \prime }(\epsilon )`$ and it is thus very important that high quality data are obtained for it.

In summary: Because of our substraction scheme we can keep the entropy increment and concomitantly the effectivity of the simulation procedure five orders of magnitude above Ref. .

We thank Michael Promberger for his advice.


text:
# Very Low Mass Stars and Brown Dwarfs of the Young Open Cluster IC23911footnote 11footnote 1Based on observations obtained at the Cerro Tololo Inter-American Observatory.

## 1 Introduction

IC2391 is one of the youngest and nearest open clusters. The Hipparcos distance modulus for the cluster is (m-M)o$``$5.82. As a compromise between various determinations of the distance modulus and Hipparcos parallaxes (Becker & Fenkart 1971; Lyngå 1987; van Leeuwen 1999; Robichon et al. 1999), we adopt (m-M)0=5.95$`\pm `$0.1, corresponding to a distance of 155 pc. The mean reddening towards the cluster is small, E(B-V)=0.06 (Patten & Simon 1996). Prior to our lithium depletion age determination, the most commonly quoted age estimate for IC2391 was $``$35 Myr (Mermilliod 1981). In Barrado y Navascués, Stauffer & Patten (1999) we derived a new age estimate of $`53\pm 5`$ Myr.

Because of its proximity and youth, IC2391 has been the target of a number of recent studies. Rotational velocities, lithium abundances and H$`\alpha `$ data for G, K and early M dwarf members can be found in Stauffer et al. (1989, 1997). X-ray fluxes, rotational periods and low resolution spectra for another set of probable cluster members are provided by Patten & Simon (1986) and Simon & Patten (1998). All these data support the conclusion that IC2391 is younger than the Alpha Per cluster ($``$85 Myr in the new lithium age scale, Stauffer et al. 1999; Barrado y Navascués, Stauffer, & Bouvier 1999; Stauffer & Barrado y Navascués 2000) but significantly older than 10 Myr (e.g., because the cluster contains no classical TTauri stars).

For an age of 35 Myr, we expected the lithium depletion boundary for IC2391 to be at IC$``$15.2 (Baraffe 1998, private communication; Stauffer, Schultz, & Kirkpatrick 1998). The ROSAT images of the cluster only supplied candidate cluster members down to IC$``$15. No appropriate deep photographic plate sequences are available for a proper motion survey to the desired depth. We therefore obtained deep, multicolor photometric imaging of the cluster in order to produce a list of candidate cluster members. Because of the low galactic latitude of the cluster (bII=-6.90), our derived candidate list may be significantly contaminated by field stars. We therefore defer extensive discussion of these candidates until detailed spectroscopic follow-up is available.

In Section 2, we provide the details of our imaging program. The method we used to identify candidate cluster members is outlined in Section 3.

## 2 Observations and Data Reduction

### 2.1 The optical photometry and the initial selection of candidates

Our optical CCD photometric data (V,RC,IC,Z) were collected at the Cerro Tololo Inter-American observatory during four different observing runs: January 22, 1998, at the 0.9m (RIC filters); January 24–25, 1998, at the Blanco 4m telescope (RIC filters); April 4, 1998, at the 1.5m (RIC filters); and January 4–6, 1999, again with the 0.9m telescope (VICZ filters). In the case of the CTIO 0.9m runs, we used the Tek 2048 #3 camera (0.43 arcsec/pix), yielding a field of view $``$14 arcmin on a side. For the CTIO 1.5m telescope, we used the CCD SITe $`2048\times 2048`$ imager at the f/7.5 focus, yielding an image scale of 0.44”/pix and a field of view of $`14.8\times 14.8`$ arcmin. For the Blanco telescope campaign, we used the Big Throughput Camera (BTC), a mosaic detector composed of four different 20482 pixel CCDs. Each CCD covers an area of 14.7$`\times `$14.7 arcmin, and the mosaic has a cross-shaped gap, 5.4 arcmin wide, between the CCDs. The total projected area on the sky of the BTC Mosaic is roughly 0.25 sq. deg., with a scale of 0.43 arcsec/pixel. The January 1998 and April 1998 runs were used to make sure we would be able to calibrate the BTC observations and to cover the gaps between the CCDs of the BTC mosaic. In total, we have covered an area close to 2.5 sq. deg. (2 sq. deg. with the BTC). Table 1 lists the coordinates of the center of each field, as well as the exposure times. Figure 1 shows the location of the observed area, together with our very low mass (VLM) stars and brown dwarf (BD) candidates (cross symbols).

The CCD images were bias-subtracted and flat-fielded using standard data reduction techniques and tools within IRAF222IRAF is distributed by the National Optical Astronomy Observatories, which is operated by the Association of Universities for Research in Astronomy, Inc., under contract with the National Science Foundation.. The APPHOT package was used to extract instrumental magnitudes for the objects of interest in each CCD field.

### 2.2 The calibration of the data

All our 0.9m and 1.5m data were obtained under photometric conditions. Unfortunately, this was not the case for the CTIO 4m run (BTC data). Therefore, we used the first two sets of data from the 0.9m and 1.5m telescopes to calibrate the BTC data, corresponding to the (RI)C filters. This calibration included several steps. First, we extracted the instrumental magnitudes for all the CCDs using small apertures ($``$2 pixels, equivalent to $``$0.9 arcsec). Then, for each image, we derived an aperture correction ($``$0.12 mag). Using this method, we minimized the photometric errors, since we did not include flux from the nearby sky. Then, for the 0.9m and 1.5m data, we corrected for the airmass. We assumed the standard CTIO extinction coefficient (0.08 mag/airmass and 0.06 mag/airmass for the RC and IC filters, respectively). Standard stars from Landolt (1992) were observed in all cases, including the fields SA98, SA104, PG1047 and PG1323. They were used to calibrate the data into the Cousins system. Finally, the calibrated photometry from the 0.9m and 1.5m runs was compared with the instrumental magnitudes of the BTC campaign. This provided a zero point shift for each filter from which we computed calibrated magnitudes for the BTC survey. We did not find any significant color term between the three cameras/filter systems.

Table 1 lists the limiting magnitudes for each field in the RC and IC filters. The internal errors should be better than 0.15 magnitudes at the limiting magnitudes; and the photon statistical errors are $``$0.01 mag for RC$``$17 and IC$``$16 for the 0.9m and 1.5m data, and 0.01 for RC$``$19 and IC$``$17.5 for the BTC data. We have estimated the completeness limits from the histogram Log Ndetections versus magnitude. These limits are defined by the magnitude where the histogram deviates from a straight line (see Wainscoat et al. 1992 and Santiago et al. 1996). In the case of the 0.9m and 1.5m data, they are RC$``$18.5 and IC$``$18.0 mag. For the BTC data, we have reached RC$``$20.5 and IC$``$20 mag.

In order to check the external accuracy of our optical data, we have carried out different comparisons. As stated before, we calibrated the BTC data using 2 different campaigns (0.9m and 1.5m telescopes) which took place immediately before and after that run. The comparison between these two runs gives dispersions of $`\sigma `$(R)=0.08, $`\sigma `$(I)=0.07 and $`\sigma `$(R-I)C=0.01. The calibration of the BTC photometry has dispersions of 0.03 and 0.04 magnitudes for the RC and IC filters, respectively (rms scatter of the difference between the raw and calibrated magnitudes of the standard stars). Finally, the optical photometry for overlapping fields of the BTC survey agrees within 0.06 and 0.04 magnitudes (one sigma) for the RC and IC filters, after the calibration.

For the last observing run (January 1999, CTIO 0.9m), we used the VICZ filters and the photometry was calibrated independently. The V-band instrumental magnitudes were corrected for the effects of atmospheric extinction and were placed on the standard system (Johnson-Kron-Cousins) using observations of photometric standard stars (Landolt 1992) and previously determined (R-I)C colors for our program objects. Because the IC and Z CCD images were taken in consecutive pairs (i.e., at the same airmass) and because at that time there were no well-established standard stars for Z filter photometry in the literature, we did not perform an absolute calibration of the IC and Z data. Instead, an (I-Z)CTIO color index was calculated using just the $`I`$ and $`Z`$ instrumental magnitudes. As in Zapatero Osorio et al. (1999), we set IC=Z for those standard stars observed with (R-I)$``$0 in order to determine a zero-point correction for our (I-Z)CTIO color index.

The internal errors of our V-band photometry are estimated to range from $``$0.02 mag for V$``$17 to $``$0.15 mag for V$``$21. The external errors on the V magnitudes may be larger since our typical program stars are much redder than any of the standard stars used to determine the transformation to the standard system. The uncertainties in the (I-Z)CTIO color are estimated from the errors in the IC and Z magnitudes as supplied by the PHOT routine in IRAF. They are computed as $`\mathrm{\Delta }`$(Z-I)2=$`\mathrm{\Delta }`$Z2$`+`$$`\mathrm{\Delta }`$I2. These should be on the order of 0.02 mag or less for the majority of the objects in our IC 2391 sample (for IC$``$17 mag).

### 2.3 Infrared photometry

The infrared data were obtained by the 2MASS survey (Skrutskie et al. 1997) at the CTIO facility during November 21-25 and November 28-29, 1998 and processed at IPAC in January, 1999. All 2MASS scans used for this analysis were observed under photometric conditions. Based on photometric error and total sources detected with magnitude we estimate that these 2MASS data still reach $`0.99`$ completeness and SNR $`10`$ at $`J=15.8`$, $`H=15.1`$, and $`Ks=14.3`$ despite any possible confusion noise in this field.

Both our optical data and the 2MASS survey have very accurate coordinates (see next subsection). Therefore, using the optical coordinates as input, we searched in the 2MASS catalog for a IR counterpart, using a 3 arcsec radii for this search, much larger than the errors in the coordinates.

### 2.4 Coordinates

To define astrometric solutions for the CCD frames we used secondary astrometric standards derived from United Kingdom Schmidt photographic plate material measured using the precision microdensitometer SuperCOSMOS (eg. Hambly et al. 1998). The global astrometric solution for the Schmidt plate was derived using the Tycho–ACT reference catalog (Urban, Corbin & Wycoff 1998), and includes correction for non-linear systematic effects caused by the mechanical deformation of the plates during exposure (eg. Irwin et al. 1998). We used the ”short red” survey plate R6843 (epoch 1981.3, field number 165) for IC2391. These exposures, taken at low galactic latitudes, are far less crowded than the sky limited survey plates and reach R$``$20 (as opposed to R$`22`$ for the deep survey plates). They are ideal for accurate astrometry of secondary standards as faint as R=20 which overlaps with unsaturated objects on the CCD frames. The RMS residual per ACT star in the global astrometric plate solution was $`0.1`$ arcsec in both coordinates; we estimate that there will be no systematic errors in the CCD astrometric solutions larger than this value. The positions of our IC2391 candidates are listed in Table 2.

## 3 Discussion

### 3.1 Initial selection of candidate members of IC2391

The initial selection of candidate members of IC2391 was carried out using the location of the detected stars on the IC versus (R-I)C color-magnitude diagram. Figure 2 displays the data obtained at the 4m CTIO/BTC telescope, where detections are shown as dots. An empirical Zero Age Main Sequence (ZAMS), based on data from Leggett (1992), is included as a solid line. The ZAMS is plotted for our assumed IC2391 distance modulus and a reddening of E(R-I)=0.007 \[David: is this E(R-I) compatible with the E(B-V) mentioned earlier in the text?\] Completeness and limiting magnitudes are depicted as dashed lines. We selected all the detected objects in a wide band well above the ZAMS. The lower envelope for the strip follows a 50 Myr isochrone (D’Antona & Mazzitelli 1997) and the upper envelope is displaced 0.75 mag brighter to include equal mass binaries. The boundaries of the strip were also adjusted to take into account the photometric errors and errors in the distance, age and the reddening. Hence, this band is wider than the 0.75 mag. Then, all candidates were visually inspected on the original images, in order to avoid the presence of false detections, and to verify the stellar-like shape of the detections and the lack of nearby bright stars or cosmic rays which could modify the photometry of the candidate. Since there is no clear separation between the field stars and the location of the cluster isochrone, we expect a strong contamination by spurious members. Specifically, those candidates located only slightly above the ZAMS should be considered with some caution, since they are well below the 50 Myr isochrone. In total, we have selected 206 candidate members, 94 from the BTC survey and another 112 from the 0.9m and 1.5m data. Figure 3 shows all our candidate members as circles. Solid circles indicate the position of the candidates whose membership has been confirmed spectroscopically, via radial velocity, spectral type and several spectroscopic features (H$`\alpha `$, NaI8200Å). These results are discussed in Barrado y Navascués, Stauffer, & Patten (1999), where we provide a lithium depletion boundary (see Barrado y Navascués, Stauffer & Bouvier 1999) age for the cluster (53$`\pm `$5 Myr). Figure 3 also displays the position of a 50 Myr isochrone (short-dashed line), adapted from D’Antona & Mazzitelli (1997). Based on our spectroscopically confirmed members and on Simon & Patten (1996) data, we have created an empirical IC2391 isochrone, shown as a long-dashed line.

Table 2 lists the positions, optical and infrared photometry, the separation between the optical detection and the IR source and the identification with stars from the literature, for each candidate.

### 3.2 Color-Magnitude and Color-Color Diagrams of IC2391

The merger of our RIC data with the 2MASS photometry, and the additional VICZ photometry collected in January 1999, allows us to create a large database of optical-infrared broadband photometry in 7 different filters. Therefore, we have been able to construct several color-magnitude and color-color diagrams. We have used these diagrams to estimate the membership status of each candidate.

Figure 4a depicts the \[V,(V-IC)\] diagram. The empirical main sequence (MS) for young disk stars by Leggett (1992) is included as a solid line. This MS was built for M0-M9 dwarfs, in the ranges 6.65 $``$ M(IC) $``$ 14.67, 4.77 $``$ M(Kcit) $``$ 10.17, 0.75 $``$ (R-I)C $``$ 2.30, and 0.17 $``$ (H-K) $``$ 0.48. Note that our 2MASS data were taken in Ks. However, there are no important differences between the Ks and Kcit systems (see Persson et al 1998). For comparison purposes, the shift of the photometry corresponding to an interstellar absorption of Av=2 is plotted as an arrow. Proposed members of IC2391 from the literature (Patten and Simon 1996; Patten and Pavlovski 1999) appear as crosses. Figure 4b displays the IC magnitudes against the (IC-Ks) color, whereas Figure 4c shows the Ks magnitudes against the (J-Ks). In all these figures, solid triangles represent initial candidate members whose membership has been rejected based on Figures 4 and 5 (probable non-members), whereas open triangles show the possible non–members. Possible and probable members are represented by open and solid circles, respectively (see next subsection).

Several color-color diagrams can be found in Figure 5. Symbols are as in Figure 4. Panel a and b depict \[(V-RC),(R-I)C\], and \[(IC-Ks),(R-I)C\], respectively. It is clear from examination of these diagrams that our original sample included both stars that are plausible cluster members and objects that are instead very likely to be heavily reddened, background stars.

### 3.3 A final list of photometric candidate members

In order to remove the spurious members present in our candidate member list, we have used several color-color diagrams and color-magnitude diagrams. This selection was carried out in a hierarchical way, stressing the spectroscopic and IR data. The scheme we followed is described below:

1. A fraction of our candidates has intermediate resolution spectroscopy (Barrado y Navascués et al. 1999), which indicates whether they are probable members or non-members.
2. We have removed the objects having a large interstellar reddening (see Figure 4b,c and Figure 5b). They were classified as probable non-members.
3. Color-color diagrams were used to remove additional probable non-members (Figures 5a,b).
4. Objects fainter than IC=17 mag, and without IR data (they are too faint to be detected by the 2MASS survey) were classified as possible members, whereas stars brighter than this value, with no IR data (they should have been detected by 2MASS) are listed as possible non-members.
5. Objects bluer than the Leggett’s (1992) IR main-sequence appear in Table 2 as possible members.
6. Objects bluer than our empirical optical isochrone or at the upper (bright) edge of our CM diagram selection strip were also classified as possible members.
7. Finally, the objects that remain are considered to be bona fide members of the cluster (“probable members”).

As a summary, we have classified our initial IC2391 candidate members in four different categories:

* Probable members. Objects located in all CM and CC diagrams with positions which indicate membership. They are identified with the flag “MEM” in the last column of Table 2 (50 objects, including 16 having spectroscopy).
* Possible members (82 objects). They appear identified with the flag “MEM?” in Table 2.
* Possible non-members (10). (“NM?” flag).
* Probable non-members (64). All of them are flagged with “NM” in Table 2.

If we only consider our probable members, we have detected objects $``$2 magnitudes fainter than the previous surveys (Rolleston & Byrne 1997; Simon & Patten 1998; Patten & Pavlovsky 1999). Our faintest possible IC2391 member has IC=20.9, which is 5 magnitudes fainter than the least massive candidate discovered to date. In total, we list 132 objects as probable or possible members of IC2391.

We have not computed the luminosity function (LF) and mass function (MF). Due to its low galactic latitude of the cluster (bII=-6.90), even when using the probable/possible member list, pollution by spurious members is likely to be significant. Visual inspection of the I,(R-I)C color-magnitude diagram (Figure 3) reveals two important characteristics: the relative high number of candidate members with IC in the range 19.5–20.9, and the apparent gap just before this clustering occurs. only additional photometry or spectroscopy would allow us to establish if these faint candidates are real. However, we expect a strong contamination by field stars in this range. If most of these objects are, indeed, not members, the gap would be illusory and the low number of candidate members below IC=18 could be a consequence of our completeness limit (dotted line in Figure 3). The question is still open, until a follow-up spectroscopic study is carried out.

### 3.4 The contamination by field stars

In a deep optical survey of the Pleiades, Bouvier et al. (1998) estimated that the contamination due to field stars was 25%. This value has been confirmed by subsequent spectroscopic follow-up of the Pleiades candidates. Since IC2391 is closer to the Galactic Plane than the Pleiades (bII=-6.90 and bII=-23.52, respectively), the contamination should be much stronger. In fact, Figure 2 of Bouvier et al. (1998) shows a clear discontinuity between field stars and the Pleiades population. This is not the case for IC2391, whose color-magnitude digram (Figure 3) depicts a smooth transition between the field and the locus of the main sequence of the cluster, indicating that contamination by spurious members should be stronger than in the case of the Pleiades. To estimate the degree of contamination by field stars, we have constructed histograms of the number of stars per magnitude bin at various (R-I)C color intervals. In Figure 6 we show the histogram for the (R-I)C=1.9-2.0 range, representative of the distribution of stars down to approximately the completeness limit of the BTC data. The locus of the cluster is indicated. Similar diagrams for adjoining color ranges provide analogous results. The last bin before the cluster contains 3 stars, whereas the brightest bin, well above the main sequence for equal mass binary cluster stars (0.75 brighter than the single star MS), has 5 stars. The bins in between have 17 stars. Assuming an average contamination of 4 stars per bin, the pollution rate would be $``$50% for the initial list of IC2391 candidates (the RIC survey). Our optical and infrared color-color and color-magnitude diagrams have allowed us to remove about 33% of the candidates from the initial 206 objects. Therefore, we estimate that objects cataloged as members in Table 2 (probable and possible members) still have a $``$25% or greater probability of being spurious.

The faintest objects, lacking VZJHK data, may be more strongly polluted by field stars.

### 3.5 Brown Dwarfs in the Cluster

Based on the position of the lithium depletion boundary in members of IC2391, Barrado y Navascués, Stauffer, & Patten (1999) have recently determined an age for the cluster of 53$`\pm `$5 Myr. For this age, the border between the stellar and substellar regimes would be at M(IC)=11.06 (Baraffe 2000, private communication). Taking into account the IC2391 distance and reddening, this leads to IC=17.03. Therefore, all our IC2391 candidates fainter than that value should be brown dwarfs if they indeed belong to the cluster. Our list of final candidates contains ten objects which have been cataloged as probable and possible members based on both optical and infrared data and which have IC in the range 17.06–17.62. Using the Baraffe (2000) models, their mass range is 0.070–0.055 M. Actually, two of these objects (CTIO-061 and CTIO-113) have been observed using low S/N, intermediate resolution spectroscopy (Barrado y Navascués, Stauffer, & Patten 1999). Several spectral characteristics, such as H$`\alpha `$ and NaI8200Å equivalent widths, spectral type, rough radial velocities, indicate that they are real members of the cluster and, therefore, brown dwarfs.

Another 50 objects in Table 2 are fainter than than the IC magnitude of the stellar/substellar limit and lack infrared data. All of them are listed as possible members in Table 2, with magnitudes down to IC=20.9, M(IC)=14.93, and M $``$ 0.025 M.

### 3.6 Comparison with previous surveys and other clusters

A comparison between our final candidate members (circles) and previous surveys (triangles) of IC2391 is shown in Figure 7a. Probable candidate members are displayed with solid circles, whereas open circles represent objects which do not have infrared counterparts (normally, because they are too faint to be detected by 2MASS). The stars from Simon & Patten (1996) are shown as solid triangle, and the objects from Patten & Pavlovsky (1999) appear as open triangles. Clearly, all these samples merge smoothly, describing a good cluster isochrone. Most of our probable candidate members are above the empirical isochrone, which was obtained based on confirmed members of the cluster (Figure 3). In fact, they are located in a wide band, slightly larger than $``$0.75 mag (the maximum shift due to binarity).

We have compared our optical data of IC2391 candidates with data from two other very well known clusters. Figure 7b contains photometry from the Pleiades (asterisks), alpha Per (plus symbols) and IC2391 (solid and open circles for probable and possible members). These open clusters have lithium ages of 125, 85 and 53 Myr, respectively (Stauffer, Schultz, & Kirkpatrick 1998; Stauffer et al. 1999; Barrado y Navascués, Stauffer, & Patten 1999). The large width of the IC2391 main sequence, which includes probable and possible members, can be appreciated. However, this scatter is greatly reduced if we only take into account probable members. Then, the scatter is similar to the value present in the other two clusters. This is consistent with our previous estimate of the contamination for the possible members which are located at the blue side of our empirical IC2391 isochrone. Only further spectroscopy (or, eventually, proper motions) will allow us to verify the membership status of these objects. The comparison between the MS lower envelopes of these two clusters with the lower MS of IC2391 probable members suggests that, indeed, IC2391 is slightly younger than Alpha Per, and that both of them are considerably younger than the Pleiades.

## 4 Summary

Using optical and infrared data, and based on the location on CM and CC diagrams, we have identified a large sample of very low mass stars and brown dwarf candidates of the young cluster IC2391. Accurate coordinates, derived from the SuperCOSMOS microdensitometer, are provided for all of them. We have established a total of 50 probable cluster members, based on their position in the cluster loci in all CC and CM diagrams. Another 82 objects have been cataloged by us as possible members. We have identified two candidate sub-stellar IC2391 members for which we have a full set of multicolor photometry, and whose location in all the CM and CC diagrams supports cluster membership. We have an additional 50 candidate substellar mass members in IC2391, but because we have less or poorer data for these objects we expect that many of them will instead be low mass field stars.

DBN thanks the “Instituto Astrofísico de Canarias” (Spain) and the “Deutsche Forschungsgemeinschaft” (Germany) for their fellowship. JRS acknowledges support from NASA Grant NAGW-2698 and 3690. This work has been partially suported by Spanish “Plan Nacional del Espacio”, under grant ESP98–1339-CO2. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center, funded by the National Aeronautics and Space Administration and the National Science Foundation.

Figure Captions