text:
# Phase ordering and roughening on growing films

## Abstract

We study the interplay between surface roughening and phase separation during the growth of binary films. Already in 1+1 dimensions, we find a variety of different scaling behaviors, depending on how the two phenomena are coupled. In the most interesting case, related to the advection of a passive scalar in a velocity field, nontrivial scaling exponents are obtained in simulations.

PACS numbers: 68.35.Rh, 05.70.Jk, 05.70.Ln, 64.60.Cn

Thin solid films are grown for a variety of technological applications, using molecular beam epitaxy (MBE) or vapor deposition. In order to create materials with specific electronic, optical, or mechanical properties, often more than one type of particle is deposited. When the particle mobility in the bulk is small, surface configurations become frozen in the bulk, leading to anisotropic structures that reflect the growth history, and are different from bulk equilibrium phases. Characterizing structures generated during composite film growth is not only of technological importance, but represents also an interesting and challenging problem in statistical physics.

In this paper, we examine the growth of binary films through vapor deposition, and study some of the rich phenomena resulting from the interplay of phase separation and surface roughening. Simple models for layer by layer growth assume either that the probability that an incoming atom sticks to a given surface site depends on the state of the neighboring sites in the layer below , or that the top layer is fully thermally equilibrated . Assuming that the bulk mobility is zero, once a site is occupied, its state does not change any more. If the growth rules are invariant under the exchange of the two particle types, the phase separation is in the universality class of an equilibrium Ising model. Correlations perpendicular to the growth direction are characterized by the critical exponent $`\nu `$ of the Ising model, and those parallel to the growth direction by the exponent $`\nu z_m`$, with $`z_m`$ being the dynamical critical exponent of the Ising model.

However, the layer by layer growth mode underlying these simple models is unstable, and the growing surface becomes rough. In many cases the fluctuations in the height $`h(𝐱,t)`$, at position $`𝐱`$ and time $`t`$ are self-affine, with correlations

$$\left[h(𝐱,t)h(𝐱^{},t^{})\right]^2|𝐱𝐱^{}|^{2\chi }g\left(t/|𝐱𝐱^{}|^{z_h}\right),$$
(1)
where $`\chi `$ is the roughness exponent, and $`z_h`$ is a dynamical scaling exponent. A computer model with local sticking probabilities that allows for a rough surface was introduced in . In 1+1 dimensions, the authors find phase separation into domains (with sizes consistent with the Ising model), and a very rough surface profile with sharp minima at the domain boundaries. We may ask the following questions: (1) Are the roughness exponents different at the phase transition point? (2) Are the critical exponents modified on a rough surface? We shall demonstrate that the coupling of roughening and phase separation leads to a rich phase diagram, and to nontrivial critical exponents already in 1+1 dimensions.

To characterize phase separation, we introduce an order parameter $`m(𝐱,t)`$, which is the difference in the densities of the two particle types at the surface at position $`𝐱`$ and time $`t`$. The interplay between the fluctuations in $`m`$, and the height $`h`$ is captured phenomenologically by the coupled Langevin equations,

$`_th`$ $`=`$ $`\nu ^2h+{\displaystyle \frac{\lambda }{2}}(h)^2+{\displaystyle \frac{\alpha }{2}}m^2+\zeta _h,`$ (2)
$`_tm`$ $`=`$ $`K(^2m+rmum^3)+ahm+bm^2h`$ (4)
$`+{\displaystyle \frac{c}{2}}m(h)^2+\zeta _m.`$
Here, we have included the lowest order (potentially relevant) terms allowed by the symmetry $`mm`$. Equation (2) is the Kardar-Parisi-Zhang (KPZ) equation for surface growth, plus a coupling to the order parameter. Equation (4) is the time dependent Landau–Ginzburg equation for a (non-conserved) Ising model, with three different couplings to the height fluctuations. The Gaussian, delta-correlated noise terms, $`\zeta _h`$ and $`\zeta _m`$, mimic the effects of faster degrees of freedom. A different set of equations was proposed by Léonard and Desai for phase separation during MBE. Their equations reflect the MBE conditions of random particle deposition (in contrast to sticking probabilities that depend on the local environment), and a conserved order parameter which evolves by surface diffusion. They do not include the KPZ nonlinearity. Computer simulations of corresponding 1+1 dimensional systems are presented in .

Dimensional analysis indicates that the couplings appearing in Eqs. (2-4) are relevant, and may lead to new universality classes. We shall leave the renormalization group analysis of these equations to a more technical paper, and focus here instead on computer simulations in 1+1 dimensions. The quantities evaluated in the computer simulations are the height correlation function in Eq. (1), and the order parameter correlation functions perpendicular and parallel to the growth direction. Allowing for the possibility of different dynamic exponents, $`z_m`$ and $`z_h`$, for the order parameter and the height variables, we fit to the scaling forms

$`G_m^{(x)}(xx^{})`$ $``$ $`m(x,t)m(x^{},t)`$ (5)
$`=`$ $`|xx^{}|^{\eta 1}g_m^{}(|xx^{}|/\xi )`$ (6)
$`G_m^{(t)}(tt^{})`$ $``$ $`m(x,t)m(x,t^{})`$ (7)
$`=`$ $`|tt^{}|^{(\eta 1)/z_m}g_m^{}(|tt^{}|/\xi ^{z_m}).`$ (8)
Our simulations were done using a “brick wall” restricted solid-on-solid model (see Fig. 1). Starting from a flat surface, particles are added such that no overhangs are formed, and with the center of each particle atop the edge of two particles in the layer below. We use two types of particles, $`A`$ and $`B`$ (black and grey in the figures). The probability for adding a particle to a given surface site, and the rule for choosing its color, depend on the local neighborhood. When $`A`$ particles are more likely to be added to $`A`$ dominated regions, and vice versa, the particles tend to phase separate and form domains. In this case, the order parameter correlation length $`\xi `$ is of the order of the average domain width. By changing the growth rules, it is possible to study cases in which some (or all) of the couplings $`a`$, $`b`$, $`c`$, and $`\alpha `$ vanish, and thus to gain a complete picture of the different ways in which the height and the order parameter influence each other.

The decoupled case, $`\alpha =a=b=c=0`$, is implemented using the following updating rules: A surface site is chosen at random, and a particle is added if it does not generate overhangs. Its color is then chosen depending on the colors of its two neighbors in the layer below. If both neighbors have the same color, the newly added particle takes this color with probability $`1p`$, and the other color with probability $`p`$ (where $`p`$ is much smaller than 1). If the two neighbors have different colors, the new particle takes either color with probability 1/2. Neighbors within the same layer are not considered.

Since the probability of adding a particle to a given surface site does not depend on its color, the surface grows exactly as with only one particle type, and is characterized by the KPZ exponents $`\chi =1/2`$, and $`z_h=3/2`$. Similarly, the choice of particle color at a given site is not affected by the height profile. The height profile determines only the moment at which a site is added, since the no-overhang condition requires both neighbors in the previous layer to be occupied. If we equate layer number with time, domain walls move to the right or left with probability 1/2 during one time unit, and a pair of new domain walls is created with probability $`p`$. This is identical to the Glauber model for a one-dimensional Ising chain with coupling $`J`$ and at temperature $`T`$, with $`p=\mathrm{exp}(4J/kT)`$. The correlation length $`\xi `$ perpendicular to the growth direction is consequently $`\xi =\mathrm{exp}(2J/kT)=1/\sqrt{p}`$, and the correlation time is $`\tau =\mathrm{exp}(4J/kT)=1/p`$. The dynamical critical exponent for the order parameter is thus $`z_m=2`$. Note that the “time” used for the order parameter (namely layer number) is different from real time, which is for each particle the moment when it is added to the growing surface. However, this difference becomes negligible for sufficiently small $`p`$ since the thickness of the surface over the correlation length, $`\sqrt{\xi }`$, is much smaller than the characteristic time, $`\xi ^2`$, for order parameter fluctuations. Simulations indeed confirm that the order parameter and height evolve completely independently. A typical profile is shown in Fig. 2a; the corresponding scaling analysis conforms to expectations, and is not presented here.

The situation $`\alpha >0`$ with $`a=b=c=0`$ can be implemented by updating sites on top of particles of different colors less often by a factor $`r<1`$ compared to sites above particles of the same color. As the order parameter is not affected by the height variable, its dynamics is still the same as that of an Ising model, with $`z_m=2`$. The height profile now has domain boundaries sitting preferentially at its local minima, with mounds forming over domains (see Fig. 2b). This leads to a surface roughness exponent of $`\chi =1`$ on length scales $`\xi `$, which is the case studied in. At these scales, changes in the height profile are slaved to domain wall motion, and the dynamic exponent is $`z_h=2`$. However, on length scales much larger than $`\xi `$, the KPZ exponents of $`\chi =1/2`$ and $`z_h=3/2`$ are regained. The crossover in the roughness can be described by a scaling form

$$[h(x,t)h(x^{},t)]^2=|xx^{}|^2g(|xx^{}|/\xi ),$$

with a constant $`g(y)`$ for $`y1`$, and $`g(y)1/y`$ for $`y1`$.

To mimic the influence of surface roughness on the order parameter (nonzero $`a`$, $`b`$, or $`c`$ in Eqs.(4)), the color of a newly added particle is made dependent not only on those of its two neighbors in the layer below, but also on the colors of its two nearest neighbors on the same layer, if these sites are already occupied. With probability $`1p`$, the newly added particle takes the color of the majority of its 2, 3, or 4 neighbors, and with probability $`p`$ it assumes the opposite color. If there is a tie, the color is chosen at random with equal probability. The height variable now affects the order parameter in two ways: (1) Domain walls are driven downhill. The reason is that the neighbor on the hillside of a site being updated is more likely to be occupied than the one on the valley side. The newly added particle is thus more likely to have the color on the hillside. (This corresponds to $`a>0`$ in Eq. (4).) (2) New domains are predominantly formed on hilltops. This is because domains on hilltops can expand more easily than those on slopes or in valleys, indicating $`b>0`$ in Eq. (4). Another consequence is that for the same $`p`$, the correlation length $`\xi `$ is much larger than in the decoupled case, as is apparent in Figs.2c,d.

For the fully coupled case depicted in Fig.2c we find essentially the same scaling behavior as in Fig.2b, i.e. a height profile slaved to the Glauber dynamics of the domains. The most interesting case, shown in Fig.2d, is when the height profile is independent of the domains ($`\alpha =0`$), evolving with KPZ dynamics, while the order parameter is influenced by the roughness. The dynamic exponent $`z_m`$ for the order parameter was first obtained by collapsing the correlation functions using Eqs. (8), as shown in Fig.3. These curves imply that $`\eta =1`$, $`\xi p^{0.542}`$, and $`\tau \xi ^{z_m}1/p`$, giving $`z_m1/0.5421.85`$.

The same non-trivial value for $`z_m`$ is obtained by a completely independent measurement of the dynamics of domain coarsening following a quench from a “high temperature” ($`p`$ close to 0.5) to zero temperature ($`p`$=0). Fig. 4 shows the domain density as function of time for a system of size $`L=16384`$. The resulting $`z_m1.85`$, is in agreement with the value from the scaling collapse.

The following simple argument fails to provide the exponent $`z_m1.85`$. Consider a Langevin equation, $`\dot{x}=\eta (t)`$, for the position $`x`$ of a single domain wall at time $`t`$. Since the motion of the domain wall is strongly influenced by the height profile, the noise $`\eta (t)`$ must have long-range correlations $`\eta (t)\eta (t^{})=D|tt^{}|^\alpha ,`$ reflecting the dynamics of surface. This choice leads to $`z_m=2`$ for $`\alpha >1`$, and $`z_m=2/(2\alpha )`$ for $`\alpha <1`$. For a colored noise dominated by the slope fluctuations, $`\alpha =2/3`$ and $`z_m=3/2`$, i.e. the height imposes its characteristic time scale on the order parameter. This would presumably be the case if the domain walls were uniformly distributed along the surface. However, due to their tendency to move downhill, they are preferentially found near valleys. A different scaling of the slope fluctuations in the valleys may be the reason for the nontrivial value of $`z_m`$. Indeed, for short times, before the domain walls have moved to their preferred positions, the exponent $`3/2`$ is seen.

The dynamics of domain walls on a growing KPZ surface bears some resemblance to the advection of a passive scalar in a turbulent velocity field, which is characterized by nontrivial exponents and multiscaling . If we neglect interactions between domain walls, and treat them as independent “dust particles” floating on the KPZ surface, the Langevin equation for the particle density $`\rho `$ is

$$_t\rho =K^2\rho +a(h\rho +\rho ^2h)+\zeta _\rho .$$
(9)
The second term describes the advection of particles along a velocity field $`\stackrel{}{v}=h`$. Indeed this transformation maps the KPZ equation into the Burgers equation for a vorticity-free, compressible fluid flow . Equation (9) is a special case of Eq. (4) for $`m`$, with $`r=u=c=0`$, $`b=a`$, and with a conserved noise $`\zeta _\rho `$. (Together with Eq. (2) for the height profile, it is also a special case of the equations used to describe the dynamic relaxation of drifting polymers.) In the remainder, we give the results of computer simulations for this case. The rules for the motion of “dust particles” are identical to those for domain walls. However, each particle is treated as if the others were not present. This means in particular that there is no creation or annihilation of particles.

Fig. 5 shows the mean square displacement of a single “dust particle” in a system of size $`L=4096`$. To obtain good statistics, we averaged over 512 independent and noninteracting particles, and used more than 40 runs. The best fit is obtained for $`z_\rho 1.74`$, distinct from the previous $`z_m1.85`$, implying that the exponents depend on whether or not the domain walls (or “dust particles”) are conserved. In contrast to the advection of a passive scalar in a turbulent velocity field, we find no sign of multiscaling. Fig. 6 shows the positions of 1024 independent “dust particles” in a system of length $`L=512`$. While there is some correlation between minima of the surface profile and wall positions, there are also clusters of particles at higher elevations, indicating that particle diffusion is not sufficiently fast to fully adjust the density to the faster changing height profile. A fit of the density-density correlation function to $`\rho (x)\rho (0)1/x^{2(1\chi _\rho )}`$, gives an exponent $`\chi _\rho 0.85`$.

In summary, the interplay between surface roughening and phase separation leads to a variety of novel critical scaling behaviors. At one extreme, the height profile adapts to the dynamics of critical domain ordering. At the other, the dynamics of domain wall motion is influenced by the roughness, exhibiting new and nontrivial scaling behaviors.

This work was supported by EPSRC (grant No. GR/K79307, for BD), and the National Science Foundation (Grant No. DMR-98-05833, for MK).


format: text

text:
# Relativistic Effects in the Pulse Profile of the 2.5 msec X-Ray Pulsar SAX J1808.4-3658

## 1 Introduction

Strong X-ray pulsations with a 2.49 msec period were discovered from SAX J1808.4-3658 in an April 1998 observation with the Rossi X-ray Timing Explorer (Wijnands & van der Klis 1998). The pulsar is a member of an accreting binary system with an orbital period of 2.01 hour and a low-mass companion (Chakrabarty & Morgan 1998). Though similar to other low-mass X-ray binaries in its timing and spectral properties (e.g. Wijnands & van der Klis 1999; Heindl & Smith 1998), SAX J1808.4-3658 is unique for its X-ray pulsations. No other such binary has shown coherent pulsations in its persistent flux despite of careful searches (Vaughan et al. 1994b, and references therein).

As such, SAX J1808.4-3658 is the fastest rotating accreting neutron star. If the pulsations are due to modulated emission from one hot spot on the neutron star surface, the 2.49 msec period corresponds to an equatorial speed of approximately 0.1c. With these high speeds, SAX J1808.4-3658 offers an excellent system for studying relativistic effects.

One such effect may be the observed lag of low-energy photons relative to high-energy photons in the pulse discovered by Cui et al. (1998). Cui et al. (1998) suggest that the lags are due to Comptonization in a relatively cool surrounding medium. Alternatively, the lags may be the result of a relativistic effect: the high-energy photons are preferentially emitted at earlier phases due to Doppler boosting along the line of sight. This possibility was suggested for the similar lags in the 549 Hz oscillations in an X-ray burst of Aql X-1 (Ford 1999), where a simple model showed that the delays roughly match those expected.

In the following, we present new measurements of the pulsed emission from SAX J1808.4-3658. We show that the energy-dependent phase lags are equivalent to a hardening pulse profile. We model this behavior in terms of a hot spot on the neutron star, including relativistic effects.

## 2 Observations & Analysis

We have used publicly available data from the proportional counter array (PCA) on board RXTE in an ‘event’ mode with high time resolution (122$`\mu `$sec) and high energy resolution (64 channels). The observations occurred from April 10 1998 to May 7 1998, when the source was in outburst.

We generate folded lightcurves in each PCA channel. This is accomplished with the fasebin tool in FTOOLS version 4.2, which applies all known XTE clock corrections and corrects photon arrival times to the solar system barycenter using the JPL DE-200 ephemeris, yielding a timing accuracy of several $`\mu `$sec (much less than the phase binning used here). As a check, we have applied this method to Crab pulsar data and the results are identical to Pravdo et al. (1997). To produce pulse profiles in the neutron star rest frame, we use the SAX J1808.4-3658 orbital ephemeris found by Chakrabarty & Morgan (1998). An example folded lightcurve is shown in Figure 1 (top) for the observation of April 18 1998 14:05:40 to April 19 1998 00:51:44 UTC.

To study the energy spectra at each phase bin, we take the rates at pulse minimum and subtract it from the rest of the data at other phases. This effectively accomplishes background subtraction and eliminates the unpulsed emission which we do not wish to consider. Note, however, that the pulsed emission may have some contribution even at the pulse minimum and this subtraction scheme represents only a best approximation to the true pulse emission. We generate detector response matrices appropriate to the observation date and data mode using pcarsp v2.38, and use XSPEC v.10.0 to fit model spectra.

We fit the spectrum here with a simple powerlaw function. Though the function itself is not meant to be a physical description, the powerlaw index provides a good measure of the spectral hardness. Fits in each phase bin have reduced $`\chi ^2`$ of 0.7 to 1.6. Including an interstellar absorption does not substantially affect these results. The powerlaw index clearly increases through the pulse phase (Figure 1, bottom), i.e. spectrum evolves from hard to soft.

We also fit the profiles of the folded lightcurves in each channel using Fourier functions at the fundamental frequency and its harmonics. From these fits we determine the phase lag in each channel relative to the fits in some baseline channel range. Results for the 18 April 1998 observation are shown in Figure 2 as solid symbols. Note that negative numbers indicate that high-energy photons precede low-energy photons. We are also able to measure lags in the first harmonic, and find that they are opposite in sign to the fundamental, i.e. low-energy photons precede high-energy photons in the first harmonic. No lags are measurable in the other harmonics.

Another way to measure energy-dependent phase lags is by Fourier cross-correlation analysis. This is the method used for SAX J1808.4-3658 by Cui et al. (1998) and for other timing signals as well, e.g. kilohertz QPOs (e.g. Kaaret et al. 1999). For a description of cross correlation analysis see e.g. Vaughan et al. (1994a). From the PCA event mode data we calculate Fourier spectra in various channel ranges with Nyquist frequency of 2048 Hz and resolution of 0.25 Hz. We then calculate cross spectra defined as $`C(j)=X_1^{}(j)X_2(j)`$, where $`X`$ are the complex Fourier coefficients for two energy bands at the pulsar frequency $`\nu _j`$. The phase lag between the two energy bands is given by the argument of $`C`$. We measure all phase delays relative to the unbinned channels range 5 to 8, i.e. 1.83 to 3.27 keV for 5 detector units in PCA gain epoch 3 (April 15 1996 to March 22 1999). The results for the 18 April 1998 observation, are shown in Figure 2 by the open symbols. The phase lags are consistent with those calculated from the lightcurve fitting. From the cross-correlation spectra we are not able to measure lags in the much weaker harmonics.

We have calculated phase lag spectra also for other RXTE observations during the April 1998 outburst. These spectra are similar to that in Figure 2. To quantify the trends, we compute an average phase delay, $`\varphi _{avg}`$, over all energies for each observation. We also fit a broken powerlaw function to each phase delay spectrum: $`\varphi =E^\alpha `$ below a break energy, $`E_b`$, and $`\varphi =\varphi _{max}`$ above this energy. Figure 3 shows the quantities $`\varphi _{avg}`$, $`E_b`$ and $`\varphi _{max}`$ versus the time of each observation.

There is a clear connection between the results of the two analyses presented here. The phase resolved spectroscopy shows that the spectrum softens and correspondingly the peak of the pulse profile appears slightly earlier in phase for higher energies (Figure 1). The method of measuring phase delays shows the same behavior: higher-energy photons preferentially lead lower-energy photons in the fundamental and the magnitude of this phase delay increases with energy (Figure 2). In the following we discuss a model that can account for the phase delays/spectral softening measured here.

## 3 Model

We calculate the expected luminosity as a function of phase in a manner similar to Pechenick et al. (1983); Strohmayer (1992) but including Doppler effects (Chen & Shaham 1989) and time of flight delays (Ftaclas et al. 1986). This treatment is based on a Schwarzschild metric, where the photon trajectories are completely determined by the compactness, $`R/M`$. The predicted luminosity as a function of phase from our code matches the results of Pechenick et al. (1983) and Chen & Shaham (1989) for the various choices of parameters. Braje et al. (2000) recently developed a model for pulse profiles using a slightly different approach.

The parameters in the model are the speed at the equator of the neutron star, $`v`$, the mass of the neutron star, $`M`$, the compactness, $`R/M`$, the angular size of the cap, $`\alpha `$, and the viewing angles, $`\beta `$ (the angle between the rotation axis and the cap center) and $`\gamma `$ (rotation axis to line of sight). Another ingredient is the emission from the spot, which we take as isotropic and isothermal. The spectrum of energy emitted from the spot is another important input. Note that if the emitted spectrum is power law like, the observed spectrum will not evolve with phase, since Doppler transformation preserves the shape (see Chen & Shaham 1989). The intrinsic spectrum must therefore have some shape that transforms to match the observed hardening and phase lags; we use a blackbody spectrum with temperature $`kT_0`$.

A fit from the model is shown in Figure 2. This fit uses the following model parameters: $`R/M=5`$, $`M=1.8\mathrm{M}_{\mathrm{}}`$, $`kT_0=0.6`$ keV, $`v=0.1`$, $`\beta =\gamma =10^\mathrm{o}`$, $`\alpha =10^\mathrm{o}`$. To derive count rates, we use table models in XSPEC and the appropriate response files as discussed above. The fit for this single set of parameters is good; we find $`\chi ^2=5.3`$ from the Fourier cross-correlation data or a reduced $`\chi ^2`$ of 0.8 for all the parameters fixed.

A full exploration of the parameter space of the model is beyond the scope of this letter. We note the following trends, however. The magnitude of the phase lags depends sensitively on $`v`$ and $`kT_0`$. Larger delays result from higher speeds, because the pulses become increasingly asymmetric due to Doppler boosting. This asymmetry is also energy dependent, so there is a dependence on $`kT_0`$ as well (see Chen & Shaham 1989). The lags also depend on $`\beta `$ and $`\gamma `$, especially at higher energies where there is the turn-over noted above. The phase lags are less sensitive to $`R/M`$ and $`\alpha `$. We have tested the assumption of isotropic emission and found that the phase lags depend only weakly on the beaming.

## 4 Discussion

The pulsed emission from SAX J1808.4-3658 evolves through its phase from a relatively hard to a soft spectrum, as shown by our phase resolved spectroscopy (Figure 1). This evolution can also be thought of, and measured as, an energy-dependent phase lag in the fundamental (Figure 2), i.e. higher-energy photons emerging earlier in phase than lower-energy photons.

We have applied a model to the data which consists of a hot spot on the rotating neutron star under a general relativistic treatment. The dominant effect accounting for energy-dependent delays is Doppler boosting, the larger boosting factors at earlier phases giving harder spectra. The model fits the data quite well. The model also provides a stable mechanism for generating the phase delays, which meshes with the fact that the characteristic delays remain stable in time to within 25% (Figure 3) even as there is a factor of two decrease in the X-ray flux, a possible tracer of accretion rate. As noted in Ford (1999), in addition to explaining the energy-dependent phase lags in SAX J1808.4-3658, this mechanism may also account for the lags in the burst oscillations of Aql X-1 (Ford 1999) and kilohertz QPOs (Vaughan et al. 1998; Kaaret et al. 1999). Phase resolved spectroscopy has not yet been possible in these signals.

The model offers a new means of measuring the neutron star mass and radius, a notoriously difficult problem in accreting binary systems. The radius is directly related to $`v`$, the equatorial speed ($`v=\mathrm{\Omega }_{\mathrm{spin}}R`$). The fits also depend on the model parameters $`M`$ and $`R/M`$, though the lag spectra are less sensitive to these quantities. The model could be independently constrained by future optical spectroscopy of the companion which could provide information on $`M`$ and the viewing angles, $`\alpha `$ and $`\beta `$.

We acknowledge stimulating discussions with the participants of the August 1999 Aspen workshop on Relativity where an early version of this work was presented. We thank Michiel van der Klis, Mariano Méndez and Luigi Stella for helpful discussions. This work was supported by NWO Spinoza grant 08-0 to E.P.J.van den Heuvel, by the Netherlands Organization for Scientific Research (NWO) under contract number 614-51-002, and by the Netherlands Researchschool for Astronomy (NOVA). This research has made use of data obtained through the High Energy Astrophysics Science Archive Research Center Online Service, provided by the NASA/Goddard Space Flight Center.


text:
# Introduction

## Introduction

In 1765, Euler proved that several important centers of a triangle are collinear; the line containing these points is named after him. The incenter $`I`$ of a triangle, however, is generally not on this line. Less than twenty years ago, Andrew P. Guinand discovered that although $`I`$ is not necessarily on the Euler line, it is always fairly close to it. Guinand’s theorem  states that for any non-equilateral triangle, the incenter lies inside, and the excenters lie outside, the circle whose diameter joins the centroid to the orthocenter; henceforth the orthocentroidal circle. Furthermore, Guinand constructed a family of locus curves for $`I`$ which cover the interior of this circle twice, showing that there are no other restrictions on the positions of the incenter with respect to the Euler line.

Here we show that the Fermat point also lies inside the orthocentroidal circle; this suggests that the neighborhood of the Euler line may harbor more secrets than was previously known. We also construct a simpler family of curves for $`I`$, covering the interior once only except for the nine-point center $`N`$, which corresponds to the limiting case of an equilateral triangle.

Triangle geometry is often believed to be exhausted, although both Davis  and Oldknow  have expressed the hope that the use of computers may revive it. New results do appear occasionally, such as Eppstein’s recent construction of two new triangle centers . This article establishes some relations among special points of the triangle, which were indeed found by using computer software.

## The locus of the incenter

To any plane triangle $`ABC`$ there are associated several special points, called *centers*. A few of these, in the standard notation , are: the centroid $`G`$, where the medians intersect; the orthocenter $`H`$, where the altitudes meet; the centers of the inscribed and circumscribed circles, called the incenter $`I`$ and the circumcenter $`O`$; and the nine-point center $`N`$, half-way between $`O`$ and $`H`$. The radii of the circumcircle and the incircle are called $`R`$ and $`r`$, respectively.

If equilateral triangles $`BPC`$, $`AQC`$ and $`ARB`$ are constructed externally on the sides of the triangle $`ABC`$, then the lines $`AP`$, $`BQ`$ and $`CR`$ are concurrent and meet at the Fermat point $`T`$. This point minimizes the distance $`TA+TB+TC`$ for triangles whose largest angle is $`120^{}`$ .

The points $`O`$, $`G`$, $`N`$, and $`H`$ lie (in that order) on a line called the Euler line, and $`OG:GN:NH=2:1:3`$. They are distinct unless $`ABC`$ is equilateral. The circle whose diameter is $`GH`$ is called the *orthocentroidal circle*.

Guinand noticed that Euler’s relation $`OI^2=R(R2r)`$ \[1, p. 85\] and Feuerbach’s theorem $`IN=\frac{1}{2}Rr`$ \[1, p. 105\] together imply that

$$OI^24IN^2=R(R2r)(R2r)^2=2r(R2r)=\frac{2r}{R}OI^2>0.$$

Therefore, $`OI>2IN`$. The locus of points $`P`$ for which $`OP=2PN`$ is a circle of Apollonius; since $`OG=2GN`$ and $`OH=2HN`$, this is the orthocentroidal circle. The inequality $`OI>2IN`$ shows that $`I`$ lies in the interior of the circle .

Guinand also showed that the angle cosines $`\mathrm{cos}A`$, $`\mathrm{cos}B`$, $`\mathrm{cos}C`$ of the triangle satisfy the following cubic equation:

$$\rho ^4(12x)^3+8\rho ^2\sigma ^2x(32x)16\sigma ^4x4\sigma ^2\kappa ^2(1x)=0,$$
(1)
where $`OI=\rho `$, $`IN=\sigma `$ and $`OH=\kappa `$. We exploit this relationship below.

The relation $`OI>2IN`$ can be observed on a computer with the software *The Geometer’s Sketchpad*$`^\text{®}`$, that allows tracking of relative positions of objects as *one* of them is moved around the screen. Let us fix the Euler line by using a Cartesian coordinate system with $`O`$ at the origin and $`H`$ at $`(3,0)`$. Consequently, $`G=(1,0)`$ and $`N=(1.5,0)`$. To construct a triangle with this Euler line, we first describe the circumcircle $`(O,R)`$, centered at $`O`$ with radius $`R>1`$ —in order that $`G`$ lie in the interior— and choose a point $`A`$ on this circle. If $`AA^{}`$ is the median passing through $`A`$, we can determine $`A^{}`$ from the relation $`AG:GA^{}=2:1`$; then $`BC`$ is the chord of the circumcircle that is bisected perpendicularly by the ray $`OA^{}`$.

It is not always possible to construct $`ABC`$ given a fixed Euler line and a circumradius $`R`$. If $`1<R<3`$ then there is an arc on $`(O,R)`$ on which an arbitrary point $`A`$ yields an $`A^{}`$ outside the circumcircle, which is absurd. If $`U`$ and $`V`$ are the intersections of $`(O,R)`$ with the orthocentroidal circle, and if $`UY`$ and $`VZ`$ are the chords of $`(O,R)`$ passing through $`G`$, then $`A`$ cannot lie on the arc $`ZY`$ of $`(O,R)`$, opposite the orthocentroidal circle (Figure 1). $`OYG`$ and $`NUG`$ are similar triangles. Indeed, $`YO=UO=2UN`$ and $`OG=2GN`$, so $`YG:GU=2:1`$; in the same way, $`ZG:GV=2:1`$. Hence, if $`A=Y`$ or $`A=Z`$, then $`A^{}=U`$ or $`A^{}=V`$ respectively, and $`ABC`$ degenerates into the chord $`UY`$ or $`VZ`$. If $`A`$ lies on the arc $`ZY`$, opposite the orthocentroidal circle, $`A^{}`$ will be outside $`(O,R)`$.

Another formula of Guinand ,

$$OH^2=R^2(18\mathrm{cos}A\mathrm{cos}B\mathrm{cos}C),$$
(2)
shows that the forbidden arc appears if and only if $`\mathrm{cos}A\mathrm{cos}B\mathrm{cos}C<0`$, that is, if and only if $`1<R<3=OH`$. The triangle $`ABC`$ is obtuse-angled in this case.

Once $`A`$ is chosen and the triangle constructed, we can find $`I`$ by drawing the angle bisectors of $`ABC`$ and marking their intersection. *The Geometer’s Sketchpad* will do this automatically, and we can also ask it to draw the locus which $`I`$ traces as we move the point $`A`$ around the circumcircle. The idea of parameterizing the loci with $`R`$, instead of an angle of the triangle (as Guinand did) was inspired by the drawing tools available in *The Geometer’s Sketchpad*. This locus turns out to be a quartic curve, as follows.

Proposition 1. In the coordinate system described above, the incenter of $`ABC`$ is on the curve

$$(x^2+y^2)^2=R^2[(2x3)^2+4y^2].$$
(3)
Proof. From equation (2) we see that, once we fix $`R`$, the product $`\mathrm{cos}A\mathrm{cos}B\mathrm{cos}C`$ is fixed. Using Viète’s formulas, this product is obtained from the constant term of Guinand’s equation (1):

$$\mathrm{cos}A\mathrm{cos}B\mathrm{cos}C=\frac{1}{8}\left(1\frac{4\sigma ^2\kappa ^2}{\rho ^4}\right),$$

so that

$$\rho ^4(18\mathrm{cos}A\mathrm{cos}B\mathrm{cos}C)=4\sigma ^2OH^2.$$
(4)
Now, $`I`$ is a point of intersection of the two circles

$$\rho ^2=x^2+y^2\text{and}\sigma ^2=(x\frac{3}{2})^2+y^2.$$

We get (3) by substituting these and (2) in (4), and dividing through by the common factor $`18\mathrm{cos}A\mathrm{cos}B\mathrm{cos}C`$, which is positive since $`OH`$ (the equilateral case is excluded).

Figure 2 is a *Mathematica* plot of the curves (3) with the orthocentroidal circle.

Every point on a locus inside the orthocentroidal circle is the incenter of a triangle. When $`R>3`$, the locus is a lobe entirely inside the circle; if the point $`A`$ travels around the circumcircle once, then $`I`$ travels around the lobe three times, since $`A`$ will pass through the three vertices of each triangle with circumradius $`R`$. When $`1<R<3`$, the interior portion of (3) isshaped like a bell (Figure 2). Let $`A`$ travel along the allowable arc from $`Y`$ to $`Z`$, passing through $`V`$ and $`U`$; then $`I`$ travels along the bell from $`U`$ to $`V`$, back from $`V`$ to $`U`$, and then from $`U`$ to $`V`$ again. While $`A`$ moves from $`V`$ to $`U`$, the orientation of the triangle $`ABC`$ is reversed. When $`R=3`$, the locus closes at $`H`$ and one vertex $`B`$ or $`C`$ also coincides with $`H`$; the triangle is right-angled. If $`A`$ moves once around the circumcircle, starting and ending at $`H`$, $`I`$ travels twice around the lobe.

Proposition 2. The curves (3), for different values of $`R`$, do not cut each other inside the orthocentroidal circle, and they fill the interior except for the point $`N`$.

Proof. Let $`(a,b)`$ be inside the orthocentroidal circle, that is,

$$(a2)^2+b^2<1.$$
(5)
If $`(a,b)`$ also lies on one of the curves (3), then

$$R=\sqrt{\frac{(a^2+b^2)^2}{(2a3)^2+4b^2}}.$$

There is only one positive value of $`R`$ and thus *at most* one curve of the type (3) on which $`(a,b)`$ can lie. Now we show $`(a,b)`$ lies on *at least* one curve of type (3); to do that, we need to show that given (5), $`R>1`$. We need only prove

$$(a^2+b^2)^2>(2a3)^2+4b^2.$$
(6)
Indeed, $`(2a3)^2+4b^2=0`$ only if $`(a,b)=(\frac{3}{2},0)`$; this point is $`N`$. It cannot lie on a locus (3), in fact it corresponds to the limiting case of an equilateral triangle as $`R\mathrm{}`$.

The inequality (5) can be restated as

$$a^2+b^2<4a3,$$
(7)
and (6) as

$$(a^2+b^2)^24a^2+3(4a3)4b^2>0.$$
(8)
From (7) it follows that

$`(a^2+b^2)^24a^2+3(4a3)4b^2`$ $`>`$ $`(a^2+b^2)^24a^2+3(a^2+b^2)4b^2`$
$`=`$ $`(a^2+b^2)(a^2+b^21).`$
But $`a^2+b^2>1`$ since $`(a,b)`$ is inside the orthocentroidal circle; therefore, (6) is true.

## The whereabouts of the Fermat point

The same set-up, a variable triangle with fixed Euler line and circumcircle, allows us to examine the loci of other triangle centers. Further experimentation with *The Geometer’s Sketchpad* suggests that the Fermat point $`T`$ also lies inside the orthocentroidal circle in all cases.

Theorem 1. The Fermat point of any non-equilateral triangle lies inside the orthocentroidal circle.

Proof. Given a triangle $`ABC`$ with the largest angle at $`A`$, we can take a coordinate system with $`BC`$ as the $`x`$-axis and $`A`$ on the $`y`$-axis. Then $`A=(0,a)`$, $`B=(b,0)`$, $`C=(c,0)`$ where $`a`$, $`b`$ and $`c`$ are all positive. Let $`BPC`$ and $`AQC`$ be the equilateral triangles constructed externally over $`BC`$ and $`AC`$ respectively (Figure 4). Then $`P=(\frac{1}{2}(cb),\frac{1}{2}\sqrt{3}(b+c))`$ and $`Q=(\frac{1}{2}(\sqrt{3}a+c),\frac{1}{2}(a+\sqrt{3}c))`$.

The coordinates of $`T`$ can be found by writing down the equations for the lines $`AP`$ and $`BQ`$ and solving them simultaneously. After a little work, we get $`T=({\displaystyle \frac{u}{d}},{\displaystyle \frac{v}{d}})`$, where

$`u`$ $`=`$ $`(\sqrt{3}bc\sqrt{3}a^2acab)(bc),`$
$`v`$ $`=`$ $`(a^2+\sqrt{3}ab+\sqrt{3}ac+3bc)(b+c),`$
$`d`$ $`=`$ $`2\sqrt{3}(a^2+b^2+c^2)+6ac+6ab+2\sqrt{3}bc.`$ (9)
The perpendicular bisectors of $`BC`$ and $`AC`$ intersect at the circumcenter $`O=(\frac{1}{2}(cb),(a^2bc)/2a)`$. The nine-point center $`N`$ is the circumcenter of the triangle whose vertices are the midpoints of the sides of $`ABC`$; we can deduce that $`N=(\frac{1}{4}(cb),(a^2+bc)/4a)`$.

To show that $`T`$ lies inside the orthocentroidal circle, we must show that $`OT>2NT`$, or

$$OT^2>4NT^2.$$
(10)
In coordinates, this inequality takes the equivalent form:

$$\left[\frac{u}{d}\left(\frac{cb}{2}\right)\right]^2+\left[\frac{v}{d}\left(\frac{a^2bc}{2a}\right)\right]^2>4\left[\frac{u}{d}\left(\frac{cb}{4}\right)\right]^2+4\left[\frac{v}{d}\left(\frac{a^2+bc}{4a}\right)\right]^2,$$

or, multiplying by $`(2ad)^2`$,

$$[2auad(cb)]^2+[2avd(a^2bc)]^2>[4auad(cb)]^2+[4avd(a^2+bc)]^2.$$

After expanding and canceling terms, this simplifies to

$$4a^2du(cb)+4adv[(2a^2+2bc)(a^2bc)]4a^2d^2bc>12a^2u^2+12a^2v^2,$$

or better,

$$adu(cb)+dv(a^2+3bc)abcd^23au^23av^2>0.$$
(11)
One way to verify this inequality is to feed the equations (9) into *Mathematica*, which expands and factors the left hand side of (11) as

$$2(b+c)(\sqrt{3}a^2+\sqrt{3}b^2+\sqrt{3}c^2+\sqrt{3}bc+3ab+3ac)(a^4+a^2b^28a^2bc+a^2c^2+9b^2c^2).$$

The first three factors are positive. The fourth factor can be expressed as the sum of two squares,

$$(a^23bc)^2+a^2(bc)^2,$$

and could be zero only if $`a^2=3bc`$ and $`b=c`$, so that $`a=\sqrt{3}b`$. This gives an equilateral triangle with side $`2b`$. Since the equilateral case is excluded, all the factors are positive, which shows that (11) is true, and therefore (10) holds.

Varying the circumradius $`R`$ and the position of the vertex $`A`$ with *The Geometer’s Sketchpad* reveals a striking parallel between the behavior of the loci of $`T`$ and those of $`I`$. It appears that the loci of $`T`$ also foliate the orthocentroidal disc, never cutting each other, in a similar manner to the loci of $`I`$ (Figure 2). The locus of $`T`$ becomes a lobe when $`R=3`$, as is the case with the locus of $`I`$. Furthermore, the loci of $`T`$ close in on the center of the orthocentroidal circle as $`R\mathrm{}`$, just as the loci of $`I`$ close in on $`N`$ (Figure 2).

It is difficult to prove these assertions with the same tools used to characterize the loci of $`I`$, because we lack an equation analogous to (1) involving $`T`$ instead of $`I`$. A quick calculation for non-equilateral isosceles triangles, however, shows that $`T`$ can be anywhere on the segment $`GH`$ except for its midpoint. This is consistent with the observation that the loci of $`T`$ close in on the center of the orthocentroidal circle.

Consider a system of coordinates like those of Figure 4. Let $`b=c`$ so that $`ABC`$ is an isosceles triangle. In this case, by virtue of (9), $`T=(0,b/\sqrt{3})`$. For this choice of coordinates, $`G=(0,a/3)`$ and $`H=(0,b^2/a)`$. $`T`$ lies on the Euler line, which can be parametrized by $`(1t)G+tH`$ for real $`t`$. This requires that

$$(1t)\frac{a}{3}+t\frac{b^2}{a}=\frac{b}{\sqrt{3}},$$

for some real $`t`$. Solving for $`t`$, this becomes

$$t=\frac{a^2\sqrt{3}ab}{a^23b^2}=\frac{a}{a+\sqrt{3}b},$$

unless $`a=\sqrt{3}b`$. This case is excluded since $`ABC`$ is not equilateral. Note that $`t\frac{1}{2}`$ as $`a\sqrt{3}b`$. Thus, $`t`$ takes real values between $`0`$ and $`1`$ except for $`\frac{1}{2}`$, so $`T`$ can be anywhere on the segment $`GH`$ except for its midpoint.

## Acknowledgments

I am grateful to Les Bryant and Mark Villarino for helpful discussions, to Joseph C. Várilly for advice and nical help, and to Julio González Cabillón for good suggestions on the use of computer software. I also thank the referees for helpful comments.


text:
# Not enough stellar Mass Machos in the Galactic Halo Based on observations made at the European Southern Observatory, La Silla, Chile.

## 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 tool to probe the dark matter content of the Galactic halo (Paczyński (1986)). The first microlensing candidates were reported in 1993, towards the lmc (Aubourg et al. (1993); Alcock et al. 1993) and the Galactic Centre (Udalski et al. (1993)) 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 a substantial fraction of the Galactic dark matter resides in planet-sized objects (Aubourg et al. (1995); Alcock et al. (1996); Renault et al. 1997; Renault et al. (1998); Alcock et al. (1998)).

However a few events were detected with longer timescales. From 6-8 candidate events towards the lmc, the macho group estimated an optical depth of order half that required to account for the dynamical mass of the standard spherical dark halo; the typical Einstein radius crossing time of the events, $`t_E`$, implied an average mass of about 0.5 M for the lenses (Alcock et al. 1997a). Based on two candidates, eros1 set an upper limit on the halo mass fraction in objects of similar masses (Ansari et al. 1996), that is below that required to explain the rotation curve of our Galaxy111 Assuming the 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 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 detection of one microlensing event (Palanque-Delabrouille et al. 1998; see also Alcock et al., 1997b). This single event, out of 5.3 million 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 (Afonso et al. 1999). In contrast, an optical detection of a halo white dwarf population was reported (Ibata et al. 1999).

In this letter, we describe the analysis of the two-year light curves from 17.5 million lmc stars. We then combine these results with our previous limits, 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 Bauer et al. (1997) and Palanque-Delabrouille et al. (1998). 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 1998 from 25 square-degrees spread over 43 fields have been analysed. In this period, about 70-110 images of each field were taken, with exposure times ranging from 3 min in the lmc center to 12 min on the periphery; the average sampling is once every 6 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, Ansari et al. (1996)) and avoid losing rare microlensing events. Preliminary results of the present analysis were reported in Lasserre (1999). We only give here a list of the various steps, as well as a short description of the principal new features; details will be provided in Lasserre et al. (2000).

We first select the 8% “most variable” light curves, a sample much larger than the number of detectable variable 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 ask for a minimum 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 curve differs 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.

After this second set of cuts, stars selected in at least one passband 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. Here we find that selected stars, known as blue bumpers (Alcock et al. 1997a ), display variations that are always smaller than 60% and lower in the visible passband than in the red one. These cannot correspond to simple microlensing, which is achromatic; they cannot correspond to microlensing plus blending with another unmagnified star either, as it would imply blending by even bluer stars, which is very unlikely. We thus reject all candidates from the second zone exhibiting these two features.

The fourth set of cuts tests for compatibility between the light curves in both passbands. We retain candidates selected in only one passband if they have no conflicting data in the other passband. For candidates selected independently in the two passbands, we require that their largest variations coincide in time.

The tuning of each cut and the calculation of the microlensing detection efficiency are done with simulated simple microlensing light curves, as described in Palanque-Delabrouille et al. (1998). For the efficiency calculation, microlensing parameters are drawn uniformly in the following intervals: time of maximum magnification $`t_0`$ within the observing period $`\pm 150`$ days, impact parameter normalised to the Einstein radius $`u_0[0,2]`$ and timescale $`t_E[5,300]`$ days. All cuts on the data were also applied to the simulated light curves.

Only two candidates remain after all cuts. Their light curves are shown in Fig. 1; microlensing fit parameters are given in Table 1. Although the candidates pass all cuts, agreement with simple microlensing is not excellent.

The efficiency of the analysis, normalised to events with an impact parameter $`u_0<1`$ and to an observing period $`T_{\mathrm{obs}}`$ of two 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 (Palanque-Delabrouille (1997)). Our final efficiency is within 10% of the naive efficiency.

## 4 EROS1 results revisited

The two eros1 microlensing candidates have been monitored by eros2. The source star in event eros-lmc-2 had been found to be variable (Ansari et al. (1995)), but microlensing fits taking into account the observed periodicity gave a good description of the measurements. Its follow-up by eros2 revealed a new bump in March 1999, eight years after the first one222We thank the macho group for communication about their data on this star.. This new variation, of about a factor two, was not well sampled but is significant. Therefore, eros-lmc-2 is no longer a candidate and we do not include it in the limit computation.

## 5 Limits on Galactic halo MACHOs

eros has observed four microlensing candidates towards the Magellanic Clouds, one from eros1 and two from eros2 towards the lmc, and one towards the smc. As discussed in Palanque-Delabrouille et al. (1998), 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 (Alcock et al. 1997b ), indicates that it is most likely due to a lens in the smc. For that reason, the limit derived below uses the three lmc candidates; for completeness, we also give the limit corresponding to all four candidates.

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 Palanque-Delabrouille et al. (1998) as model 1. The model predictions are computed for each eros data set in turn, taking into account the corresponding detection efficiencies (Ansari et al. (1996); Renault et al. 1998; Afonso et al. 1999; 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 Ansari et al. (1996). We consider two ranges of timescale $`t_E`$, smaller or larger than $`t_E^{\mathrm{lim}}=10`$ days. (All candidates have $`t_E>t_E^{\mathrm{lim}}`$.) 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 with $`t_E<t_E^{\mathrm{lim}}`$ and three or less (alt. four or less) with $`t_E>t_E^{\mathrm{lim}}`$. 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 $`t_E^{\mathrm{lim}}`$, the difference ($`5\%`$) is noticeable only for masses around $`0.02\mathrm{M}_{}`$. Furthermore, we consider 10 days to be a conservative choice.

Figure 2 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 three lmc candidates; it is the main result of this letter. (The dashed line includes the smc candidate in addition.) This limit rules out a standard spherical halo model fully comprised of objects with any mass function inside the range $`[10^74]M_{}`$. In the region of stellar mass objects, where this result improves most on previous ones, the new lmc data contribute about 60% to our total sensitivity (the smc and eros1 lmc data contribute 15% and 25% respectively). The total sensitivity, that is proportional to the sum of $`N_{}T_{\mathrm{obs}}ϵ(t_E)`$ over the four eros data sets, is 2.4 times larger than that of Alcock et al. (1997a). We observe that a large fraction of the domain previously allowed by Alcock et al. (1997a) is excluded by the limit in Fig. 2.

## 6 Discussion and conclusion

After eight years of monitoring the Magellanic Clouds, eros has a meager crop of three microlensing candidates towards the lmc and one towards the smc, whereas 27 events are expected for a spherical halo fully comprised of $`0.5\mathrm{M}_{}`$ objects. These were obtained from four different data sets analysed by four independent, cross-validated programs. So, the small number of observed events is unlikely to be due to bad 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}_{}`$, $`4\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 Alcock et al. (1997a), $`f=0.5`$ and $`0.5\mathrm{M}_{}`$, is incompatible with the limits in Fig. 2 at the 99.7% C.L. (but see the note added below).

What are possible reasons for such a difference? Apart from a potential bias in the detection efficiencies, several differences 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. Second, eros covers a larger solid angle (43 deg2 in the lmc and 10 deg2 in the smc) than macho, which monitors primarily the central 11 deg2 of the lmc. The eros rate should thus be less contaminated by self-lensing that is more common in the central regions - the importance of self-lensing was first stressed by Wu (1994) and Sahu (1994). Third, the macho data have a more frequent time sampling. Finally, while the eros limit uses both Clouds, the macho result is based only on the lmc. For halo lensing, the timescales towards the two Clouds should be nearly identical and the optical depths comparable. In this regard, we remark that the smc event is longer than all lmc candidates from macho and eros.

Finally, given the scarcity of our candidates and the possibility that some observed microlenses actually lie in the Magellanic Clouds, eros is not willing to quote at present 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.

Note added. While the writing of this letter was being finalised, the analysis of 5.7 yrs of lmc observations by the macho group was made public (Alcock et al. (2000)). The new favoured estimate of the halo mass fraction in the form of compact objects, $`f=0.20`$, is 2.5 times lower than that of Alcock et al. (1997a) and is compatible with the limit presented here. None of the conclusions in this article have to be reconsidered. A detailed comparison of our results with those of Alcock et al. (2000) will be available in our forthcoming publication (Lasserre et al. (2000)).

###### 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.


text:
# On The Reddening in X-ray Absorbed Seyfert 1 Galaxies

## 1 Introduction

The presence of absorption edges of O VII and O VIII in the X-ray (Reynolds 1997; George et al. 1998a) indicates that there is a significant amount of intrinsic ionized material along our line-of-sight to the nucleus in a large fraction ($``$ 0.5) of Seyfert 1 galaxies. In addition to highly ionized gas (referred to as an X-ray or “warm” absorber), X-ray spectra often show evidence for a less-ionized absorber. This component has been modeled using neutral gas (cf. George et al. 1998b), and its relationship to the “warm” absorber is unclear. Interestingly, there are several instances in which this additional neutral column is too small by as much as an order of magnitude to explain the reddening of the continuum and emission lines, assuming typical Galactic dust/gas ratios (cf. Shull & van Steenburg 1985). This inconsistency was first noted in regard to the absence of high ionization emission lines in the IUE spectra of MCG -6-30-15 (Reynolds & Fabian 1995). The first quantitative comparison of the neutral columns inferred from the X-ray data to that derived from the reddening was for the QSO IRAS 13349+2438 by Brandt, Fabian, & Pounds (1996), who suggested that the dust exists within the highly ionized X-ray absorber (ergo, a dusty warm absorber). It has been suggested that dusty warm absorbers are present in several other Seyferts (NGC 3227: Komossa & Fink 1997a; NGC 3786: Komossa & Fink 1997b; IRAS 17020+4544: Leighly et al. 1997, Komossa & Bade 1998; MCG -6-30-15: Reynolds et al. 1997). Since it is unlikely that dust could form within the highly ionized gas responsible for the O VII and O VIII absorption, it has been suggested that the dust is evaporated off the putative molecular torus (at $``$ 1 pc) and, subsequently, swept up in an radially outflowing wind (cf. Reynolds 1997).

In this paper, we present an alternative explanation. It is possible that there is a component of dusty gas (which we will refer to as the “lukewarm” absorber), with an ionization state such that hydrogen is nearly completely ionized but the O VII and O VIII columns are negligible, which has a sufficient total column to account for the reddening. Such a possibility has been mentioned by Brandt et al. (1996), while Reynolds et al. (1997) have suggested that the dusty warm absorber in MCG -6-30-15 may have multiple zones. Here we suggest that the lukewarm absorber lies far into the narrow-line region (NLR). Such a component has been detected in the Seyfert galaxy NGC 4151 (Kraemer et al. 1999), and it lies at sufficient radial distance to cover much of the NLR. We will demonstrate that the combination of a dusty lukewarm absorber and a more highly ionized (O VII and O VIII) absorber is consistent with the observed X-ray data and with the reddening of the narrow emission lines in the Seyfert 1 galaxy NGC 3227.

## 2 Absorption and Reddening in NGC 3227

NGC 3227 ($`z`$ $`=`$ 0.003) is a well studied Sb galaxy with an active nucleus, usually classified as a Seyfert 1.5 (Osterbrock & Martel 1993). X-ray observations of NGC 3227 with ROSAT and ASCA reveal the presence of ionized gas along the line-of-sight to the nucleus (Ptak et al. 1994; Reynolds 1997; George et al. 1998b). Using the combined ASCA and ROSAT dataset obtained in 1993, George et al. characterized the absorber with an ionization parameter (number of photons with energies $``$ 1 Ryd per hydrogen atom at the ionized face of the cloud) U $``$ 2.4, and a column density NH $``$ 3 x 1021 cm-2.

The UV and optical emission lines and continuum in NGC 3227 are heavily reddened. Cohen (1983) measured a narrow H$`\alpha `$/H$`\beta `$ ratio $``$ 4.68, and derived a reddening of EB-V $`=`$ 0.51 $`\pm `$ 0.04, assuming the intrinsic decrement to be equal to the Case B value (Osterbrock 1989). Cohen derived a somewhat larger reddening from the \[S II\] lines, EB-V $`=`$ 0.94 $`\pm `$0.23, which may be less reliable due to the weakness of \[S II\] $`\lambda `$4072 (see Wampler 1968). The ratio of broad H$`\alpha `$/H$`\beta `$ $``$ 5.1, indicating that the broad and narrow lines are similarly reddened. Winge et al. (1995) used the total (narrow $`+`$ broad) H$`\alpha `$/H$`\beta `$ ratio to derive a somewhat smaller reddening, EB-V $``$ 0.28. IUE spectra show the UV continuum of NGC 3227 is also heavily reddened (Komossa & Fink 1997a). Based on the Balmer lines, assuming Galactic dust properties and dust/gas ratio, the derived reddening requires a hydrogen (H I and H II combined) column density $``$ 2 x 1021 cm-2 (cf. Shull & Van Steenberg 1985), which is much greater than the estimated neutral column, but similar to that of the ionized gas detected in X-rays in 1993 ($`N_\mathrm{H}`$ $``$ 3 x 1021 cm-2; George et al. 1998b).

Several workers have suggested that NGC 3227 contains a screen of neutral material (in addition to that in the Galaxy) along the line-of-sight to the nucleus. Besides the ionized absorber, both Komossa & Fink (1997a) and George et al. (1998b) found a column density of $``$ 3 x 1020 cm-2 of neutral material, in addition to the Galactic column ($``$ 2.1 x 1020 cm-2, cf. Murphy et al. 1996), is required to model the X-ray spectrum below $``$500 eV. A higher column density ($``$ 6 x 1020 cm-2) has been suggested based on 21-cm VLA observations (Mundell et al. 1997). However the angular resolution of the VLA data is poor (12′′, or $``$ 850 pc) and Mundell et al. did not make a direct detection of H I absorption against the radio continuum source in the inner nucleus. In this paper we argue that based on the current data, there is no reason to include a significant column of completely neutral material. We show that a dusty lukewarm absorber lying outside the NLR is consistent with both the reddening of the optical continuum and narrow lines, and with the attenuation of the X-ray spectrum below $``$500 eV.

## 3 Modeling The Absorber

### 3.1 The Lukewarm Component

The photoionization code we use has been described in previous publications (cf. Kraemer et al. 1994). For the sake of simplicity, we assume that the lukewarm absorber can be represented as a single zone, described by one set of initial conditions (i.e., density, ionization parameter, elemental abundances, and dust fraction). The gas is ionized by the continuum radiation emitted by the central source in the active nucleus of NGC 3227.

In order to fit the SED, we first determined the intrinsic luminosity at the Lyman limit. Since NGC 3227 is heavily reddened, we fit the value at the Lyman limit based on the optical continuum flux. From the average fluxes measured by Winge et al. (1995), after correcting for a reddening of EB-V $`=`$ 0.4 (the average of the reddening quoted by Cohen (1983) and Winge et al.), assuming the reddening curve of Savage & Mathis (1979), we find that the intrinsic flux at 5525 Å is Fλ $``$ 2.5 x 10-14 ergs s-1 cm-2 Å-1. Interestingly, the optical luminosities of NGC 3227 and NGC 4151 are roughly equal and, therefore, we have made the assumption that the two galaxies have similar optical-UV SEDs. Using the same ratio of optical to UV flux for NGC 3227 as in NGC 4151 (Nelson et al. 1999), we determine Fν at the Lyman limit to be $``$ 6.2 x 10-26 ergs s-1 cm-2 Hz-1. The X-ray continuum from 2 – 10 keV can be fit with an index, $`\alpha `$ $``$ 0.6 (George et al. 1998b), from which we derive a flux at 2 keV of $``$ 2.20 x 10-29 ergs s-1 cm-2 Hz-1. Since an extrapolation of the X-ray continuum underpredicts the flux at the Lyman limit by more than two orders of magnitude, the continuum must steepen below 2 keV. Hence, we have modeled the EUV to X-ray SED as a series of power-laws of the form Fν $``$ $`\nu ^\alpha `$, with $`\alpha =1`$ below 13.6 eV, $`\alpha =2`$ over the range 13.6 eV $`h\nu <`$ 500 eV, and $`\alpha =0.6`$ above 500 eV.

Given this simple parameterization, the steepening of the continuum cannot occur at a much lower energy, otherwise the EUV continuum would be too soft222It is possible to have a lower break energy and a sufficient number of He II ionizing photons if the EUV continuum has a significant “Big Blue Bump”, as suggested by Mathews & Ferland (1987). Although assuming such a continuum does not appreciably affect our predictions, a full exploration of parameter space is beyond the scope of this paper. to produce the observed He II $`\lambda `$4686/H$`\beta `$ ratio ($``$ 0.23), specifically since the strong \[O I\] $`\lambda `$6300 line indicates that much of the NLR gas is optically thick (for the relative narrow emission-line strengths, see Cohen 1983). The luminosity in ionizing photons, from 13.6 – 10,000 eV, is $``$ 1.5 x 1053 photons s-1.

We have assumed roughly solar element abundances (cf. Grevesse & Anders 1989), which are, by number relative to H, as follows: He $`=`$ 0.1, C $`=`$ 3.4 x 10-4, N $`=`$ 1.2 x 10-4, O $`=`$ 6.8 x 10-4, Ne $`=`$ 1.1 x 10-4, Mg $`=`$ 3.3 x 10-5, Si $`=`$ 3.1 x 10-5, S $`=`$ 1.5 x 10-5, and Fe $`=`$ 4.0 x 10-5. We assume that both silicate and carbon dust grains are present in the gas, with a power-law distribution in sizes (see Mathis, Rumpl, & Nordsieck 1977). Thus, we have modified the abundances listed above by depletion of elements from gas-phase onto dust grains, as follows (cf. Snow & Witt 1996): C, 20%; O, 15%; Si, Mg and Fe, 50%.

For our model, we require that 1) the absorber lies outside the majority of the NLR emission, and 2) the column of gas is fixed to obtain the observed reddening. Based on the WFPC2 narrow-band \[O III\] $`\lambda `$5007 imaging (Schmitt & Kinney 1996), we have placed the lukewarm absorber at least 100 pc from the central source, and truncated the model at a hydrogen column density NH $`=`$ 2 x 1021 cm-2. We adjusted the ionization parameter such that the model produced a reasonable match to the absorption in the observed soft X-ray continuum (see below).

### 3.2 Comparison to the X-ray data

To compare our model predictions to the X-ray data, we used the 1993 ASCA (0.6–10 keV) and ROSAT PSPC (0.1–2.5 keV) data described in George et al (1998b), excluding the flare (“t3” in Fig 4 of George et al). Following standard practice, the normalizations of each of the 4 ASCA instruments and of the ROSAT dataset were allowed to vary independently. The 5-7 keV band was also excluded from the analysis due to the presence of the strong, broad Fe emission line. We assumed the continuum described above, except that the spectral index above 500 eV was allowed to vary during the analysis. In addition to the lukewarm absorber (and Galactic absorption), any highly-ionized gas was modeled by a series of edges of fixed energy. This is somewhat problematic since neither the ASCA nor ROSAT instruments have sufficient spectral resolution and/or sensitivity to resolve all the possible edges. We have therefore tested for the edges most likely to be visible in highly ionized gas (e.g., O VII and O VIII).

An acceptable fit to the data was obtained ($`\chi ^2=1195`$ for 1176 degrees of freedom; $`\chi _\nu ^2=1.02`$) with the following parameters for the lukewarm absorber: U $`=`$ 0.13, nH $`=`$ 20 cm-3, and the distance of the cloud from the ionizing source is $``$ 120 pc. The predicted electron temperature at the ionized face of this component is $``$ 18,000K and, therefore, it is thermally stable (cf. Krolik, McKee, & Tarter 1981). The best-fitting value for the spectral index above 500 eV was $`0.58\pm 0.03`$, consistent with our initial assumptions. We find evidence of absorption by several ions, with the following column densities: C V, 6.5 x 1017 cm-2; O VI, 2.9 x 1017 cm-2; O VII, 7.8 x 1017 cm-2; O VIII, 1.0 x 1018 cm-2; and, Ne IX, 4.9 x 1017 cm-2. The O VII and O VIII edges translate to an effective hydrogen column density of highly-ionized gas of $``$2 x 1021 cm-2. The data/model ratios from this fit are shown in Fig. 1. The slight underprediction of the absorption below 300 eV in the ROSAT band is easily rectified by a small ($``$ 20%) increase in the column density of the lukewarm absorber.

The ionic column densities for the lukewarm absorber are listed in Table 1 and, as expected, the column densities of O VII and are too small to make detectable contributions to the X-ray absorption edges ($`\tau _{\mathrm{OVII}}`$ $`<`$ 0.01, as opposed to $``$ 0.19 for the highly-ionized gas). Therefore, this component does not resemble the X-ray absorbers most frequently discussed to date (eg Reynolds 1997; George et al 1998a). On the other hand, the model predicts substantial columns for H I, N V, Si IV, C IV, and Mg II, which would result in strong and, for the most part, saturated, UV resonance absorption lines.

ASCA observed NGC 3227 again in 1995. George et al (1998b) have shown that during this epoch the observed spectrum was significantly different, consistent with a thick, highly-ionized cloud moving into and attenuating $``$85% of the line-of-sight to X-ray source. Under our hypothesis that the lukewarm absorber is located $``$100 pc from the nucleus, hence, we do not expect it to have varied between these two epochs. We have therefore checked and found that indeed the 1995 ASCA dataset are consistent with the soft X-ray attenuation from our lukewarm absorber. (Although the lack of simultaneous ROSAT PSPC data during 1995 prevents a stringent test.)

## 4 Discussion

We have shown that the X-ray spectrum of NGC 3227 is consistent with attenuation by the sum of a highly-ionized absorber and a lukewarm absorber. We suggest that these are physically different components of the circumnuclear material surrounding NGC 3227. The characteristics of the highly-ionized absorber are similar to those previously suggested for NGC 3227 (George et al. 1998b); this is in the range of and generally similar to those in other Seyfert 1s (Reynolds 1997; George et al. 1998a). Such absorbers have been observed to vary on timescales $``$3 yr (and much faster in some cases), and are probably due to gas well within the NLR.

The main result of this paper is that the second component, our “lukewarm absorber”, has the appropriate physical conditions to simultaneously explain the absorption seen below 0.5 keV in the X-ray band (previously modeled as completely neutral gas) and the reddening seen in the optical/UV. Agreement with the soft X-ray data is the result of the lukewarm gas containing significant opacity due to He II. Agreement with the reddening of the narrow emission lines places the component outside the NLR.

Although the lukewarm absorber has the appropriate physical conditions and radial distance to redden the NLR, it must also have a sufficiently high covering fraction to be detected. For example, Reynolds (1997) found that 4/20 of radio-quiet active galaxies show both intrinsic X-ray absorption and reddening. Thus, the global covering factor of the dusty ionized absorber must be 20%, within the solid angle that we see these objects (cf., Antonucci 1993). Kraemer et al. (1999) have shown that the covering factor for optically thin gas in NGC 4151, similar to our lukewarm model, can be quite large ($``$ 30%). In addition, Crenshaw et al. (1999) find that $``$ 60% of Seyfert 1 galaxies have UV absorbers with a global covering factor $``$ 50%, and an ionization parameter similar to the lukewarm absorber, but with lower columns on average (cf. Crenshaw & Kraemer 1999). Therefore, it is entirely plausible that there would be optically thin NLR gas along our line-of-sight to the nucleus in a fraction of Seyfert 1s.

The lukewarm model predicts a column of Mg II of 3.3 x 1014 cm-2, which would produce strong Mg II $`\lambda `$2800 absorption. It is interesting that NGC 3227 is one of the few Seyfert 1s to show evidence of Mg II $`\lambda `$2800 in absorption (Ulrich 1988). While Kriss (1998) has shown that Mg II absorption can arise in clouds characterized by small column density and low ionization parameter (NH $``$ 1019.5 cm-2, U $``$ 10-2.5), our results predict that it may also arise in a large column of highly ionized NLR gas, even if a substantial fraction of Mg is depleted onto dust grains.

The lukewarm model predicts average grain temperatures of 30K – 60K, for grains with radii from 0.25 $`\mu `$m – 0.005 $`\mu `$m, respectively. The reradiated IR continuum, which is produced primarily by the silicate grains (cf. Mezger, Mathis, & Panagia 1982), peaks near 60 $`\mu `$m. Assuming a covering factor of unity, this component only accounts for $``$ 1% of the observed IR flux from NGC 3227, (which is $``$ 7.98 Jy at 60 $`\mu `$m; The IRAS Point Source Catalog ). It is likely that most of the thermal IR emission in NGC 3227 arises in the dense (nH $``$ 103 cm-3) NLR gas in which the narrow emission lines are formed, as is the case for the Seyfert 2 galaxy, Mrk 3 (Kraemer & Harrington 1986).

## 5 Summary

Using the combined 1993 ROSAT/ASCA dataset for NGC 3227, and photoionization model predictions, we have demonstrated that the observed reddening may occur in dusty, photoionized gas which is in a much lower ionization state than X-ray absorbers detected (by their O VII and O VIII edges) in Seyfert 1 galaxies (cf. Reynolds 1997; George et al. 1998a) and the dusty warm absorbers that have been proposed (cf. Komossa & Fink 1997a). This component (the lukewarm absorber) is $``$ 120 pc from the central ionizing source and its physical conditions are similar to those in optically thin gas present in the NLR of the Seyfert 1 galaxy NGC 4151. If this model is correct, we predict that strong UV resonance absorption lines with high column densities from the lukewarm absorber will be observed in NGC 3227.

We have confirmed earlier results regarding the presence of an X-ray absorber within NGC 3227 with O VII and O VIII optical depths similar to those determined by Reynolds (1997). This component lies closer to the central source than the lukewarm absorber, but is essentially transparent to EUV and soft X-ray radiation and, hence, does not effectively screen the lukewarm gas. We find no requirement for neutral gas in addition to the Galactic column.

These results illustrate that a moderately large ($``$ 1021 cm-2) column of ionized gas can produce significant soft X-ray absorption if much of the helium is in the singly ionized state. Since clouds with large He II columns may be a common feature of the NLR of Seyfert galaxies, such a component should be included in modeling the X-ray absorption. If such a component is present along our line-of-sight in an active galaxy, it is likely that the intrinsic neutral column has been overestimated.

S.B.K. and D.M.C. acknowledge support from NASA grant NAG5-4103. T.J.T. acknowledges support from UMBC and NASA/LTSA grant NAG5-7835.