text:
# An Optical and Near Infrared Search for Brown Dwarfs in the Pleiades Cluster

## 1 Introduction

Many brown dwarf (BD) searches have targeted young open clusters or star forming regions rather than general fields, because here BDs are still relatively bright and warm during their early phase. The Pleiades region is one of the best clusters for this purpose; it is a fairly young, nearby, rich and compact cluster. It has an age of 125 Myr (Stauffer, Schultz, & Kirkpatrick, 1998) at a distance of 130 pc (Crawford & Perry, 1976; Pinfield et al., 2000). Twelve hundred member stars are located within 2.5 degrees of the cluster centre (Pinfield, Jameson, & Hodgkin, 1998).

A number of optical surveys searching for BDs in the Pleiades cluster have been done to date and identified a numerous population of bona-fide BDs and BD candidates; “CFHT-PL” objects found by an $`RI`$ survey (Bouvier et al., 1998); “NPL” objects by an $`RIJK`$ survey (Festin, 1998); “MHObd” objects by a $`VI`$ survey (Stauffer et al., 1998b); “IPMBD” objects by an $`RI`$ survey (Hambly et al., 1999); “Roque” objects by an $`IZ`$ survey (Zapatero Osorio et al., 1999); “BPL” objects by an $`IZ`$ survey (Pinfield et al., 2000); “int-pl-IZ” objects by an $`IZ`$ survey (Dobbie et al., 2002b). About 20 of these objects are confirmed as Pleiades brown dwarfs on the basis of observed lithium abundance, radial velocity and proper motion (e.g., Rebolo, Zapatero-Osorio, & Martin, 1995; Rebolo et al., 1996; Basri, Marcy, & Graham, 1996; Zapatero Osorio et al., 1997a; Stauffer, Schultz, & Kirkpatrick, 1998; Martín et al., 2000). However, to date, only one L-type BD candidate (Roque25; Martín et al., 1998) has been spectroscopically identified in the Pleiades cluster. Detection of L-type BDs is difficult at optical wavelengths because they should emit most of their bolometric flux in the infrared region. We have carried out a brown dwarf survey in both the optical and near infrared of a small region of the central region of the Pleiades cluster or southern portion of Merope. Most existing surveys have ignored this region because of its proximity to a small molecular cloud which results in large and variable extinction (Breger, 1987).

## 2 Observations

We surveyed an area of $`14\mathrm{}\times 23\mathrm{}`$ centred at $`\mathrm{R}.\mathrm{A}.3^\mathrm{h}46^\mathrm{m}59^\mathrm{s}.6,\mathrm{Dec}.+23\mathrm{°}44\mathrm{}59\mathrm{}(\mathrm{J2000}.0)`$ both in the optical and near infrared bands. The observed area is shown as the shaded rectangle in Fig. 1. The four open rectangles show the area covered at optical wavelengths.

### 2.1 Near Infrared Survey

Our near infrared survey was made on 2000 October 12, 17 and 18 with SIRIUS (Simultaneous InfraRed Imager for Unbiased Survey) on the University of Hawaii 2.2 m telescope, atop Mauna Kea. SIRIUS employs three 1024 $`\times `$ 1024 pixel HgCdTe arrays which provides $`J`$ (1.25$`\mu `$m), $`H`$ (1.65$`\mu `$m), and $`Ks`$ (2.15$`\mu `$m) band images simultaneously (Nagashima et al., 1999; Nagayama et al., 2002). The field of view in the each band is $`4\stackrel{}{.}9\times 4\stackrel{}{.}9`$ with a pixel scale of $`0\stackrel{}{.}28`$ at the Cassegrain focus of f/10. We mapped an area of $`14\mathrm{}\times 23\mathrm{}`$ with 3 $`\times `$ 5 tiles of the field of view. We obtained each of the tiles by dithering 18 frames. The exposure time was 60 sec for each frame. The limiting magnitudes at 10 $`\sigma `$ are 20.0, 19.2, and 18.5 at $`J`$, $`H`$, and $`Ks`$, respectively. Typical seeing in the $`Ks`$ band was $`1\stackrel{}{.}0`$.

Dark and dome-flat frames were taken at the beginning or the end of each night. We observed standard stars, Nos. 9105, 9107, 9108, 9116 and 9188, in the faint NIR standard star catalog of Persson et al. (1998) on the same nights for photometric calibration, and the red standard stars, LDN 547, BRI B0021-0214 and BRI B2202-1119, in Table 3 of Persson et al. (1998) on 2001 August 31 for determination of the colour conversion formula from ours to the CIT system. We obtained the following equations as colour-conversion from ours to the CIT system.

$`J_{\mathrm{CIT}}=J_{\mathrm{obs}}(0.013\pm 0.017)\times (JKs)_{\mathrm{obs}}+(0.012\pm 0.030)`$
$`H_{\mathrm{CIT}}=H_{\mathrm{obs}}(0.029\pm 0.029)\times (JH)_{\mathrm{obs}}+(0.007\pm 0.033)`$
$`K_{\mathrm{CIT}}=Ks_{\mathrm{obs}}(0.024\pm 0.020)\times (JKs)_{\mathrm{obs}}+(0.018\pm 0.035)`$
We used pipeline software based on NOAO’s IRAF (Imaging Reduction & Analysis Facility)111 IRAF is distributed by the National Optical Astronomy Observatory, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation. package to reduce the data. We have applied standard NIR image reduction procedures, including dark current subtraction, sky subtraction and flat fielding. Source detection, morphological classification and photometry were performed with the APPHOT package in IRAF.

### 2.2 Optical Survey

The $`i^{}`$ (0.77 $`\mu `$m, SDSS system) and $`Z`$ (0.9 $`\mu `$m, RGO) band image was taken on 2000 October 20 with the Wide Field Camera (WFC) on the Isaac Newton Telescope (INT) at La Palma, as part of the INT Wide Angle Survey (McMahon et al., 2001). The WFC accommodates four 2048 $`\times `$ 4196 pixel EEV CCDs. The field of view for each CCD is $`11\stackrel{}{.}4\times 22\stackrel{}{.}8`$ with a pixel scale of $`0\stackrel{}{.}33`$ at the prime focus of the INT and the mosaic layout of the four CCDs enables us to encompass the whole field of $`14\mathrm{}\times 23\mathrm{}`$ taken in the near infrared bands (Fig. 1). We obtained one frame for each filter; the exposure time of each frame was 600 s. The 10 $`\sigma `$ detection limits are 22.0 and 20.5 at $`i^{}`$ and $`Z`$, respectively, where $`i^{}`$ and $`Z`$ are on natural system of WFC.

The observed data were reduced at the Institute of Astronomy, Cambridge, using the WFC data reduction pipeline software (Irwin & Lewis, 2001). All the frames were bias-subtracted and then corrected for linearity before the science fields were flatfielded and defringed. The pipeline was also used to generate photometrically and astrometrically calibrated object catalogues. On the night in question, the $`i^{}`$\- and $`Z`$-band photometric zeropoints were measured from three standard star fields (observed once each) in the catalogue of Landolt (1992) and were found to be stable to about 1 %. The astrometric solution takes account of the radial distortion across the WFC field of view, and resulted in residuals of $`0\stackrel{}{.}22`$.

## 3 Selection of Pleiades BD Candidates

### 3.1 The colour-colour diagram

We constructed the colour-colour ($`JK,i^{}J`$) diagram for all the objects, 846 sources in total, detected with S/N ratios higher than 10 in the near infrared and optical bands, differentiating them as extended (149; open circles) and point (697; filled circles) sources in Fig. 2. Seven sources detected with S/N ratios higher than 10 in the near infrared but lower than 5 in the optical are also plotted with large error bars. “Extended objects” were sorted with a criterion of the FWHMs larger by at least 20 % than the seeing size of the night in the $`Ks`$ images.

The thick solid and dashed curves from the bottom-left corner represents the loci of main sequence and red giant stars (Bessell, 1991, 1979); spectral types of K7 dwarf at the cusp at around (1,0.8) and of M7 at the tip (3.5,1.2) for ($`i^{}J,JK`$), while the arrow at the upper-right corner shows the reddening vector of $`\mathrm{A}_\mathrm{V}`$ = 3 mag (Schultz & Wiemer, 1975; Rieke & Lebofsky, 1985). Most of the point sources (filled circles) are concentrated above the locus from K7 to M7, along the extinction vector up to 3 mag from the main sequence (reddening band), indicating they are K- or M-type dwarfs lying intrinsically on the main sequence. On the other hand, most of the extended sources and some point sources spread upward of the reddening band. About 20 of the extended sources have galaxy-like shapes. The K-correction curves from z = 0 - 0.1 of elliptical and spiral galaxies (Furusawa, private communication) are plotted as dotted lines. Thus, objects spread upward of the reddening band including the point sources (182 objects in total) are likely to be elliptical or spiral type galaxies reddened with K-corrections of z $``$ 0 to 0.1 and $`\mathrm{A}_\mathrm{V}`$ of 0 to several magnitudes. A major problem in this kind of survey is to distinguish stars from distant galaxies at faint magnitudes (e.g. Festin, 1997). The optical and near-infrared colour-colour diagram allows us to easily distinguish them. Note that the near infrared colour-colour diagram ($`JH,HK`$) alone cannot distinguish low-mass stars from these galaxies, due to merging of them in the diagram (e.g. Kirkpatrick et al., 1999; Jarrett et al., 2000).

To highlight the expected location of Pleiades BDs in the colour-colour diagram, two theoretical isochrones for the Pleiades cluster are also plotted in Fig. 2; the dashed line is the 125 Myr isochrone of the NextGen model (Baraffe et al., 1998), while the dot-dashed line is the 120 Myr isochrone of the Dusty model (Chabrier et al., 2000, ; Baraffe, private communication for 120 Myr model). The $`I_\mathrm{C}`$ magnitudes of these models have been transformed on to our system using the following relation, $`i^{}=I_\mathrm{C}+0.211(R_\mathrm{C}I_\mathrm{C})+0.011(VI_\mathrm{C})^20.003(VI_\mathrm{C})`$. This transformation is derived via Landolt system ($`I_{\mathrm{Lan}}`$); the transformation between $`i^{}`$ and $`I_{\mathrm{Lan}}`$ is from an INT web page 222http://www.ast.cam.ac.uk/~wfcsur/photom.php, and that between $`I_{\mathrm{Lan}}`$ and $`I_\mathrm{C}`$ is from Bessell (1990). The NextGen model is based on the nongray dust-free atmosphere (Hauschildt, Allard, & Baron, 1999), and describes successfully various observed properties of early-mid M dwarfs, but predicts near infrared colours too blue compared to the observed ones for cooler objects. The Dusty model includes the effect of grain formation, that is, 1) the photospheric depletion of dust-forming elements, and 2) scattering and absorption by dust, explaining successfully the near infrared colours of the late M and L types with 1,800 $`\mathrm{T}_{\mathrm{eff}}2,200`$ K (e.g. Jameson et al., 2002). Because of the inadequacy of the NextGen model in the cooler BD region, we have adopted the Dusty model in the following discussion. Using near infrared photometry and lithium tests of low mass cluster members, Zapatero Osorio et al. (1997b) estimated the stellar / substellar boundary to be located at an absolute magnitude of ($``$ 0.075 $`\mathrm{M}_{\mathrm{}}`$) for Pleiades ($``$120 Myr) to be at $`\mathrm{M}_\mathrm{I}`$ = 12.4, $`\mathrm{M}_\mathrm{J}`$ = 10.0 and $`\mathrm{M}_\mathrm{K}`$ = 9.0 mag. This corresponds to the point ($`i^{}J,JK`$) = (2.9, 1.0) in our colour-colour diagram, highlighted by a cross, at approximately spectral type M6. Twelve late-type field dwarfs from Dobbie et al. (2002a) were also plotted as asterisks. We adopt the redder side of the line running from the cross in parallel with the extinction vector as the location of Pleiades BDs in Fig. 2. Note that this location coincides with that of late-M and L-type field dwarfs. There remain 12 sources on the redder side of the line.

### 3.2 The colour-magnitude diagram

To separate Pleiades BDs from field dwarfs, we construct the colour-magnitude diagram ($`J,i^{}K`$). The 12 red objects selected by Fig. 2 are plotted in Fig. 3. The theoretical tracks representative of the low mass stellar and substellar members of the Pleiades are also plotted; NextGen model (Baraffe et al., 1998) and Dusty model (Chabrier et al., 2000). We choose the $`J`$ magnitude as the ordinate of the diagram because $`i^{}`$ magnitudes of the red objects have larger errors and $`K`$ magnitudes are less sensitive to the mass difference on the theoretical tracks. The arrow at the upper-right corner indicates the reddening vector of $`\mathrm{A}_\mathrm{V}`$ = 3 mag, which runs nearly in parallel with the tracks.

We adopt the distance modulus used by Dobbie et al. (2002b), $`(mM)_0=5.53`$, to estimate the ordinate displacement ($`J`$ magnitude) of the isochrone. Different age estimates (70 – 150 Myr) result in a displacement of $``$ -0.3 and $``$ +0.1 mag. Different distance estimates result in a displacement of $``$ -0.2 mag (Hipparcos) and $``$ +0.1 mag (photometry). The effect of the cluster depth results in a displacement of $`\pm `$ 0.2 mag. Additionally, location of unresolved binaries is by 0.75 mag above the single-star sequence (e.g. Pinfield et al., 2000). Accounting for all of these uncertainties and allowing for a small degree of error in both of the theoretical models and photometry, BDs belonging to the Pleiades cluster should be located between 0.3 mag below and 1.0 mag above the relevant single-star isochrone.

With these criteria, we identified four Pleiades BD candidates from 12 red objects. The remaining eight objects lie more than 1 mag below the isochrone in the colour-magnitude diagram, although they fall in the Pleiades BD domain in the colour-colour diagram. They could be late-M and L-type dwarfs lying background of Pleiades. Finally, we list the four sources (PL- # of Fig. 3) as BD candidates from our five-colour survey. All of them are previously unpublished. The positions and magnitudes of four BD candidates are provided in Table 1, and that of the remaining eight objects are in Table 2. Finding charts of four BD candidates are presented in Fig. 4.

## 4 Discussion

### 4.1 Estimation of the BD candidate masses

To estimate masses of BD candidates from their observed magnitudes, we should determine the level of interstellar extinction towards them. The distribution of the background field stars in Fig. 2 indicates that the interstellar extinction in the surveyed portion of Pleiades ranges from 0 to 3 mag in $`\mathrm{A}_\mathrm{V}`$. The foreground extinction toward Pleiades is relatively uniform with the value of $`\mathrm{A}_\mathrm{V}`$ 0.12 mag, except for a region south of Merope, which contains a small molecular cloud (Crawford & Perry, 1976). Because the cloud is assigned to be inside the cluster from polarimetry and reddening of the member stars and the field stars (Breger, 1987), the extinctions estimated from background stars should be regarded as upper limit for the BD candidates.

To estimate the maximum interstellar extinction for the individual BD candidates, we constructed an extinction map in the surveyed area from all the field stars (small dots between the two straight lines in Fig. 2), assuming them to be K- or M-type main sequence stars (the thick line in Fig. 2). Visual extinction, $`\mathrm{A}_\mathrm{V}`$, toward each star was calculated from its $`i^{}J`$ and $`JK`$ colours with the reddening law of Schultz & Wiemer (1975) and Rieke & Lebofsky (1985), i.e. $`\mathrm{E}(JK)/\mathrm{E}(i^{}J)=0.56,\mathrm{A}_\mathrm{V}=6.10\times [(JK)_{\mathrm{obs}}(JK)_{\mathrm{intrinsic}}]`$. We meshed the area with 17 $`\times `$ 10 grids of $`1\mathrm{}`$ square and adopted the maximum $`\mathrm{A}_\mathrm{V}`$ within the grid, and averaged them over each mesh with Gaussian weights. The final $`\mathrm{A}_\mathrm{V}`$ map is shown in Fig. 5 and the values toward the BD candidates are listed in Table 3.

The estimated maximum extinctions for individual candidates are shown in Fig. 3 with dotted arrows. If extinction is neglected, PL-2 and -3 are estimated to have M $``$ 0.035 $`\mathrm{M}_{\mathrm{}}`$, similar to Roque 25 (Martín et al., 1998) and int-pl-IZ-69, -81, -84 (Dobbie et al., 2002b), the coolest candidate members so far detected. On the other hand, assuming the maximum likely level of extinction, their masses are estimated to be M $``$ 0.04 $`\mathrm{M}_{\mathrm{}}`$. Similarly we estimate the mass of the faintest candidate to lie in the range of 0.033 – 0.028 $`\mathrm{M}_{\mathrm{}}`$. Therefore, PL-4 maybe the lowest mass member of the Pleiades identified to date. The estimated masses from $`J`$ magnitudes for the case of no extinction and that of maximum possible extinction toward individual candidates are shown in Table 3. Note that uncertainties in the age and the distance of the cluster and the effect of the cluster depth contribute an additional error of 10 – 20 % in the mass estimation; e.g. the mass of PL-4 for the case of no extinction is $`0.028_{0.001}^{+0.002}\mathrm{M}_{\mathrm{}}`$. We note that our mass estimates would increase slightly if we were to take into account the “missing M dwarf gap” (Dobbie et al., 2002c).

### 4.2 Broad-band energy distributions

In Fig. 6 we compare the four SEDs with those of cool field dwarfs, determined in previous surveys (M6 – M8; Kirkpatrick & McCarthy 1994, L0 – L5; Kirkpatrick et al. 1999,2000). Solid and dashed lines represent SEDs of the four BD candidates for the cases of no extinction and maximum extinction discussed in § 4.1, respectively. We have assumed the distance modulus for Pleiades of $`(\mathrm{m}\mathrm{M})_0`$ = 5.53 (Crawford & Perry, 1976; Pinfield et al., 2000) to convert the observed magnitudes to absolute magnitudes for the cluster member. The zero-magnitude fluxes are taken from Bessell, Castelli, & Plez (1998).

All the candidates show a steep decline towards the shorter wavelengths and a broad hump between 1 and 2 $`\mu `$m. Within the uncertainties arising from the correction for interstellar extinction, the four SEDs of the candidate BDs are well matched by the energy distributions of the late-M and L-type field dwarfs: PL-1 is consistent with a mid to late-M spectral type, PL-2 and PL-3 with late-M to early L spectral type and PL-4, an early-mid L spectral type. We summarize these results in Table 3.”

### 4.3 Contamination of field stars

As discussed in section 3.1, the optical and near-infrared colour-colour diagram allows us to separate distant galaxies and giants from BD candidates. However, field dwarfs are still possible sources of contamination.

The Dusty model indicates that Pleiades BDs are overluminous with respect to field stars of similar spectral type by $``$ 1 mag. The criterion of selecting objects between 0.3 mag below and 1.0 mag above the isochrone means that we are sensitive to single and binary dwarfs at distance of 71 – 130 and 96 – 184 pc, respectively. This corresponds to space volumes of 15.5 and 45.1 $`\mathrm{pc}^3`$ for single and binary field stars, respectively. Assuming a space density of M8 – L4.5 dwarfs of 0.0066 $`\mathrm{pc}^3`$ (Gizis et al., 2000) and a binary fraction of 50 % (e.g. Steele & Jameson, 1995), we estimate the contamination of $``$ 0.2 field star in our BD candidates. Thus, the level of the contamination of the BD candidates by non-members is negligibly low.

### 4.4 The cluster luminosity function and mass function

We have derived the cluster luminosity function for $`J`$ magnitudes between 15.5 to 20.0 which correspond to the BD limit of Zapatero Osorio et al. (1997b) and the 90 % completeness limit of our survey (almost the same as 10 $`\sigma `$ detection limit), respectively. In our calculation we assume that the interstellar extinction toward Pleiades in our survey area is spatially uniform with the value of 1.5 mag at $`V`$ or 0.4 mag at $`J`$, we subdivided this magnitude range into 2-mag bins, given in Table 4.

Using the Dusty model, we have calculated the masses corresponding to each of our luminosity bins. We estimated the total number, $`N_{\mathrm{annul}}`$, of member of the whole cluster from the number of our sample. First we calculate the number, $`N_{\mathrm{annul}}`$, of members inside an annulus. Our survey area covers 7.7 % of the area of the annulus of an inner and an outer radius of $`r_i=0.26\mathrm{°}`$ and $`r_o=0.64\mathrm{°}`$, respectively (Fig 1). $`N_{\mathrm{annul}}`$ is calculated simply from the area ratio. Next we calculate $`N_{\mathrm{total}}`$. Pinfield, Jameson, & Hodgkin (1998) and Raboud & Mermilliod (1998) show that the surface density distribution of the Pleiades members can be well fitted by a King distribution (King, 1962) whose core radius increases as the stellar mass decreases. The following equation, which is obtained by integrating the function of King distribution, provides the number of cluster members inside the circle of a radius $`r`$,

$$n(r)=\pi r_c^2k\left[\mathrm{ln}(1+x)4\frac{\sqrt{1+x}1}{\sqrt{1+x_t}}+\frac{x}{1+x_t}\right],$$
(1)
where $`x=(r/r_c)^2`$, $`x_t=(r_t/r_c)^2`$, $`k`$ is the surface density of the cluster, $`r`$ the radius from the cluster centre, $`r_c`$ the core radius at which the surface density of members falls to half its central value, and $`r_t`$ the tidal radius at which the external gravitational potential equals the cluster potential. The total number of cluster members, $`N_{\mathrm{total}}`$, can be determined independently of $`k`$ by taking the ratio $`n(r_t)/[n(r_o)n(r_i)]`$, where $`n(r_o)n(r_i)`$ is equal to $`N_{\mathrm{annul}}`$ and $`n(r_t)`$ is equal to $`N_{\mathrm{total}}`$. We adopt a value of $`r_c=2.2\mathrm{°}`$ from Jameson et al. (2002) for low-mass members of Pleiades, and a value of $`r_t=5.54\mathrm{°}`$ from Pinfield, Jameson, & Hodgkin (1998).

$$\frac{N_{\mathrm{total}}}{N_{\mathrm{annul}}}=\frac{n(r_t)}{[n(r_o)n(r_i)]}=12.78.$$
(2)
Finally, we calculate the values of substellar mass function of the cluster, $`N_{\mathrm{total}}`$ per unit mass. The values are given in Table 4 and are plotted in Fig. 7 with several previous estimates (Dobbie et al., 2002b; Hodgkin & Jameson, 2000; Hambly et al., 1999). The shape of mass function of the Pleiades across and below the stellar / substellar boundary has been reported to be well represented by a power-law with an index $`\alpha `$ lying between 0.4 – 1.0, dependent on investigator (e.g. Jameson et al., 2002; Bouvier et al., 1998; Hodgkin & Jameson, 2000; Zapatero Osorio et al., 1997a). Our new points seem to support the larger values of $`\alpha `$. However, as the present survey covers only a small fraction of the cluster (i.e., $`<0.5`$ %), this result must be treated with a measure of caution. Nevertheless, we note that Béjar et al. (2001) have recently determined the substellar mass function of the much younger $`\sigma `$ Orionis cluster to be well matched by a power-law with index $`\alpha `$ = 0.8.

## 5 Conclusions

We have undertaken an optical ($`i^{},Z`$) and near infrared ($`J,H,Ks`$) survey of a $`14\mathrm{}\times 23\mathrm{}`$ area near the centre of the Pleiades cluster to search for brown dwarf members. We have unearthed 4 BD candidates (PL-1, -2, -3, and -4), out of 853 sources. PL-1 is likely an M-type BD, PL-2 and -3 are likely late-M or early-L type BDs, while PL-4 is probably an L-type BD. The estimated mass of our faintest candidate, PL-4, is in the range of 0.028 – 0.033 $`\mathrm{M}_{\mathrm{}}`$, lower than that of Roque 25 (Martín et al., 1998) and can be lower than that of int-pl-IZ-69 (Dobbie et al., 2002b), possibly making it lowest mass Pleiad identified to date.

## Acknowledgments

We thank the staffs of the UH 2.2 m telescope and the Subaru telescope for supporting the observation of SIRIUS. This work is financially supported by Sumitomo foundation, Grant-in-Aid for Scientific Researches on Priority Areas (A), and International Scientific Research, the Ministry of Education, Culture, Sports, Science, and Technology. The Isaac Newton Telescope is operated by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofisica de Canarias. Nagashima, Nagayama, and Y. Nakajima are financially supported by the Japan Society for the Promotion of Science. PDD is a PPARC postdoctral research associate. We also thank an anonymous referee for his/her useful comments and advice.


format: text

text:
# THE X-RAY HALO OF G21.5–0.9

## INTRODUCTION

In various respects supernova remnant G21.5–0.9 looks like a standard plerion. In fact it shows: i. a filled-center morphology, both in radio and in X rays; ii. a flat radio synchrotron spectrum; iii. a non-thermal X-ray spectrum, with photon index $`2`$. There is no direct evidence of the pulsar powering the nebula (no pulsed emission has been detected as yet), but G21.5–0.9 harbours an almost point-like central X-ray source, that has been interpreted as the termination shock of the pulsar wind (Slane et al., 2000).

However a non-thermal X-ray halo has been recently detected, extending well beyond the boundaries of the radio emitting plerion (Warwick et al., 2001). This is the opposite of what observed in other plerions, as well as of what standard models of plerions would predict. Our main goal is to investigate the origin of this halo.

Moreover, a result of past X-ray spectral analyses of the halo region (Warwick et al., 2001; Safi-Harb et al., 2001) has been that the overall spectrum of the halo does not show any thermal signature. No thermal emission has been even detected in “the North Spur”, a feature in the northern part of the halo, that morphologically resembles a partial shell. In this paper we shall also discuss this point.

## THE NATURE OF THE X-RAY HALO

### Difficulties with the standard interpretation

An X-ray synchrotron nebula with a size larger than the radio emitting nebula represents a challenge to the existing models of plerion. In fact the synchrotron emitting relativistic electrons are expected to be injected in the immediate surroundings of the associated pulsar, and then to be advected and/or diffused to fill the whole nebula; in the meanwhile synchrotron losses are selectively quenching higher energy particles. Therefore X-ray emitting electrons should typically disappear before reaching the boundary of the magnetic bubble, which instead should be rather well outlined at radio wavelengths. The above scenario nicely explains (at least at a qualitative level) the behavior of most plerions, in which the apparent size of the nebula is shrinking with increasing photon energy.

In order to account for the opposite behavior, as observed in G21.5–0.9, one may for instance invoke higher diffusion for higher energy electrons. But then, in order to effectively get synchrotron emission, strong magnetic fields are required: therefore the magnetic bubble should extend to the whole X-ray halo; but in this case it is hard to explain why the radio emitting electrons do not fill the whole magnetic bubble. More detailed modelling may be devised, but the above reasoning should have made clear the theoretical embarrassment rising after the G21.5–0.9 case.

### An alternative explanation: dust scattering

An alternative interpretation is that the plerionic nebula is not intrinsically more extended in X rays than in radio; but that the observed X-ray halo is just an effect of dust scattering in the foreground medium.

Following the empirical relation derived by Predehl and Schmitt (1995), the measured column density $`N_\mathrm{H}=2.3\times 10^{22}\mathrm{cm}^2`$ (as from Warwick et al., 2001) corresponds to $`\tau _{\mathrm{sca},1}=1.07`$ (where $`\tau _{\mathrm{sca},1}`$ indicates the dust scattering optical depth at 1 keV). For such a high value the expected fractional halo intensity is compatible with what measured (38%) in the 0.1–2.4 keV band. However the Predehl and Schmitt statistical relation allows a considerable dispersion, and therefore it is not sufficient by itself to prove that dust scattering is the main cause of the halo observed in G21.5–0.9.

We should add that fits to the halo spectrum imply $`N_\mathrm{H}`$ values significantly lower than in the inner part of the plerion (Warwick et al., 2001; Safi-Harb et al., 2001). This decrease of $`N_\mathrm{H}`$ with increasing radius cannot be intrinsic; it may be regarded, instead, as a sign that dust scattering distorts to some extent the local spectra, biasing the spectral fits. Scattering may partly affect even the observed radial softening in the X-ray halo. The analysis described below is intended to explore in depth the dust scattering hypothesis, by investigating the halo properties in various X-ray energy bands, taking advantage of energy scaling relations intrinsic of the scattering physics.

## MODELLING A DUST SCATTERING HALO

### The data

We have reanalyzed EPIC/PN data obtained in a XMM-Newton Cal/PV observation of G21.5–0.9. In Figure 1 the source is imaged in three different bands, 0.5–2, 2–5 and 5–8 keV. Overlaid to the soft image are the selection areas for the spectral analysis described in the last section: an annulus with inner and outer radius of respectively $`70^{\prime \prime }`$ and $`110^{\prime \prime }`$ has been subdivided into 4 sectors, labelled as halo1 (P.A. from $`30^{}`$ to $`60^{}`$), halo2 (P.A. from $`60^{}`$ to $`150^{}`$), halo3 (P.A. from $`150^{}`$ to $`240^{}`$), halo4 (P.A. from $`240^{}`$ to $`330^{}`$); furthermore, a $`45^{\prime \prime }`$ diameter circle centered on the North Spur has been labelled Nspur.

### The model assumptions

Our analysis is based on some assumptions, which are listed below.

* We assume that the scattering is adequately described by the Rayleigh-Gans Theory. This implies well defined energy dependences, both for the halo scale length ($`E^1`$) and for the scattering optical depth ($`E^2`$; see e.g. Mauche and Gorenstein, 1986). A deviation from these power laws is expected in soft X rays, and is actually observed in our data (see below).
* The dust properties (namely dust density and grain size distribution) are taken to be homogeneous along the line of sight. In the case of variations in the distribution, a sort of average is obtained. An uneven density distribution along the line of sight, instead, does not affect the profile of the halo tail (on which we base our analysis).
* We model the intrinsic X-ray size of the source as only weakly dependent on energy. In fact in other plerions the X-ray size may slightly shrink with increasing energy. Anyway this would not affect the dust halo tail, where the original source can be hardly distiguished from a point-like source.
* We assume EPIC/PN Point Spread Function (PSF) to be energy independent. Looking at the radial profiles of a point-like source (LMC X–3; Figure 2) it can be easily verified to be a valid approximation.

To summarize, all above assumptions look reasonable. Moreover the final conclusions do not strictly require them to be valid. Therefore we believe our results are rather robust.

### Analysis procedure

This procedure is based on the comparison of radial profiles in various X-ray bands. In order to obtain profiles of sufficient quality to detect even the outermost tail, we have limited our analysis to 4 spectral bands, namely 1.3–1.7 keV, 2–3 keV, 3–5 keV, and 5–8 keV.

Figure 2 shows the 4 profiles, for G21.5–0.9 as well as for LMC X–3 (a bright point-like source with low $`N_\mathrm{H}`$, used to determine the instrumental PSF). The profiles are background subtracted and vignetting corrected, and out-of-time events have been taken away. In the case of G21.5–0.9 we have also excluded regions of the image containing the North Spur and the star SS 397. All profiles are shown normalized to their peak value. In the case of G21.5–0.9 the radial profile is getting broader in softer X rays: this is not an instrumental effect, since from Figure 2 it is apparent that the PSF is energy independent till the faintest parts of its wings.

In order to fit the dust scattering halo of G21.5–0.9 we adopt a model with 3 parameters: $`\tau _{\mathrm{sca},1}`$ (optical depth at 1 keV), $`a_{\mathrm{max}}`$ (maximum grain size) and $`q`$ (power-law index of the grain distribution; Predehl and Schmitt, 1995). Taking the 5–8 keV profile G21.5–0.9 as a good guess of the halo-free source profile, a fit to the 2–3 keV source profile (see the relative curve in Figure 3) gives $`\tau _{\mathrm{sca},1}=1.50`$, $`a_{\mathrm{max}}=0.90\mu \mathrm{m}`$, $`q=3.3`$. Then, keeping all the best-fit parameters fixed, we have used the energy scaling laws for dust scattering in Rayleigh-Gans approximation (size $`E^1`$, $`a_{\mathrm{max}}E^2`$) to successfully reproduce also the profile in the 3–5 keV X-ray band, while the data at 1.3–1.7 keV look slightly overestimated (Figure 3).

## DISCUSSION

The main results of the above analysis are the following:

* By comparing our best-fit $`a_{\mathrm{max}}`$ ($`0.90\mu \mathrm{m}`$) and $`q`$ (3.3) with the values averaged over a large sample of sources ($`a_{\mathrm{max}}=0.17\mu \mathrm{m}`$ and $`q=3.9`$; Predehl and Schmitt, 1995), the distribution of grain sizes turns out to be flatter and more extended to larger sizes than the average; while, owing to the uncertainties involved in the Predehl and Schmitt empirical relation, the best-fit $`\tau _{\mathrm{sca},1}`$ (1.50) is in reasonable agreement with the value (1.07) derived from $`N_\mathrm{H}`$,
* The detected slight discrepancy from the 1.3–1.7 keV data (about 30%, with the model overestimating the data) is naturally explained in terms of the invalidity of the Rayleigh-Gans theory in soft X rays. Simulations with a similar $`a_{\mathrm{max}}/E`$ ratio (Smith and Dwek, 1998) give an offset in the same direction as observed, ranging from 25% (for graphite grains) to 50% (for silicate grains).
* The radial profile in the 5–8 keV band, that on the basis of the above analysis should be only slightly affected by dust scattering, shows instead a tail broader than the PSF (Figure 4). This effect, although involving very weak surface brightnesses ($`10^4`$ times the peak), is significant. It may either be an instrumental bias that we have not identified yet (and that is not documented), or the evidence of an intrinsic extended component, even though much less prominent than previously thought.
* We have tried to reproduce the data in the 5–8 keV band convolving various simple source models with the PSF. No source smaller in size than the radio plerion has been found to reproduce correctly the X-ray wings. A good fit is obtained modelling the intrinsic halo as a uniform sphere of radius $`200^{\prime \prime }`$ and surface brightness $`10^4`$ times the inner component (for a total flux which is only 2% of that of the inner source, Figure 4), although this solution is not unique.

## THE NORTH SPUR

A spur has been detected in the northern sector of the halo (see Figure 1, as well as Warwick et al., 2001; Safi-Harb et al., 2001). The presence of this intriguing feature cannot be explained by our dust model. A pure power-law fit to the spur EPIC/PN spectrum gives acceptable $`\chi ^2`$ values. However a pure power-law fit to the spur gives $`N_\mathrm{H}=1.2\times 10^{22}\mathrm{cm}^2`$, namely a factor 2 less than for the other spectra: this discrepancy may indicate that a pure power-law fit is not appropriate.

If a thermal component is included in the fit model, it is significatively detected in the spur, while it is not detected in other halo regions. This is shown by the (1-$`\sigma `$, 2-$`\sigma `$, 3-$`\sigma `$) confidence contours plotted (Figure 5) in the parameters plane for the thermal component: the close contours refer to the area labelled Nspur (see Figure 1 for the areas definition), while fits to the areas labelled halo1, halo2, and halo3 give only upper limits to the thermal fraction (open contours). In fact, adding a thermal component in the fit to the Nspur region, the estimated $`N_\mathrm{H}`$ is now between 1.7 and $`4.0\times 10^{22}\mathrm{cm}^2`$, values consistent with the other halo regions.

## ACKNOWLEDGEMENTS

This work has been partly supported by the Italian Ministry for University and Research (MIUR) under Grant Cofin 2001–02–10.


text:
# 1 Introduction

## 1 Introduction

Recently computer simulations on cellular automata type models have been suggested which can be applied to simulate the development of a surface layer under different physical conditions. A randomly occupied square lattice is taken and particles are allowed to stick at vacant sites or evaporate from occupied sites according to a prescribed set of rules. Appropriate probabilities for the ‘growing’ and ‘culling’ are set according to temperature and pressure conditions required.We can also introduce an effective interaction between the particles by making the probabilities a function of the surroundings of a site, e.g. the number of nearest neighbours.

In our earlier work  we started with a randomly occupied square lattice and let particles evaporate with a probability which decreased with the number of occupied neighbouring sites. This amounted to an attractive interaction between the particles. The condensation was a function of only the superincumbent pressure. The present paper is an extension of the earlier work, here we introduce an effective repulsive interaction between particles, by allowing the condensation probability at an empty site to decrease with the number of occupied neighbours. The evaporation probability is a function of temperature only. The results are very interesting. We study as in the previous work  evolution of a mobile layer by analytical and numerical methods and an immobile layer by computer simulation. In both cases we find a steady final coverage under some conditions, but a two-point limit cycle showing oscillations in the coverage when temperature is high and pressure low.

In the next section we describe the model. The results are presented in the third section and discussion and conclusions in the last.

## 2 The Model

The process of condensation/evaporation or adsorption/desorption is described traditionally by two different sets of models – one for a mobile adsorption layer and one for an immobile layer. Probabilities for sticking and evaporation on a two-dimensional monolayer are specified according to the physics behind the model.These are functions of the temperature, superincumbent pressure and the existing coverage. At equilibrium, the sticking and evaporation probabilities are set equal and the resulting equation is solved to get the equilibrium coverage at that temperature and pressure.

The simplest mobile layer model is a two-dimensional ideal gas, and the improved versions include interaction between particles, similar to a two-dimensional Van der Waal gas. The so-called immobile layer models introduce a sticking probability, depending on how long a molecule in the vapor above the surface is in contact with a surface site. The simplest ‘immobile model’ is the Langmuir equation derived as follows.

For vapor condensation the particle flux, i.e. the number of particles deposited per unit time per unit surface area is equal to:

$$c=\frac{P\lambda }{h}$$
(1)
where $`P`$ is the pressure,

$$\lambda =\sqrt{h^2/(2\pi mkT)}$$

is de Broglie length, $`h`$ is Planck constant. The probability of condensation per unit time is obtained by dividing $`c`$ by the total number of sites available for condensation. And the same may be done for the evaporation probability. Flux of evaporation from saturated surface (at $`p=1`$) may be approximated by  :

$$d=\frac{kT}{h\lambda ^2}\mathrm{exp}(E_e/kT)$$
(2)
Condensation probability is set equal to evaporation probability, giving equilibrium:

$$\frac{P\lambda ^3}{kT}(1p)=\mathrm{exp}(E_e/kT)p$$
(3)
and we have the simple Langmuir equation

$$Pb=\frac{p}{1p}$$
(4)
where

$$b=\frac{\lambda ^3}{kT}\mathrm{exp}(E_e/kT),$$

$`E_e`$ is activation energy for evaporation.

In this approximation it is assumed that evaporation energy is the same for all configurations.

In the present work we introduce a repulsive interaction between particles, as follows \- condensing particles are less likely to stick at sites surrounded by more neighbors. This mimics a situation where the adsorbing particle has to overcome the potential barrier due the repelling neighbors to reach a vacant site on the substrate. It is implied that the substrate attracts the adsorbing particle.

We define a probability for condensation on a vacant site $`f_c`$ determined by the pressure as

$$f_c=P/P_s$$
(5)
where $`P_s`$ is a saturation pressure and $`0<f_c<1`$. The growth probability at a vacant site with $`n`$ occupied neighbors is proportional to $`f_{c}^{}{}_{}{}^{n}`$. We have

$$P_{gr}=(1p)^5f_{c}^{}{}_{}{}^{4}+4(1p)^4pf_{c}^{}{}_{}{}^{3}+6(1p)^3p^2f_{c}^{}{}_{}{}^{2}+4(1p)^2p^3f_c+(1p)p^4$$
(6)
The successive terms in the RHS above represent the probability of a vacant site having respectively 4,3,2,1 and 0 vacant near neighbor sites. $`P_{gr}`$ can be written in the compact form

$$P_{gr}=(1p)(pf_c+1p)^4$$
(7)
The culling probability is assumed to be a function of temperature only

$$P_{cull}=p(kT/P_s\lambda ^3)exp(E_e/kT)$$
(8)
$`P_s`$ is temperature dependent, and this must be taken into account. But we can approximate it by the simple relation that follows from equilibrium relation on the surface of fluid

$$P_s=\frac{kT}{\lambda ^3}\mathrm{exp}(E_{vap}/kT)$$
(9)
where $`E_{vap}`$ is the enthalpy of vaporization of a fluid that supports vapour.

In this formulation the simplest condensation evaporation equilibrium condition is

$$\frac{P}{P_s}(1p)=\frac{kT}{P_s\lambda ^3}\mathrm{exp}(E_e/kT)p=\mathrm{exp}\mathrm{\Delta }E/kTp$$
(10)
where $`\mathrm{\Delta }E=E_eE_{vap}`$

Here too condensation probability $`f_c`$ (or reduced pressure $`P^{}`$) is

$$f_c=P/P_s$$
(11)
and the evaporation probability is

$$f_e=\mathrm{exp}(\mathrm{\Delta }E/kT)$$
(12)
That means that the molecules evaporate from the our surface with probability 1 when $`E_e=E_{vap}`$, when energy of evaporation from surface $`E_e`$ equals to energy of evaporation from surface of liquid that support vapor $`E_{vap}`$.

So, we can rewrite $`P_{cull}`$ as

$$P_{cull}=p\mathrm{exp}(1/T^{})$$
(13)
where $`f_c=P/P_s`$ and $`T^{}=kT/\mathrm{\Delta }E`$

Our model is complementary to the Fowler-Guggenheim model , where the condensation probability is dependent only on pressure, but the evaporation probability decreases with the number of occupied near neighbors. The neighbors are assumed to attract the particle and create a barrier to evaporation. This type of interaction was incorporated in Dutta et al .

Now we study the behavior of this system for different values of the temperature and pressure parameters. We look at the two different situations the mobile surface layer and the immobile layer.

### 2.1 Mobile surface layer

Let us suppose that the energy barrier for motion of the deposited particles along the surface is very low so a continuous rearrangement is going on. In this case the distribution of particles always becomes random and the probabilities in eq(6) are always valid. So one can equate the growth and culling probabilities in the steady state and solve for the coverage $`p`$. We confine ourselves to ranges of $`T^{}`$ and $`P^{}`$, which give growth and culling probabilities in the range 0 – 1. We find that meaningful positive values for the coverage between 0 and 1 are obtained for a very narrow range of temperature and pressure. The results are shown in the next section. We may also see the how the coverage evolves using the following iterative procedure

$$p_f=p_i+P_{gr}P_{cl}$$
(14)
and

$$p_i=p_f$$
(15)
here $`p_i`$ is the initial coverage and $`p_f`$ the final. If there is a steady solution the final coverage $`p_f`$ converges to a steady value. We find however that in some combinations of the pressure and temperature parameters the coverage converges to a stable limit cycle, oscillating between two positive and realistic values. The solution for p in eq.(6) now becomes an unstable fixed point. This is very interesting and reminiscent of the discovery of oscillatory behaviour in real chemical systems .

The results obtained here can be verified from a computer simulation on a finite system. This was done as follows.

In the mobile case, the adsorbed molecules can move laterally on the surface. The physical situation simulated is described by eq.(6) . A two-dimensional square lattice of unit spacing and size $`300\times 300`$, is occupied randomly with an initial coverage $`p_{initial}`$. Every occupied site is assigned the value $`1`$, empty sites are assigned the value $`0`$. The occupied sites are then culled parallally with a probability which is the same for each occupied site

$$g=exp(1/T^{})$$

. For very large $`T^{}`$ the culling probability approaches 1. The vacant sites are filled with a probability dependent on the number of vacant near neighbor sites n as $`f_{c}^{}{}_{}{}^{n}`$. After one round of growth and culling is complete, the concentration of the occupied sites $`p_{final}`$ is calculated.

In the next time-step, the $`p_{final}`$ of the previous time-step becomes the new $`p_{initial}`$. The square lattice is then randomly occupied afresh with this $`p_{initial}`$. A complete time-step begins with the random occupation of all sites with a $`p_{initial}`$ and ends with the assignment of the $`p_{final}`$ to the $`p_{initial}`$ of the next time-step. This iterative process stops when $`p_{final}`$ saturates with increasing time to a definite value. In the simulation, we checked up to 50,000 time-steps. The ‘mobility’ of the molecules is simulated by the randomization of the concentration $`p_{initial}`$ in the beginning of every time-step.

### 2.2 Immobile Surface Layer

If we consider the energy barrier to surface movement of the particles to be large enough the arguments presented in the last subsection are not valid. We start with a certain concentration of particles randomly distributed on a square lattice and allow growing and culling according to our rules. But after one time-step the distribution is no longer random, so we cannot find out the number of new growth or evaporation sites as in eq.(6). We can however simulate the system as above without the randomization procedure after each round of growth and culling. This mimics the development of a surface layer where the particles are not free to move laterally on the surface. The results in this case are quite different .

## 3 Results

We present the results for the mobile and immobile layer case. The mobile layer results are done analytically or by numerical iteration and verified by computer simulation. The immobile layer results are done by computer simulation as calculation is not possible here.

### 3.1 Results for the mobile layer

The results are shown in figure 1. The final coverage $`p`$ for each $`P^{}`$, on increasing $`T^{}`$, from a low value shows a bifurcation at a certain point to a two state limit cycle. The temperature and coverage at this point are denoted respectively by $`T_{}^{}{}_{b}{}^{}`$ and $`p_b`$, the inset shows the variation of $`T_{}^{}{}_{b}{}^{}`$ and $`p_b`$ with $`P^{}`$.

These results are independent of the initial coverage. For $`T^{}`$ very large and $`P`$=1 there is an oscillation between $`p_{initial}`$ and $`1p_{initial}`$, as is also evident from eq.(6).

### 3.2 Results for the immobile surface layer

The immobile situation can be studied only by computer simulation. Starting with an initial random distribution of particles the system is allowed to evolve and the final coverage for different $`P`$ and $`T^{}`$ is studied. In this case as there is no randomization, the system takes a very long time to reach the steady state, and we have observed up to 30,000 time steps. Steady state solutions of $`p_{final}`$ are obtained for $`P_{cull}<1`$. The steady state isobars and isotherms are shown in Fig.2. and Fig.3. respectively. At $`P_{cull}`$=1.0, oscillations between a two-point limit sets in for all values of $`P`$. However for every $`P`$, there is a certain initial coverage, $`p_\alpha `$ that yields a final coverage oscillating in a narrower confine between $`0.4`$ and $`0.5`$ for all $`P`$. Even a departure of 0.05 from $`p_\alpha `$, shifts the final coverage to the fixed point limits of the wider limit cycle. The rate of convergence to $`p_{final}`$ from different $`p_{initial}`$s other than $`p_\alpha `$, is slower for higher values of $`P`$. For $`p_{initial}`$ very close to $`p_\alpha `$ in the case of $`P=0.9`$, we checked for 50,000 time steps to detect whether the final coverage was approaching the fixed points of the wider limit cycle. The $`p_\alpha `$ increases slowly with $`P`$ as displayed in Table.I.

## 4 Discussion

We have demonstrated that the condensation/evaporation problem seen as a modification of the percolation problem on a square lattice shows very interesting possibilities. There is firstly, the question of mobile and immobile layer models of chemical physics . From our point of view equating the growth and culling probabilities to find the steady state coverage is valid only when the particles can rearrange very fast to preserve the random distribution. And we demonstrate that making the particles truly immobile,on the two-dimensional surface gives a quite different final coverage. This is also dependent on the initial coverage. This was discussed in the previous work  with an attractive interaction between particles. Careful experiments are needed to test the justification of these arguments.

The most striking result of this work is the observation of oscillations in the condensation/evaporation problem, using such an extremely simple model, without taking any details of the system into account. Earlier well known work on non-linear mathematical equations  show the appearance of bifurcations, repeated period doubling sometimes leading to chaos. And such situations are found in real life in problems of physics, chemistry and biology. Chemical systems are known to exhibit oscillations in concentration of reactants and products in chemical reactions . There is also some study of such a phenomenon in condensation/evaporation problems . Experiments with repelling particles adsorbing on a substrate are reported , but in this experiment two different types of particles are involved.

One possibility of having a repulsive interaction is with similarly charged particles condensing on an oppositely charged substrate.

So in conclusion the simple stochastic model with growth and evaporation on a square lattice gives rise to a host of interesting phenomena. With an attractive interaction between particles phase transitions and hysteresis was observed , While a repulsive interaction gives oscillations in coverage under certain conditions.

| $`P`$ | 0.1 | 0.2 | 0.3 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- |
| $`p_\alpha `$ | 0.23 | 0.25 | 0.265 | 0.305 | 0.33 | 0.35 | 0.39 | 0.43 |

Table I : The variation of $`p_\alpha `$ with $`P`$

Figure captions:

Figure 1 : Isobars for the final coverage for the mobile layer. At $`T_{}^{}{}_{b}{}^{}`$, the steady solution for $`p_b`$ shows the bifurcation to a limit cycle. The inset shows variation of $`T_{}^{}{}_{b}{}^{}`$ and $`p_b`$ with $`P^{}`$.

Figure 2 : Isobars for the immobile case showing variation of coverage with $`exp(1/T^{})`$. The values of $`P`$ for each curve are shown in the figure.

Figure 3 : Isotherms for the immobile case showing variation of coverage with pressure $`P`$. Stable $`p`$ are obtained for $`exp(1/T^{})`$ upto 0.9, for $`exp(1/T^{})`$=1.0 oscillations are obtained.


text:
# Goodness-of-fit of the Heston model

## 1 Assessed model

Standard models of stock-market fluctuations predict a normal (Gaussian) distribution for stock price log-returns . However, the empirical distribution exhibits significant kurtosis with a greater probability mass in its tails and centre of the distribution . In Quantitative Finance 2002 (), Drăgulescu and Yakovenko (DY) proposed an analytical formula for the probability density function (*pdf*) of stock price log-returns, based on the Heston model  — a geometric Brownian motion for the log-returns time series coupled with a stochastic mean-reverting volatility. The resulting *pdf* is claimed to outperform the Gaussian on a large scale of time lags ($`t=1,5,20,40`$ and $`250`$ days).

### 1.1 Heston model and DY formula

This model starts with the usual assumption that the price $`S_t`$ follows a geometric Brownian motion described by the following stochastic differential equation

$$dS_t=\mu S_tdt+\sigma _tS_tdW_t^{(1)}$$
(1)
| where | $`\mu `$ is the trend of the market, |
| --- | --- |
|  | $`\sigma _t`$ is the volatility, |
|  | $`W_t^{(1)}`$ is a standard Wiener process. |

Log-returns $`r_t=Log\frac{S_t}{S_0}`$ and centred log-returns $`x_t=r_t\mu t`$ are then introduced:

$$dr_t=(\mu \frac{v_t}{2})dt+\sqrt{v_t}dW_t^{(1)}\text{since}\sigma _t=\sqrt{v_t}$$
(2)
and

$$dx_t=\frac{v_t}{2}dt+\sqrt{v_t}dW_t^{(1)}$$
(3)
Then, instead of having a constant volatility $`\sigma _t=\sigma `$ as in the Bachelier-Osborne model , the Heston model assumes the variance $`v_t=\sigma _t^2`$ obeys the following mean-reverting stochastic differential equation

$$dv_t=\gamma (v_t\theta )dt+k\sqrt{v_t}dW_t^{(2)}$$
(4)
| where | $`\theta `$ is the long time mean of $`v_t`$, |
| --- | --- |
|  | $`\gamma `$ is the relaxation rate of this mean, |
|  | $`k`$ is a constant parameter called the variance noise, |
|  | $`dW_t^{(2)}`$ is another standard Wiener process, not necessarily correlated with $`dW_t^{(1)}`$. |

DY solve the forward Kolmogorov equation that governs the time evolution of the joint probability $`P_t(x,v|v_i)`$ of having the log-return $`x`$ and the variance $`v`$ for the time lag $`t`$, given the initial value $`v_i`$ of the variance

$`{\displaystyle \frac{\delta }{\delta t}}P`$ $`=`$ $`\gamma {\displaystyle \frac{\delta }{\delta v}}[(v\theta )P]+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\delta }{\delta x}}(vP)`$ (5)
$`+\rho k{\displaystyle \frac{\delta ^2}{\delta x\delta v}}(vP)+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\delta ^2}{\delta x^2}}(vP)+{\displaystyle \frac{k^2}{2}}{\displaystyle \frac{\delta ^2}{\delta v^2}}(vP)`$
They introduce a Fourier transform to solve analytically this equation, and obtain the following expression for the probability distribution of centred log-returns $`x`$, given a time lag $`t`$:

$$P_t(x)=\frac{1}{2\pi }_{\mathrm{}}^+\mathrm{}𝑑p_xe^{ip_x+F_t(p_x)}$$
(6)
with

$$F_t(p_x)=\frac{\gamma \mathrm{\Gamma }\theta t}{k^2}\frac{2\gamma \theta }{k^2}ln[\frac{cosh\frac{\mathrm{\Omega }t}{2}+\frac{\gamma }{\mathrm{\Omega }}sinh\frac{\mathrm{\Omega }t}{2}}{cosh\frac{\mathrm{\Omega }t}{2}+\frac{\mathrm{\Gamma }}{\mathrm{\Omega }}sinh\frac{\mathrm{\Omega }t}{2}}]\frac{2\gamma \theta }{k^2}ln[cosh\frac{\mathrm{\Omega }t}{2}+\frac{\mathrm{\Omega }^2\mathrm{\Gamma }^2+2\gamma \mathrm{\Gamma }}{2\gamma \mathrm{\Omega }}sinh\frac{\mathrm{\Omega }t}{2}]$$
(7)
| where | $`\mathrm{\Gamma }=\gamma +i\rho kp_x`$, |
| --- | --- |
|  | $`\rho `$ is the correlation coefficient between the two Wieners $`W_t^{(1)}`$ and $`W_t^{(2)}`$, |
|  | $`\mathrm{\Omega }=\sqrt{\mathrm{\Gamma }^2+k^2(p_x^2ip_x)}`$, |
|  | $`\gamma `$, $`\theta `$, $`k`$ and $`\mu `$ are the parameters of the Heston model. |

Eqn. (6), hereafter designated as the “DY formula”, is the central result of DY’s paper . It gives the expected probability distribution of centred log-returns $`x`$, given a time lag $`t`$. An asymptotic analysis111See  Part VI of $`P_t(x)`$ shows that it predicts a Gaussian distribution for small values of $`|x|`$, and exponential, time dependent tails for large values of $`|x|`$.

To confront their formula with observed log-returns, DY take the Dow-Jones Industrial Average from January 04, 1982 to December 31, 2001, and train the four parameters of the Heston model, $`\gamma `$, $`\theta `$, $`k`$ and $`\mu `$, to fit the empirical distribution by minimising the following square-mean deviation error

$$E=\underset{x,t}{}|logP_t^{}(x)logP_t(x)|$$

for all available values of log-returns $`x`$, and time lags $`t=1,5,20,40`$ and $`250`$ days, where $`P_t^{}(x)`$ is the empirical probability mass and $`P_t(x)`$ the one predicted by the DY formula. In their results, they set the correlation coefficient $`\rho `$ to zero, since (i) their trained parameter $`\rho ^{trained}`$ is almost null ($`\rho ^{trained}0`$), and (ii) they do not observe any significant difference, when fitting empirical data, between taking $`\rho ^{trained}0`$ or $`\rho =0`$. Hence, they reduce the complexity of their formula. Minimising the deviation of the log instead of the absolute difference $`|P_t^{}(x)P_t(x)|`$ forces the parameters to fit the fat tails instead of the middle of the distribution, where the probability mass is very high.

A replication of their results, using the same dataset and the same method, is shown in Figure 1. The *pdfs* for different time lags from 1 to 250 are shown each shifted upwards by a factor of 10 for clarity. Theoretically, this result is brilliant, since the authors obtain an analytic expression for the *pdf* $`P_t(x)`$. Furthermore, their model (plain line) seems to fit the empirical data (dots) far better than the Gaussian (dash line), especially if we look at the fat tails.

But we argue that the method they use to evaluate the goodness-of-fit of their model suffers from two major drawbacks when pre-processing the data.

### 1.2 Pre-processing the data

Firstly, DY trimmed the log-returns time series, rejecting any value out of the boundaries presented in Table 1.222This step is not mentioned in DY’s paper. Before applying this trimming method, strange points used to appear in our results. Then we contacted the authors who informed us they had trimmed the log-returns time series using the boundaries specified in Table 1. By doing so, they remove most of the leptokurtosis of the original dataset (the positive peakiness of the *pdf* in the centre), which is precisely one of the discrepancies from the Gaussian that they should try to fit.

We visualise in Figure 2 the effect of trimming the data: all of the log-returns outside the boundaries, represented here by the two horizontal lines, were trimmed. We believe this way of trimming the data is unfair, because it removes information from the dataset. Given that the model is supposed to outperform the Bachelier-Osborne model, and specially to fit the kurtosis and the fat tails, removing extreme values (that belong to the fat tails, and produce kurtosis) is counter-productive. Even the normal distribution could fit the data quite well in these conditions . To prove this, we compare the kurtosis of trimmed and untrimmed empirical data in Table 2.

Obviously, the kurtosis disappears when data are trimmed, which makes the dataset considerably smoother. Back in 1965, Fama  had already made similar criticism of Kendall’s experiments , in which the latter considered outliers so extreme that he just dropped them.

A second major concern is that DY use a single log-return time series of overlapping returns. For a given index *I* at a given period, let us say the Dow Jones Industrial Average from January 04, 1982, to December 31, 2001, and a given time lag $`\tau `$, let us say $`\tau =5`$ days, the raw close price dataset *closePrice* is composed of *n* close prices, here $`n=5050`$. When DY compute the log-returns dataset *logReturns* starting from *closePrice*, they obtain the following time series:

$$logReturns=\{r_t|t[1,n\tau ]\}$$

where $`r_t=log\frac{P_{t+\tau }}{P_t},t[1,n\tau ]`$.
In our example, we would have

$`logReturns`$ $`=`$ $`\{r_1,r_2,\mathrm{},r_{n\tau }\}`$
$`=`$ $`\{log{\displaystyle \frac{P_{1+\tau }}{P_1}},log{\displaystyle \frac{P_{2+\tau }}{P_2}},\mathrm{},log{\displaystyle \frac{P_n}{P_{n\tau }}}\}`$
$`=`$ $`\{log{\displaystyle \frac{P_6}{P_1}},log{\displaystyle \frac{P_7}{P_2}},\mathrm{},log{\displaystyle \frac{P_{5050}}{P_{5045}}}\}`$
Thus, they obtain a single dataset of $`n\tau `$ log-returns. We believe this way of computing the log-returns time series is biased, because it “re-uses” the data. Indeed, let us assume that a crash occurs at time $`t^{}`$. Then they will take into account this specific event $`\tau `$ times exactly in their dataset in log-returns $`\{r_{t^{}\tau },r_{t^{}\tau +1},\mathrm{},r_{t^{}1}\}`$. Working on this time series will definitely fatten the tails, since every shock in the original price time series will appear exactly $`\tau `$ times in the log-returns. Hence, the resulting log-returns time series is not an accurate representation of the market moves anymore; this gives a decisive advantage to their model compared with the Gaussian, since they train their parameters to fit specifically the fat tails.

Finally, DY do not statistically test the goodness-of-fit of their formula, but rather rely exclusively on plots that truly look good.

## 2 Goodness-of-fit tests

For these reasons, we decided to keep the outliers in our dataset and to use multiple non-overlapping log-return time series. We perform our tests on the Dow-Jones Industrial Average, from January 04, 1982 to December 31, 2001. From the raw price time series, we derive $`\tau `$ log-return time series $`logReturns_j`$ ($`j[1,\tau ]`$) of cardinality $`m=[\frac{n}{\tau }]`$, called “paths”, instead of a single log-return time series $`logReturn`$ of cardinality $`n\tau `$. This method sounds more reasonable for very different reasons: firstly, when economists or traders talk about weekly returns ($`\tau =5`$ days, holidays excluded), they mean returns from Mondays to Mondays (for instance), and not from Monday to Monday, Tuesday to Tuesday, etc. Secondly, this method gives a better view of the real market, since tails are not artificially fattened. And finally, this method will enable us to assess the statistical accuracy of the estimators we will compute.
We can now perform our goodness-of-fit tests on our different models: the Gaussian (*normPDF*), the curve resulting from the DY formula (*draguPDF*), and a Neural Network333Even if *draguPDF* is supposed to fit the empirical distribution, *empPDF*, better than *normPDF*, we want to compare it with the best fit possible, the one obtained with a Neural Network. This Neural Network must be as simple as possible, but should fit the main characteristics of the empirical time series, fat tail and kurtosis. The structure chosen was the following: it is a feed-forward back-propagation network, with a five node hidden layer and a single node output layer. The transfer functions are respectively $`tansig`$ and $`purelin`$, where $`tansig(n)=\frac{2}{1+e^{2n}}1`$ and $`purelin(n)=n`$. This structure appears to be a good trade-off between the complexity and the goodness of fit. intentionally overfitting the data, that will be used as a benchmark (*nnPDF*).

### 2.1 Kolmogorov-Smirnov Statistic

DY claim that their model fits the empirical data of the concerned dataset better than the Gaussian for any time lag. To check this, we use the Kolmogorov-Smirnov Statistic, based on the maximal discrepancy between the expected and the observed cumulative distributions, for any log-return $`x`$. This statistic is suitable for testing only a simple hypothesis, for instance a Gaussian with known parameters $`\mu `$ and $`\sigma `$, but not a composite hypothesis (a class of Gaussians, or a Gaussian with $`\mu `$ and $`\sigma `$ derivated from the sample dataset being tested). Unfortunately, whatever the model, we always derive the parameters ($`\mu `$ and $`\sigma `$ for *normPDF*, $`\gamma `$, $`\theta `$, $`k`$ and $`\mu `$ for *draguPDF*, the weights and biases for *nnPDF*) from the initial dataset. By performing this test with parameters derived from the dataset being tested, we expect the statistic to be large enough to reject the simple hypothesis, and *a fortiori* the composite hypothesis (). But if the value of the statistic $`Z`$ is small enough to accept the simple hypothesis, it does not mean that we can accept the composite hypothesis.

Methodology - For each time lag, we compute the log-returns dataset, and we divide it into paths. For each path and each model, we build the empirical cumulative density function *empCDF* and the expected CDF *modelCDF* (*normCDF*, *draguCDF* or *nnCDF*), and we compute the KS-statistic $`Z`$ (see Figure 3 for an example). We present in Tables 3, 4 and 5 the mean $`\stackrel{~}{Z}`$ and standard deviation $`\sigma _Z`$ of $`Z`$ over the different paths, and the associated p-value444The p-value is the probability of observing the given sample result under the assumption that the null hypothesis (the tested model) is true. If the p-value is less than the level of significance $`\alpha `$, then you reject the null hypothesis. For example, if $`\alpha =0.05`$ and the p-value is 0.03, then you reject the null hypothesis. interval $`p(\stackrel{~}{Z}+\sigma _Z)p(\stackrel{~}{Z})p(\stackrel{~}{Z}\sigma _Z)`$.

Results - Firstly, we observe an important variance, over the different paths, in the statistic $`Z`$: the standard deviation $`\sigma _Z`$ is not negligible in comparison with the mean $`\stackrel{~}{Z}`$. This is an evidence that the paths are not equivalent, which legitimates, *a posteriori*, our sampling method. It comes from the high heterogeneity of the dataset, which makes our tests less robust. But any test performed on this heterogeneous dataset would suffer from the same problem. Even though this apparent lack of consistency prevents us from drawing any strong and global conclusion, the knowledge of the mean and standard deviation $`\stackrel{~}{Z}\pm \sigma _Z`$ provides us with a fair overview of the statistic $`Z`$.
On plots, the DY formula seems to fit the empirical cumulative distribution better than the Gaussian. But in fact, on average, both models are rejected for high frequencies (for $`\tau =1\text{and}5`$ days) at the 0.01 level of significance. Even the Neural Network is rejected for a one day time lag. This rejection of the three models may come from the fact that this test is based on the maximum discrepancy between the empirical and the theoretical cumulative distributions, *for any* $`x`$. To pass this test, a model must fit the observed data sufficiently well everywhere, i.e. in the tails (problem of fat tails) and in the middle (problem of high kurtosis for high frequencies) of the distribution.
We point out that even if the DY formula is rejected for a one day time lag, the statistic $`Z`$ is smaller than the Gaussian one (0.109 vs 0.131), which is an indication that the model fits the data a bit better. For other time lags, the p-value are equivalent: both models are systematically rejected for 5 days ($`p0.01`$), sometimes rejected for 20 days ($`p(\stackrel{~}{Z}+\sigma _Z)0.01`$, but $`p(\stackrel{~}{Z}\sigma _Z)0.05`$) and never rejected for higher frequencies ($`p0.05`$). For medium and low frequencies, the fact that the simple hypothesis is not rejected does not guarantee that the composite hypothesis can be accepted.

Conclusion - The Kolmogorov-Smirnov goodness-of-fit test rejects both the Gaussian and the DY formula for high frequencies ($`\tau =1\text{and}5`$ days). For medium and low frequencies, we cannot come to a firm conclusion because of the theoretical limits of this test. To continue with the investigation, we need a more powerful statistical test that can be used even when the parameters of the model are derived from the dataset being tested. The $`\chi ^2`$ statistic is suitable in those conditions.

### 2.2 $`\chi ^2`$ Statistic

The $`\chi ^2`$ goodness-of-fit test, based on binned data, is a powerful statistical tool to test if an empirical sample comes from a given distribution. Contrary to the Kolmogorov-Smirnov test, it is designed to evaluate a composite hypothesis, i.e. the parameters of the model can be derivated from the empirical dataset tested. This test is a good trade-off between the goodness-of-fit of a model (the better fit, the smaller the $`\chi ^2`$ statistic) and its complexity (the more complex, the smaller p-value). Indeed, even if a model fits the empirical data very well, a too large complexity may penalise its p-value, so that it can still be rejected. Finally, to be meaningful, this test must be performed using relatively large bins, and a critical value of 5 expected observations per bin is regarded as a minimum.

Methodology - If we perform this test with equal size bins, then the fats tails will be trimmed (there are less than 5 expected log-returns per bin in the tails) and will not participate in the value of the statistic, making the test irelevant. Instead, we split the log-return axis into *equal expected frequency bins*, so that all of the log-returns participate in the value of the statistic. We use an expected frequency of 5 log-returns per bin. Unfortunately, this test cannot be performed for large time lags, because of the lack of data. Indeed, for the Dow-Jones index from 1982 to 2001 for instance, we have initially around 5000 close prices, which means that for a time lag of 250 days, each path will have only about 20 log-returns. In those conditions, because of the critical value of 5 log-returns per bin, we will have at best 4 bins, which is too small to perform a relevant test.

Results - We present our results of the $`\chi ^2`$ goodness-of-fit test in Tables 6, 7 and 8. The degree of freedom is given by $`df=noBins1m`$, where $`m`$ is the number of parameters of the model ($`m=2`$ for the Gaussian, $`m=4`$ in the DY formula, and $`m=11`$ for the Neural Networks if we count the weights and the biases). For large time lags, $`df`$ becomes smaller and smaller because $`noBins`$ decreases, as explained above.
Concerning the Neural Network, we cannot perform this test for time lags higher than 40 days, or else the degree of freedom decreases to zero. This is due to the relatively high number of parameters ($`m=11`$). With a structure even more complicated, we could not have performed the test at all, except for high frequencies.
First we notice that the Neural Network’s $`\chi ^2`$ statistic is slightly smaller than the DY’s one, itself smaller than the Gaussian’s one, for all paths with a time lag from $`\tau =1`$ to $`\tau =80`$ days. It means that the Neural Networks fits empirical data better than the DY formula, which itself has a better fit than the Gaussian. But there is a price to pay, in terms of complexity: due to too many parameters (and then a lower degree of freedom), the p-value of the Neural Networks and DY formula are not systematically higher than the p-value of the Gaussian. And it is precisely the p-value that is used to accept or reject a model, not directly the $`\chi ^2`$ statistic.
If we look at the p-value in detail, we observe that

* For $`\tau =1`$, the three models are rejected at a 0.05 level of significance
* For $`\tau =5`$, only the Neural Network is systematically accepted. The Gaussian and DY formula are only accepted in the best situation ($`p(\stackrel{~}{\chi ^2})<0.05<p(\stackrel{~}{\chi ^2}\sigma _{\chi ^2})`$)
* For $`\tau =20`$, the three models are accepted and the DY formula is better than the Neural Network and the Gaussian
* For $`\tau =40`$ and $`\tau =80`$, the DY formula is still accepted, but its p-value is smaller than the one of the Gaussian

Conclusion - Thanks to the $`\chi ^2`$ goodness-of-fit test, we can assert that the DY formula fits empirical data slightly better than the Gaussian, for high and medium frequencies. Nevertheless, both models are rejected for high frequencies ($`\tau =1\text{and}5`$ days), at a 0.05 level of significance. In this sense, these results are consistent with the Kolmogorov-Smirnov goodness-of-fit test.
We also observe a clear shift in the goodness-of-fit of the models around $`\tau =40`$ days: the probability of accepting the Gaussian becomes larger than the probability to accept the DY formula (and even the Neural Network) due to the lower complexity of the Gaussian (two parameters instead of four and eleven respectively).
To put it in a nutshell, using a complex model, such as the DY formula or a Neural Network, is only worthwhile for $`\tau =1,5\text{and}20`$ days. For lower frequencies ($`\tau 40`$ days), the Gaussian is preferable because it is simpler. Given that for these frequencies, we had observed neither fat tails nor kurtosis in the empirical datasets, the Gaussian represents the best trade-off between goodness-of-fit and complexity.

## 3 Conclusion

The DY analytical formula is quite exciting, and would validate the Heston model, if the resulting expected distribution was similar to the empirical one. Our method of pre-processing the data is different from the one suggested by Drăgulescu and Yakovenko, and, we believe, more truly captures the market. We find that the DY formula consistently outperforms the Gaussian in terms of best fit (the error between the model and the observed data), but a higher complexity is the price to pay: the number of parameters is greater in the Heston model, which enables the fitted *pdf* to exhibit kurtosis and fat tails when the empirical dataset does, but also penalizes it when it doesn’t. Hence, in terms of goodness-of-fit (the trade-off between best fit and complexity), the Gaussian is preferable for low frequencies ($`40t250`$ days), since the empirical dataset is then quite smooth and does not exhibit kurtosis. For medium frequencies ($`t=20`$ days), both models are accepted at a .05 confidence level. Finally, for high frequencies, although it performs better than the Gaussian, the DY formula is still rejected ($`t=1`$ and $`5`$ days), mainly because of the extremely high kurtosis of the observed data.


text:
# LSS 1135: an O-Type Spectroscopic Binary in the galactic OB association Bochum 7

## 1 Introduction

LSS 1135 ($`\alpha _{2000}`$ = 8h 43m 46s and $`\delta _{2000}`$ = -46 07$`^{^{}}`$ 14$`^{^{\prime \prime }}`$) was assigned a spectral type OB and a photographic magnitude mpg = 11.2 in the catalog of Luminous Stars in The Southern Milky Way (LSS) (Stephenson & Sanduleak). Moffat & Vogt (1975) proposed from photoelectric photometry that LSS 1135 together with the Wolf-Rayet star LSS 1145 and other 7 OB stars listed in the LSS catalogue, form an OB association with galactic coordinates l = 265.20; b = -2.1. They called this star group Bochum 7 (Bo 7).

Recently Sung et al. (1999) from UBVI CCD photometry of the Bo 7 region have suggested that Bo 7 is part of the larger OB association Vela OB 3.

In this paper we present optical spectroscopic data of LSS 1135, which show that this star is a close binary system. We also present spectral types and radial velocities for the other 7 OB stars included by Moffat & Vogt (1975) in Bo 7, namely LSS 1131, 1132, 1137, 1140, 1144, 1146, 1147.

The paper is organized as follows: in Sect. 2 we describe the observations. In Sect. 3 we discuss the results and we present radial velocities and spectral classifications for stars in Bo 7, and the orbital parameters for LSS 1135. In Sect. 4 we summarize our main results.

## 2 Observations

The observational material consists of 21 photographic spectrograms and 29 digital CCD spectrograms of LSS 1135. We also obtained 13 photographic and 12 digital spectra of the other stars listed by Moffat & Vogt (1975) as members of Bo 7. The instrumental configurations used are detailed in Table 1.

### 2.1 Photographic spectrograms

The photographic spectrograms were all obtained by VSN between February 1982 and March 1985, at the Cerro Tololo Interamerican Observatory (CTIO), Chile. These spectra were secured with the Carnegie Image Tube Spectrograph (CITS) attached to the 1m Yale reflector telescope. All the exposures were made on Kodak III a-J emulsion and widened to 1mm. A He-Ar lamp was used as comparison source.

The photographic spectrograms were measured by VSN for the determination of radial velocities with an oscilloscope microdensitometer (GRANT engine) at Instituto de Astronomía y Física del Espacio, Buenos Aires, Argentina. We also digitized some photographic spectrograms with a GRANT engine at La Plata Observatory, Argentina.

### 2.2 Digital spectra

Digital spectral images of LSS 1135 were obtained between February 1997 and 1999, with the Boller & Chivens (B & C) and the REOSC Cassegrain echelle spectrographs attached to the 2.15 - m telescope at Complejo Astronómico El Leoncito (CASLEO)111Operated under agreement between CONICET, SeCyT, and the Universities of La Plata, Córdoba and San Juan, Argentina. in San Juan, Argentina.

Fourteen spectra of LSS 1135, and eleven of other members of Bo 7, were secured with the B & C spectrograph using a PM 512$`\times `$512 pixels CCD detector with pixel size of 20$`\mu `$m. We used a 600 l mm-1 grating and the slit width was set to 200$`\mu `$.

Fifteen spectra of LSS 1135 were obtained with the REOSC Cassegrain echelle spectrograph using as detector a TEK 1024x1024 pixels CCD, with pixel size of 24$`\mu `$m. We used a 400 l mm-1 grating as cross disperser and the slit width was set to 250$`\mu `$ and 300$`\mu `$. Ten échelle spectra of LSS 1135 were obtained in February 1998 binning the CCD by a factor 2.

He-Ar (or Th-Ar with REOSC spectrograph) comparison arc images were observed at the same telescope position as the stellar images immediately after or before the stellar exposures. Also bias and flat-field frames were obtained every night, as well as spectra of flux and radial velocity standard stars.

All digital, and digitized photographic spectra, were processed and analysed with IRAF222IRAF is distributed by NOAO, operated by AURA, Inc., under agreement with NSF. routines at La Plata Observatory. Radial velocities were determined by fitting Gaussian profiles to the spectral lines. Typical errors in our the radial velocities are approximately 20 km s-1(standard error of the mean) for the instrumental configurations I and III, and 15 km s-1 for the instrumental configuration II. The heliocentric radial velocities of the interestellar absorption lines measured in the high dispersion echelle spectra of LSS 1135 with their respective standard error are: Caii K $`\lambda `$3933 Å = 25 $`\pm `$ 5 km s-1, Caii H $`\lambda `$3968 Å = 24 $`\pm `$ 4 km s-1, Nai $`\lambda `$5890 Å = 25 $`\pm `$ 4 km s-1 and Nai $`\lambda `$5896 Å = 27 $`\pm `$ 4 km s-1. Only one component of each interstellar line is observed in our spectra.

## 3 Results and their discussion

### 3.1 Spectral types and radial velocities of stars in the region of Bo 7

Our photographic spectrograms of stars listed as members of Bo 7 by Moffat & Vogt (1975) at first glance confirm them as OB stars, as their spectra show absorption lines of hydrogen and helium. Approximate spectral types were determined initially by eye estimates of relative absorption line strengths in the photographic spectrograms, and then comparing the digitized and digital spectra with the digital spectral atlas of OB stars (Walborn & Fitzpatrick, 1990). Spectral classifications are indicated in Fig. 1 depicting the spectra of stars observed in the region of Bo 7.

For determination of radial velocities, we measured all absorption lines visible in the spectra. The radial velocity of LSS 1135 was found to be variable from night to night, and this star was therefore included in our program of studies of O type spectroscopic binaries. The results for LSS 1135 are presented separately below. Radial velocities for stars in the region of Bo 7 are listed in Table2.

### 3.2 The spectrum of LSS 1135

Fig. 2 depicts one of the higher S/N digital spectra of LSS 1135 obtained at CASLEO with the instrumental configuration III. This spectrum shows the absorption lines of Heii, Hei and H characteristic of early O spectral type. In addition, several of our spectra show faint emission of Niii $`\lambda `$4634–40 Å. The diffuse interestellar band (DIB) $`\lambda `$4428 Å and other interestellar absorptions are also present in the spectrum of LSS 1135, as well as in the spectra of other stars in Bo 7, except LSS 1132 which is a foreground object (cf. Fig1).

For spectral classification of LSS 1135 in our medium resolution spectra we compared our spectra with the Digital Atlas of OB Stars published by Walborn & Fitzpatrick (1990) in the blue optical spectral region. This comparison yields a spectral type O6.5((f)). In our higher resolution spectra we also used the ratio of Hei $`\lambda `$5875 to Heii$`\lambda `$5411 as described by Walborn (1980) for classification of O–type spectra in the yellow–red optical spectral region. In the blue spectral region of our high resolution spectra we used the equivalent width ratios of the quantitative spectral clasification for O-type stars of Conti & Alschuler (1971) and Conti & Frost (1977). All these confirm for LSS 1135 the spectral type O6.5V((f)).

As an O type binary with a short period (see below), LSS 1135 appears as a potential candidate for showing the effects of colliding stellar winds. For this reason we searched for the Struve-Sahade (S-S) effect, which is produced by colliding stellar winds according to the model proposed by Gies et al. (1997). Examining our spectra in detail, we compared carefully one quadrature with the other in order to detect systematic differences in the relative strength of the He lines observed in the approaching and the receding phases. We did not detect in our data any spectral changes beyond the errors; therefore the S-S effect does not seem to be present in this binary system.

The contribution of the secondary component to the spectral lines of LSS 1135 is not evident in our spectra, except perhaps in the Hydrogen Balmer lines whose radial velocity variations appeared of somewhat lower amplitude than those of He lines, when phased with the binary period (see below). We may therefore assume that the visual absolute magnitude of the secondary component is at least 2 magnitudes fainter than the O6.5((f)) star, which would correspond to an early B spectral type.

### 3.3 The radial velocity orbit of LSS 1135

Table 3 lists the journal of observations and the mean heliocentric radial velocities, obtained averaging Hei $`\lambda `$4471, $`\lambda `$5875 and Heii $`\lambda `$4199, $`\lambda `$4541 and $`\lambda `$4685 absorption lines in the spectra of LSS 1135. The numbers following the mean velocities of the absorption lines indicate how many lines were included in each average.

As is evident from Table 3, the radial velocities of the absorption lines in the spectrum of LSS 1135 show large variations from one night to the other, but the radial velocities obtained during the same night do not show appreciable differences, implying a binary period of a few days. A period search algorithm (Marraco & Muzzio, 1980) was applied to the radial velocities of Table 3.

The best period found was P = 2.75318 $`\pm `$ 0.00002 days and no alias periods with similar probabilities were present. This period was entered as an initial value to calculate the orbital elements of LSS 1135. These were calculated with an improved version of the program initially published by Bertiau & Grobben (1969). In the determination of the orbital elements we assigned weight 10 to the spectra obtained with the instrumental configuration II, and weight 1 for all others. The orbital elements are listed in Table 4. The errors quoted in this table are those calculated by the above mentioned program. The radial velocity variations and orbit of LSS 1135 are illustrated in Fig. 3.

### 3.4 Distance

In Table 5 we list data for the stars in the region of Bo 7. The succesive columns in this table indicate for each star: the number in the LSS catalogue (Stephenson & Sanduleak, 1971); our estimate of the spectral type; V and B-V values from Moffat & Vogt (1975); E(B-V) calculated from intrinsic colors corresponding to the spectral types; the spectrophotometric distance modulus determined adopting the absolute magnitudes corresponding to the spectral types, according the calibration by Schmidt–Kaler (1982); stellar radial velocity (in km s-1) referred to the LSR.

The mean distance modulus for the stars in Bo 7 ( excluding the foreground object LSS 1132, and the emission line star LSS 1147) is 13.50, which corresponds to a distance of 5.0 kpc. This is in good agreement with the distance found by Sung et al. (1999) from CCD photometry. We also note that the photometric distance of LSS 1145 = WR 12 quoted by van der Hucht (2001), is 5 kpc, confirming that this star belongs to Bo 7 star group. The mean (LSR) radial velocities of the cluster members are between 55 and 40 km s-1, (excluding LSS 1146, cf. note in Table 5) implying a kinematical distance between 5.0 and 6.0 kpc in circular galactic rotation model. Thus Bo 7 appears to be located at the extreme of the Perseus spiral arm in our Galaxy (cf. Fig.5 in Russeil, (2003)).

## 4 Summary

A radial velocity study of the spectral lines of LSS 1135, a member of the galactic OB association Bo 7, shows that it is a single–lined O-type binary, which we classify as O6.5V((f)). We find an orbital period of 2.75320 days for this binary. The value of the semi-amplitude (K) of the radial velocity variations of the He lines is 122 km s-1 and the mass function (F(M)) of the binary system is 0.52 M. The secondary component is not detected in our spectra, which means that it should be at least 2 magnitudes fainter, corresponding to an early B type star.

We also have observed 7 other OB stars in the region of Bo 7, for which we present spectral types and radial velocities. From these data we derive a spectroscopic and kinematical distance of about 5.0 kpc for Bo 7, and conclude that the Wolf-Rayet star LSS 1145 = WR 12 is also member of this association.

Acknowledgements We would like to thank the Directors and staff of CTIO and CASLEO for the use of their facilities. We also acknowledge the use at CASLEO of CCD and data acquisition system partly financed by U.S. NSF Grant AST-90-15827 to Dr. R.M. Rich.
This research has received financial support from IALP, CONICET, Argentina.

References Bertiau, F.C. & Grobben, J., 1969, Ric. Astron. Sp. Vaticano 8, 1
Conti, P. & Alschuler, W., 1971, ApJ 170, 325
Conti, P. & Frost, S., 1977, ApJ 212, 728
Gies, D.R., Bagnuolo, W.G. & Penny, L.R., 1997, ApJ 479, 408
Lundstrom, I. & Stenholm, B., 1984, A&AS 58, 163
Marraco, H. & Muzzio, J.C., 1980, PASP 92, 700
Moffat, A.F.J. & Vogt, N., 1975, A&AS 20, 85
Niemela, V.S., 1982, In Wolf-Rayet Stars: Observations, Physics, Evolution, Proc. IAU Sump. 99, eds. C. de Loore & A.J. Willis, (Reidel - Dordrecht)
Russeil, D., 2003, A&A 397, 133
Schaller, G., Schaerer, D., Meynet, G. & Maeder, A., 1992, A&AS 96, 269
Schmidt-Kaler, Th. 1982, In Landolt-Bornstein & H.H. Voigt, (Springer-Verlag, Berlin)
Stephenson, C. & Sanduleak, N., 1971, Publ. Warner and Swasey Obs. 1, 1
Sung, H., Bessell, M.S., Park, B.G. & Kang, Y.H., 1999, JKAS 32, 109
van der Hucht, K.A., 2001, New Astronomy Review 45, 135
Walborn, N., 1980, ApJS 44, 535
Walborn, N. & Fitzpatrick, E., 1990, PASP 102, 379