format: text

text:
# Scattered Light Models of Protostellar Envelopes: Multiple Outflow Cavities and Misaligned Circumstellar Disks

## 1. Introduction

In the standard picture of low mass star formation (e.g., Shu, Adams, & Lizano 1987), Class I protostars are young sources (age $`10^5`$ years) that are still surrounded by massive ($`0.01M_{}`$) circumstellar disks and partially obscured by their natal envelopes. These stars appear highly extincted in the optical, with spectral energy distributions that peak in the mid- to far-infrared (e.g., Lada & Wilking 1984; Beichman et al. 1986; Lada 1987; Myers et al. 1987; Wilking, Lada, & Young 1989; Kenyon et al. 1990; Kenyon et al. 1993a), and exhibit highly polarized reflection nebulae (e.g., Whitney & Hartmann 1993; Whitney, Kenyon, & Gomez 1997; Lucas & Roche 1997, 1998). From observations of Class I protostars has emerged a coherent picture of the protostellar environment. The basic model consists of a dusty infalling envelope with cavities carved by the energetic outflows associated with this phase of star formation. The envelope absorbs the optical light and emits this radiation at longer wavelengths (e.g., Adams & Shu 1986; Adams, Lada, & Shu 1987; Kenyon et al. 1993a); the reflection nebulae are formed by starlight scattering off the cavity walls (Bastien & Menard 1990; Tamura et al. 1991; Whitney & Hartmann 1993; Kenyon et al. 1993b; Whitney et al. 1997; Lucas & Roche 1997, 1998; Lucas, Blundell, & Roche 2000; Reipurth et al. 2000).

Using this basic picture as a foundation, models have been constructed including disks, envelopes, and bipolar cavities that can reproduce in gross detail the observed spectral energy distributions and low-resolution ground-based imaging of many Class I sources (Whitney & Hartmann 1992, 1993; Whitney et al. 1997; Fischer, Henning, & Yorke 1996; Lucas & Roche 1997, 1998; Stapelfeldt et al. 1998b). However, detailed agreement between these models and the observed scattered light patterns of many Class I sources has been less forthcoming, in large part because the existing 2-D models have necessarily assumed axisymmetric circumstellar geometries. High-resolution HST images show that many sources exhibit very complicated scattered light patterns indicative of non-axisymmetric geometries and/or illumination (Burrows et al. 1996; Padgett et al. 1999; Stapelfeldt et al. 1995, 1998a,b, 1999; Terebey et al. 1998; Hartmann et al. 1999; Krist et al. 1998, 1999; Reipurth et al. 2000). Recent imaging polarimetry has provided additional compelling evidence that some sources exhibit departures from an axisymmetric disk-envelope-cavity system (Whitney et al. 1997; Lucas & Roche 1997, 1998; Lucas et al. 2000).

Two particularly good examples of sources with non-axisymmetric circumstellar geometries are T Tauri and IRAS 04016+2610. Ground based near-IR imaging polarimetry of IRAS 04016+2610 indicates two sets of bipolar outflow cavities (Lucas & Roche 1997) and HST/NICMOS images show the complex scattered light pattern at higher resolution (Padgett et al. 1999). Ground based $`J`$, $`H`$, and $`K`$ images of T Tau (Weintraub et al. 1992) were modeled by Whitney & Hartmann (1993) with a standard envelope + bipolar cavity system, but the WFPC2 images (Stapelfeldt et al. 1998b) and near-IR adaptive optics images (Roddier et al. 1999, 2000) clearly show non-axisymmetric nebulosity.

Motivated by these observed departures from axisymmetry, in this paper we develop scattered light models for IRAS 04016+2610 and T Tauri that for the first time include 3-D geometries and illumination by multiple sources. We demonstrate that including the effects of multiple illumination sources, multiple outflow cavities, and mis-aligned circumstellar disks can explain the observed scattered light patterns much more fully than has been possible with 2-D modeling efforts hitherto. We have not undertaken an exhaustive search of parameter space when modeling the two objects as such investigations have been presented previously for axisymmetric models (Whitney & Hartmann 1993; Kenyon et al. 1993b; Whitney et al. 1997). For our models we adopt protostellar envelope parameters derived from previous 2-D modeling of scattered light images and spectral energy distributions of T Tau and IRAS 04016+2610. The effects of changing the evnelope parameters (e.g., accretion rate, centrifugal radius, inclination) have been described in the aforementioned papers and our focus is on the departures from axisymmetry required to obtain a better match between the scattered light models and data. Furthermore, since we view our models as a first step towards unravelling the complexities of the scattered light patterns, we have not attempted to quantify any figure of merit for model fits to the data. Rather, the models we present are the first numerical simulations of published suggestions as to the origin of the asymmetric scattered light patterns of the two sources. In § 2 we describe the scattered light images that form the basis for our modeling, in § 3 we describe the ingredients of our models, in § 4 we present the models for the individual sources, and we conclude in § 5 with a discussion of our findings and further observational tests of our suggested geometries.

## 2. Scattered Light Images

The NICMOS data (F110W, F160W, F187W, F205W) for IRAS 04016+2610 were presented by Padgett et al. (1999), and the WFPC2 data (F555W) on T Tau by Stapelfeldt et al. (1998).

The reduced and PSF-subtracted F555W image of T Tau was provided to us by Karl Stapelfeldt. For IRAS 04016+2610 the NICMOS data were obtained from the HST archive and were re-reduced using the IRAF STSDAS package with the most recent set of reference files available. We created bad pixel masks using the data-quality frames associated with each image, manually masking transient bad pixels, and we also manually masked the bad column 128 (Bergeron & Skinner 1997). The subsets of each image association were registered and combined using the IRAF STSDAS pipeline calnicb task.

In order to subtract the point-spread function (PSF) from each of the re-calibrated HST images, we re-binned the images by a factor of 2 and deconvolved them using PSFs computed with TINYTIM version 4.4 (Krist 1995), which were also oversampled by a factor of 2. The deconvolution was performed within IRAF using the STSDAS lucy task, which uses a Richardson-Lucy based algorithm. The deconvolution process is an iterative one, and we stopped the process after 5 iterations in the F110W filter, 10 iterations in F160W, 15 iterations in F187W, 20 iterations in F205W.

For IRAS 04016+2610 we also obtained an $`I`$ band image using the WIYN 111The WIYN Observatory is a joint facility of the University of Wisconsin-Madison, Indiana University, Yale University, and the National Optical Astronomy Observatories. telescope’s CCD imager. The WIYN observations were obtained on 13/14 December 1995 using the WIYN imager camera, as part of a feasibility study conducted under the WIYN queue observing program. Conditions during the night were good, with some light cirrus and seeing of 0.8 - 1.4 arcsec. Observations were made in each of the VRI bands, consisting of 3 - 300s exposures in each band. Offsets of 10 arcsec were made between each exposure. The data were reduced using the standard IRAF 222IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Associations of Universities for Research in Astronomy, Inc., under cooperative agreement with the NSF. data reduction packages for CCD images. After reduction, the individual images for each exposure were combined into one image (900s total exposure) for each band. Here we report only I-band results (the morphology seen at V and R is similar). The combined I-band image is shown in Figure 4a.

## 3. Model Ingredients

In this section we present the basic ingredients of the 3-D radiation transfer code used to model T Tau and IRAS 04016+2610. The code is comprehensive, and necessarily involves a large number of parameters that in principle can be tuned to match observations. Wherever possible, we constrain our model parameters based on known properties of the sources.

The basic geometry we consider consists of an infalling envelope with bipolar cavities, illuminated from within by a single source or a binary system. For binary illumination, the two components possess circumstellar disks that are not necessarily aligned with one another or with the envelope’s rotation axis.

### 3.1. Envelope Geometry

Of the two sources we model in this paper, T Tau is a known binary (Dyck, Simon, & Zuckerman 1982) and IRAS 04016+2610 shows evidence for binarity with possibly two bipolar outflow cavities (Lucas & Roche 1997). Despite the binary nature of the sources, we adopt the rotational collapse geometry of Terebey, Shu, & Cassen (1984, hereafter TSC) developed for single star systems. We thus in effect assume that the circumbinary density is given by the TSC solution and neglect the effects of binary motion on the envelope structure. The focus of our present study is the 3-D nature of the cavities and envelope illumination; future studies may include circumbinary density distributions predicted from binary star formation models.

The TSC envelope density used in our simulations is given by

$$\rho =\frac{\dot{M}}{4\pi }\left(\frac{GM}{r_c^3}\right)^{1/2}\left(\frac{r}{r_c}\right)^{3/2}\left(1+\frac{\mu }{\mu _0}\right)^{1/2}\left(\frac{\mu }{\mu _0}+\frac{2\mu _0^2r_c}{r}\right)^1,$$
(1)
where $`\dot{M}`$ is the mass infall rate, $`r_c`$ is the centrifugal radius, $`\mu =\mathrm{cos}\theta `$ and $`\mu _0=\mu (r\mathrm{})`$ is determined by

$$\mu _0^3+\mu _0(r/r_c1)\mu (r/r_c)=0.$$
(2)
For our models, we vary $`\dot{M}`$ and $`r_c`$, in the range $`\dot{M}=(210)\times 10^6M_{}`$yr-1 and $`r_c=50100`$AU (Strom 1994; Kenyon et al. 1993b; Whitney, et al. 1997). Our radiation transfer code uses a linear Cartesian grid (§ 3.5) and in order to resolve the TSC density structure we adopt an outer radius for the envelope of 1500AU. We note that in the models of Whitney & Hartmann (1993) and Whitney et al. (1997), outer radii of around 10000AU were adopted. We have found that while our smaller envelopes can reproduce the overall morphology of the scattered light images, larger envelopes produce redder colors (see the IRAS 04016+2610 models in §4.2).

### 3.2. Cavity Geometries

We explore three dimensional geometries by carving cavities in the axisymmetric TSC envelope. By virtue of being 2-D, previous modeling (e.g., Whitney & Hartmann 1993; Whitney et al. 1997; Lucas & Roche 1997; Stapelfeldt et al. 1998) used bipolar cavities that were aligned with the envelope’s rotational axis. In our models we adopt the same cavity shapes as used previously, but now allow the cavities to be offset from the envelope’s symmetry axis and also allow for the existence of multiple sets of cavities, as might result in a binary system.

We consider three cavity shapes in our modeling: streamline, cylindrical, and curved. In the streamline cavity shape, the cavity walls trace out the trajectory of an infalling particle as given in the TSC envelope solution. This cavity shape is approximately conical at large radii, but at small radii becomes more curved. The curvature increases the cross-section of the cavity walls exposed to direct starlight and results in more scattering than a purely conical shape and consequently a brighter cavity (see Whitney & Hartmann 1993). Cylindrical cavities have straight walls parallel to the outflow axis. Curved cavities take the form $`z=a\omega ^b`$ (e.g., Stapelfeldt et al. 1998b), where $`\omega =\sqrt{x^2+y^2}`$ and $`a`$ and $`b`$ are parameters that can be adjusted to provide a good match to the data or to approximate an opening angle of a molecular outflow. We also allow for the presence of dust within the cavities.

### 3.3. Illumination

As stated earlier, T Tau is a binary and there are indications that IRAS 04016+2610 may also be a binary. Therefore we assume the illumination is by binary star + disk systems as used by Wood, Crosas, & Ghez (1999) in their scattered light model of GG Tau’s circumbinary disk. In this model we assumed two point sources of illumination, but with the radiation angular distribution of a star + disk system. In this approximation we do not explicitly treat the radiation transfer through the circumstellar disks, but assign an angular distribution for the flux emerging from the star+disk systems. The angular dependence of the illumination of these “star + disk” systems is displayed in Figure 1, showing the small emergent flux level at high inclinations due to the dense equatorial disk (Whitney & Hartmann 1992). The angular distribution of the illuminating sources was set up so that the sources would be at the center of the cavities, with the polar direction of the illumination pattern aligned with the cavity axis.

### 3.4. Dust Opacity and Scattering Parameters

As with previous modeling of Taurus protostars (Kenyon et al. 1993b; Whitney et al. 1997; Wood et al. 1998) we adopt the albedo, opacity, and scattering parameters for a Kim, Martin, & Hendry (1994, hereafter KMH) grain mixture, typical of grains in the interstellar medium (ISM). The scattering phase function is assumed to be a Henyey-Greenstein phase function (Henyey & Greenstein 1941) with asymmetry parameter $`g`$. In more evolved disk systems the spectral index and level of the long wavelength continuum emission provides strong evidence for dust grains that are larger than ISM grains (Beckwith et al. 1990; Beckwith & Sargent 1991). However, Whitney et al. (1997) found that the KMH parameters gave a good match to the colors and sizes of images Class I sources, but that they underestimated the polarization level. The KMH wavelength dependent opacity and scattering parameters we adopt in this paper are presented in Table 1.

### 3.5. Radiation Transfer

We perform the radiation transfer with a Monte Carlo code that allows for arbitrary sources of emission and multiple scattering within an arbitrary geometry. The scattering code is based on that described by Code & Whitney (1995), but has now been modified to run on a three dimensional linear Cartesian grid (Wood & Reynolds 1999) and includes forced first scattering (Witt 1977) and a “peeling off” procedure (Yusef-Zadeh, Morris, & White 1984). These modifications enable us to construct model images of three dimensional systems from specified viewing angles very efficiently. We have tested our 3D code against results from the axisymmetric code of Whitney & Hartmann (1993), finding excellent agreement between the two codes.

The output of each Monte Carlo simulation is a two dimensional image including the spatially resolved polarization (Stokes $`I`$, $`Q`$, and $`U`$ images). In order to compare our models with the observations we convolve our high resolution models with appropriate PSFs using the TINYTIM software (Krist 1995). From these images we calculate the integrated polarization and by combining simulations at different wavelengths we calculate colors for comparison with observations. The colors we compute are subject to a normalization set by the intrinsic color of the illuminating source(s). The observed color is $`JH=(J^{}H^{})+(J_0H_0)`$, where $`J_0H_0`$ is the source color and $`J^{}H^{}`$ is the color of the Monte Carlo simulation. In their simulations, Whitney et al. (1997) adopted source colors for a typical M0 T Tauri star plus reprocessing disk from Kenyon et al. (1994), $`J_0H_0=0.77`$ and $`H_0K_0=0.43`$. In our models of IRAS 04016+2610 we adopt the same source colors. However, for T Tau the spectral type is K0 or earlier and since we are modeling only the WFPC2 F555W image the intrinsic source colors are not required.

## 4. Results

In this section we present our models for T Tau and IRAS 04016+2610. As our aim in this paper is to demonstrate the usefulness of our 3-D modeling method, we report here on those models that best reproduce the observed scattered light images, and discuss the salient properties of those models. As we discuss in detail for each source, we constrain the parameters of our models based on the available data from the literature. Nonetheless, the uniqueness of the models presented is not demonstrated. These models should be regarded as existence proofs and as suggestions for more detailed 3-D investigations of these and other Class I sources in the future.

### 4.1. T Tau

#### 4.1.1 Previous 2-D models and evidence for source asymmetry

T Tau is a binary star with one component seen directly and the secondary ($`0.7`$ arcsec distant) only detected at IR and longer wavelengths (Dyck, Simon, & Zuckerman 1982; Ghez et al. 1991). The bolometric luminosities of the sources are estimated to be $`9L_{}`$ for the primary (Ghez et al. 1991) and up to $`14L_{}`$ for the secondary (Cohen, Emerson, & Beichman 1989) . Therefore, although undetected at optical wavelengths, the secondary may be emitting a lot of radiation which we only observe when it is scattered into our sightline. The optical light of the secondary may be obscured by either its location within a dense circumbinary envelope, or occulted by its own circumstellar disk edge-on to our line of sight. The WFPC2 scattered light pattern shows departures from axisymmetry with one side of the cavity appearing brighter than the other (Figure 2a; Stapelfeldt et al. 1998b). A similar asymmetric scattered light pattern is seen in the near-IR adaptive optics images of Roddier et al. (1999).

The radial distribution of the near-IR scattered light surrounding T Tau (Weintraub et al. 1992) was modeled by Whitney & Hartmann (1993) adopting a near pole-on orientation of a TSC envelope ($`\dot{M}=2\times 10^6M_{}`$yr-1, $`r_c=100`$AU) and evacuated bipolar cavities. They found that cavities with streamline shapes give too steep a gradient in the radial distribution of the scattered light and suggested that curved cavities could provide a better fit to the data. Stapelfeldt et al. (1998b) modeled T Tau’s scattered light pattern assuming single scattering in a TSC envelope and investigated what viewing angle and cavity curvature is required to match their HST images. They estimated the inclination of the system to be in the range $`i=25^{}`$ to $`45^{}`$ with curved cavities ($`z=ar^b`$, $`a=5\times 10^3`$, $`b=1.95`$). Figure 2b shows the axisymmetric single source model of Stapelfeldt et al. (1998b). The scattered light distribution in the model is necessarily symmetric, and does not reproduce the strong asymmetry observed in the WFPC2 image.

#### 4.1.2 New 3-D model

A binary system with shadowing caused by misaligned circumstellar disks may yield a scattered light pattern such as that observed in T Tau. Indeed, Stapelfeldt et al. (1998b) suggested that the unseen secondary may play an important role in forming the scattered light pattern. We have therefore constructed a scattered light model for T Tau incorporating illumination from a binary star system and allowed the alignment of the circumstellar disks to depart from coplanarity.

Our modeling of T Tau begins with the axisymmetric TSC envelope model of Whitney & Hartmann (1993) with the curved cavities derived by Stapelfeldt et al. (1998b), and we adopt $`i=40^{}`$ which results in a visual extinction through the envelope to the primary of $`A_V2`$, close to the observed value. However, in our model the axisymmetry is broken by the illuminating sources. The illumination is by a binary star + disk system where the primary source is aligned with the rotational axis of the envelope and the secondary (larger bolometric luminosity) has its disk aligned so that its polar axis points towards one side of the cavity wall. The direction of the secondary’s axis is chosen to agree with the large outflow discovered by Reipurth, Bally, & Devine (1997), which is believed to originate from the secondary.

We assume the secondary emits 1.5 times more photons than the primary, based on their estimated bolometric luminosity ratio. The secondary is located 90AU from the primary in the envelope’s equatorial plane with its disk rotational axis at $`\varphi =120^{}`$, $`\theta =70^{}`$ (where $`\theta `$ and $`\varphi `$ are the usual polar coordinates). For our adopted viewing angle the secondary disk is edge-on to our line of sight and consequently its direct light is highly extincted at optical wavelengths, consistent with its very red spectral energy diistribution. Furthermore the disk orientation blocks the light from the secondary from reaching one side of the cavity. In this model we assume the primary and secondary are within the same cavity, created by a combination of outflows from both sources.

In Figures 2c and 2d we show new 3-D scattered light models that better reproduce the observed scattered light image of T Tau. We rotated these model images by $`20^{}`$ to allow for comparison with the WFPC2 data. The model in Fig. 2c shows the scattered light pattern when the illumination is by misaligned circumstellar disks where the angular dependence of the illumination for each star + disk system is that of the WH92 model in Fig. 1. The non-axisymmetric illumination yields an asymmetric nebula with one side brighter than the other, similar to what is observed. Direct light from the secondary is blocked due to our adopted edge-on orientation of its circumstellar disk, consistent with the lack of observed optical light from the secondary.

In Fig. 2d we allowed the secondary’s illumination pattern to depart from that of the WH92 star + disk of Fig. 1 in order to investigate whether we could increase the asymmetry in our scattered light models. In this model, the illumination pattern from the secondary is more concentrated towards the polar regions (see Fig. 1). Such an illumination may arise if the secondary’s circumstellar material is arranged in a highly flared disk (e.g., Chiang & Goldreich 1999) or a small envelope (Roddier et al. 2000). With the more extreme illumination pattern, the asymmetry of the model scattered light nebula is enhanced.

#### 4.1.3 Efficacy and limitations of 3-D model

Our model for T Tau resembles the WFPC2 image in that it reproduces the asymmetry in the scattered light nebula. Support for our model may be found in observations which show a very complex pattern of outflows in CO (Schuster, Harris, & Russell 1997; Leverault 1998) and molecular hydrogen (Herbst, et al. 1996; van Langevelde et al. 1994). As with IRAS 04016+2610 (see below), the plethora of emission knots and filaments is suggestive of a precessing jet and/or multiple outflow sources. The most compelling support for the orientation of the secondary in our models is that its disk axis is aligned with the direction of the outflow discovered by Reipurth et al. (1997) and attributed to the secondary.

In our model the circumstellar disk axes are offset from one another, the primary viewed at a low inclination and the secondary oriented close to edge-on to our line of sight. This is not an unreasonable assumption, given the fact that the secondary in the T Tau system is undetected at optical wavelengths. Also, a system with such a geometry has been observed in the binary HK Tau where the primary is seen directly and the secondary’s scattered light nebula indicates a highly inclined circumstellar disk (Stapelfeldt et al. 1998a; Koresko 1998).

The inclusion of a secondary source not only breaks the asymmetry in the scattered light image, but also increases the contrast between the bright cavity walls and the region interior to this. In their axisymmetric modeling, Stapelfeldt et al. (1998b) produced a contrast no greater than about two between the nebular arc and the inner region, whereas the WFPC2 image shows a contrast in excess of four. Our models can produce contrasts in excess of four as shown in Fig. 3. This figure shows horizontal intensity cuts at various offsets across the axisymmetric model (Fig. 2b) and our new model (Fig. 2d). The larger brightness contrast of our non-axisymmetric model (solid line) compared to the axisymmetric model (dashed line) is clearly evident.

While the overall agreement between our $`V`$ model and the WFPC2 image is very good, we cannot reproduce the thickness of T Tau’s nebular arc. This discrepancy is also present in the model of Stapelfeldt et al. (1998b). In our model we have assumed that T Tau’s IR companion is a disk source seen almost edge-on. However, as it is not detected in the optical, but bright in the IR, we may be viewing at grazing incidence through the upper layers of a flared disk or envelope. As such, detailed 3D modeling of the spectral energy distribution will help in refining our model. We are currently developing techniques for such 3D radiative equilibrium modeling within our Monte Carlo codes (Bjorkman & Wood 2001).

### 4.2. IRAS 04016+2610

#### 4.2.1 Previous 2-D models and evidence for source asymmetry

IRAS 04016+2610 is a Taurus Class I source with a bolometric luminosity of $`3.7L_{}`$ (Kenyon & Hartmann 1995) and displays a very complex scattered light pattern (Whitney et al. 1997; Lucas & Roche 1997; Padgett et al. 1999; Lucas et al. 2000). The HH knots mapped by Gomez et al. (1997) appear scattered throughout the environs, indicative of a precessing jet or multiple outflows. Evidence in favor of multiple outflows is found in the C18O molecular line profiles of Zhou et al. (1994) that show multiple peaks both red and blueshifted. Perhaps the most compelling evidence of multiple outflows is the near-IR imaging polarimetry presented by Lucas & Roche (1997, Fig. 11) and Whitney et al. (1997, Fig. 3). Lucas & Roche presented contours of polarized flux (essentially mapping the scattered light) which they interpreted as arising from scattering in an envelope with two sets of dust-filled bipolar outflow cavities with projected axes that are almost perpendicular to one another. Our WIYN $`I`$ band image (shown in Fig. 4a) displays bright elongated limbs that are aligned with the directions of the suspected outflow cavities mapped by Lucas & Roche (1997).

In previous axisymmetric modeling of IRAS 04016+2610, Kenyon et al. (1993a,b) obtained a good match to the spectral energy distribution and low resolution near-IR imagery adopting a TSC envelope with bipolar streamline cavities viewed at $`i=60^{}`$. Whitney et al. (1997) adopted a similar model ($`\dot{M}=5\times 10^6M_{}`$yr-1, $`r_c=50`$AU) for their near-IR imaging polarimetry and found that while it yielded colors close to those observed, to match the observed nebular pattern would require a fully 3D model (as further demonstrated by Lucas & Roche 1997; Lucas et al. 2000). A single wide outflow cavity cannot reproduce the distinct V-shape in the WIYN image and imaging polarimetry. A single cavity results in a conical scattered light image, but cannot reproduce the large brightness contrast between the arms of the V-shape and the interior region, e.g., see the axisymmetric T Tau model (Fig. 2b) and also Whitney & Hartmann (1993) and Stapelfeldt et al. (1998b).

#### 4.2.2 New 3-D model

In our 3D modeling of IRAS 04016+2610 we adopt the TSC envelope geometry and test the multiple cavity scenario proposed by Lucas & Roche (1997). The TSC envelope has $`\dot{M}=5\times 10^6M_{}`$yr-1, $`r_c=50`$AU (Whitney et al. 1997). Our scattered light model shown in Fig. 4 has two sets of bipolar cylindrical, dust-filled cavities carved in the TSC envelope. The primary cavity (aligned along the rotation axis of the TSC envelope) has a radius of 50AU and the secondary cavity has a radius of 80AU. The viewing angle is $`i=50^{}`$ relative to the primary cavity. Both cavities have a $`J`$ band optical depth of 2 measured from the center of the envelope along the cavity axis. To match the bright, elongated limbs in the WIYN image, we rotated the axis of the second bipolar cavity to lie along $`\varphi =85^{}`$, $`\theta =80^{}`$ (see Fig. 5). In this model the two bipolar cavities are almost perpendicular as noted by Lucas & Roche (1997). The scattered light images have been rotated clockwise by $`220^{}`$ so that the scattered light pattern formed by the primary cavity (aligned with the TSC rotation axis) is aligned with the primary outflow cavity identified by Lucas & Roche (1997), see the labeling in Lucas et al. (2000, Fig. 2).

Lucas & Roche (1997) interpreted their data as possible evidence for a binary carving two outflow cavities. As such, we have adopted an unresolved binary for the illuminating sources. The NICMOS images of Padgett et al. (1999) and centimeter observations of Lucas et al. (2000) did not resolve a binary companion. We also constucted a model with the two outflow cavities described above, but illumination from a single source. Because the disk sources depart from isotropy only along equatorial directions (Fig. 1), we found that the scattered light pattern from a single source is not distinguished from that of the unresolved binary model of Fig. 4. Therefore, unlike our T Tau model, we infer that for IRAS 04016+2610 it is the envelope structure that is most important in producing the image morphology and not asymmetric illumination.

#### 4.2.3 Efficacy and limitations of 3-D model

Our model reproduces the overall size and morphology of the WIYN and NICMOS images very well. In the $`I`$ band model, the envelope is optically thick and it is the cavities that dominate the image morphology. However, at optically thinner longer wavelengths the envelope structure is important in shaping the resulting image. Our model has a single TSC envelope with the primary cavity along the rotation axis and the second cavity almost perpendicular to this. If the source is indeed a binary then the envelope geometry is likely different from TSC. Indeed, recent modeling of imaging spectroscopy by Hogerheijde (2001) adopted a 2000AU radius axisymmetric rotating disklike structure. Such a “disk” should not be confused with the smaller circumstellar disks we adopt and are seen in other Classical T Tauri stars (e.g., Burrows et al. 1996; Padgett et al. 1999). The TSC rotational collapse geometry produces a large equatorial density enhancement which may correspond to Hogerheijde’s parameterized equatorial structure. As such, the equatorial enhancement of TSC and Hogerheijde’s “disk” represent a different evolutionary stage from disks around Classical T Tauri stars (Shu et al. 1987). A more detailed model of IRAS 04016+2610 than we have undertaken should simultaneously address the envelope geometry and velocity structure. Nevertheless, our models which incorporate the TSC parameters determined by Kenyon et al. (1993b) and Whitney et al. (1997) provide support for Lucas & Roche’s (1997) dual bipolar cavity interpretation of their imaging polarimetry.

Our model colors are $`J^{}H^{}=1.34`$, $`H^{}K^{}=0.79`$ giving actual colors of $`JH=2.11`$, $`HK=1.22`$ when combined with the adopted intrinsic colors (§ 2.5). The observed colors for IRAS 04016+2610 are $`JH=1.95`$, $`HK=2.02`$ (Kenyon & Hartman 1995; Whitney et al. 1997). The failure of both our models and those of Whitney et al. (1997) to reproduce these colors (both models predict $`[JH]>[HK]`$) may indicate that the intrinsic source colors are different than those we have adopted.

While we have been succesful at reproducing the image morphology, our models suffer from the same problem as those of Whitney et al. (1997) in that they predict lower polarizations than observed. Whitney et al. attribute this to the low intrinsic polarization of the KMH grains. In our case, the resolved polarization is lower than observed due to the dust grains, but we also find much lower levels for the total unresolved polarization in our models. We typically predict an integrated polarization of around 3% in the near-IR, whereas the observed polarization ranges from $`P_J=21\%`$ to $`P_K=8\%`$ (Whitney et al. 1997). Another problem is that neither our models or those of Whitney et al. (1997) can reproduce the observed pattern of polarization vectors close to the source(s). So, while the dual cavity scenario proposed by Lucas & Roche (1997) gives a good match to the overall image morphology, its failure to reproduce the observed polarization level and pattern in IRAS 04016+2610 indicates that the details of the model need to be refined. In particular, polarization due to scattering off aligned grains can yield larger polarization values and we are investigating grain alignment in protostars to reconcile models with the large observed polarizations (Whitney & Wolff 2001, in preparation).

## 5. Summary of Results

We have constructed 3D models for the scattered light nebulae associated with IRAS 04016+2610, and T Tau. Our models assumed an axisymmetric infalling circumstellar envelope, but the axisymmetry was broken by multiple outflow cavities and non-axisymmetric illumination from binary stars.

The asymmetry seen in the nebula surrounding T Tau was modeled by introducing a secondary star+disk system with its rotation axis misaligned from the axis of the TSC envelope. The orientation of the secondary was chosen so that its disk was oriented almost edge-on to our line of sight, and was selected to be perpendicular to the Herbig-Haro outflow attributed to the secondary (Reipurth et al. 1997). This orientation results in the secondary being undetected in our scattered light models as its disk extincts the direct optical light.

We reproduced the scattered light pattern associated with IRAS 04016+2610 with a model which had two sets of bipolar outflow cavities in the TSC envelope. This model affirms the notion of Lucas & Roche (1997), who interpreted their near-IR imaging polarimetry as arising from scattering in an envelope with orthogonal bipolar cavities. Further evidence for multiple cavities in the envelope is found in the complex pattern of outflow emission knots surrounding this source (Gomez et al. 1997). Our model fails, however, to reproduce the level and pattern of the near-IR polarization and we suspect that a more accurate treatment of polarization from scattering off aligned grains may provide a better match to the observations (Whitney & Wolff 2001, in preparation).

In the models presented in this paper, we adopted the axisymmetric TSC solution for the circumstellar envelopes. This may not be appropriate as T Tau is a known binary and IRAS 04016+2610 shows strong evidence for binarity. Also, other studies indicate that protostellar environments have material arranged with a shallower radial gradient than TSC (e.g., Barsony & Chandler 1993). The extension of our Monte Carlo scattering codes to three dimensions now enables us to incorporate more complex circumstellar envelope structures, such as those predicted from binary star formation models. However, this is among the first modeling investigations to explore the effects of three dimensional cavity structures and illuminations and has proved successful in reproducing the scattered light morphologies and in testing several scenarios proposed by others to explain the observations. While we have been fairly succesful in reproducing the scattered light morphologies, spectral energy distribution modeling incorporating the 3D structures will provide a crucial test of our models. This requires current axisymmetric radiation transfer radiative equilibrium codes to be extended to three dimensions and we are developing such techniques (Bjorkman & Wood 2001).

We acknowledge financial support from NASA’s Long Term Space Astrophysics Research Program, NAG5 6039 (KW), NAG5 8412 (BW); the National Science Foundation, AST 9909966 (BW and KW); the HST Archival Research Program AR-08367.01-97A, (BW, KW, KS); NASA’s Astrophysical Data Program, NAG5 3904 (MJW). KSB is a Cottrell Scholar of the Research Corporation, and gratefully acknowledges their support. DSS was funded by the NSF Research Experience for Undergraduates program at the Harvard-Smithsonian Center for Astrophysics and thanks Christine Jones, Jonathan McDowell, and Tania Ruiz for the generous donation of time and self to make the SAO Intern Program as rewarding as it has been. We thank Karl Stapelfeldt for providing the WFPC2 image of T Tau and Ted von Hippel for assistance with the WIYN queue observations.


text:
# 1 𝜒² contours for 𝑚₁₆=1500 GeV (Left) and 𝑚₁₆=2000 GeV (Right). The shaded region is excluded by the chargino mass limit 𝑚_𝜒̃⁺>103 GeV.

CERN-TH/2001-176

July 2001

Predictions for Higgs and SUSY spectra from SO(10) Yukawa Unification with $`\mu >0`$

T. Blažek, R. Dermíšek and S. Raby

Department of Physics, University of Southampton, Southampton, UK
and Faculty of Mathematics and Physics, Comenius Univ., Bratislava, Slovakia
Department of Physics, The Ohio State University, 174 W. 18th Ave., Columbus, Ohio 43210

Abstract

> We use $`t,b,\tau `$ Yukawa unification to constrain SUSY parameter space. We find a narrow region survives for $`\mu >0`$ (suggested by $`bs\gamma `$ and the anomalous magnetic moment of the muon) with $`A_01.9m_{16}`$, $`m_{10}1.4m_{16}`$, $`m_{16}12003000`$  GeV and $`\mu ,M_{1/2}100500`$  GeV . Demanding Yukawa unification thus makes definite predictions for Higgs and sparticle masses.

Minimal supersymmetric \[SUSY\] SO(10) grand unified theories \[GUTs\] have many profound features : all fermions in one family sit in one 16 dimensional spinor representation; the two Higgs doublets of the minimal supersymmetric standard model sit in one 10 dimensional fundamental representation, and gauge coupling unification at a GUT scale $`M_G3\times 10^{16}`$  GeV fits well with the low energy data . In addition in the simplest version of SO(10) the third generation Yukawa couplings are given by a single term in the superpotential $`W=\lambda \mathbf{\hspace{0.33em}16\hspace{0.33em}10\hspace{0.33em}16}`$ resulting in Yukawa unification $`\lambda _t=\lambda _b=\lambda _\tau =\lambda _{\nu _\tau }\lambda `$ and a prediction for $`M_t`$ with large $`\mathrm{tan}\beta 50`$ . 111Note, GUT scale threshold corrections to this Yukawa unification boundary condition are naturally small ( $`<`$ 1% ), since they only come at one loop from the SO(10) gauge sector and the third generation - Higgs Yukawa coupling . This is in contrast to GUT scale threshold corrections to gauge coupling unification which may be significant, coming from doublet/triplet splitting in the Higgs sector and, even more importantly, the SO(10) breaking sector which typically has many degrees of freedom. The data requires $`ϵ_3`$ $`=\frac{\alpha _3(M_G)\alpha _G(M_G)}{\alpha _G(M_G)}`$ -4% . This beautiful result is however marred by potentially large weak scale threshold corrections

$$m_b(M_Z)=\lambda _b(M_Z)\frac{v}{\sqrt{2}}cos\beta (1+\mathrm{\Delta }m_b^{\stackrel{~}{g}}+\mathrm{\Delta }m_b^{\stackrel{~}{\chi }^+}+\mathrm{\Delta }m_b^{\stackrel{~}{\chi }^0}+\mathrm{\Delta }m_b^{\mathrm{log}}).$$

For $`\mu >0`$the gluino term is positive and in most regions of SUSY parameter space it is the dominant contribution to $`\mathrm{\Delta }m_b`$. Reasonable fits prefer $`\mathrm{\Delta }m_b`$$`<0`$; hence Yukawa unification is easy to satisfy with $`\mu <0`$.

The decay $`bs\gamma `$ and the muon anomalous magnetic moment also get significant corrections proportional to $`\mathrm{tan}\beta `$. The SUSY contribution to $`bs\gamma `$ comes from one loop diagrams similar to those contributing to the bottom mass. The chargino term typically dominates and has opposite sign to the SM and charged Higgs contributions, thus reducing the branching ratio for $`\mu >0`$. This is necessary to fit the data since the SM contribution is somewhat too big. $`\mu <0`$ would on the other hand constructively add to the branching ratio and is problematic. In addition, the recent measurement of the anomalous magnetic moment of the muon $`a_\mu ^{NEW}=(g2)/2=43(16)\times 10^{10}`$ also favors $`\mu >0`$ . Thus it is important to confirm that Yukawa unification can work consistently with $`\mu >0`$.

In this paper we assume exact Yukawa unification and search, using a $`\chi ^2`$ analysis, for regions of SUSY parameter space with $`\mu >0`$providing good fits to the low energy data. We show that Yukawa unification dramatically constrains the Higgs and SUSY spectra. These results are sensitive to the SUSY breaking mechanism.

It is much easier to obtain EWSB with large $`\mathrm{tan}\beta `$ when the Higgs up/down masses are split ($`m_{H_u}^2<m_{H_d}^2`$. In our analysis we consider two particular schemes we refer to as universal and D term splitting. In the first case the third generation squark and slepton soft masses are given by the universal mass parameter $`m_{16}`$ , and only Higgs masses are split: $`m_{(H_u,H_d)}^2=m_{10}^2(1\mathrm{\Delta }m_H^2)`$. In the second case we assume D term splitting, i.e. that the D term for $`U(1)_X`$ is non-zero, where $`U(1)_X`$ is obtained in the decomposition of $`SO(10)SU(5)\times U(1)_X`$. In this second case, we have $`m_{(H_u,H_d)}^2=m_{10}^22D_X`$, $`m_{(Q,\overline{u},\overline{e})}^2=m_{16}^2+D_X`$, $`m_{(\overline{d},L)}^2=m_{16}^23D_X`$. The universal case does not at first sight appear to be similarly well motivated. It is quite clear however that in any SUSY model the Higgs bosons are very special. R parity is used to distinguish Higgs superfields from quarks and leptons. In addition, a supersymmetric mass term $`\mu `$ with value of order the weak scale is needed for an acceptable low energy phenomenology. Since $`\mu `$ is naturally of order $`M_G`$, one needs some symmetry argument why it is suppressed. Of course, if the Higgs are special, then it is reasonable to assume splitting of Higgs, while maintaining universal squark and slepton masses. This may be achieved by GUT scale threshold corrections to the soft SUSY breaking scalar masses . Here we present the most compelling mechanism . In $`SO(10)`$, neutrinos necessarily have a Yukawa term coupling active neutrinos to the “sterile” neutrinos present in the 16. In fact for $`\nu _\tau `$ we have $`\lambda _{\nu _\tau }\overline{\nu }_\tau LH_u`$ with $`\lambda _{\nu _\tau }=\lambda _t=\lambda _b=\lambda _\tau \lambda `$. In order to obtain a tau neutrino with mass $`m_{\nu _\tau }0.05`$ eV (consistent with atmospheric neutrino oscillations), the “sterile” $`\overline{\nu }_\tau `$ must obtain a Majorana mass $`M_{\overline{\nu }_\tau }10^{13}`$ GeV. Moreover, since neutrinos couple to $`H_u`$ (and not to $`H_d`$) with a fairly large Yukawa coupling (of order 0.7), they naturally distinguish the two Higgs multiplets. With $`\lambda =0.7`$ and $`M_{\overline{\nu }_\tau }=10^{13}`$ GeV, we obtain a significant GUT scale threshold correction with $`\mathrm{\Delta }m_H^210`$%, remarkably close to the value needed to fit the data. At the same time, we obtain a small threshold correction to Yukawa unification $`2.5`$% (for more details see ).

Our analysis is a top-down approach with 11 input parameters, defined at $`M_G`$, varied to minimize a $`\chi ^2`$ function composed of 9 low energy observables. The 11 input parameters are: $`M_G,\alpha _G(M_G),`$ $`ϵ_3`$; the Yukawa coupling $`\lambda `$, and the 7 soft SUSY breaking parameters $`\mu ,M_{1/2},A_0,\mathrm{tan}\beta `$, $`m_{16}^2,m_{10}^2`$, $`\mathrm{\Delta }m_H^2`$($`D_X`$) for universal (D term) case. We use two (one)loop renormalization group \[RG\] running for dimensionless (dimensionful) parameters from $`M_G`$ to $`M_Z`$ and complete one loop threshold corrections at $`M_Z`$ . We require electroweak symmetry breaking using an improved Higgs potential, including $`m_t^4`$ and $`m_b^4`$ corrections in an effective 2 Higgs doublet model below $`M_{stop}`$ . The $`\chi ^2`$ function includes the 9 observables; 6 precision electroweak data $`\alpha _{EM},G_\mu ,\alpha _s(M_Z),M_Z,M_W,\rho _{NEW}`$ and the 3 fermion masses $`M_{top},m_b(m_b),M_\tau `$. The experimental values used for the low energy observables are given in the table.

Fig. 1 shows the constant $`\chi ^2`$ contours for $`m_{16}=1500`$ and $`2000`$  GeV in the case of universal squark and slepton masses. We find acceptable fits ($`\chi ^2<3`$) for $`A_01.9m_{16}`$, $`m_{10}1.4m_{16}`$ and $`m_{16}1.2`$ TeV. The best fit is for $`m_{16}2000`$  GeV with $`\chi ^2<1`$. Note, electroweak symmetry breaking in this region of parameter space requires splitting Higgs up/down masses, $`\mathrm{\Delta }m_H^2`$$``$ O(13%). This range of soft SUSY parameters is consistent with solution (B) of Olechowski and Pokorski . In the table we present the input parameters, the fits and the predicted Higgs and SUSY spectra for two representative points with universal squark and slepton masses and the best fit value for D term splitting. We have not presented the contour plots for D term splitting since as can be seen from the best fit point in the table, the bottom quark mass is poorly fit in this case and $`\chi ^2>5`$. Recall, since we have 11 input parameters and only 9 observables, we consider such poor fits unacceptable.

Fig. 2 gives the constant $`m_b(m_b)`$ and $`\mathrm{\Delta }m_b`$contours for $`m_{16}=2000`$  GeV . We see that the best fits, near its central value, are found with $`\mathrm{\Delta }m_b`$$`2`$%. Why does Yukawa unification only work in this narrow region of SUSY parameter space? The log corrections $`\mathrm{\Delta }m_b^{\mathrm{log}}46`$% (total contribution from gluino, neutralino, chargino and electroweak loops) are positive and they must be cancelled in order to obtain $`\mathrm{\Delta }m_b`$$`2`$ %. The leading mass insertion corrections proportional to $`\mathrm{tan}\beta `$ are approximately given by

$$\mathrm{\Delta }m_b^{\stackrel{~}{g}}\frac{2\alpha _3}{3\pi }\frac{\mu m_{\stackrel{~}{g}}}{m_{\stackrel{~}{b}}^2}tan\beta \mathrm{and}\mathrm{\Delta }m_b^{\stackrel{~}{\chi }^+}\frac{\lambda _t^2}{16\pi ^2}\frac{\mu A_t}{m_{\stackrel{~}{t}}^2}tan\beta .$$

They can naturally be as large as 40%. The chargino contribution is typically opposite in sign to the gluino, since $`A_t`$ runs to an infrared fixed point $`M_{1/2}`$(see for example, Carena et al. ). Hence in order to cancel the positive contribution of both the log and gluino contributions, a large negative chargino contribution is needed. This can be accomplished for $`A_t>m_{\stackrel{~}{g}}`$ and $`m_{\stackrel{~}{t}_1}<<m_{\stackrel{~}{b}_1}`$. The first condition can be satisfied for $`A_0`$ large and negative, which helps pull $`A_t`$ away from its infrared fixed point. The second condition is also aided by large $`A_t`$. However in order to obtain a large enough splitting between $`m_{\stackrel{~}{t}_1}`$ and $`m_{\stackrel{~}{b}_1}`$, large values of $`m_{16}`$ are needed. Note, that for universal scalar masses, the lightest stop is typically lighter than the sbottom. We typically find $`m_{\stackrel{~}{b}_1}3m_{\stackrel{~}{t}_1}`$. On the other hand, D term splitting with $`D_X>0`$ gives $`m_{\stackrel{~}{b}_1}m_{\stackrel{~}{t}_1}`$. As a result in the case of universal boundary conditions excellent fits are obtained for top, bottom and tau masses; while for D term splitting the best fits give $`m_b(m_b)4.59`$ GeV .

Finally in Fig. 3 we show the constant light Higgs mass contours for $`m_{16}=1500`$ and $`2000`$  GeV (solid lines) with the constant $`\chi ^2`$ contours overlayed (dotted lines). Yukawa unification for $`\chi ^21`$ clearly prefers light Higgs mass in a narrow range, 112 - 118 GeV . In this region the CP odd, the heavy CP even Higgs and the charged Higgs bosons are also quite light (see fit 2 in the table).222 It would be interesting to see how sensitive our results, for Higgs masses, may be to alternative electroweak symmetry breaking approximations. In this paper we have used the effective 2 Higgs doublet analysis of  with an estimated 3  GeV uncertainty in Higgs masses. This approximation may be particularly well suited to the light Higgs spectrum we obtain in our analysis. The alternative scheme, in which the Higgs tadpoles are evaluated at a scale of order $`M_{stop}`$  is however more frequently used in the literature. In addition we find the mass of $`\stackrel{~}{t}_1(150250)`$ GeV, $`\stackrel{~}{b}_1(450650)`$ GeV, $`\stackrel{~}{\tau }_1(200500)`$ GeV, $`\stackrel{~}{g}(6001200)`$ GeV, $`\stackrel{~}{\chi }^+(100250)`$ GeV, and $`\stackrel{~}{\chi }^0(80170)`$ GeV. All first and second generation squarks and sleptons have mass of order $`m_{16}`$. The light stop and chargino may be visible at the Tevatron. With this spectrum we expect $`\stackrel{~}{t}_1\stackrel{~}{\chi }^+b`$ with $`\stackrel{~}{\chi }^+\stackrel{~}{\chi }_1^0\overline{l}\nu `$ to be dominant. Lastly $`\chi _1^0`$ is the LSP and possibly a good dark matter candidate .

The region of SUSY parameter space preferred by Yukawa unification may be consistent with a supergravity mechanism for SUSY breaking at $`M_{Pl}`$ with RG running from $`M_{Pl}`$ to $`M_G`$ (see for example Murayama et al. ). It however cannot be obtained with gauge mediated or gaugino mediated SUSY breaking mechanisms where $`A_0=0`$ at zeroth order. It may also be obtained in anomally mediated schemes but in this case one still has to worry about slepton masses squared and also the fact that in this case, since the gluino and chargino masses have opposite sign, it is difficult to fit both $`bs\gamma `$ and $`a_\mu `$.

In a future paper  we present the sparticle spectrum in more detail and consequences for Tevatron searches. We discuss the sensitivity of our results to small GUT scale threshold corrections to Yukawa unification with both universal and D term Higgs up/down splitting. We also check the robustness of the Higgs spectrum by artificially adjusting the CP odd Higgs mass using a penalty in $`\chi ^2`$. We find that $`\chi ^2`$ increases by at most 40% for any $`m_{A^0}`$ less than $`350`$ GeV. The light Higgs mass $`m_h`$ is rather insensitive to the value of $`m_{A^0}`$; whereas $`m_H,m_{H^+}`$ are linearly dependent on $`m_{A^0}`$. We also consider constraints resulting from the processes $`bs\gamma `$, $`B_s\mu ^+\mu ^{}`$, $`a_\mu ^{NEW}`$ and the proton lifetime in a semi model independent way. We shall only make a few short comments here. In order to fit $`bs\gamma `$ we find that the coefficient $`C_7^{MSSM}C_7^{SM}`$ (see for example, Eqn. 9 in Ref. ) with the chargino term dominating by a factor of order 5 over all other contributions. This is due to the light stop $`\stackrel{~}{t}_1`$. In fact, $`bs\gamma `$ is more sensitive to $`m_{\stackrel{~}{t}_1}`$ than $`m_b(m_b)`$. Fitting the central value $`B(bs\gamma )=2.96\times 10^4`$  requires a heavier $`\stackrel{~}{t}_1`$ with $`(m_{\stackrel{~}{t}_1})_{MIN}500`$ GeV; significantly larger than the range which provides the best fits to $`m_b`$. We now find $`m_b(m_b)_{MIN}4.3`$. Moreover no other sparticle masses are affected. The process $`B_s\mu ^+\mu ^{}`$ provides a lower bound on $`m_{A^0}200`$ GeV (see recent work of ). 333We thank K.S. Babu and C. Kolda for discussions. However this has only a minor impact on $`\chi ^2`$ as discussed above. We recall that proton decay experiments prefer values of $`m_{16}>2000`$ GeV and $`m_{16}>>M_{1/2}`$ (see ref. ). This is in accord with the range of SUSY parameters found consistent with third generation Yukawa unification. There is however one experimental result which is not consistent with either Yukawa unification or proton decay and that is the anomalous magnetic moment of the muon. Large values of $`m_{16}1200`$ GeV lead to very small values for $`a_\mu ^{NEW}16\times 10^{10}`$. Hence a necessary condition for Yukawa unification is that forthcoming BNL data  and/or a reanalysis of the strong interaction contributions to $`a_\mu ^{SM}`$ will significantly decrease the discrepancy between the data and the standard model value of $`a_\mu `$.

In summary, most of the results of our analysis including only third generation fermions remain intact when incorporating flavor mixing. The light Higgs mass and most sparticle masses receive only small corrections. The lightest stop mass increases, due to $`bs\gamma `$. Nevertheless there is still a significant $`\stackrel{~}{t}_1\stackrel{~}{t}_2`$ splitting and $`m_{\stackrel{~}{t}_1}<<m_{\stackrel{~}{b}_1}`$. The $`A^0,H,H^+`$ masses are necessarily larger in order to be consistent with $`B_s\mu ^+\mu ^{}`$ , which suggests that this process should be observed soon; possibly at Run II of the Tevatron. Finally, the central value for $`a_\mu ^{NEW}`$ must significantly decrease. The “smoking guns” of SO(10) Yukawa unification, presented in this letter, should be observable at Run III of the Tevatron or at LHC. Also, in less than a year we should have more information on $`a_\mu ^{NEW}`$.

In previous works Yukawa unification with $`\mu >0`$was not possible.444 In ref.  the corrections $`\mathrm{\Delta }m_b^{\mathrm{log}}`$ were neglected and Yukawa unification with $`\mu >0`$was obtained in a large region of parameter space with no constraints on the Higgs and SUSY spectra. Pierce et al.  assume $`\mathrm{\Delta }m_H^2`$$`=0`$ and, as a result, they are not able to enter the region of SUSY parameter space consistent with both EWSB and Yukawa unification. Baer et al.  also cannot obtain Yukawa unification with $`\mu >0`$. This is because they use D term splitting for Higgs up/down which as discussed typically leads to sbottom lighter than stop.

While completing this article, the paper by Baer and Ferrandis  appeared which confirmed our results  on the existence of a preferred region of SUSY parameter space consistent with Yukawa unification and $`\mu >0`$. Their results however require significant GUT threshold corrections to $`\lambda _t=\lambda _b=\lambda _\tau `$ of order 8 - 15% which helps them obtain $`m_{\stackrel{~}{t}_1}<m_{\stackrel{~}{b}_1}`$. They also claim that better fits are obtained with D term splitting than with the universal splitting case. We believe the latter is only true because the authors do not allow their SUSY parameters, in particular $`m_{16}`$ and $`A_0`$, to explore the region of parameter space discussed in  and this paper. 555In version 3 of this paper, the authors now cover a larger region of soft SUSY breaking parameter space. Now solutions with D term splitting are only found with $`28`$% GUT threshold corrections to Yukawa unification. Moreover, they still claim that D term splitting works better than the universal case. We have no explanation for this discrepancy.

Acknowledgements

S.R. and R.D. are partially supported by DOE grant DOE/ER/01545-816. S.R. graciously acknowledges the Alexander von Humboldt Foundation, the Universität Bonn and the Theory Division at CERN for partial support. R.D. and S.R. thank the Physics Department at the Universität Bonn and the Theory Division at CERN for their kind hospitality while working on this project. R.D. and S.R. thank A. Dedes for his help and collaboration in the early stages of this analysis.


text:
# Variable gamma-ray emission from the Be/X-ray transient A0535+26?

## 1 Introduction

The EGRET instrument of the late Compton Gamma Ray Observatory detected 271 point-like gamma-ray sources at energies $`E>100`$ MeV during its lifetime (Hartman et al. 1999). The majority of these sources still remain unidentified. Excluding 6 possible artifacts, there are 75 sources at low galactic latitudes $`|b|<10^{}`$ still not conclusively identified with objects seen at lower energies (Romero et al. 1999). Several of these low-latitude, presumably galactic sources seem to display significant levels of variability over timescales of months (Tompkins 1999, Torres et al. 2001a, b). A few particular cases of variable gamma-ray sources at low latitudes have been recently studied by Tavani et al. (1997, 1998), Paredes et al. (2000) and Punsly et al. (2000), but there is not yet consensus on their nature. The variety of behaviours displayed by the sources seems to suggest that there exist more than a single class of gamma-ray emitting objects in the Galaxy. To establish the nature of these peculiar objects is one of the most urgent problems of high-energy astrophysics.

In this paper we present a study of the gamma-ray source 3EG J0542+2610. We shall show that the only known object within the 95% confidence location contour of the source capable of generating the observed gamma-ray emission is the Be/X-ray transient A0535+26. We shall argue that the gamma-rays are produced during the accretion disk formation and loss phases in each orbit, through hadronic interactions between relativistic protons accelerated in an electrostatic gap in the pulsar magnetosphere and the matter in the disk. The structure of the paper is as follows. In the next section we summarize the main characteristics of 3EG J0542+2610 and the results of our analysis of the radio field within the EGRET location error contours. We then present the peculiarities of the system A0535+26. In Section 4 we present a model for the gamma-ray production based on the accelerator gap model by Cheng & Ruderman (1989, 1991). The main difference between our treatment and the accreting pulsar gamma-ray model of Cheng et al. (1991) lies in the introduction of time variability through the variation of the matter content of the disk, which is a transient structure according to recent observations.

## 2 The gamma-ray source 3EG J0542+2610 and the surrounding radio field

The best estimated position of the gamma-ray source 3EG J0542+2610 is at ($`l`$, $`b`$)$``$(182.02,$`1.99`$). Its 95% confidence location contour overlaps with the shell-type supernova remnant (SNR) G180.0-1.7 (Romero et al. 1999, Torres et al. 2001c). The gamma-ray flux for the combined EGRET viewing periods is $`(14.7\pm 3.2)\mathrm{\hspace{0.33em}10}^8`$ ph cm-2 s-1. Since the flux is highly variable on timescales of months, it should be originated in a compact object and not in the extended SNR. In Figure 1 we show the EGRET light curve. The source switches between periods of clear detections and periods when only upper bounds to the flux can be determined. The variability analysis of Torres et al. (2001b) assigns to 3EG J0542+2610 a variability index $`I=3.16`$, which means that the source variability level is $`4.32\sigma `$ above the average (spurious) variability of all known gamma-ray pulsars. Tompkins’ (1999) variability index for this source ($`\tau =0.7`$) also indicates that the source is variable.

If a pulsar origin is discarded due to the high variability and the steep spectral index ($`\mathrm{\Gamma }=2.67\pm 0.22`$), we are left with two main possibilities: 1) the source is a background, unnoticed gamma-ray blazar seen through the galactic plane, or 2) it is a galactic compact object with an energy budget high enough as to generate significant gamma-ray emission and, at the same time, has not the stable properties usually associated with isolated pulsars.

In the standard model for gamma-ray blazars, the high energy photons are produced by ultra-relativistic electrons through inverse Compton interactions with seed lower energy photons that can come from external sources (e.g the accretion disk) or from the jet itself (see Krolik 1999 and references therein). The model is supported by the fact that all identified gamma-ray AGNs are strong flat-spectrum radio sources (Mattox et al. 1997).

In Figure 2, lower panel, we show a 1.4-GHz VLA map made with data from the NVSS Sky Survey (Condon et al. 1998), where all point-like radio sources within the 68% confidence contour of 3EG J0542+2610 can be seen. The measured characteristics of these 29 sources are listed in Table 1. Most of them have no entry in any existing point source catalog. In those cases where we were able to find positional counterparts at other frequencies we have estimated the spectral indices, which are also shown in the table. Most of these sources are very weak, at the level of a few mJy. No strong (at Jy level), flat or nonthermal source is within the location error box of the gamma-ray source. The strongest radio source (No. 23 in our table) has a rather steep spectrum and is a factor $`10`$ below the minimum flux density of firm gamma-ray blazar identifications given by Mattox et al. (1997). The nature of this source is not clear at present; it could be a background weak radio quasar. The fact that it is not seen at X-rays seems to argue against a galactic microquasar or any other kind of accreting source.

At X-ray energies the most significant source within the EGRET error box is the X-ray transient A0535+26, which is discussed in the next section. This source does not present significant radio emission and consequently it cannot be seen in our maps. We have indicated its position with a star symbol in Fig. 2, middle panel. We also show in this figure the direction of the proper motion of the system, as determined by Lee Clark & Dolan (1999). The upper panel shows the entire radio field as determined from radio observations with Effelsberg 100-m single dish telescope at 1.408 GHz (data from Reich et al. 1997). The middle panel present an enhanced image obtained at 2.695 GHz with the same telescope (Fuerst et al. 1990), where the gamma-ray location probability contours have been superposed (Hartman et al. 1999). We have processed these large-scale images using the background filtering techniques described by Combi et al. (1998).

## 3 The X-ray transient A0535+26

A0536+26 is a Be/X-ray transient where the compact object is a 104s pulsar in an eccentric orbit around the B0III star HDE 245770 (Giovanelli & Sabau Graziati 1992). Be stars are rapidly rotating objects which eject mass irregularly forming gaseous disks on their equatorial planes. If there is a compact companion in a close orbit, accretion from the star can result in strong X-ray emission. In the case of A0535+26, strong and recurrent X-ray outbursts are observed with a period of 111 days, which has been identified with the orbital period (Giovanelli & Sabau Graziati 1992). It is generally agreed that these outbursts occur when the accretion onto the neutron star increases at the periastron passage. The average ratio of the X-ray luminosity at the periastron to that of the apoastron is $`100`$ (Janot-Pacheco et al. 1987). In Table 2 we list the main characteristics of the A0535+26 system.

During a major outburst of A0535+26 in 1994, the BATSE instrument of the Compton Gamma Ray Observatory detected a broad quasi-periodic oscillation (QPO) in the power spectra of the X-ray flux (Finger et al. 1996). The QPO component was detected during 33 days, with a central frequency that was well correlated with both the hard X-ray flux and neutron star spin-up rate inferred from pulse timing. Finger et al.’s (1996) observations are the first clear evidence that an accretion disk is formed during giant outbursts.

Using the simultaneous variations of the spin and the QPO frequencies in the context of the beat frequency model (Alpar & Shaham 1985), Li (1997) has determined the evolution of the ratio $`\xi `$ of the inner accretion disk radius to the Alfvén radius. He found that during the initial rise of the outbursts $`\xi `$ quickly increased from $`0.6`$ to $`1`$, indicating a transition of the accretion process from spherical accretion before the outbursts to disk accretion during the high X-ray luminosity phase.

Prior to the direct evidence for a transient accretion disk in A0535+26, Motch et al. (1991) have already suggested, on the basis of an analysis of the long-term X-ray, UV, and optical history of the system, that two different types of interactions exist between the Be star and the pulsar, accordingly with the dynamical state of the highly variable circumstellar envelope. “Normal” outbursts would occur for high equatorial wind velocities ($`200`$ km s-1), whereas “giant” outburst would result from lower equatorial wind velocities ($`2080`$ km s-1) which allow the formation of a transient accretion disk.

It is interesting to notice that, during the initial stage of the 1994 outburst, when QPOs were found in A0535+26, the gamma-ray source 3EG J0542+2610 was detected by EGRET with a flux of $`(39.1\pm 12.5)\mathrm{\hspace{0.33em}10}^8`$ ph cm-2 s-1 (viewing period 321.1: February 8 - 15, 1994). But when A0535+26 was at the peak of its X-ray luminosity, on Ferbuary 18, the gamma-ray source was not detected (viewing period 321.5). The gamma-ray emission seems to have been quenched precisely when the accretion disk was well-formed and maximally rotating. In the next section we present a model that can account for the gamma-ray production in A0535+26 necessary to explain the EGRET source 3EG J0542+2610 and that is in agreement with our present knowledge of the Be/X-ray transient source.

## 4 Model

Our purpose in this section is to show that there exist a plausible mechanism that could explain the gamma-ray emission of the source 3EG J0542+2610 as originated in the X-ray transient A0535+26. This mechanism should be capable of predicting the observed gamma-ray flux, the variability in the lightcurve, the fact that no gamma-ray emission was observed on February 18, when A0535+26 was at the peak of the X-ray outburst, but also that it was positively detected in the previous days when the X-ray flux was rising, and finally the fact that A0535+26 is not a non-thermal radio source. This latter restriction seems to suggest a hadronic origin for the gamma rays. Otherwise, relativistic electrons should also produce synchrotron radio emission.

Cheng & Ruderman (1989, 1991) have studied the disturbances produced in the magnetosphere of an accreting pulsar when the Keplerian disk rotates more rapidly than the star. As in the case of equal angular velocities, when $`\mathrm{\Omega }_{}<\mathrm{\Omega }_\mathrm{d}`$, inertial effects of electrons ($``$) and ions ($`+`$) lead to a complete charge separation around the “null surface” $`𝛀_{}𝐁=0`$. However, since now the equatorial plasma between the inner accretion disk radius $`r_0`$ and the Alfvén radius $`r_\mathrm{A}`$ co-rotates with the disk, whereas the rest of plasma co-rotates with the star, an electrostatic gap with no charge at all is created around the “null surface” (see Figure 3). In this gap $`𝐄𝐁0`$ and a strong potential drop is established (see Cheng & Ruderman 1991 for details).

Cheng & Ruderman (1989) have shown that the potential drop along the magnetic field lines through the gap is:

$`\mathrm{\Delta }V_{\mathrm{max}}`$ $``$ $`{\displaystyle \frac{B_\mathrm{s}R^3\mathrm{\Omega }_\mathrm{d}(r_0)}{r_0c}}`$ (1)
$``$ $`\mathrm{4\hspace{0.33em}10}^{14}\beta ^{5/2}\left({\displaystyle \frac{M}{M_{}}}\right)^{1/7}R_6^{4/7}L_{37}^{5/7}B_{12}^{3/7}\mathrm{V},`$
where $`B_\mathrm{s}`$ is the neutron star’s surface dipole magnetic field, $`R`$ is its radius, $`r_0`$ is the inner accretion disk radius, $`M`$ is the compact star mass, and $`L_{37}`$ is the X-ray luminosity in units of $`10^{37}`$ erg s-1. The radius and magnetic field in the second expression are in units of $`10^6`$ cm and $`10^{12}`$ Gauss, respectively. The parameter $`\beta 2r_0/r_\mathrm{A}`$ is twice the ratio of the inner accretion disk radius to the Alfvén radius. For the particular case of A0535+26 we can adopt $`\beta 1`$, according to the estimates by Li (1997).

Protons entering into the gap are accelerated up to energies above $`E_pe\mathrm{\Delta }V_{\mathrm{max}}`$ whereas inverse Compton and curvature losses would limit the energy gain of electrons and positrons to lower values. The maximum current that can flow through the gap can be determined from the requirement that the azimuthal magnetic field induced by the current cannot exceed that of the initial field $`𝐁`$ (Cheng & Ruderman 1989):

$`J_{\mathrm{max}}`$ $``$ $`cB_\mathrm{s}R^3r_0^2`$
$``$ $`\mathrm{1.5\hspace{0.33em}10}^{24}\beta ^2\left({\displaystyle \frac{M}{M_{}}}\right)^{2/7}R_6^{1/7}L_{37}^{4/7}B_{12}^{1/7}\mathrm{esu}\mathrm{s}^1.`$
As it can be seen from Fig. 3, the proton current flux is directed from the polar cap of the star, where the accreting material impacts producing strong X-ray emission and abundant ions to be captured by the gap, towards the accretion disk. Electrons move in the opposite sense in order to cancel any net charge current flow, keeping in this way $`𝐄.𝐁0`$ in the gap.

The collision of the relativistic proton beam into the disk will produce hadronic interactions with copious $`\pi ^0`$ production (Cheng et al. 1990). These $`\pi ^0`$ will quickly decay into $`\gamma `$-rays that could escape only if the disk density is sufficiently low. Otherwise, they will be absorbed by the matter and then re-emitted as X-rays. The interaction of relativistic protons with a thin hydrogen layer has been studied by Cheng et al. (1990) in the context of their model for the Crab pulsar. The $`\pi ^0`$-decay $`\gamma `$-rays will only escape insofar as the column density of the disk would be $`\mathrm{\Sigma }100`$ g cm-2. From Novikov & Thorne (1972), we have:

$`\mathrm{\Sigma }`$ $``$ $`\mathrm{7\hspace{0.33em}10}^2\beta ^{3/5}\left({\displaystyle \frac{M}{M_{}}}\right)^{3/35}R_6^{30/35}L_{37}^{27/35}B_{12}^{12/35}\mathrm{g}\mathrm{cm}^2.`$
In a system like A0535+26, the column density of the disk will evolve with time following the variations in the X-ray luminosity. All other physical parameters in Eq. (4) remain constant along the orbital period. If at some point near the periastron, when the X-ray luminosity is close to its maximum, the column density exceeds the value $`100`$ g cm-2, the medium will be no longer transparent to gamma-rays and the gamma-emission will be quenched (we are assuming that the disk is seen from the opposite side to that where the current impacts). This could explain the fact that A0535+26 has not been detected by EGRET at the peak of its X-ray outbursts on February 18, 1994. The detection, instead, was clear a few days before, during the previous viewing period.

The expected hadronic $`\gamma `$-ray flux from A0535+26 on Earth will be:

$$F(E_\gamma )=\frac{1}{\mathrm{\Delta }\mathrm{\Omega }D^2}_{E_\gamma }^{\mathrm{}}\frac{d^2N_p(E_{\pi ^0})}{dtdE_{\pi ^0}}𝑑E_p,$$
(4)
where $`D`$ is the distance to the source, $`\mathrm{\Delta }\mathrm{\Omega }`$ is the beaming solid angle and:

$$\frac{d^2N_p(E_{\pi ^0})}{dtdE_{\pi ^0}}=\frac{2J_{\mathrm{max}}}{e}\frac{\mathrm{\Sigma }<m_{\pi ^0}>}{m_p}\frac{d\sigma (E_p,E_{\pi ^0})}{dE_{\pi ^0}}$$
(5)
is the differential production rate of neutral pions from $`pp`$ interactions. In this latter expression, $`d\sigma (E_p,E_{\pi ^0})/dE_{\pi ^0}`$ is the differential cross-section for the production of $`\pi ^0`$-mesons of energy $`E_{\pi ^0}`$ by a proton of energy $`E_p`$ in a $`pp`$ collision, and $`<m_{\pi ^0}>`$ is the mean multiplicity of $`\pi ^0`$.

For calculation purposes we shall adopt the cross section given by Dermer (1986a, b) and a beaming factor $`\mathrm{\Delta }\mathrm{\Omega }/4\pi 0.3`$ as used by Cheng et al. (1991). We shall evaluate the total gamma-ray luminosity between 100 MeV and 20 GeV of A0535+2610 at an early epoch $`t`$ of the disk formation, when the X-ray luminosity ratio is $`L_X^{\mathrm{max}}/L_X^t10`$. If the absorption feature observed near 110 keV is a cyclotron line, the polar magnetic field of the neutron star results $`\mathrm{9.5\hspace{0.33em}10}^{12}`$ G (Finger et al. 1996). In our calculations we adopt this value along with $`\beta 1`$, $`R_6=1`$ and $`M=1.4M_{}`$ (Janot-Pacheco et al. 1987). We then obtain that the column density of the disk is $`42.5`$ g cm-2 at this stage and, consequently, gamma-rays are not absorbed in the disk material. The potential drop in the electrostatic gap results $`\mathrm{\Delta }V_{\mathrm{max}}\mathrm{2.5\hspace{0.33em}10}^{13}`$V whereas the proton current deposited into the disk is $`N_p=J_{\mathrm{max}}/e\mathrm{4.8\hspace{0.33em}10}^{32}`$ s-1. Using Eq. (4) for photons in EGRET’s energy range we obtain:

$$F(100\mathrm{MeV}<E_\gamma <20\mathrm{GeV})\mathrm{5.1\hspace{0.33em}10}^8\mathrm{ph}\mathrm{cm}^2\mathrm{s}^1.$$
(6)
This flux is consistent with the lower EGRET detections. As the accretion onto the neutron star increases the gamma-ray flux also increases, reaching a maximum when $`\mathrm{\Sigma }100`$ g cm-2. At this point the flux has grown by a factor $`4`$ (notice the dependence on $`J_{\mathrm{max}}`$ in addition to $`\mathrm{\Sigma }`$) and photon absorption into the disk becomes important quenching the radiation. The gamma-ray source, then, is not detected when the peak of the X-ray luminosity occurs. The additional gamma-ray flux expected from secondary electrons and positrons radiating in the the strong magnetic fields anchored in the disk can explain even higher EGRET fluxes (Cheng et al. 1991). A more detailed analysis in this sense will be presented elsewhere, including spectral considerations. The observed spectrum depends on the proton injection spectrum at the gap. Strong shocks near the polar gap can accelerate the protons up to a power law which will be preserved through the propagation across the gap and imprinted into the $`\pi ^0`$ gamma-ray spectrum. Additionally, $`e^+e^{}`$ cascades induced in the disk material can result in a relativistic bremsstrahlung contribution. Electromagnetic shower simulations are in progress in order to determine this contribution (Romero et al., in preparation).

## 5 Discussion

The model here outlined does not imply strict periodicity because of the highly chaotic nature of the Be stellar winds, which can significantly vary on short timescales producing strong changes in the accretion rate. Notwithstanding, the general prediction that the peak of the gamma-ray emission should not be coincident with the maximum X-ray luminosity during a given outburst can be used to test the general scenario proposed here. At present, the poor time resolution of the gamma-ray lightcurve does not allow correlation studies. EGRET data for the best sampled X-ray outburst (in February 1994) consist only of two viewing periods, as it was mentioned, and in one of them the source was not detected. In the future, however, new instruments like GLAST could provide the tools for these kind of investigations.

TeV emission should be produced in the accretion disk according to Eq. (4), although degradation effects during the propagation in the strong magnetic and photon fields around the accretion disk could suppress much of it. The magnetic field at the inner accretion disk radius is (e.g. Cheng et al. 1991):

$$B(r_0)=\mathrm{3\hspace{0.33em}10}^5\beta ^3B_{12}^{5/7}R_6^{9/7}L_{37}^{6/7}(M/M_{})^{3/7}\mathrm{G}$$
(7)
which for the typical values of A0535+26 yields $`B10^4`$ G. TeV gamma-rays can be absorbed in the field through one-photon pair production above the threshold given by (e.g. Bednarek 1993):

$$\chi \frac{E_\gamma B\mathrm{sin}\theta }{2m_ec^2B_{\mathrm{cr}}}=\frac{1}{15},$$
(8)
where $`\theta `$ is the angle between the photon momentum and the magnetic field and $`B_{\mathrm{cr}}`$ is the critical magnetic field given by $`B_{\mathrm{cr}}=m^2c^3/e\mathrm{}\mathrm{4.4\hspace{0.33em}10}^{13}`$ G. This process, then, suppress gamma-ray photons with energies higher than $`300`$ TeV in A0535+26. Since this value is well above the energy of the protons that impact on the disk, we find that absorption in the magnetic field should not occur.

However, TeV photons with lower energies should be absorbed by two-photon pair production in the accretion disk X-ray photosphere. Calculations by Bednarek (1993) show that the optical depth quickly goes to values above 1 for disks with luminosities $`L_{37}1`$. The opacity effects can reach even GeV energies leading to a steepening in the spectrum respect to what is expected from a pure pion-decay mechanism. The fact that the observed spectrum in 3EG J0542+4610 has an index $`\mathrm{\Gamma }2.7`$ seems to support the idea that important absorption is occurring in the X-ray photosphere of this source. Additional spectral modifications should be produced by the emission of secondary pairs in the magnetic field close to the disk (Cheng et al. 1991, Romero et al., in preparation).

In their original model, Cheng & Ruderman (1989) suggested that if the disk is sufficiently dense then the gamma-rays could be produced only in a moving low-density “window”. This window would collimate a pencil beam of gamma-rays aligned with the magnetic axis. Consequently, a pulsed emission could exist with the same period of the X-ray source (104 sec in the case of A0235+26). Small changes in the gamma-ray period might be produced by the radial motion of the “window”. It would be interesting to test whether this pulses are present in A0235+26. However, since the disk is a transient structure and only would excess the critical density during a few days per orbit, the number of photon counts in EGRET data are too low to allow a periodicity analysis as in the case of isolated and stable gamma-ray pulsars. Instruments with higher sensitivity like GLAST could sum up over several viewing periods in order to look for these features.

## 6 Summary

We have shown that the Be/X-ray transient system A0535+26 can be also a transient gamma-ray source under very reasonable assumptions. The existence of a QPO phenomenon detected by BATSE during the 1994 X-ray outburst provided direct evidence of the formation of a transient accretion disk near the periastron passage. In the beat frequency model for QPO the Keplerian orbital frequency of the material at the inner edge of the accretion disk should exceed the spin frequency of the neutron star. Cheng & Ruderman (1989) have shown that in such a circumstance an electrostatic gap is open in the magnetosphere around the null surface determined by $`𝛀_{}𝐁=\mathrm{𝟎}`$. This gap can accelerate protons up to energies of tens of TeV, producing a hadronic current that impacts into the accretion disk generating $`\pi ^0`$ gamma-rays. The transient character of the disk makes the high-energy gamma radiation highly variable, as was found in the observed gamma-ray flux evolution. A specific prediction of the model is the suppression of gamma-ray emission when the column density of the disk exceeds a critical value, near the peak of the X-ray luminosity. Future GeV and TeV observations of this source with instruments of high temporal resolution, like GLAST or 5@5 (see Aharonian et al. 2001), could be used to test the proposed model and, if it is basically correct, to probe the evolution of the matter content on the accretion disk in this extraordinary X-ray binary.

###### Acknowledgements.

GER thanks valuable comments by Dr. Pablo H. Posta. An anonymous referee made interesting suggestions that lead to a substantial improvement of the manuscript. This work was supported by CONICET (PIP 0430/98), ANPCT (PICT 03-04881) and Fundación Antorchas.